[GRASS-SVN] r66760 - grass-addons/grass7/imagery/i.spec.sam

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Nov 6 06:21:23 PST 2015


Author: ychemin
Date: 2015-11-06 06:21:23 -0800 (Fri, 06 Nov 2015)
New Revision: 66760

Modified:
   grass-addons/grass7/imagery/i.spec.sam/global.h
   grass-addons/grass7/imagery/i.spec.sam/main.c
   grass-addons/grass7/imagery/i.spec.sam/open.c
   grass-addons/grass7/imagery/i.spec.sam/spec_angle.c
Log:
Fixed CELL to DCELL throughout, fixed spec_angle tmp vect

Modified: grass-addons/grass7/imagery/i.spec.sam/global.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/global.h	2015-11-06 13:59:48 UTC (rev 66759)
+++ grass-addons/grass7/imagery/i.spec.sam/global.h	2015-11-06 14:21:23 UTC (rev 66760)
@@ -20,13 +20,13 @@
 extern char *group;
 extern struct Ref Ref;
 
-extern CELL **cell;
+extern DCELL **cell;
 extern int *cellfd;
 
-extern CELL **result_cell;
+extern DCELL **result_cell;
 extern int *resultfd;
 
-extern CELL **error_cell;
+extern DCELL **error_cell;
 extern int  error_fd;
 
 extern char result_name[80];

Modified: grass-addons/grass7/imagery/i.spec.sam/main.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/main.c	2015-11-06 13:59:48 UTC (rev 66759)
+++ grass-addons/grass7/imagery/i.spec.sam/main.c	2015-11-06 14:21:23 UTC (rev 66760)
@@ -45,21 +45,21 @@
 
 struct Ref Ref;
 
-CELL **cell;
+DCELL **cell;
 int *cellfd;
 
-CELL **result_cell;
+DCELL **result_cell;
 int *resultfd;
 
-CELL **error_cell;
+DCELL **error_cell;
 int  error_fd;
 
 char result_name[80];
 char *result_prefix;
 
 mat_struct  *open_files(char * matrixfile, char *img_grp);
-float spectral_angle();
-CELL myround(double x);
+float spectral_angle(vec_struct * Avector1, vec_struct * Avector2, int vtype);
+DCELL myround(double x);
 
 int main(int argc,char * argv[])
 {
@@ -125,7 +125,7 @@
 	 if (j !=i)
 	    {
 	     b = G_matvect_get_row(A, j);      /* compare with next col in A */
-	     spectangle = spectral_angle(Avector, b);
+	     spectangle = spectral_angle(Avector, b, RVEC);
 	     anglefield[i][j]= spectangle;
 	     G_vector_free(b);
 	    }
@@ -161,9 +161,9 @@
      *        using singular value decomposition.  IGARSS 1989: 12th Canadian
      *           symposium on Remote Sensing. Vol.4 pp.2069-2072
      */
-    G_message("Singular values ouput vector init");
+    G_verbose_message("Singular values ouput vector init");
     double *svdvalues = G_alloc_vector(A->cols);
-    G_message("Singular values of Matrix A: copy A into double **");
+    G_verbose_message("Singular values of Matrix A: copy A into double **");
     double **Avals = G_alloc_matrix(A->cols,A->rows);
     int count=0;
     for (i = 0; i < A->cols ; i++){
@@ -172,35 +172,34 @@
       count++;
      }
     }
-    G_message("Singular values of Matrix A: run svdval");
-    if(G_math_svdval( svdvalues, Avals, A->cols, A->rows))
+    G_verbose_message("Singular values of Matrix A: run svdval");
+    if(G_math_svdval(svdvalues, Avals, A->cols, A->rows))
         G_fatal_error("Error in singular value decomposition, exiting...\n");
-    G_message("/*Experimental: display values (replace v_output() in original version)*/");
-    /*Experimental: display values (replace v_output() in original version)*/
+    G_verbose_message("/*display values (replace v_output() in original version)*/");
+    /*display values (replace v_output() in original version)*/
     for(i=0;i<A->cols;i++)
         G_message("%f", svdvalues[i]);
 
-    G_message("/* alright, start Spectral angle mapping */");
+    G_verbose_message("/* alright, start Spectral angle mapping */");
     /* alright, start Spectral angle mapping */
     nrows = Rast_window_rows();
     ncols = Rast_window_cols();
     
     G_verbose_message("Calculating for %i x %i = %i pixels:\n",nrows,ncols, (ncols * ncols));
     G_verbose_message("%s ... ", G_program_name());
-    G_message("Calculating for %i x %i = %i pixels:\n",nrows,ncols, (ncols * ncols));
     G_message("%s ... ", G_program_name());
 
     for (row = 0; row < nrows; row++)                 /* rows loop*/
     {
 	G_percent(row, nrows, 2);
 	for (band = 0; band < Ref.nfiles; band++)     /* get row for all bands*/
-	    Rast_get_c_row (cellfd[band], cell[band], row);
+	    Rast_get_d_row (cellfd[band], cell[band], row);
 
 	for (col = 0; col < ncols; col++)             /* cols loop, work pixelwise for all bands */
 	{
 	    /* get pixel values of each band and store in b vector: */
 	     /*b = v_get(A->m);*/                   /* m=rows; dimension of vector = matrix size = Ref.nfiles*/
-            b = G_vector_init(A->rows,A->rows,RVEC);
+            b = G_vector_init(A->cols,A->cols,CVEC);
 	    for (band = 0; band < Ref.nfiles-1; band++)
                  b->vals[band] = cell[band][col];  /* read input vector */
 	   
@@ -209,8 +208,9 @@
              * and write result in full degree */
              for (i = 0; i < Ref.nfiles; i++) /* Ref.nfiles = matrixsize*/
              {
-              Avector = G_matvect_get_row(A, i);  /* go row-wise through matrix*/
-	      spectangle = spectral_angle(Avector, b);
+              Avector = G_matvect_get_column(A, i);  /* go row-wise through matrix*/
+              G_verbose_message("Av: %f %f %f %f",Avector->vals[0],Avector->vals[1],Avector->vals[2],Avector->vals[3]);
+	      spectangle = spectral_angle(Avector, b, CVEC);
 	      result_cell[i][col] = myround (spectangle);
 	      G_vector_free(Avector);
              }
@@ -220,7 +220,8 @@
 
 	/* write the resulting rows: */
         for (i = 0; i < Ref.nfiles; i++)
-          Rast_put_c_row (resultfd[i], result_cell[i]);
+        for (i = 0; i < Ref.nfiles; i++)
+          Rast_put_d_row (resultfd[i], result_cell[i]);
 
     } /* rows loop */
     G_percent(row, nrows, 2);
@@ -242,9 +243,9 @@
 } /* main*/
 
 
