r.texture: code optimization

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@47698 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.texture/h_measure.h

@@ -23,25 +23,9 @@
  *
  *****************************************************************************/
 
-float h_measure(int **, int, int, int, int);
-float f1_asm(float **, int);
-float f2_contrast(float **, int);
-float f3_corr(float **, int);
-float f4_var(float **, int);
-float f5_idm(float **, int);
-float f6_savg(float **, int);
-float f7_svar(float **, int, double);
-float f8_sentropy(float **, int);
-float f9_entropy(float **, int);
-float f10_dvar(float **, int);
-float f11_dentropy(float **, int);
-float f12_icorr(float **, int);
-float f13_icorr(float **, int);
-float f14_maxcorr(float **, int);
-float *vector(int, int);
-float **matrix(int, int, int, int);
-void results(char *, float *);
-void simplesrt(int, float[]);
-void mkbalanced(float **, int);
-void reduction(float **, int);
-void hessenberg(float **, int, float[], float[]);
+float h_measure(int);
+void alloc_vars(int, int);
+int set_vars(int **grays, int curr_rrow, int curr_col,
+                int size, int offset, int t_d);
+int set_angle_vars(int angle, int have_px, int have_py, int have_sentr,
+                   int have_pxpys, int have_pxpyd);

raster/r.texture/main.c

@@ -29,46 +29,71 @@
 #include <grass/glocale.h>
 #include "h_measure.h"
 
