[GRASS-SVN] r47698 - grass/trunk/raster/r.texture
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Aug 17 08:45:12 EDT 2011
Author: mmetz
Date: 2011-08-17 05:45:12 -0700 (Wed, 17 Aug 2011)
New Revision: 47698
Modified:
grass/trunk/raster/r.texture/h_measure.c
grass/trunk/raster/r.texture/h_measure.h
grass/trunk/raster/r.texture/main.c
Log:
r.texture: code optimization
Modified: grass/trunk/raster/r.texture/h_measure.c
===================================================================
--- grass/trunk/raster/r.texture/h_measure.c 2011-08-17 12:36:58 UTC (rev 47697)
+++ grass/trunk/raster/r.texture/h_measure.c 2011-08-17 12:45:12 UTC (rev 47698)
@@ -25,6 +25,7 @@
#include <stdlib.h>
#include <math.h>
#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
#define RADIX 2.0
@@ -44,316 +45,303 @@
#define F11 "Difference Entropy "
#define F12 "Measure of Correlation-1 "
#define F13 "Measure of Correlation-2 "
-#define F14 "Max Correlation Coeff "
#define SIGN(x,y) ((y)<0 ? -fabs(x) : fabs(x))
#define SWAP(a,b) {y=(a);(a)=(b);(b)=y;}
#define PGM_MAXMAXVAL 255
+#define MAX_MATRIX_SIZE 512
-void results(), hessenberg(), mkbalanced(), reduction(), simplesrt();
-float f1_asm(), f2_contrast(), f3_corr(), f4_var(), f5_idm(),
-f6_savg(), f7_svar(), f8_sentropy(), f9_entropy(), f10_dvar(),
-f11_dentropy(), f12_icorr(), f13_icorr(), f14_maxcorr(), *vector(),
-**matrix();
+float **matrix(int nr, int nc);
+float *vector(int n);
+float f1_asm(float **P, int Ng);
+float f2_contrast(float **P, int Ng);
+float f3_corr(float **P, int Ng);
+float f4_var(float **P, int Ng);
+float f5_idm(float **P, int Ng);
+float f6_savg(float **P, int Ng);
+float f7_svar(float **P, int Ng, double S);
+float f8_sentropy(float **P, int Ng);
+float f9_entropy(float **P, int Ng);
+float f10_dvar(float **P, int Ng);
+float f11_dentropy(float **P, int Ng);
+float f12_icorr(float **P, int Ng);
+float f13_icorr(float **P, int Ng);
+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 tones = 0;
+static float sentropy = 0.0;
+static float *px, *py;
+static float Pxpys[2 * PGM_MAXMAXVAL + 2];
+static float Pxpyd[2 * PGM_MAXMAXVAL + 2];
-float h_measure(int **grays, int rows, int cols, int t_m, int t_d)
+
+void alloc_vars(int size, int dist)
{
- int R0, R45, R90, angle, d, x, y;
- int row, col;
- int itone, jtone, tones;
- static int memory_all = 0;
- float sentropy[4];
+ int Ng;
- d = t_d;
+ /* 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)
+ Ng = size * size;
+ else
+ Ng = 256;
+
+ px = vector(Ng + 1);
+ py = vector(Ng + 1);
+}
+
+int set_vars(int **grays, int curr_row, int curr_col,
+ int size, int offset, int t_d)
+{
+ int R0, R45, R90, R135, x, y;
+ int row, col, row2, col2, rows, cols;
+ int itone, jtone;
+
+ rows = cols = size;
+
/* Determine the number of different gray scales (not maxval) */
- for (row = PGM_MAXMAXVAL; row >= 0; --row)
+ for (row = 0; row <= PGM_MAXMAXVAL; row++)
tone[row] = -1;
- for (row = rows - 1; row >= 0; --row)
- for (col = 0; col < cols; ++col) {
- if (grays[row][col] < 0) /* No data pixel found */
- G_fatal_error(_("Negative or no data pixel found. "
- "This module is not yet able to process no data holes in a map, "
- "please fill with r.fillnulls or other algorithms"));
+ for (row = curr_row - offset; row <= curr_row + offset; row++) {
+ for (col = curr_col - offset; col <= curr_col + offset; col++) {
+ if (grays[row][col] < 0) { /* No data pixel found */
+ return 0;
+ }
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];
}
- for (row = PGM_MAXMAXVAL, tones = 0; row >= 0; --row)
- if (tone[row] != -1)
- tones++;
+ }
/* Collapse array, taking out all zero values */
- for (row = 0, itone = 0; row <= PGM_MAXMAXVAL; row++)
+ tones = 0;
+ for (row = 0; row <= PGM_MAXMAXVAL; row++) {
if (tone[row] != -1)
- tone[itone++] = tone[row];
+ tone[tones++] = tone[row];
+ }
+
/* Now array contains only the gray levels present (in ascending order) */
- /* Allocate memory for gray-tone spatial dependence matrix */
-
- if (!memory_all) {
- P_matrix0 = matrix(0, 512, 0, 512);
- P_matrix45 = matrix(0, 512, 0, 512);
- P_matrix90 = matrix(0, 512, 0, 512);
- P_matrix135 = matrix(0, 512, 0, 512);
- memory_all = 1;
- }
- for (row = 0; row < tones; ++row)
- for (col = 0; col < tones; ++col) {
+ for (row = 0; row < tones; row++)
+ for (col = 0; col < tones; col++) {
P_matrix0[row][col] = P_matrix45[row][col] = 0;
P_matrix90[row][col] = P_matrix135[row][col] = 0;
}
/* Find gray-tone spatial dependence matrix */
- for (row = 0; row < rows; ++row)
- for (col = 0; col < cols; ++col)
- for (x = 0, angle = 0; angle <= 135; angle += 45) {
- while (tone[x] != grays[row][col])
- x++;
- if (angle == 0 && col + d < cols) {
- y = 0;
- while (tone[y] != grays[row][col + d])
- y++;
- P_matrix0[x][y]++;
- P_matrix0[y][x]++;
- }
- if (angle == 90 && row + d < rows) {
- y = 0;
- while (tone[y] != grays[row + d][col])
- y++;
- P_matrix90[x][y]++;
- P_matrix90[y][x]++;
- }
- if (angle == 45 && row + d < rows && col - d >= 0) {
- y = 0;
- while (tone[y] != grays[row + d][col - d])
- y++;
- P_matrix45[x][y]++;
- P_matrix45[y][x]++;
- }
- if (angle == 135 && row + d < rows && col + d < cols) {
- y = 0;
- while (tone[y] != grays[row + d][col + d])
- y++;
- P_matrix135[x][y]++;
- P_matrix135[y][x]++;
- }
+ for (row = 0; row < rows; row++) {
+ row2 = curr_row - offset + row;
+ for (col = 0; col < cols; col++) {
+ col2 = curr_col - offset + col;
+ x = 0;
+ while (tone[x] != grays[row2][col2])
+ x++;
+ if (col + t_d < cols) {
+ y = 0;
+ while (tone[y] != grays[row2][col2 + t_d])
+ y++;
+ P_matrix0[x][y]++;
+ P_matrix0[y][x]++;
}
+ if (row + t_d < rows) {
+ y = 0;
+ while (tone[y] != grays[row2 + t_d][col2])
+ y++;
+ P_matrix90[x][y]++;
+ P_matrix90[y][x]++;
+ }
+ if (row + t_d < rows && col - t_d >= 0) {
+ y = 0;
+ while (tone[y] != grays[row2 + t_d][col2 - t_d])
+ y++;
+ P_matrix45[x][y]++;
+ P_matrix45[y][x]++;
+ }
+ if (row + t_d < rows && col + t_d < cols) {
+ y = 0;
+ while (tone[y] != grays[row2 + t_d][col2 + t_d])
+ y++;
+ P_matrix135[x][y]++;
+ P_matrix135[y][x]++;
+ }
+ }
+ }
/* Gray-tone spatial dependence matrices are complete */
/* Find normalizing constants */
R0 = 2 * rows * (cols - 1);
R45 = 2 * (rows - 1) * (cols - 1);
R90 = 2 * (rows - 1) * cols;
+ R135 = R45;
/* Normalize gray-tone spatial dependence matrix */
- for (itone = 0; itone < tones; ++itone)
- for (jtone = 0; jtone < tones; ++jtone) {
+ for (itone = 0; itone < tones; itone++) {
+ for (jtone = 0; jtone < tones; jtone++) {
P_matrix0[itone][jtone] /= R0;
P_matrix45[itone][jtone] /= R45;
P_matrix90[itone][jtone] /= R90;
- P_matrix135[itone][jtone] /= R45;
+ P_matrix135[itone][jtone] /= R135;
}
+ }
+ return 1;
+}
+
+int set_angle_vars(int angle, int have_px, int have_py, int have_sentr,
+ int have_pxpys, int have_pxpyd)
+{
+ int i, j, Ng;
+ float **P;
+
+ 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;
+ }
+
+ if (have_sentr)
+ sentropy = f8_sentropy(P_matrix, tones);
+
+ Ng = tones;
+ 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)
+{
switch (t_m) {
+ /* Angular Second Moment */
case 0:
- return (f1_asm(P_matrix0, tones));
+ return (f1_asm(P_matrix, tones));
break;
+
+ /* Contrast */
case 1:
- return (f1_asm(P_matrix45, tones));
+ return (f2_contrast(P_matrix, tones));
break;
+
+ /* Correlation */
case 2:
- return (f1_asm(P_matrix90, tones));
+ return (f3_corr(P_matrix, tones));
break;
+
+ /* Variance */
case 3:
- return (f1_asm(P_matrix135, tones));
+ return (f4_var(P_matrix, tones));
break;
+ /* Inverse Diff Moment */
case 4:
- return (f2_contrast(P_matrix0, tones));
+ return (f5_idm(P_matrix, tones));
break;
+
+ /* Sum Average */
case 5:
- return (f2_contrast(P_matrix45, tones));
+ return (f6_savg(P_matrix, tones));
break;
+
+ /* Sum Entropy */
case 6:
- return (f2_contrast(P_matrix90, tones));
+ return (sentropy);
break;
+
+ /* Sum Variance */
case 7:
- return (f2_contrast(P_matrix135, tones));
+ return (f7_svar(P_matrix, tones, sentropy));
break;
+ /* Entropy */
case 8:
- return (f3_corr(P_matrix0, tones));
+ return (f9_entropy(P_matrix, tones));
break;
+
+ /* Difference Variance */
case 9:
- return (f3_corr(P_matrix45, tones));
+ return (f10_dvar(P_matrix, tones));
break;
+
+ /* Difference Entropy */
case 10:
- return (f3_corr(P_matrix90, tones));
+ return (f11_dentropy(P_matrix, tones));
break;
+
+ /* Measure of Correlation-1 */
case 11:
- return (f3_corr(P_matrix135, tones));
+ return (f12_icorr(P_matrix, tones));
break;
+ /* Measure of Correlation-2 */
case 12:
- return (f4_var(P_matrix0, tones));
+ return (f13_icorr(P_matrix, tones));
break;
- case 13:
- return (f4_var(P_matrix45, tones));
- break;
- case 14:
- return (f4_var(P_matrix90, tones));
- break;
- case 15:
- return (f4_var(P_matrix135, tones));
- break;
-
- case 16:
- return (f5_idm(P_matrix0, tones));
- break;
- case 17:
- return (f5_idm(P_matrix45, tones));
- break;
- case 18:
- return (f5_idm(P_matrix90, tones));
- break;
- case 19:
- return (f5_idm(P_matrix135, tones));
- break;
-
- case 20:
- return (f6_savg(P_matrix0, tones));
- break;
- case 21:
- return (f6_savg(P_matrix45, tones));
- break;
- case 22:
- return (f6_savg(P_matrix90, tones));
- break;
- case 23:
- return (f6_savg(P_matrix135, tones));
- break;
-
- case 24:
- sentropy[0] = f8_sentropy(P_matrix0, tones);
- return (sentropy[0]);
- break;
- case 25:
- sentropy[1] = f8_sentropy(P_matrix45, tones);
- return (sentropy[1]);
- break;
- case 26:
- sentropy[2] = f8_sentropy(P_matrix90, tones);
- return (sentropy[2]);
- break;
- case 27:
- sentropy[3] = f8_sentropy(P_matrix135, tones);
- return (sentropy[3]);
- break;
- case 28:
- return (f7_svar(P_matrix0, tones, sentropy[0]));
- break;
- case 29:
- return (f7_svar(P_matrix45, tones, sentropy[1]));
- break;
- case 30:
- return (f7_svar(P_matrix90, tones, sentropy[2]));
- break;
- case 31:
- return (f7_svar(P_matrix135, tones, sentropy[3]));
- break;
-
- case 32:
- return (f9_entropy(P_matrix0, tones));
- break;
- case 33:
- return (f9_entropy(P_matrix45, tones));
- break;
- case 34:
- return (f9_entropy(P_matrix90, tones));
- break;
- case 35:
- return (f9_entropy(P_matrix135, tones));
- break;
-
- case 36:
- return (f10_dvar(P_matrix0, tones));
- break;
- case 37:
- return (f10_dvar(P_matrix45, tones));
- break;
- case 38:
- return (f10_dvar(P_matrix90, tones));
- break;
- case 39:
- return (f10_dvar(P_matrix135, tones));
- break;
-
- case 40:
- return (f11_dentropy(P_matrix0, tones));
- break;
- case 41:
- return (f11_dentropy(P_matrix45, tones));
- break;
- case 42:
- return (f11_dentropy(P_matrix90, tones));
- break;
- case 43:
- return (f11_dentropy(P_matrix135, tones));
- break;
-
- case 44:
- return (f12_icorr(P_matrix0, tones));
- break;
- case 45:
- return (f12_icorr(P_matrix45, tones));
- break;
- case 46:
- return (f12_icorr(P_matrix90, tones));
- break;
- case 47:
- return (f12_icorr(P_matrix135, tones));
- break;
-
- case 48:
- return (f13_icorr(P_matrix0, tones));
- break;
- case 49:
- return (f13_icorr(P_matrix45, tones));
- break;
- case 50:
- return (f13_icorr(P_matrix90, tones));
- break;
- case 51:
- return (f13_icorr(P_matrix135, tones));
- break;
-
- case 52:
- return (f14_maxcorr(P_matrix0, tones));
- break;
- case 53:
- return (f14_maxcorr(P_matrix45, tones));
- break;
- case 54:
- return (f14_maxcorr(P_matrix90, tones));
- break;
- case 55:
- return (f14_maxcorr(P_matrix135, tones));
- break;
}
- exit(0);
+
+ return 0;
}
void MatrixDealloc(float **A, int N)
@@ -361,78 +349,67 @@
/*A is NxN */
int i;
- for (i = 0; i < N; ++i)
+ for (i = 0; i < N; i++)
G_free(A[i]);
G_free(A);
}
-float f1_asm(float **P, int Ng)
/* Angular Second Moment */
+/*
+ * The angular second-moment feature (ASM) f1 is a measure of homogeneity
+ * of the image. In a homogeneous image, there are very few dominant
+ * gray-tone transitions. Hence the P matrix for such an image will have
+ * fewer entries of large magnitude.
+ */
+float f1_asm(float **P, int Ng)
{
int i, j;
float sum = 0;
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
+ for (i = 0; i < Ng; i++)
+ for (j = 0; j < Ng; j++)
sum += P[i][j] * P[i][j];
return sum;
-
- /*
- * The angular second-moment feature (ASM) f1 is a measure of homogeneity
- * of the image. In a homogeneous image, there are very few dominant
- * gray-tone transitions. Hence the P matrix for such an image will have
- * fewer entries of large magnitude.
- */
}
-
-float f2_contrast(float **P, int Ng)
-
/* Contrast */
+/*
+ * The contrast feature is a difference moment of the P matrix and is a
+ * measure of the contrast or the amount of local variations present in an
+ * image.
+ */
+float f2_contrast(float **P, int Ng)
{
int i, j, n;
- float sum = 0, bigsum = 0;
+ float sum, bigsum = 0;
- for (n = 0; n < Ng; ++n) {
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- if ((i - j) == n || (j - i) == n)
+ for (n = 0; n < Ng; n++) {
+ sum = 0;
+ for (i = 0; i < Ng; i++) {
+ for (j = 0; j < Ng; j++) {
+ if ((i - j) == n || (j - i) == n) {
sum += P[i][j];
+ }
+ }
+ }
bigsum += n * n * sum;
-
- sum = 0;
}
return bigsum;
-
- /*
- * The contrast feature is a difference moment of the P matrix and is a
- * measure of the contrast or the amount of local variations present in an
- * image.
- */
}
-float f3_corr(float **P, int Ng)
-
/* Correlation */
+/*
+ * This correlation feature is a measure of gray-tone linear-dependencies
+ * in the image.
+ */
+float f3_corr(float **P, int Ng)
{
int i, j;
- float sum_sqrx = 0, sum_sqry = 0, tmp, *px;
+ float sum_sqrx = 0, sum_sqry = 0, tmp = 0;
float meanx = 0, meany = 0, stddevx, stddevy;
- px = vector(0, Ng);
- for (i = 0; i < Ng; ++i)
- px[i] = 0;
- /*
- * px[i] is the (i-1)th entry in the marginal probability matrix obtained
- * by summing the rows of p[i][j]
- */
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- px[i] += P[i][j];
-
-
/* Now calculate the means and standard deviations of px and py */
/*- fix supplied by J. Michael Christensen, 21 Jun 1991 */
@@ -440,33 +417,23 @@
/*- further modified by James Darrell McCauley, 16 Aug 1991
* after realizing that meanx=meany and stddevx=stddevy
*/
- for (i = 0; i < Ng; ++i) {
+ for (i = 0; i < Ng; i++) {
meanx += px[i] * i;
sum_sqrx += px[i] * i * i;
+
+ for (j = 0; j < Ng; j++)
+ tmp += i * j * P[i][j];
}
meany = meanx;
sum_sqry = sum_sqrx;
stddevx = sqrt(sum_sqrx - (meanx * meanx));
stddevy = stddevx;
- /* Finally, the correlation ... */
- for (tmp = 0, i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- tmp += i * j * P[i][j];
-
- G_free(px);
-
return (tmp - meanx * meany) / (stddevx * stddevy);
- /*
- * This correlation feature is a measure of gray-tone linear-dependencies
- * in the image.
- */
}
-
-float f4_var(float **P, int Ng)
-
/* Sum of Squares: Variance */
+float f4_var(float **P, int Ng)
{
int i, j;
float mean = 0, var = 0;
@@ -475,131 +442,91 @@
* calculates the mean intensity level instead of the mean of
* cooccurrence matrix elements
*/
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
+ for (i = 0; i < Ng; i++)
+ for (j = 0; j < Ng; j++)
mean += i * P[i][j];
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
+ for (i = 0; i < Ng; i++)
+ for (j = 0; j < Ng; j++)
var += (i + 1 - mean) * (i + 1 - mean) * P[i][j];
return var;
}
-float f5_idm(float **P, int Ng)
-
/* Inverse Difference Moment */
+float f5_idm(float **P, int Ng)
{
int i, j;
float idm = 0;
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
+ for (i = 0; i < Ng; i++)
+ for (j = 0; j < Ng; j++)
idm += P[i][j] / (1 + (i - j) * (i - j));
return idm;
}
-float Pxpy[2 * PGM_MAXMAXVAL];
-
-float f6_savg(float **P, int Ng)
-
/* Sum Average */
+float f6_savg(float **P, int Ng)
{
- int i, j;
- extern float Pxpy[2 * PGM_MAXMAXVAL];
+ int i;
float savg = 0;
- for (i = 0; i <= 2 * Ng; ++i)
- Pxpy[i] = 0;
+ for (i = 0; i < 2 * Ng - 1; i++)
+ savg += (i + 2) * Pxpys[i];
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- Pxpy[i + j + 2] += P[i][j];
- for (i = 2; i <= 2 * Ng; ++i)
- savg += i * Pxpy[i];
-
return savg;
}
-
-float f7_svar(float **P, int Ng, double S)
-
/* Sum Variance */
+float f7_svar(float **P, int Ng, double S)
{
- int i, j;
- extern float Pxpy[2 * PGM_MAXMAXVAL];
+ int i;
float var = 0;
- for (i = 0; i <= 2 * Ng; ++i)
- Pxpy[i] = 0;
+ for (i = 0; i < 2 * Ng - 1; i++)
+ var += (i + 2 - S) * (i + 2 - S) * Pxpys[i];
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- Pxpy[i + j + 2] += P[i][j];
-
- for (i = 2; i <= 2 * Ng; ++i)
- var += (i - S) * (i - S) * Pxpy[i];
-
return var;
}
-float f8_sentropy(float **P, int Ng)
-
/* Sum Entropy */
+float f8_sentropy(float **P, int Ng)
{
- int i, j;
- extern float Pxpy[2 * PGM_MAXMAXVAL];
- float sentropy = 0;
+ int i;
+ float sentr = 0;
- for (i = 0; i <= 2 * Ng; ++i)
- Pxpy[i] = 0;
+ for (i = 0; i < 2 * Ng - 1; i++)
+ sentr -= Pxpys[i] * log10(Pxpys[i] + EPSILON);
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- Pxpy[i + j + 2] += P[i][j];
-
- for (i = 2; i <= 2 * Ng; ++i)
- sentropy -= Pxpy[i] * log10(Pxpy[i] + EPSILON);
-
- return sentropy;
+ return sentr;
}
-
-float f9_entropy(float **P, int Ng)
-
/* Entropy */
+float f9_entropy(float **P, int Ng)
{
int i, j;
float entropy = 0;
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
+ for (i = 0; i < Ng; i++) {
+ for (j = 0; j < Ng; j++) {
entropy += P[i][j] * log10(P[i][j] + EPSILON);
+ }
+ }
return -entropy;
}
-
-float f10_dvar(float **P, int Ng)
-
/* Difference Variance */
+float f10_dvar(float **P, int Ng)
{
- int i, j, tmp;
- extern float Pxpy[2 * PGM_MAXMAXVAL];
+ int i, tmp;
float sum = 0, sum_sqr = 0, var = 0;
- for (i = 0; i <= 2 * Ng; ++i)
- Pxpy[i] = 0;
-
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- Pxpy[abs(i - j)] += P[i][j];
-
/* Now calculate the variance of Pxpy (Px-y) */
- for (i = 0; i < Ng; ++i) {
- sum += Pxpy[i];
- sum_sqr += Pxpy[i] * Pxpy[i];
+ for (i = 0; i < Ng; i++) {
+ sum += Pxpyd[i];
+ sum_sqr += Pxpyd[i] * Pxpyd[i];
}
tmp = Ng * Ng;
var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
@@ -607,416 +534,80 @@
return var;
}
-float f11_dentropy(float **P, int Ng)
-
/* Difference Entropy */
+float f11_dentropy(float **P, int Ng)
{
- int i, j;
- extern float Pxpy[2 * PGM_MAXMAXVAL];
+ int i;
float sum = 0;
- for (i = 0; i <= 2 * Ng; ++i)
- Pxpy[i] = 0;
+ for (i = 0; i < Ng; i++)
+ sum += Pxpyd[i] * log10(Pxpyd[i] + EPSILON);
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j)
- Pxpy[abs(i - j)] += P[i][j];
-
- for (i = 0; i < Ng; ++i)
- sum += Pxpy[i] * log10(Pxpy[i] + EPSILON);
-
return -sum;
}
-float f12_icorr(float **P, int Ng)
-
/* Information Measures of Correlation */
+float f12_icorr(float **P, int Ng)
{
int i, j;
- float *px, *py;
- float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
+ float hx = 0, hy = 0, hxy = 0, hxy1 = 0;
- px = vector(0, Ng);
- py = vector(0, Ng);
-
- /*
- * px[i] is the (i-1)th entry in the marginal probability matrix obtained
- * by summing the rows of p[i][j]
- */
- for (i = 0; i < Ng; ++i) {
- for (j = 0; j < Ng; ++j) {
- px[i] += P[i][j];
- py[j] += P[i][j];
- }
- }
-
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j) {
+ for (i = 0; i < Ng; i++)
+ for (j = 0; j < Ng; j++) {
hxy1 -= P[i][j] * log10(px[i] * py[j] + EPSILON);
- hxy2 -= px[i] * py[j] * log10(px[i] * py[j] + EPSILON);
hxy -= P[i][j] * log10(P[i][j] + EPSILON);
}
/* Calculate entropies of px and py - is this right? */
- for (i = 0; i < Ng; ++i) {
+ for (i = 0; i < Ng; i++) {
hx -= px[i] * log10(px[i] + EPSILON);
hy -= py[i] * log10(py[i] + EPSILON);
}
- G_free(px);
- G_free(py);
+
/* fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */
return ((hxy - hxy1) / (hx > hy ? hx : hy));
}
-float f13_icorr(float **P, int Ng)
-
/* Information Measures of Correlation */
+float f13_icorr(float **P, int Ng)
{
int i, j;
- float *px, *py;
- float hx = 0, hy = 0, hxy = 0, hxy1 = 0, hxy2 = 0;
+ float hxy = 0, hxy2 = 0;
- px = vector(0, Ng);
- py = vector(0, Ng);
-
- /*
- * px[i] is the (i-1)th entry in the marginal probability matrix obtained
- * by summing the rows of p[i][j]
- */
- for (i = 0; i < Ng; ++i) {
- for (j = 0; j < Ng; ++j) {
- px[i] += P[i][j];
- py[j] += P[i][j];
- }
- }
-
- for (i = 0; i < Ng; ++i)
- for (j = 0; j < Ng; ++j) {
- hxy1 -= P[i][j] * log10(px[i] * py[j] + EPSILON);
+ for (i = 0; i < Ng; i++) {
+ for (j = 0; j < Ng; j++) {
hxy2 -= px[i] * py[j] * log10(px[i] * py[j] + EPSILON);
hxy -= P[i][j] * log10(P[i][j] + EPSILON);
}
-
- /* Calculate entropies of px and py */
- for (i = 0; i < Ng; ++i) {
- hx -= px[i] * log10(px[i] + EPSILON);
- hy -= py[i] * log10(py[i] + EPSILON);
}
- G_free(px);
- G_free(py);
+
/* fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */
return (sqrt(abs(1 - exp(-2.0 * (hxy2 - hxy)))));
}
-float f14_maxcorr(float **P, int Ng)
-
-/* Returns the Maximal Correlation Coefficient */
+float *vector(int n)
{
- int i, j, k;
- float *px, *py, **Q;
- float *x, *iy, tmp;
-
- px = vector(0, Ng);
- py = vector(0, Ng);
- Q = matrix(1, Ng + 1, 1, Ng + 1);
- x = vector(1, Ng);
- iy = vector(1, Ng);
-
- /*
- * px[i] is the (i-1)th entry in the marginal probability matrix obtained
- * by summing the rows of p[i][j]
- */
- for (i = 0; i < Ng; ++i) {
- for (j = 0; j < Ng; ++j) {
- px[i] += P[i][j];
- py[j] += P[i][j];
- }
- }
-
- /* Find the Q matrix */
- for (i = 0; i < Ng; ++i) {
- for (j = 0; j < Ng; ++j) {
- Q[i + 1][j + 1] = 0;
- for (k = 0; k < Ng; ++k)
- Q[i + 1][j + 1] += P[i][k] * P[j][k] / px[i] / py[k];
- }
- }
-
- /* Balance the matrix */
- mkbalanced(Q, Ng);
- /* Reduction to Hessenberg Form */
- reduction(Q, Ng);
- /* Finding eigenvalue for nonsymetric matrix using QR algorithm */
- hessenberg(Q, Ng, x, iy);
- /* simplesrt(Ng,x); */
- /* Returns the sqrt of the second largest eigenvalue of Q */
- for (i = 2, tmp = x[1]; i <= Ng; ++i)
- tmp = (tmp > x[i]) ? tmp : x[i];
-
- MatrixDealloc(Q, Ng);
- G_free(px);
- G_free(py);
- G_free(x);
- G_free(iy);
- return sqrt(x[Ng - 1]);
-}
-
-float *vector(int nl, int nh)
-{
float *v;
- v = (float *)G_malloc((unsigned)(nh - nl + 1) * sizeof(float));
+ v = (float *)G_malloc(n * sizeof(float));
if (!v)
G_fatal_error(_("Unable to allocate memory")), exit(EXIT_FAILURE);
return v;
}
-
-float **matrix(int nrl, int nrh, int ncl, int nch)
-
/* 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((unsigned)(nrh - nrl + 1) * sizeof(float *));
+ m = (float **)G_malloc(nr * sizeof(float *));
/* allocate rows */
- for (i = 0; i < (nrh - nrl + 1); i++) {
- m[i] = (float *)G_malloc((unsigned)(nch - ncl + 1) * sizeof(float));
+ 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 simplesrt(int n, float arr[])
-{
- int i, j;
- float a;
-
- for (j = 2; j <= n; j++) {
- a = arr[j];
- i = j - 1;
- while (i > 0 && arr[i] > a) {
- arr[i + 1] = arr[i];
- i--;
- }
- arr[i + 1] = a;
- }
-}
-
-void mkbalanced(float **a, int n)
-{
- int last, j, i;
- float s, r, g, f, c, sqrdx;
-
- sqrdx = RADIX * RADIX;
- last = 0;
- while (last == 0) {
- last = 1;
- for (i = 1; i <= n; i++) {
- r = c = 0.0;
- for (j = 1; j <= n; j++)
- if (j != i) {
- c += fabs(a[j][i]);
- r += fabs(a[i][j]);
- }
- if (c && r) {
- g = r / RADIX;
- f = 1.0;
- s = c + r;
- while (c < g) {
- f *= RADIX;
- c *= sqrdx;
- }
- g = r * RADIX;
- while (c > g) {
- f /= RADIX;
- c /= sqrdx;
- }
- if ((c + r) / f < 0.95 * s) {
- last = 0;
- g = 1.0 / f;
- for (j = 1; j <= n; j++)
- a[i][j] *= g;
- for (j = 1; j <= n; j++)
- a[j][i] *= f;
- }
- }
- }
- }
-}
-
-
-void reduction(float **a, int n)
-{
- int m, j, i;
- float y, x;
-
- for (m = 2; m < n; m++) {
- x = 0.0;
- i = m;
- for (j = m; j <= n; j++) {
- if (fabs(a[j][m - 1]) > fabs(x)) {
- x = a[j][m - 1];
- i = j;
- }
- }
- if (i != m) {
- for (j = m - 1; j <= n; j++)
- SWAP(a[i][j], a[m][j])
- for (j = 1; j <= n; j++)
- SWAP(a[j][i], a[j][m])
- }
- if (x) {
- for (i = m + 1; i <= n; i++) {
- if ((y = a[i][m - 1])) {
- y /= x;
- a[i][m - 1] = y;
- for (j = m; j <= n; j++)
- a[i][j] -= y * a[m][j];
- for (j = 1; j <= n; j++)
- a[j][m] += y * a[j][i];
- }
- }
- }
- }
-}
-
-void hessenberg(float **a, int n, float wr[], float wi[])
-{
- int nn, m, l, k, j, its, i, mmin;
- float z, y, x, w, v, u, t, s, r, q, p, anorm;
-
- anorm = fabs(a[1][1]);
- for (i = 2; i <= n; i++)
- for (j = (i - 1); j <= n; j++)
- anorm += fabs(a[i][j]);
- nn = n;
- t = 0.0;
- while (nn >= 1) {
- its = 0;
- do {
- for (l = nn; l >= 2; l--) {
- s = fabs(a[l - 1][l - 1]) + fabs(a[l][l]);
- if (s == 0.0)
- s = anorm;
- if ((float)(fabs(a[l][l - 1]) + s) == s)
- break;
- }
- x = a[nn][nn];
- if (l == nn) {
- wr[nn] = x + t;
- wi[nn--] = 0.0;
- }
- else {
- y = a[nn - 1][nn - 1];
- w = a[nn][nn - 1] * a[nn - 1][nn];
- if (l == (nn - 1)) {
- p = 0.5 * (y - x);
- q = p * p + w;
- z = sqrt(fabs(q));
- x += t;
- if (q >= 0.0) {
- z = p + SIGN(z, p);
- wr[nn - 1] = wr[nn] = x + z;
- if (z)
- wr[nn] = x - w / z;
- wi[nn - 1] = wi[nn] = 0.0;
- }
- else {
- wr[nn - 1] = wr[nn] = x + p;
- wi[nn - 1] = -(wi[nn] = z);
- }
- nn -= 2;
- }
- else {
- if (its == 30)
- G_fatal_error(_("Too many iterations to required to find %s - giving up"),
- F14), exit(1);
- if (its == 10 || its == 20) {
- t += x;
- for (i = 1; i <= nn; i++)
- a[i][i] -= x;
- s = fabs(a[nn][nn - 1]) + fabs(a[nn - 1][nn - 2]);
- y = x = 0.75 * s;
- w = -0.4375 * s * s;
- }
- ++its;
- for (m = (nn - 2); m >= l; m--) {
- z = a[m][m];
- r = x - z;
- s = y - z;
- p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
- q = a[m + 1][m + 1] - z - r - s;
- r = a[m + 2][m + 1];
- s = fabs(p) + fabs(q) + fabs(r);
- p /= s;
- q /= s;
- r /= s;
- if (m == l)
- break;
- u = fabs(a[m][m - 1]) * (fabs(q) + fabs(r));
- v = fabs(p) * (fabs(a[m - 1][m - 1]) + fabs(z) +
- fabs(a[m + 1][m + 1]));
- if ((float)(u + v) == v)
- break;
- }
- for (i = m + 2; i <= nn; i++) {
- a[i][i - 2] = 0.0;
- if (i != (m + 2))
- a[i][i - 3] = 0.0;
- }
- for (k = m; k <= nn - 1; k++) {
- if (k != m) {
- p = a[k][k - 1];
- q = a[k + 1][k - 1];
- r = 0.0;
- if (k != (nn - 1))
- r = a[k + 2][k - 1];
- if ((x = fabs(p) + fabs(q) + fabs(r))) {
- p /= x;
- q /= x;
- r /= x;
- }
- }
- if ((s = SIGN(sqrt(p * p + q * q + r * r), p))) {
- if (k == m) {
- if (l != m)
- a[k][k - 1] = -a[k][k - 1];
- }
- else
- a[k][k - 1] = -s * x;
- p += s;
- x = p / s;
- y = q / s;
- z = r / s;
- q /= p;
- r /= p;
- for (j = k; j <= nn; j++) {
- p = a[k][j] + q * a[k + 1][j];
- if (k != (nn - 1)) {
- p += r * a[k + 2][j];
- a[k + 2][j] -= p * z;
- }
- a[k + 1][j] -= p * y;
- a[k][j] -= p * x;
- }
- mmin = nn < k + 3 ? nn : k + 3;
- for (i = l; i <= mmin; i++) {
- p = x * a[i][k] + y * a[i][k + 1];
- if (k != (nn - 1)) {
- p += z * a[i][k + 2];
- a[i][k + 2] -= p * r;
- }
- a[i][k + 1] -= p * q;
- a[i][k] -= p;
- }
- }
- }
- }
- }
- } while (l < nn - 1);
- }
-}
Modified: grass/trunk/raster/r.texture/h_measure.h
===================================================================
--- grass/trunk/raster/r.texture/h_measure.h 2011-08-17 12:36:58 UTC (rev 47697)
+++ grass/trunk/raster/r.texture/h_measure.h 2011-08-17 12:45:12 UTC (rev 47698)
@@ -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);
Modified: grass/trunk/raster/r.texture/main.c
===================================================================
--- grass/trunk/raster/r.texture/main.c 2011-08-17 12:36:58 UTC (rev 47697)
+++ grass/trunk/raster/r.texture/main.c 2011-08-17 12:45:12 UTC (rev 47698)
@@ -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,144 +101,128 @@
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";
+ 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";
- /* Define the different flags */
+ 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");
- /* "Normalized" unused in the code ???
- flag0 = G_define_flag() ;
- flag0->key = 'N' ;
- flag0->description = _("Normalized") ;
- flag0->guisection = _("Features");
- */
+ flag_ind = G_define_flag();
+ flag_ind->key = 's';
+ flag_ind->description = _("Separate output for each angle (0, 45, 90, 135)");
- flag2 = G_define_flag();
- flag2->key = 'a';
- flag2->description = _("Angular Second Moment");
- flag2->guisection = _("Features");
+ flag_all = G_define_flag();
+ flag_all->key = 'a';
+ flag_all->description = _("Calculate all textural measurements");
- flag3 = G_define_flag();
- flag3->key = 'c';
- flag3->description = _("Contrast");
- flag3->guisection = _("Features");
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
- flag4 = G_define_flag();
- flag4->key = 'k';
- flag4->description = _("Correlation");
- flag4->guisection = _("Features");
+ 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"));
- flag5 = G_define_flag();
- flag5->key = 'v';
- flag5->description = _("Variance");
- flag5->guisection = _("Features");
+ 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);
- flag6 = G_define_flag();
- flag6->key = 'i';
- flag6->description = _("Inverse Diff Moment");
- flag6->guisection = _("Features");
+ 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++;
+ }
+ }
- flag7 = G_define_flag();
- flag7->key = 's';
- flag7->description = _("Sum Average");
- flag7->guisection = _("Features");
+ /* 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;
- 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");
-
-
- 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."));
-
infd = Rast_open_old(name, "");
/* determine the inputmap type (CELL/FCELL/DCELL) */
@@ -222,9 +231,35 @@
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 @@
}
/* 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,15 +307,8 @@
Rast_close(infd);
G_free(dcell_row);
- /* Now raster map is into memory. */
+ /* Now raster map is loaded to 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++;
-
/* *************************************************************************************************
*
* Compute of the matrix S.G.L.D. (Spatial Gray-Level Dependence Matrices) or co-occurrence matrix.
@@ -288,90 +317,82 @@
*
***************************************************************************************************/
- 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';
+ 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);
- strcpy(mapname, filename);
- strcat(mapname, suffixes[t_measure]);
- G_important_message(_("Calculating measure #%d <%s> (56 measures available)"),
- (t_measure + 1), mapname);
+ for (i = 0; i < n_outputs; i++)
+ Rast_set_f_null_value(fbuf[i], ncols);
- for (row = 0; row < nrows - (size - 1); row++) {
- G_percent(row, nrows, 2);
+ /*process the data */
+ for (col = first_col; col < last_col; col++) {
- /*process the data */
- for (col = 0; col < ncols - (size - 1); 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;
+ }
- /* Pointer for the s.w. */
- data = data + row;
- for (j = 0; j < size; j++)
- *(data + j) += col;
+ /* 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++) {
-/***********************************************************************************************
- *
- * 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(measure_idx[j]);
- 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;
+ 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);
More information about the grass-commit
mailing list