-CELL myround (double x)
+DCELL myround (double x)
   {
-    CELL n;
+    DCELL n;
     
     if (x >= 0.0)
         n = x + .5;

Modified: grass-addons/grass7/imagery/i.spec.sam/open.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/open.c	2015-11-06 13:59:48 UTC (rev 66759)
+++ grass-addons/grass7/imagery/i.spec.sam/open.c	2015-11-06 14:21:23 UTC (rev 66760)
@@ -32,11 +32,11 @@
 
     G_message("matrixfile read");
 
-    A = G_matrix_init(A_input.rows, A_input.cols, A_input.rows);
+    A = G_matrix_init(A_input.cols, A_input.rows, A_input.cols);/*changed r/c*/
     if (A == NULL)
         G_fatal_error(_("Unable to allocate memory for matrix"));
 
-    A = G_matrix_copy(&A_input);
+    A = G_matrix_transpose(&A_input);/*transposed for spec angle process*/
 
     if(A->rows < A->cols)
 	G_fatal_error("Need m (rows) >= n (cols) to obtain least squares fit\n");
@@ -48,7 +48,7 @@
 	m_output(A);
     }
     */
-    matrixsize=A->cols;
+    matrixsize=A->rows; /*changed cols to rows*/
 
     G_message("/* open input files from group */");
 /* open input files from group */
@@ -75,11 +75,11 @@
 
    G_message("/* get memory for input files */");
    /* get memory for input files */
-    cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+    cell = (DCELL **) G_malloc (Ref.nfiles * sizeof (DCELL *));
     cellfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
     for (i=0; i < Ref.nfiles; i++)
     {
-	cell[i] = Rast_allocate_c_buf();
+	cell[i] = Rast_allocate_d_buf();
 	name = Ref.file[i].name;
 	mapset = Ref.file[i].mapset;
 	G_verbose_message("Opening input file no. %i [%s]\n", (i+1), name);
@@ -89,14 +89,14 @@
 
     G_message("/* open files for results*/");
   /* open files for results*/
-    result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+    result_cell = (DCELL **) G_malloc (Ref.nfiles * sizeof (DCELL *));
     resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
     for (i=0; i < Ref.nfiles; i++)
     {
 	sprintf(result_name, "%s.%d", result_prefix, (i+1));
 	G_verbose_message("Opening output file [%s]\n", result_name);	 
-	result_cell[i] = Rast_allocate_c_buf();
-	resultfd[i] = Rast_open_c_new (result_name);
+	result_cell[i] = Rast_allocate_d_buf();
+	resultfd[i] = Rast_open_new (result_name, DCELL_TYPE);
     }
 
     G_message("open.c: Returning A");

Modified: grass-addons/grass7/imagery/i.spec.sam/spec_angle.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/spec_angle.c	2015-11-06 13:59:48 UTC (rev 66759)
+++ grass-addons/grass7/imagery/i.spec.sam/spec_angle.c	2015-11-06 14:21:23 UTC (rev 66760)
@@ -30,7 +30,7 @@
  *       Geocarto International, Vol.12, no.3 (Sept.). pp. 27-40
  */
 
-float spectral_angle(vec_struct * Avector1, vec_struct * Avector2)
+float spectral_angle(vec_struct * Avector1, vec_struct * Avector2, int vtype)
 {
     vec_struct *vtmp1;
     double norm1, norm2, norm3;
@@ -38,8 +38,16 @@
     /* Measure spectral angle */
 
     /* multiply one A column with second */
-    vtmp1 = G_vector_init (Avector1->cols, Avector1->cols, RVEC);
-    vtmp1 = G_vector_product(Avector1, Avector2,vtmp1);
+    if(vtype == RVEC){
+        G_verbose_message("spec_angle.c: Using RVEC type for vtmp1");
+        /*vtmp1 = G_vector_init (Avector1->cols, Avector1->cols, RVEC);*/
+        vtmp1 = G_vector_copy(Avector1, RVEC);
+    }else{
+        G_verbose_message("spec_angle.c: Using CVEC type for vtmp1");
+        /*vtmp1 = G_vector_init (Avector1->cols, Avector1->cols, CVEC);*/
+        vtmp1 = G_vector_copy(Avector1, CVEC);
+    }
+    vtmp1 = G_vector_product(Avector1, Avector2, vtmp1);
     norm1 = G_vector_norm1(vtmp1);	/* calculate 1-norm */
     norm2 = G_vector_norm_euclid(Avector1);	/* calculate 2-norm (Euclidean) */
     norm3 = G_vector_norm_euclid(Avector2);	/* calculate 2-norm (Euclidean) */



More information about the grass-commit mailing list