-static const char *suffixes[56] = {
-    "_ASM_0", "_ASM_45", "_ASM_90", "_ASM_135",
-    "_Contr_0", "_Contr_45", "_Contr_90", "_Contr_135",
-    "_Corr_0", "_Corr_45", "_Corr_90", "_Corr_135",
-    "_Var_0", "_Var_45", "_Var_90", "_Var_135",
-    "_IDM_0", "_IDM_45", "_IDM_90", "_IDM_135",
-    "_SA_0", "_SA_45", "_SA_90", "_SA_135",
-    "_SV_0", "_SV_45", "_SV_90", "_SV_135",
-    "_SE_0", "_SE_45", "_SE_90", "_SE_135",
-    "_Entr_0", "_Entr_45", "_Entr_90", "_Entr_135",
-    "_DV_0", "_DV_45", "_DV_90", "_DV_135",
-    "_DE_0", "_DE_45", "_DE_90", "_DE_135",
-    "_MOC-l_0", "_MOC-l_45", "_MOC-l_90", "_MOC-l_135",
-    "_MOC-2_0", "_MOC-2_45", "_MOC-2_90", "_MOC-2_135",
-    "_MCC_0", "_MCC_45", "_MCC_90", "_MCC_135"
+struct menu
+{
+    char *name;			/* measure name */
+    char *desc;			/* menu display - full description */
+    char *suffix;		/* output suffix */
+    char useme;			/* calculate this measure if set */
+    int idx;			/* measure index */
+};
+
+/* modify this table to add new measures */
+static struct menu menu[] = {
+    {"asm",      "Angular Second Moment",    "_ASM",   0,  0},
+    {"contrast", "Contrast",                 "_Contr", 0,  1},
+    {"corr",     "Correlation",              "_Corr",  0,  2},
+    {"var",      "Variance",                 "_Var",   0,  3},
+    {"idm",      "Inverse Diff Moment",      "_IDM",   0,  4},
+    {"sa",       "Sum Average",              "_SA",    0,  5},
+    {"se",       "Sum Entropy",              "_SE",    0,  6},
+    {"sv",       "Sum Variance",             "_SV",    0,  7},
+    {"entr",     "Entropy",                  "_Entr",  0,  8},
+    {"dv",       "Difference Variance",      "_DV",    0,  9},
+    {"de",       "Difference Entropy",       "_DE",    0, 10},
+    {"moc1",     "Measure of Correlation-1", "_MOC-1", 0, 11},
+    {"moc2",     "Measure of Correlation-2", "_MOC-2", 0, 12},
+    {NULL, NULL, NULL, 0, -1}
 };
 
+static int find_measure(const char *measure_name)
+{
+    int i;
+
+    for (i = 0; menu[i].name; i++)
+	if (strcmp(menu[i].name, measure_name) == 0)
+	    return i;
+
+    G_fatal_error(_("Unknown measure <%s>"), measure_name);
+
+    return -1;
+}
+
 int main(int argc, char *argv[])
 {
     struct Cell_head cellhd;
-    char *name, *result, *filename;
-    unsigned char *outrast;
+    char *name, *result;
+    char **mapname;
+    FCELL **fbuf;
+    int n_measures, n_outputs, *measure_idx;
     int nrows, ncols;
-    int row, col, i, j;
+    int row, col, first_row, last_row, first_col, last_col;
+    int i, j;
     CELL **data;		/* Data structure containing image */
     DCELL *dcell_row;
     struct FPRange range;
     DCELL min, max, inscale;
     FCELL measure;		/* Containing measure done */
-    int t_measure, dist, size;	/* dist = value of distance, size = s. of sliding window */
-    int infd, outfd;
-    int a, c, corr, v, idm, sa, sv, se, e, dv, de, moc1, moc2, mcc;
+    int dist, size;	/* dist = value of distance, size = s. of sliding window */
+    int offset;
+    int have_px, have_py, have_sentr, have_pxpys, have_pxpyd;
+    int infd, *outfd;
     RASTER_MAP_TYPE data_type, out_data_type;
     struct GModule *module;
-    char mapname[GNAME_MAX];
-    struct Option *input, *output, *size_O, *dist_O;
-    struct Flag *flag2, *flag3, *flag4, *flag5,
-	*flag6, *flag7, *flag8, *flag9, *flag10, *flag11,
-	*flag12, *flag13, *flag14, *flag15;
+    struct Option *opt_input, *opt_output, *opt_size, *opt_dist, *opt_measure;
+    struct Flag *flag_ind, *flag_all;
     struct History history;
+    char p[1024];
 
     G_gisinit(argv[0]);
 
@@ -76,143 +101,127 @@ int main(int argc, char *argv[])
     G_add_keyword(_("raster"));
     G_add_keyword(_("algebra"));
     G_add_keyword(_("statistics"));
+    G_add_keyword(_("texture"));
     module->description =
 	_("Generate images with textural features from a raster map.");
 
     /* Define the different options */
 
-    input = G_define_standard_option(G_OPT_R_INPUT);
+    opt_input = G_define_standard_option(G_OPT_R_INPUT);
 
-    output = G_define_option();
-    output->key = "prefix";
-    output->type = TYPE_STRING;
-    output->required = YES;
-    output->gisprompt = "new,cell,raster";
-    output->description = _("Prefix for output raster map(s)");
+    opt_output = G_define_option();
+    opt_output->key = "prefix";
+    opt_output->type = TYPE_STRING;
+    opt_output->required = YES;
+    opt_output->gisprompt = "new,cell,raster";
+    opt_output->description = _("Prefix for output raster map(s)");
 
-    size_O = G_define_option();
-    size_O->key = "size";
-    size_O->key_desc = "value";
-    size_O->type = TYPE_INTEGER;
-    size_O->required = NO;
-    size_O->description = _("The size of sliding window (odd and >= 3)");
-    size_O->answer = "3";
+    opt_size = G_define_option();
+    opt_size->key = "size";
+    opt_size->key_desc = "value";
+    opt_size->type = TYPE_INTEGER;
+    opt_size->required = NO;
+    opt_size->description = _("The size of sliding window (odd and >= 3)");
+    opt_size->answer = "3";
 
     /* Textural character is in direct relation of the spatial size of the texture primitives. */
 
-    dist_O = G_define_option();
-    dist_O->key = "distance";
-    dist_O->key_desc = "value";
-    dist_O->type = TYPE_INTEGER;
-    dist_O->required = NO;
-    dist_O->description = _("The distance between two samples (>= 1)");
-    dist_O->answer = "1";
-
-    /* Define the different flags */
-
-    /* "Normalized" unused in the code ??? 
-       flag0 = G_define_flag() ;
-       flag0->key         = 'N' ;
-       flag0->description = _("Normalized") ;
-       flag0->guisection  = _("Features");
-     */
-
-    flag2 = G_define_flag();
-    flag2->key = 'a';
-    flag2->description = _("Angular Second Moment");
-    flag2->guisection = _("Features");
-
-    flag3 = G_define_flag();
-    flag3->key = 'c';
-    flag3->description = _("Contrast");
-    flag3->guisection = _("Features");
-
-    flag4 = G_define_flag();
-    flag4->key = 'k';
-    flag4->description = _("Correlation");
-    flag4->guisection = _("Features");
-
-    flag5 = G_define_flag();
-    flag5->key = 'v';
-    flag5->description = _("Variance");
-    flag5->guisection = _("Features");
-
-    flag6 = G_define_flag();
-    flag6->key = 'i';
-    flag6->description = _("Inverse Diff Moment");
-    flag6->guisection = _("Features");
-
-    flag7 = G_define_flag();
-    flag7->key = 's';
-    flag7->description = _("Sum Average");
-    flag7->guisection = _("Features");
-
-    flag8 = G_define_flag();
-    flag8->key = 'w';
-    flag8->description = _("Sum Variance");
-    flag8->guisection = _("Features");
-
-    flag9 = G_define_flag();
-    flag9->key = 'x';
-    flag9->description = _("Sum Entropy");
-    flag9->guisection = _("Features");
-
-    flag10 = G_define_flag();
-    flag10->key = 'e';
-    flag10->description = _("Entropy");
-    flag10->guisection = _("Features");
-
-    flag11 = G_define_flag();
-    flag11->key = 'd';
-    flag11->description = _("Difference Variance");
-    flag11->guisection = _("Features");
-
-    flag12 = G_define_flag();
-    flag12->key = 'p';
-    flag12->description = _("Difference Entropy");
-    flag12->guisection = _("Features");
-
-    flag13 = G_define_flag();
-    flag13->key = 'm';
-    flag13->description = _("Measure of Correlation-1");
-    flag13->guisection = _("Features");
-
-    flag14 = G_define_flag();
-    flag14->key = 'n';
-    flag14->description = _("Measure of Correlation-2");
-    flag14->guisection = _("Features");
-
-    flag15 = G_define_flag();
-    flag15->key = 'o';
-    flag15->description = _("Max Correlation Coeff");
-    flag15->guisection = _("Features");
-
+    opt_dist = G_define_option();
+    opt_dist->key = "distance";
+    opt_dist->key_desc = "value";
+    opt_dist->type = TYPE_INTEGER;
+    opt_dist->required = NO;
+    opt_dist->description = _("The distance between two samples (>= 1)");
+    opt_dist->answer = "1";
+
+    for (i = 0; menu[i].name; i++) {
+	if (i)
+	    strcat(p, ",");
+	else
+	    *p = 0;
+	strcat(p, menu[i].name);
+    }
+    opt_measure = G_define_option();
+    opt_measure->key = "measure";
+    opt_measure->type = TYPE_STRING;
+    opt_measure->required = NO;
+    opt_measure->multiple = YES;
+    opt_measure->options = p;
+    opt_measure->description = _("Textural measurement");
+
+    flag_ind = G_define_flag();
+    flag_ind->key = 's';
+    flag_ind->description = _("Separate output for each angle (0, 45, 90, 135)");
+
+    flag_all = G_define_flag();
+    flag_all->key = 'a';
+    flag_all->description = _("Calculate all textural measurements");
 
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-    name = input->answer;
-    result = output->answer;
-    a = (!flag2->answer);
-    c = (!flag3->answer);
-    corr = (!flag4->answer);
-    v = (!flag5->answer);
-    idm = (!flag6->answer);
-    sa = (!flag7->answer);
-    sv = (!flag8->answer);
-    se = (!flag9->answer);
-    e = (!flag10->answer);
-    dv = (!flag11->answer);
-    de = (!flag12->answer);
-    moc1 = (!flag13->answer);
-    moc2 = (!flag14->answer);
-    mcc = (!flag15->answer);
-    size = atoi(size_O->answer);
-    dist = atoi(dist_O->answer);
-
-    if (a && c && corr && v && idm && sa && sv && se && e && dv && de && moc1
-	&& moc2 && mcc)
-	G_fatal_error(_("Nothing to compute. Use at least one of the flags."));
+    name = opt_input->answer;
+    result = opt_output->answer;
+    size = atoi(opt_size->answer);
+    if (size <= 0)
+	G_fatal_error(_("Size of the moving window must be > 0"));
+    if (size % 2 != 1)
+	G_fatal_error(_("Size of the moving window must be odd"));
+    dist = atoi(opt_dist->answer);
+    if (dist <= 0)
+	G_fatal_error(_("The distance between two samples must be > 0"));
+
+    n_measures = 0;
+    if (flag_all->answer) {
+	for (i = 0; menu[i].name; i++) {
+	    menu[i].useme = 1;
+	}
+	n_measures = i;
+    }
+    else {
+	for (i = 0; opt_measure->answers[i]; i++) {
+	    if (opt_measure->answers[i]) {
+		const char *measure_name = opt_measure->answers[i];
+		int n = find_measure(measure_name);
+
+		menu[n].useme = 1;
+		n_measures++;
+	    }
+	}
+    }
+    if (!n_measures)
+	G_fatal_error(_("Nothing to compute. Use at least one textural measure."));
+	
+    measure_idx = G_malloc(n_measures * sizeof(int));
+    j = 0;
+    for (i = 0; menu[i].name; i++) {
+	if (menu[i].useme == 1) {
+	    measure_idx[j] = menu[i].idx;
+	    j++;
+	}
+    }
+
+    /* variables needed */
+    if (menu[2].useme || menu[11].useme || menu[12].useme)
+	have_px = 1;
+    else
+	have_px = 0;
+    if (menu[11].useme || menu[12].useme)
+	have_py = 1;
+    else
+	have_py = 0;
+    if (menu[6].useme || menu[7].useme)
+	have_sentr = 1;
+    else
+	have_sentr = 0;
+    if (menu[5].useme || menu[6].useme || menu[7].useme)
+	have_pxpys = 1;
+    else
+	have_pxpys = 0;
+    if (menu[9].useme || menu[10].useme)
+	have_pxpyd = 1;
+    else
+	have_pxpyd = 0;
 
     infd = Rast_open_old(name, "");
 
@@ -222,9 +231,35 @@ int main(int argc, char *argv[])
     Rast_get_cellhd(name, "", &cellhd);
 
     out_data_type = FCELL_TYPE;
-    /* Allocate output buffer, use FCELL data_type */
-    outrast = Rast_allocate_buf(out_data_type);
+    /* Allocate output buffers, use FCELL data_type */
+    n_outputs = n_measures;
+    if (flag_ind->answer) {
+	n_outputs = n_measures * 4;
+    }
 
+    fbuf = G_malloc(n_outputs * sizeof(FCELL *));
+    mapname = G_malloc(n_outputs * sizeof(char *));
+    for (i = 0; i < n_outputs; i++) {
+	mapname[i] = G_malloc(GNAME_MAX * sizeof(char));
+	fbuf[i] = Rast_allocate_buf(out_data_type);
+    }
+
+    /* open output maps */
+    outfd = G_malloc(n_outputs * sizeof(int));
+    for (i = 0; i < n_measures; i++) {
+	if (flag_ind->answer) {
+	    for (j = 0; j < 4; j++) {
+		sprintf(mapname[i * 4 + j], "%s%s_%d", result,
+		        menu[measure_idx[i]].suffix, j * 45);
+		outfd[i * 4 + j] = Rast_open_new(mapname[i * 4 + j], out_data_type);
+	    }
+	}
+	else {
+	    sprintf(mapname[i], "%s%s", result,
+	            menu[measure_idx[i]].suffix);
+	    outfd[i] = Rast_open_new(mapname[i], out_data_type);
+	}
+    }
     nrows = Rast_window_rows();
     ncols = Rast_window_cols();
 
@@ -253,6 +288,7 @@ int main(int argc, char *argv[])
     }
 
     /* Read in cell map values */
+    /* TODO: use r.proj cache */
     G_important_message(_("Reading raster map..."));
     for (j = 0; j < nrows; j++) {
 	Rast_get_row(infd, dcell_row, j, DCELL_TYPE);
@@ -271,14 +307,7 @@ int main(int argc, char *argv[])
     Rast_close(infd);
     G_free(dcell_row);
 
-    /* Now raster map is into memory. */
-
-    /* Open image's files, their names show the measure */
-    filename = G_malloc(strlen(result) + 11);
-    strcpy(filename, result);
-    result = filename;
-    while (*result != '\0')
-	result++;
+    /* Now raster map is loaded to memory. */
 
     /* *************************************************************************************************
      *
@@ -288,90 +317,82 @@ int main(int argc, char *argv[])
      *
      ***************************************************************************************************/
 
-    for (t_measure = 0; t_measure < 56; t_measure++)
-	if ((t_measure == 0 && a)
-	    || (t_measure == 4 && c)
-	    || (t_measure == 8 && corr)
-	    || (t_measure == 12 && v)
-	    || (t_measure == 16 && idm)
-	    || (t_measure == 20 && sa)
-	    || (t_measure == 24 && sv)
-	    || (t_measure == 28 && se)
-	    || (t_measure == 32 && e)
-	    || (t_measure == 36 && dv)
-	    || (t_measure == 40 && de)
-	    || (t_measure == 44 && moc1)
-	    || (t_measure == 48 && moc2)
-	    || (t_measure == 52 && mcc))
-	    t_measure += 3;
-	else {
-	    outfd = Rast_open_new(strcat(filename, suffixes[t_measure]),
-				  out_data_type);
-	    *result = '\0';
-
-	    strcpy(mapname, filename);
-	    strcat(mapname, suffixes[t_measure]);
-	    G_important_message(_("Calculating measure #%d <%s> (56 measures available)"),
-				(t_measure + 1), mapname);
-
-	    for (row = 0; row < nrows - (size - 1); row++) {
-		G_percent(row, nrows, 2);
-
-		/*process the data */
-		for (col = 0; col < ncols - (size - 1); col++) {
-
-		    /* Pointer for the s.w. */
-		    data = data + row;
-		    for (j = 0; j < size; j++)
-			*(data + j) += col;
+    offset = size / 2;
+    first_row = first_col = offset;
+    last_row = nrows - offset;
+    last_col = ncols - offset;
+    Rast_set_f_null_value(fbuf[0], ncols);
+    for (row = 0; row < first_row; row++) {
+	for (i = 0; i < n_outputs; i++) {
+	    Rast_put_row(outfd[i], fbuf[0], out_data_type);
+	}
+    }
+    if (n_measures > 1)
+	G_message(_("Calculating %d texture measures"), n_measures);
+    else
+	G_message(_("Calculating %s"), menu[measure_idx[0]].desc);
+    alloc_vars(size, dist);
+    for (row = first_row; row < last_row; row++) {
+	G_percent(row, nrows, 2);
+
+	for (i = 0; i < n_outputs; i++)
+	    Rast_set_f_null_value(fbuf[i], ncols);
+
+	/*process the data */
+	for (col = first_col; col < last_col; col++) {
+
+	    if (!set_vars(data, row, col, size, offset, dist)) {
+		for (i = 0; i < n_outputs; i++)
+		    Rast_set_f_null_value(&(fbuf[i][col]), 1);
+		continue;
+	    }
 
-/***********************************************************************************************
- *
- * Parameter of the h_measure:
- *   - data: structure containing image;
- *   - size: size of the s. w., defined by operator; 
- *           The s.w. can be rectangolar but not make sense.
- *   - t_measure: measure's type required;
- *   - dist: distance required, defined by operator.
- *
- **********************************************************************************************/
-
-		    measure =
-			(FCELL) h_measure(data, size, size, t_measure, dist);
-		    /* The early (size/2) samples take value from (size/2+1)'th sample */
-		    if (col < (size / 2))
-			for (j = 0; j < (size / 2); j++)
-			    ((FCELL *) outrast)[j] = measure;
-		    ((FCELL *) outrast)[col + ((size / 2))] = measure;
-		    /* The last few (size/2) samples take value from nrows-(size/2+1)'th sample */
-		    if (col == (ncols - size))
-			for (j = 0; j <= (size / 2); j++)
-			    ((FCELL *) outrast)[ncols - j] = measure;
-
-		    for (j = (size - 1); j >= 0; j--)
-			*(data + j) -= col;
-		    data = data - row;
+	    /* for all angles (0, 45, 90, 135) */
+	    for (i = 0; i < 4; i++) {
+		set_angle_vars(i, have_px, have_py, have_sentr, have_pxpys, have_pxpyd);
+		/* for all requested textural measures */
+		for (j = 0; j < n_measures; j++) {
+
+		    measure = (FCELL) h_measure(measure_idx[j]);
+
+		    if (flag_ind->answer) {
+			/* output for each angle separately */
+			fbuf[j * 4 + i][col] = measure;
+		    }
+		    else {
+			/* use average over all angles for each measure */
+			if (i == 0)
+			    fbuf[j][col] = measure;
+			else if (i < 3)
+			    fbuf[j][col] += measure;
+			else 
+			    fbuf[j][col] = (fbuf[j][col] + measure) / 4.0;
+		    }
 		}
-		/* The early (size/2) samples take value from (size/2+1)'th sample */
-		if (row == 0)
-		    for (j = 0; j < (size / 2); j++)
-			Rast_put_row(outfd, outrast, out_data_type);
-
-		Rast_put_row(outfd, outrast, out_data_type);
 	    }
-	    /* The last few (size/2) samples take value from nrows-(size/2+1)'th sample */
-	    if ((row >= nrows - (size - 1)) && (row < nrows))
-		for (j = 0; j < (size / 2); j++)
-		    Rast_put_row(outfd, outrast, out_data_type);
+	}
+	for (i = 0; i < n_outputs; i++) {
+	    Rast_put_row(outfd[i], fbuf[i], out_data_type);
+	}
+    }
+    Rast_set_f_null_value(fbuf[0], ncols);
+    for (row = last_row; row < nrows; row++) {
+	for (i = 0; i < n_outputs; i++) {
+	    Rast_put_row(outfd[i], fbuf[0], out_data_type);
+	}
+    }
+    G_percent(nrows, nrows, 1);
 
-	    Rast_close(outfd);
+    for (i = 0; i < n_outputs; i++) {
+	Rast_close(outfd[i]);
 
-	    Rast_short_history(mapname, "raster", &history);
-	    Rast_command_history(&history);
-	    Rast_write_history(mapname, &history);
+	Rast_short_history(mapname[i], "raster", &history);
+	Rast_command_history(&history);
+	Rast_write_history(mapname[i], &history);
+	G_free(fbuf[i]);
+    }
 
-	}
-    G_free(outrast);
+    G_free(fbuf);
     G_free(data);
 
     exit(EXIT_SUCCESS);