From ef6bae8fb597c12283ee4ddf050535e0c33ff52b Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Sat, 15 Jun 2024 04:39:22 -0400 Subject: [PATCH] r.texture: Refactor code, remove global variables (#3785) * Use local variables instead of global variables. * Separate execute part to prepare it for parallelization. * Remove unused variables and defines. * Use pairs of .c and .h files. * Use structs for parameters and flags. --- raster/r.texture/execute.c | 151 ++++++++++++ raster/r.texture/execute.h | 23 ++ raster/r.texture/h_measure.c | 436 ++++++----------------------------- raster/r.texture/h_measure.h | 32 ++- raster/r.texture/main.c | 323 +++++++++----------------- raster/r.texture/matvec.c | 307 ++++++++++++++++++++++++ raster/r.texture/matvec.h | 50 ++++ 7 files changed, 734 insertions(+), 588 deletions(-) create mode 100644 raster/r.texture/execute.c create mode 100644 raster/r.texture/execute.h create mode 100644 raster/r.texture/matvec.c create mode 100644 raster/r.texture/matvec.h diff --git a/raster/r.texture/execute.c b/raster/r.texture/execute.c new file mode 100644 index 00000000000..c8b15a269ef --- /dev/null +++ b/raster/r.texture/execute.c @@ -0,0 +1,151 @@ +#include +#include +#include +#include +#include +#include + +#include "execute.h" + +/* ************************************************************************************************* + * + * Compute of the matrix S.G.L.D. (Spatial Gray-Level Dependence Matrices) + *or co-occurrence matrix. The image is analyzed for piece, every piece is + *naming moving window (s.w.). The s.w. must be square with number of size's + *samples odd, that because we want the sample at the center of matrix. + * + ***************************************************************************************************/ +int execute_texture(CELL **data, struct dimensions *dim, + struct menu *measure_menu, int *measure_idx, + struct output_setting *out_set) + +{ + int size = dim->size; + int dist = dim->dist; + int nrows = dim->nrows; + int ncols = dim->ncols; + int n_outputs = dim->n_outputs; + int n_measures = dim->n_measures; + int *outfd = out_set->outfd; + RASTER_MAP_TYPE out_data_type = out_set->out_data_type; + struct Flag *flag_null = out_set->flag_null; + struct Flag *flag_ind = out_set->flag_ind; + + int offset = size / 2; + int i, j, row, col, first_row, first_col, last_row, last_col; + int have_px, have_py, have_pxpys, have_pxpyd; + FCELL **fbuf; + FCELL measure; /* Containing measure done */ + struct matvec *mv; + + fbuf = G_malloc(n_outputs * sizeof(FCELL *)); + for (i = 0; i < n_outputs; i++) + fbuf[i] = Rast_allocate_buf(out_data_type); + + mv = G_malloc(sizeof(struct matvec)); + alloc_vars(size, mv); + + /* variables needed */ + if (measure_menu[2].useme || measure_menu[11].useme || + measure_menu[12].useme) + have_px = 1; + else + have_px = 0; + if (measure_menu[11].useme || measure_menu[12].useme) + have_py = 1; + else + have_py = 0; + if (measure_menu[5].useme || measure_menu[6].useme || measure_menu[7].useme) + have_pxpys = 1; + else + have_pxpys = 0; + if (measure_menu[9].useme || measure_menu[10].useme) + have_pxpyd = 1; + else + have_pxpyd = 0; + + if (!flag_null->answer) { + first_row = first_col = offset; + last_row = nrows - offset; + last_col = ncols - offset; + } + else { + /* no cropping at window margins */ + first_row = first_col = 0; + last_row = nrows; + last_col = ncols; + } + + 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(n_("Calculating %d texture measure", + "Calculating %d texture measures", n_measures), + n_measures); + else + G_message(_("Calculating %s..."), measure_menu[measure_idx[0]].desc); + + 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(mv, data, row, col, size, offset, dist, + flag_null->answer)) { + for (i = 0; i < n_outputs; i++) + Rast_set_f_null_value(&(fbuf[i][col]), 1); + continue; + } + /* for all angles (0, 45, 90, 135) */ + for (i = 0; i < 4; i++) { + set_angle_vars(mv, i, have_px, have_py, have_pxpys, have_pxpyd); + /* for all requested textural measures */ + for (j = 0; j < n_measures; j++) { + + measure = + (FCELL)h_measure(measure_menu[measure_idx[j]].idx, mv); + + 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; + } + } + } + } + + 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); + + for (i = 0; i < n_outputs; i++) + G_free(fbuf[i]); + G_free(fbuf); + dealloc_vars(mv); + G_free(mv); + + return 0; +} diff --git a/raster/r.texture/execute.h b/raster/r.texture/execute.h new file mode 100644 index 00000000000..a8f07594fa2 --- /dev/null +++ b/raster/r.texture/execute.h @@ -0,0 +1,23 @@ +#include +#include +#include "h_measure.h" + +typedef struct dimensions { + int size; /* size of moving window */ + int dist; /* value of distance */ + int nrows; + int ncols; + int n_outputs; + int n_measures; +} dimensions; + +typedef struct output_setting { + int *outfd; + RASTER_MAP_TYPE out_data_type; + struct Flag *flag_null; + struct Flag *flag_ind; +} output_setting; + +int execute_texture(CELL **data, struct dimensions *dim, + struct menu *measure_menu, int *measure_idx, + struct output_setting *out_set); diff --git a/raster/r.texture/h_measure.c b/raster/r.texture/h_measure.c index e7d65bcacc9..d81104b27b3 100644 --- a/raster/r.texture/h_measure.c +++ b/raster/r.texture/h_measure.c @@ -28,380 +28,80 @@ #include #include -#define BL "Direction " -#define F1 "Angular Second Moment " -#define F2 "Contrast " -#define F3 "Correlation " -#define F4 "Variance " -#define F5 "Inverse Diff Moment " -#define F6 "Sum Average " -#define F7 "Sum Variance " -#define F8 "Sum Entropy " -#define F9 "Entropy " -#define F10 "Difference Variance " -#define F11 "Difference Entropy " -#define F12 "Measure of Correlation-1 " -#define F13 "Measure of Correlation-2 " - -#define PGM_MAXMAXVAL 255 -#define MAX_MATRIX_SIZE 512 - -float **matrix(int nr, int nc); -float *vector(int n); - -float f1_asm(void); -float f2_contrast(void); -float f3_corr(void); -float f4_var(void); -float f5_idm(void); -float f6_savg(void); -float f7_svar(void); -float f8_sentropy(void); -float f9_entropy(void); -float f10_dvar(void); -float f11_dentropy(void); -float f12_icorr(void); -float f13_icorr(void); - -static float **P_matrix = NULL; -static float **P_matrix0 = NULL; -static float **P_matrix45 = NULL; -static float **P_matrix90 = NULL; -static float **P_matrix135 = NULL; - -int tone[PGM_MAXMAXVAL + 1]; -static int Ng = 0; -static float *px, *py; -static float Pxpys[2 * PGM_MAXMAXVAL + 2]; -static float Pxpyd[2 * PGM_MAXMAXVAL + 2]; - -void alloc_vars(int size) -{ - int msize2; - - /* Allocate memory for gray-tone spatial dependence matrix */ - P_matrix0 = matrix(MAX_MATRIX_SIZE + 1, MAX_MATRIX_SIZE + 1); - P_matrix45 = matrix(MAX_MATRIX_SIZE + 1, MAX_MATRIX_SIZE + 1); - P_matrix90 = matrix(MAX_MATRIX_SIZE + 1, MAX_MATRIX_SIZE + 1); - P_matrix135 = matrix(MAX_MATRIX_SIZE + 1, MAX_MATRIX_SIZE + 1); - - if (size * size < 256) - msize2 = size * size; - else - msize2 = 256; - - px = vector(msize2 + 1); - py = vector(msize2 + 1); -} - -static int bsearch_gray(int *array, int n, int val) -{ - int lo, hi, mid; - - lo = 0; - hi = n - 1; - - while (lo <= hi) { - mid = (lo + hi) >> 1; - - if (array[mid] == val) - return mid; - - if (array[mid] > val) - hi = mid - 1; - else - lo = mid + 1; - } - - return -1; -} - -int set_vars(int **grays, int curr_row, int curr_col, int size, int offset, - int t_d, int with_nulls) -{ - int R0, R45, R90, R135, x, y; - int row, col, row2, col2, rows, cols; - int rowmin, rowmax, colmin, colmax, wrows, wcols, rowd, cold; - int itone, jtone; - int cnt; - - rows = cols = size; - wrows = Rast_window_rows(); - wcols = Rast_window_cols(); - - /* Determine the number of different gray scales (not maxval) */ - for (row = 0; row <= PGM_MAXMAXVAL; row++) - tone[row] = -1; - cnt = 0; - rowmin = curr_row - offset; - if (rowmin < 0) - rowmin = 0; - rowmax = curr_row + offset; - if (rowmax > wrows - 1) - rowmax = wrows - 1; - colmin = curr_col - offset; - if (colmin < 0) - colmin = 0; - colmax = curr_col + offset; - if (colmax > wcols - 1) - colmax = wcols - 1; - for (row = rowmin; row <= rowmax; row++) { - for (col = colmin; col <= colmax; col++) { - if (grays[row][col] < 0) { /* No data pixel found */ - continue; - } - if (grays[row][col] > PGM_MAXMAXVAL) - G_fatal_error(_("Too many categories (found: %i, max: %i). " - "Try to rescale or reclassify the map"), - grays[row][col], PGM_MAXMAXVAL); - tone[grays[row][col]] = grays[row][col]; - cnt++; - } - } - /* what is the minimum number of pixels - * to get reasonable texture measurements ? - * at the very least, any of R0, R45, R90, R135 must be > 1 */ - if (cnt < size * size / 4 || (!with_nulls && cnt < size * size)) - return 0; - - /* Collapse array, taking out all zero values */ - Ng = 0; - for (row = 0; row <= PGM_MAXMAXVAL; row++) { - if (tone[row] != -1) - tone[Ng++] = tone[row]; - } - - /* Now array contains only the gray levels present (in ascending order) */ - - for (row = 0; row < Ng; row++) - for (col = 0; col < Ng; col++) { - P_matrix0[row][col] = P_matrix45[row][col] = 0; - P_matrix90[row][col] = P_matrix135[row][col] = 0; - } - - /* Find normalizing constants */ - /* not correct in case of NULL cells: */ - /* - R0 = 2 * rows * (cols - t_d); - R45 = 2 * (rows - t_d) * (cols - t_d); - R90 = 2 * (rows - t_d) * cols; - R135 = R45; - */ - - /* count actual cooccurrences for each angle */ - R0 = R45 = R90 = R135 = 0; - - /* Find gray-tone spatial dependence matrix */ - for (row = 0; row < rows; row++) { - row2 = curr_row - offset + row; - if (row2 < 0 || row2 >= wrows) - continue; - for (col = 0; col < cols; col++) { - col2 = curr_col - offset + col; - if (col2 < 0 || col2 >= wcols) - continue; - if (grays[row2][col2] < 0) - continue; - x = bsearch_gray(tone, Ng, grays[row2][col2]); - rowd = row2; - cold = col2 + t_d; - if (col + t_d < cols && cold < wcols && grays[rowd][cold] >= 0) { - y = bsearch_gray(tone, Ng, grays[rowd][cold]); - P_matrix0[x][y]++; - P_matrix0[y][x]++; - R0 += 2; - } - rowd = row2 + t_d; - cold = col2; - if (row + t_d < rows && rowd < wrows && grays[rowd][cold] >= 0) { - y = bsearch_gray(tone, Ng, grays[rowd][cold]); - P_matrix90[x][y]++; - P_matrix90[y][x]++; - R90 += 2; - } - rowd = row2 + t_d; - cold = col2 - t_d; - if (row + t_d < rows && rowd < wrows && col - t_d >= 0 && - cold >= 0 && grays[rowd][cold] >= 0) { - y = bsearch_gray(tone, Ng, grays[rowd][cold]); - P_matrix45[x][y]++; - P_matrix45[y][x]++; - R45 += 2; - } - rowd = row2 + t_d; - cold = col2 + t_d; - if (row + t_d < rows && rowd < wrows && col + t_d < cols && - cold < wcols && grays[rowd][cold] >= 0) { - y = bsearch_gray(tone, Ng, grays[rowd][cold]); - P_matrix135[x][y]++; - P_matrix135[y][x]++; - R135 += 2; - } - } - } - /* Gray-tone spatial dependence matrices are complete */ - - /* Normalize gray-tone spatial dependence matrix */ - for (itone = 0; itone < Ng; itone++) { - for (jtone = 0; jtone < Ng; jtone++) { - P_matrix0[itone][jtone] /= R0; - P_matrix45[itone][jtone] /= R45; - P_matrix90[itone][jtone] /= R90; - P_matrix135[itone][jtone] /= R135; - } - } - - return 1; -} - -int set_angle_vars(int angle, int have_px, int have_py, int have_pxpys, - int have_pxpyd) -{ - int i, j; - float **P; +#include "h_measure.h" - switch (angle) { - case 0: - P_matrix = P_matrix0; - break; - case 1: - P_matrix = P_matrix45; - break; - case 2: - P_matrix = P_matrix90; - break; - case 3: - P_matrix = P_matrix135; - break; - } - - P = P_matrix; - - /* - * px[i] is the (i-1)th entry in the marginal probability matrix obtained - * by summing the rows of p[i][j] - */ - /* Pxpy sum and difference */ - - /* reset variabless */ - if (have_px || have_py || have_pxpys || have_pxpyd) { - for (i = 0; i < Ng; i++) { - if (have_px || have_py) { - px[i] = py[i] = 0; - } - if (have_pxpys || have_pxpyd) { - Pxpys[i] = Pxpyd[i] = 0; - } - } - if (have_pxpys) { - for (j = Ng; j < 2 * Ng; j++) { - Pxpys[j] = 0; - } - } - } - - if (have_pxpys || have_pxpyd || have_px || have_py) { - for (i = 0; i < Ng; i++) { - for (j = 0; j < Ng; j++) { - if (have_px || have_py) { - px[i] += P[i][j]; - py[j] += P[i][j]; - } - if (have_pxpys) { - Pxpys[i + j] += P[i][j]; - } - if (have_pxpyd) { - Pxpyd[abs(i - j)] += P[i][j]; - } - } - } - } - - return 1; -} - -float h_measure(int t_m) +float h_measure(int t_m, struct matvec *mv) { switch (t_m) { /* Angular Second Moment */ case 1: - return (f1_asm()); + return (f1_asm(mv)); break; /* Contrast */ case 2: - return (f2_contrast()); + return (f2_contrast(mv)); break; /* Correlation */ case 3: - return (f3_corr()); + return (f3_corr(mv)); break; /* Variance */ case 4: - return (f4_var()); + return (f4_var(mv)); break; /* Inverse Diff Moment */ case 5: - return (f5_idm()); + return (f5_idm(mv)); break; /* Sum Average */ case 6: - return (f6_savg()); + return (f6_savg(mv)); break; /* Sum Variance */ case 7: - return (f7_svar()); + return (f7_svar(mv)); break; /* Sum Entropy */ case 8: - return (f8_sentropy()); + return (f8_sentropy(mv)); break; /* Entropy */ case 9: - return (f9_entropy()); + return (f9_entropy(mv)); break; /* Difference Variance */ case 10: - return (f10_dvar()); + return (f10_dvar(mv)); break; /* Difference Entropy */ case 11: - return (f11_dentropy()); + return (f11_dentropy(mv)); break; /* Measure of Correlation-1 */ case 12: - return (f12_icorr()); + return (f12_icorr(mv)); break; /* Measure of Correlation-2 */ case 13: - return (f13_icorr()); + return (f13_icorr(mv)); break; } return 0; } -void MatrixDealloc(float **A, int N) -{ - /*A is NxN */ - int i; - - for (i = 0; i < N; i++) - G_free(A[i]); - G_free(A); -} - /* Angular Second Moment */ /* * The angular second-moment feature (ASM) f1 is a measure of homogeneity @@ -409,11 +109,12 @@ void MatrixDealloc(float **A, int N) * gray-tone transitions. Hence the P matrix for such an image will have * fewer entries of large magnitude. */ -float f1_asm(void) +float f1_asm(struct matvec *mv) { int i, j; float sum = 0; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; /* for (i = 0; i < Ng; i++) @@ -436,11 +137,13 @@ float f1_asm(void) * measure of the contrast or the amount of local variations present in an * image. */ -float f2_contrast(void) +float f2_contrast(struct matvec *mv) { int i, j /*, n */; float /* sum, */ bigsum = 0; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; + int *tone = mv->tone; /* the three-loop version does not work * when gray tones that do not occur in the current window @@ -475,12 +178,15 @@ float f2_contrast(void) * This correlation feature is a measure of gray-tone linear-dependencies * in the image. */ -float f3_corr(void) +float f3_corr(struct matvec *mv) { int i, j; float sum_sqr = 0, tmp = 0; float mean = 0, stddev; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; + int *tone = mv->tone; + float *px = mv->px; /* Now calculate the means and standard deviations of px and py */ @@ -505,11 +211,13 @@ float f3_corr(void) } /* Sum of Squares: Variance */ -float f4_var(void) +float f4_var(struct matvec *mv) { int i, j; float mean = 0, var = 0; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; + int *tone = mv->tone; /*- Corrected by James Darrell McCauley, 16 Aug 1991 * calculates the mean intensity level instead of the mean of @@ -527,11 +235,13 @@ float f4_var(void) } /* Inverse Difference Moment */ -float f5_idm(void) +float f5_idm(struct matvec *mv) { int i, j; float idm = 0; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; + int *tone = mv->tone; /* for (i = 0; i < Ng; i++) @@ -550,11 +260,13 @@ float f5_idm(void) } /* Sum Average */ -float f6_savg(void) +float f6_savg(struct matvec *mv) { int i, j, k; float savg = 0; - float *P = Pxpys; + float *P = mv->Pxpys; + int Ng = mv->Ng; + int *tone = mv->tone; /* for (i = 0; i < 2 * Ng - 1; i++) @@ -572,12 +284,14 @@ float f6_savg(void) } /* Sum Variance */ -float f7_svar(void) +float f7_svar(struct matvec *mv) { int i, j, k; float var = 0; - float *P = Pxpys; - float savg = f6_savg(); + float *P = mv->Pxpys; + int Ng = mv->Ng; + int *tone = mv->tone; + float savg = f6_savg(mv); float tmp; /* @@ -597,11 +311,12 @@ float f7_svar(void) } /* Sum Entropy */ -float f8_sentropy(void) +float f8_sentropy(struct matvec *mv) { int i; float sentr = 0; - float *P = Pxpys; + float *P = mv->Pxpys; + int Ng = mv->Ng; for (i = 0; i < 2 * Ng - 1; i++) { if (P[i] > 0) @@ -612,11 +327,12 @@ float f8_sentropy(void) } /* Entropy */ -float f9_entropy(void) +float f9_entropy(struct matvec *mv) { int i, j; float entropy = 0; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; /* for (i = 0; i < Ng; i++) { @@ -640,11 +356,13 @@ float f9_entropy(void) } /* Difference Variance */ -float f10_dvar(void) +float f10_dvar(struct matvec *mv) { int i, tmp; float sum = 0, sum_sqr = 0, var = 0; - float *P = Pxpyd; + float *P = mv->Pxpyd; + int Ng = mv->Ng; + int *tone = mv->tone; /* Now calculate the variance of Pxpy (Px-y) */ for (i = 0; i < Ng; i++) { @@ -661,11 +379,12 @@ float f10_dvar(void) } /* Difference Entropy */ -float f11_dentropy(void) +float f11_dentropy(struct matvec *mv) { int i; float sum = 0; - float *P = Pxpyd; + float *P = mv->Pxpyd; + int Ng = mv->Ng; for (i = 0; i < Ng; i++) { if (P[i] > 0) @@ -676,11 +395,14 @@ float f11_dentropy(void) } /* Information Measures of Correlation */ -float f12_icorr(void) +float f12_icorr(struct matvec *mv) { int i, j; float hx = 0, hy = 0, hxy = 0, hxy1 = 0; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; + float *px = mv->px; + float *py = mv->py; for (i = 0; i < Ng; i++) { for (j = 0; j < Ng; j++) { @@ -705,11 +427,14 @@ float f12_icorr(void) } /* Information Measures of Correlation */ -float f13_icorr(void) +float f13_icorr(struct matvec *mv) { int i, j; float hxy = 0, hxy2 = 0; - float **P = P_matrix; + float **P = mv->P_matrix; + int Ng = mv->Ng; + float *px = mv->px; + float *py = mv->py; for (i = 0; i < Ng; i++) { for (j = 0; j < Ng; j++) { @@ -723,30 +448,3 @@ float f13_icorr(void) /* fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */ return (sqrt(fabs(1 - exp(-2.0 * (hxy2 - hxy))))); } - -float *vector(int n) -{ - float *v; - - v = (float *)G_malloc(n * sizeof(float)); - if (!v) - G_fatal_error(_("Unable to allocate memory")), exit(EXIT_FAILURE); - return v; -} - -/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */ -float **matrix(int nr, int nc) -{ - int i; - float **m; - - /* allocate pointers to rows */ - m = (float **)G_malloc(nr * sizeof(float *)); - - /* allocate rows */ - for (i = 0; i < nr; i++) { - m[i] = (float *)G_malloc(nc * sizeof(float)); - } - /* return pointer to array of pointers to rows */ - return m; -} diff --git a/raster/r.texture/h_measure.h b/raster/r.texture/h_measure.h index 6f51e0cbdc0..5a759b4ec72 100644 --- a/raster/r.texture/h_measure.h +++ b/raster/r.texture/h_measure.h @@ -1,4 +1,3 @@ -/* h_measure.c */ /**************************************************************************** * * MODULE: r.texture @@ -22,9 +21,28 @@ * *****************************************************************************/ -float h_measure(int); -void alloc_vars(int); -int set_vars(int **grays, int curr_rrow, int curr_col, int size, int offset, - int t_d, int with_nulls); -int set_angle_vars(int angle, int have_px, int have_py, int have_pxpys, - int have_pxpyd); +#include "matvec.h" + +typedef 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 */ +} menu; + +float h_measure(int t_m, struct matvec *mv); + +float f1_asm(struct matvec *mv); +float f2_contrast(struct matvec *mv); +float f3_corr(struct matvec *mv); +float f4_var(struct matvec *mv); +float f5_idm(struct matvec *mv); +float f6_savg(struct matvec *mv); +float f7_svar(struct matvec *mv); +float f8_sentropy(struct matvec *mv); +float f9_entropy(struct matvec *mv); +float f10_dvar(struct matvec *mv); +float f11_dentropy(struct matvec *mv); +float f12_icorr(struct matvec *mv); +float f13_icorr(struct matvec *mv); diff --git a/raster/r.texture/main.c b/raster/r.texture/main.c index 2c3623fa611..c21bcc1cac3 100644 --- a/raster/r.texture/main.c +++ b/raster/r.texture/main.c @@ -27,18 +27,11 @@ #include #include #include -#include "h_measure.h" -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 */ -}; +#include "execute.h" /* modify this table to add new measures */ -static struct menu menu[] = { +static struct menu measure_menu[] = { {"asm", "Angular Second Moment", "_ASM", 0, 1}, {"contrast", "Contrast", "_Contr", 0, 2}, {"corr", "Correlation", "_Corr", 0, 3}, @@ -58,8 +51,8 @@ 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) + for (i = 0; measure_menu[i].name; i++) + if (strcmp(measure_menu[i].name, measure_name) == 0) return i; G_fatal_error(_("Unknown measure <%s>"), measure_name); @@ -72,25 +65,26 @@ int main(int argc, char *argv[]) struct Cell_head cellhd; char *name, *result; char **mapname; - FCELL **fbuf; - int n_measures, n_outputs, *measure_idx, overwrite; - int nrows, ncols; - int row, col, first_row, last_row, first_col, last_col; + int *measure_idx, overwrite; 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 dist, size; /* dist = value of distance, size = s. of moving window */ - int offset; - int have_px, have_py, have_pxpys, have_pxpyd; int infd, *outfd; + RASTER_MAP_TYPE out_data_type; struct GModule *module; - struct Option *opt_input, *opt_output, *opt_size, *opt_dist, *opt_measure; - struct Flag *flag_ind, *flag_all, *flag_null; + struct { + struct Option *input, *output, *size, *dist, *measure, *nproc; + } parm; + struct { + struct Flag *ind, *all, *null; + } flag; struct History history; + + struct dimensions dim; + struct output_setting out_set; char p[1024]; G_gisinit(argv[0]); @@ -105,156 +99,138 @@ int main(int argc, char *argv[]) /* Define the different options */ - opt_input = G_define_standard_option(G_OPT_R_INPUT); + parm.input = G_define_standard_option(G_OPT_R_INPUT); - opt_output = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT); + parm.output = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT); - 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 moving window (odd and >= 3)"); - opt_size->answer = "3"; + parm.nproc = G_define_standard_option(G_OPT_M_NPROCS); + + parm.size = G_define_option(); + parm.size->key = "size"; + parm.size->key_desc = "value"; + parm.size->type = TYPE_INTEGER; + parm.size->required = NO; + parm.size->description = _("The size of moving window (odd and >= 3)"); + parm.size->answer = "3"; /* Textural character is in direct relation of the spatial size of the * texture primitives. */ - 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->label = _("The distance between two samples (>= 1)"); - opt_dist->description = + parm.dist = G_define_option(); + parm.dist->key = "distance"; + parm.dist->key_desc = "value"; + parm.dist->type = TYPE_INTEGER; + parm.dist->required = NO; + parm.dist->label = _("The distance between two samples (>= 1)"); + parm.dist->description = _("The distance must be smaller than the size of the moving window"); - opt_dist->answer = "1"; + parm.dist->answer = "1"; - for (i = 0; menu[i].name; i++) { + for (i = 0; measure_menu[i].name; i++) { if (i) strcat(p, ","); else *p = 0; - strcat(p, menu[i].name); + strcat(p, measure_menu[i].name); } - opt_measure = G_define_option(); - opt_measure->key = "method"; - opt_measure->type = TYPE_STRING; - opt_measure->required = NO; - opt_measure->multiple = YES; - opt_measure->options = p; - opt_measure->description = _("Textural measurement method"); - - flag_ind = G_define_flag(); - flag_ind->key = 's'; - flag_ind->label = _("Separate output for each angle (0, 45, 90, 135)"); - flag_ind->description = + parm.measure = G_define_option(); + parm.measure->key = "method"; + parm.measure->type = TYPE_STRING; + parm.measure->required = NO; + parm.measure->multiple = YES; + parm.measure->options = p; + parm.measure->description = _("Textural measurement method"); + + flag.ind = G_define_flag(); + flag.ind->key = 's'; + flag.ind->label = _("Separate output for each angle (0, 45, 90, 135)"); + flag.ind->description = _("Angles are counterclockwise from east: " "0 is East to West, 45 is North-East to South-West"); - flag_all = G_define_flag(); - flag_all->key = 'a'; - flag_all->description = _("Calculate all textural measurements"); + flag.all = G_define_flag(); + flag.all->key = 'a'; + flag.all->description = _("Calculate all textural measurements"); - flag_null = G_define_flag(); - flag_null->key = 'n'; - flag_null->label = _("Allow NULL cells in a moving window"); - flag_null->description = + flag.null = G_define_flag(); + flag.null->key = 'n'; + flag.null->label = _("Allow NULL cells in a moving window"); + flag.null->description = _("This will also avoid cropping along edges of the current region"); if (G_parser(argc, argv)) exit(EXIT_FAILURE); - name = opt_input->answer; - result = opt_output->answer; - size = atoi(opt_size->answer); - if (size <= 0) + name = parm.input->answer; + result = parm.output->answer; + dim.size = atoi(parm.size->answer); + dim.dist = atoi(parm.dist->answer); + + if (dim.size <= 0) G_fatal_error(_("Size of the moving window must be > 0")); - if (size % 2 != 1) + if (dim.size % 2 != 1) G_fatal_error(_("Size of the moving window must be odd")); - dist = atoi(opt_dist->answer); - if (dist <= 0) + if (dim.dist <= 0) G_fatal_error(_("The distance between two samples must be > 0")); - if (dist >= size) + if (dim.dist >= dim.size) G_fatal_error(_("The distance between two samples must be smaller than " "the size of the moving window")); - n_measures = 0; - if (flag_all->answer) { - for (i = 0; menu[i].name; i++) { - menu[i].useme = 1; + dim.n_measures = 0; + if (flag.all->answer) { + for (i = 0; measure_menu[i].name; i++) { + measure_menu[i].useme = 1; } - n_measures = i; + dim.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]; + for (i = 0; parm.measure->answers[i]; i++) { + if (parm.measure->answers[i]) { + const char *measure_name = parm.measure->answers[i]; int n = find_measure(measure_name); - menu[n].useme = 1; - n_measures++; + measure_menu[n].useme = 1; + dim.n_measures++; } } } - if (!n_measures) + if (!dim.n_measures) G_fatal_error( _("Nothing to compute. Use at least one textural measure.")); - measure_idx = G_malloc(n_measures * sizeof(int)); + measure_idx = G_malloc(dim.n_measures * sizeof(int)); j = 0; - for (i = 0; menu[i].name; i++) { - if (menu[i].useme == 1) { + for (i = 0; measure_menu[i].name; i++) { + if (measure_menu[i].useme == 1) { measure_idx[j] = i; 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[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, ""); Rast_get_cellhd(name, "", &cellhd); out_data_type = FCELL_TYPE; /* Allocate output buffers, use FCELL data_type */ - n_outputs = n_measures; - if (flag_ind->answer) { - n_outputs = n_measures * 4; + dim.n_outputs = dim.n_measures; + if (flag.ind->answer) { + dim.n_outputs = dim.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 = G_malloc(dim.n_outputs * sizeof(char *)); + for (i = 0; i < dim.n_outputs; i++) mapname[i] = G_malloc(GNAME_MAX * sizeof(char)); - fbuf[i] = Rast_allocate_buf(out_data_type); - } overwrite = G_check_overwrite(argc, argv); /* open output maps */ - outfd = G_malloc(n_outputs * sizeof(int)); - for (i = 0; i < n_measures; i++) { - if (flag_ind->answer) { + outfd = G_malloc(dim.n_outputs * sizeof(int)); + for (i = 0; i < dim.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); + measure_menu[measure_idx[i]].suffix, j * 45); if (!G_find_raster(mapname[i * 4 + j], G_mapset()) || overwrite) { outfd[i * 4 + j] = @@ -267,7 +243,8 @@ int main(int argc, char *argv[]) } } else { - sprintf(mapname[i], "%s%s", result, menu[measure_idx[i]].suffix); + sprintf(mapname[i], "%s%s", result, + measure_menu[measure_idx[i]].suffix); if (!G_find_raster(mapname[i], G_mapset()) || overwrite) { outfd[i] = Rast_open_new(mapname[i], out_data_type); } @@ -277,8 +254,8 @@ int main(int argc, char *argv[]) } } } - nrows = Rast_window_rows(); - ncols = Rast_window_cols(); + dim.nrows = Rast_window_rows(); + dim.ncols = Rast_window_cols(); /* Load raster map. */ @@ -286,9 +263,9 @@ int main(int argc, char *argv[]) dcell_row = Rast_allocate_d_buf(); /* Allocate appropriate memory for the structure containing the image */ - data = (int **)G_malloc(nrows * sizeof(int *)); - for (i = 0; i < nrows; i++) { - data[i] = (int *)G_malloc(ncols * sizeof(int)); + data = (int **)G_malloc(dim.nrows * sizeof(int *)); + for (i = 0; i < dim.nrows; i++) { + data[i] = (int *)G_malloc(dim.ncols * sizeof(int)); } /* read input range */ @@ -307,9 +284,9 @@ 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++) { + for (j = 0; j < dim.nrows; j++) { Rast_get_row(infd, dcell_row, j, DCELL_TYPE); - for (i = 0; i < ncols; i++) { + for (i = 0; i < dim.ncols; i++) { if (Rast_is_d_null_value(&(dcell_row[i]))) data[j][i] = -1; else if (inscale) { @@ -324,107 +301,29 @@ int main(int argc, char *argv[]) Rast_close(infd); G_free(dcell_row); - /* Now raster map is loaded to memory. */ - - /* ************************************************************************************************* - * - * Compute of the matrix S.G.L.D. (Spatial Gray-Level Dependence Matrices) - *or co-occurrence matrix. The image is analyzed for piece, every piece is - *naming moving window (s.w.). The s.w. must be square with number of size's - *samples odd, that because we want the sample at the center of matrix. - * - ***************************************************************************************************/ - - offset = size / 2; - - if (!flag_null->answer) { - first_row = first_col = offset; - last_row = nrows - offset; - last_col = ncols - offset; - } - else { - /* no cropping at window margins */ - first_row = first_col = 0; - last_row = nrows; - last_col = ncols; - } - - 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(n_("Calculating %d texture measure", - "Calculating %d texture measures", n_measures), - n_measures); - else - G_message(_("Calculating %s..."), menu[measure_idx[0]].desc); - alloc_vars(size); - 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, - flag_null->answer)) { - for (i = 0; i < n_outputs; i++) - Rast_set_f_null_value(&(fbuf[i][col]), 1); - continue; - } + /* variables needed */ + out_set.outfd = outfd; + out_set.out_data_type = out_data_type; + out_set.flag_null = flag.null; + out_set.flag_ind = flag.ind; - /* for all angles (0, 45, 90, 135) */ - for (i = 0; i < 4; i++) { - set_angle_vars(i, have_px, have_py, have_pxpys, have_pxpyd); - /* for all requested textural measures */ - for (j = 0; j < n_measures; j++) { - - measure = (FCELL)h_measure(menu[measure_idx[j]].idx); - - 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; - } - } - } - } - 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); + execute_texture(data, &dim, measure_menu, measure_idx, &out_set); - for (i = 0; i < n_outputs; i++) { + for (i = 0; i < dim.n_outputs; i++) { Rast_close(outfd[i]); - Rast_short_history(mapname[i], "raster", &history); Rast_command_history(&history); Rast_write_history(mapname[i], &history); - G_free(fbuf[i]); } - G_free(fbuf); + /* Free allocated memory */ + for (i = 0; i < dim.n_outputs; i++) + G_free(mapname[i]); + for (i = 0; i < dim.nrows; i++) + G_free(data[i]); + + G_free(measure_idx); + G_free(mapname); G_free(data); exit(EXIT_SUCCESS); diff --git a/raster/r.texture/matvec.c b/raster/r.texture/matvec.c new file mode 100644 index 00000000000..56595543f35 --- /dev/null +++ b/raster/r.texture/matvec.c @@ -0,0 +1,307 @@ +#include +#include +#include +#include +#include +#include + +#include "matvec.h" + +int bsearch_gray(int *array, int n, int val) +{ + int lo, hi, mid; + + lo = 0; + hi = n - 1; + + while (lo <= hi) { + mid = (lo + hi) >> 1; + + if (array[mid] == val) + return mid; + + if (array[mid] > val) + hi = mid - 1; + else + lo = mid + 1; + } + + return -1; +} + +float *vector(int n) +{ + float *v; + + v = (float *)G_malloc(n * sizeof(float)); + if (!v) + G_fatal_error(_("Unable to allocate memory")), exit(EXIT_FAILURE); + return v; +} + +/* Allocates a float matrix with range [nrl..nrh][ncl..nch] */ +float **matrix(int nr, int nc) +{ + int i; + float **m; + + /* allocate pointers to rows */ + m = (float **)G_malloc(nr * sizeof(float *)); + + /* allocate rows */ + for (i = 0; i < nr; i++) { + m[i] = (float *)G_malloc(nc * sizeof(float)); + } + /* return pointer to array of pointers to rows */ + return m; +} + +void MatrixDealloc(float **A, int N) +{ + /*A is NxN */ + int i; + + for (i = 0; i < N; i++) + G_free(A[i]); + G_free(A); +} + +void alloc_vars(int size, struct matvec *mv) +{ + int msize2 = (size * size < 256) ? size * size : 256; + size = PGM_MAXMAXVAL; + + /* Allocate memory for gray-tone spatial dependence matrix */ + mv->P_matrix0 = matrix(size + 1, size + 1); + mv->P_matrix45 = matrix(size + 1, size + 1); + mv->P_matrix90 = matrix(size + 1, size + 1); + mv->P_matrix135 = matrix(size + 1, size + 1); + + mv->px = vector(msize2 + 1); + mv->py = vector(msize2 + 1); + mv->Pxpys = vector(2 * msize2 + 2); + mv->Pxpyd = vector(2 * msize2 + 2); + + mv->tone = (int *)G_malloc((PGM_MAXMAXVAL + 1) * sizeof(int)); + mv->Ng = 0; +} + +void dealloc_vars(struct matvec *mv) +{ + /* Deallocate memory for gray-tone spatial dependence matrix */ + MatrixDealloc(mv->P_matrix0, PGM_MAXMAXVAL + 1); + MatrixDealloc(mv->P_matrix45, PGM_MAXMAXVAL + 1); + MatrixDealloc(mv->P_matrix90, PGM_MAXMAXVAL + 1); + MatrixDealloc(mv->P_matrix135, PGM_MAXMAXVAL + 1); + + G_free(mv->px); + G_free(mv->py); + G_free(mv->Pxpys); + G_free(mv->Pxpyd); + + G_free(mv->tone); +} + +int set_vars(struct matvec *mv, int **grays, int curr_row, int curr_col, + int size, int offset, int t_d, int with_nulls) +{ + int R0, R45, R90, R135, x, y; + int row, col, row2, col2, rows, cols; + int rowmin, rowmax, colmin, colmax, wrows, wcols, rowd, cold; + int itone, jtone; + int cnt; + + rows = cols = size; + wrows = Rast_window_rows(); + wcols = Rast_window_cols(); + + /* Determine the number of different gray scales (not maxval) */ + for (row = 0; row <= PGM_MAXMAXVAL; row++) + mv->tone[row] = -1; + cnt = 0; + rowmin = curr_row - offset; + if (rowmin < 0) + rowmin = 0; + rowmax = curr_row + offset; + if (rowmax > wrows - 1) + rowmax = wrows - 1; + colmin = curr_col - offset; + if (colmin < 0) + colmin = 0; + colmax = curr_col + offset; + if (colmax > wcols - 1) + colmax = wcols - 1; + for (row = rowmin; row <= rowmax; row++) { + for (col = colmin; col <= colmax; col++) { + if (grays[row][col] < 0) { /* No data pixel found */ + continue; + } + if (grays[row][col] > PGM_MAXMAXVAL) + G_fatal_error(_("Too many categories (found: %i, max: %i). " + "Try to rescale or reclassify the map"), + grays[row][col], PGM_MAXMAXVAL); + mv->tone[grays[row][col]] = grays[row][col]; + cnt++; + } + } + /* what is the minimum number of pixels + * to get reasonable texture measurements ? + * at the very least, any of R0, R45, R90, R135 must be > 1 */ + if (cnt < size * size / 4 || (!with_nulls && cnt < size * size)) + return 0; + + /* Collapse array, taking out all zero values */ + mv->Ng = 0; + for (row = 0; row <= PGM_MAXMAXVAL; row++) { + if (mv->tone[row] != -1) + mv->tone[mv->Ng++] = mv->tone[row]; + } + + /* Now array contains only the gray levels present (in ascending order) */ + + for (row = 0; row < mv->Ng; row++) + for (col = 0; col < mv->Ng; col++) { + mv->P_matrix0[row][col] = mv->P_matrix45[row][col] = 0; + mv->P_matrix90[row][col] = mv->P_matrix135[row][col] = 0; + } + + /* Find normalizing constants */ + /* not correct in case of NULL cells: */ + /* + R0 = 2 * rows * (cols - t_d); + R45 = 2 * (rows - t_d) * (cols - t_d); + R90 = 2 * (rows - t_d) * cols; + R135 = R45; + */ + + /* count actual cooccurrences for each angle */ + R0 = R45 = R90 = R135 = 0; + + /* Find gray-tone spatial dependence matrix */ + for (row = 0; row < rows; row++) { + row2 = curr_row - offset + row; + if (row2 < 0 || row2 >= wrows) + continue; + for (col = 0; col < cols; col++) { + col2 = curr_col - offset + col; + if (col2 < 0 || col2 >= wcols) + continue; + if (grays[row2][col2] < 0) + continue; + x = bsearch_gray(mv->tone, mv->Ng, grays[row2][col2]); + rowd = row2; + cold = col2 + t_d; + if (col + t_d < cols && cold < wcols && grays[rowd][cold] >= 0) { + y = bsearch_gray(mv->tone, mv->Ng, grays[rowd][cold]); + mv->P_matrix0[x][y]++; + mv->P_matrix0[y][x]++; + R0 += 2; + } + rowd = row2 + t_d; + cold = col2; + if (row + t_d < rows && rowd < wrows && grays[rowd][cold] >= 0) { + y = bsearch_gray(mv->tone, mv->Ng, grays[rowd][cold]); + mv->P_matrix90[x][y]++; + mv->P_matrix90[y][x]++; + R90 += 2; + } + rowd = row2 + t_d; + cold = col2 - t_d; + if (row + t_d < rows && rowd < wrows && col - t_d >= 0 && + cold >= 0 && grays[rowd][cold] >= 0) { + y = bsearch_gray(mv->tone, mv->Ng, grays[rowd][cold]); + mv->P_matrix45[x][y]++; + mv->P_matrix45[y][x]++; + R45 += 2; + } + rowd = row2 + t_d; + cold = col2 + t_d; + if (row + t_d < rows && rowd < wrows && col + t_d < cols && + cold < wcols && grays[rowd][cold] >= 0) { + y = bsearch_gray(mv->tone, mv->Ng, grays[rowd][cold]); + mv->P_matrix135[x][y]++; + mv->P_matrix135[y][x]++; + R135 += 2; + } + } + } + /* Gray-tone spatial dependence matrices are complete */ + + /* Normalize gray-tone spatial dependence matrix */ + for (itone = 0; itone < mv->Ng; itone++) { + for (jtone = 0; jtone < mv->Ng; jtone++) { + mv->P_matrix0[itone][jtone] /= R0; + mv->P_matrix45[itone][jtone] /= R45; + mv->P_matrix90[itone][jtone] /= R90; + mv->P_matrix135[itone][jtone] /= R135; + } + } + + return 1; +} + +int set_angle_vars(struct matvec *mv, int angle, int have_px, int have_py, + int have_pxpys, int have_pxpyd) +{ + int i, j; + float **P; + + switch (angle) { + case 0: + mv->P_matrix = mv->P_matrix0; + break; + case 1: + mv->P_matrix = mv->P_matrix45; + break; + case 2: + mv->P_matrix = mv->P_matrix90; + break; + case 3: + mv->P_matrix = mv->P_matrix135; + break; + } + + P = mv->P_matrix; + + /* + * px[i] is the (i-1)th entry in the marginal probability matrix obtained + * by summing the rows of p[i][j] + */ + /* Pxpy sum and difference */ + + /* reset variabless */ + if (have_px || have_py || have_pxpys || have_pxpyd) { + for (i = 0; i < mv->Ng; i++) { + if (have_px || have_py) { + mv->px[i] = mv->py[i] = 0; + } + if (have_pxpys || have_pxpyd) { + mv->Pxpys[i] = mv->Pxpyd[i] = 0; + } + } + if (have_pxpys) { + for (j = mv->Ng; j < 2 * mv->Ng; j++) { + mv->Pxpys[j] = 0; + } + } + } + + if (have_pxpys || have_pxpyd || have_px || have_py) { + for (i = 0; i < mv->Ng; i++) { + for (j = 0; j < mv->Ng; j++) { + if (have_px || have_py) { + mv->px[i] += P[i][j]; + mv->py[j] += P[i][j]; + } + if (have_pxpys) { + mv->Pxpys[i + j] += P[i][j]; + } + if (have_pxpyd) { + mv->Pxpyd[abs(i - j)] += P[i][j]; + } + } + } + } + + return 1; +} diff --git a/raster/r.texture/matvec.h b/raster/r.texture/matvec.h new file mode 100644 index 00000000000..07c3b9376ef --- /dev/null +++ b/raster/r.texture/matvec.h @@ -0,0 +1,50 @@ +/**************************************************************************** + * + * MODULE: r.texture + * AUTHOR(S): Carmine Basco - basco@unisannio.it + * with hints from: + * prof. Giulio Antoniol - antoniol@ieee.org + * prof. Michele Ceccarelli - ceccarelli@unisannio.it + + * PURPOSE: Intended to explain GRASS raster programming. + * Create map raster with textural features. + * + * COPYRIGHT: (C) 2002 by University of Sannio (BN) - Italy + * + * This program is free software under the GNU General Public + * License (>=v2). Read the file COPYING that comes with GRASS + * for details. + * + * Permission to use, copy, modify, and distribute this software and its + * documentation for any purpose and without fee is hereby granted. This + * software is provided "as is" without express or implied warranty. + * + *****************************************************************************/ + +#define PGM_MAXMAXVAL 255 + +typedef struct matvec { + float **P_matrix; + float **P_matrix0; + float **P_matrix45; + float **P_matrix90; + float **P_matrix135; + + float *px, *py; + float *Pxpys, *Pxpyd; + + int *tone; + int Ng; +} matvec; + +int bsearch_gray(int *array, int n, int val); + +float **matrix(int nr, int nc); +float *vector(int n); + +void alloc_vars(int, struct matvec *); +void dealloc_vars(struct matvec *); +int set_vars(struct matvec *, int **grays, int curr_row, int curr_col, int size, + int offset, int t_d, int with_nulls); +int set_angle_vars(struct matvec *mv, int angle, int have_px, int have_py, + int have_pxpys, int have_pxpyd);