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

svn_grass at osgeo.org svn_grass at osgeo.org
Wed Dec 2 02:38:57 PST 2015


Author: ychemin
Date: 2015-12-02 02:38:57 -0800 (Wed, 02 Dec 2015)
New Revision: 67003

Added:
   grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py
Modified:
   grass-addons/grass7/imagery/i.spec.sam/main.c
   grass-addons/grass7/imagery/i.spec.sam/open.c
Log:
Added a manual test for specangle in Python

Modified: grass-addons/grass7/imagery/i.spec.sam/main.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/main.c	2015-12-02 09:49:26 UTC (rev 67002)
+++ grass-addons/grass7/imagery/i.spec.sam/main.c	2015-12-02 10:38:57 UTC (rev 67003)
@@ -57,7 +57,7 @@
 char result_name[80];
 char *result_prefix;
 
-mat_struct  *open_files(char * matrixfile, char *img_grp);
+mat_struct  *open_files(char * matrixfile, char *img_grp, int type);
 float spectral_angle(vec_struct * Avector1, vec_struct * Avector2, int vtype);
 DCELL myround(double x);
 
@@ -75,6 +75,7 @@
     } parm;
 
     mat_struct *A; /*first use in open.c G_matrix_set()*/
+    mat_struct *AT; /*first use in open.c G_matrix_set()*/
     char *group;
     float spectangle; /*numerical value of the spectral angle*/
 
@@ -88,7 +89,7 @@
 
 
     parm.group = G_define_standard_option(G_OPT_I_GROUP);
-    parm.group->description = "Imagery group to target for Spectral Mixture Analyis";
+    parm.group->description = "Imagery group to target for Spectral Mixture Analysis";
 
     parm.matrixfile = G_define_standard_option(G_OPT_F_INPUT);
     parm.matrixfile->description = "Matrix file containing spectral signatures";
@@ -107,7 +108,7 @@
     G_message("%s",result_prefix);
 
     /*Creating A, the spectral signature matrix here*/
-    A = open_files(parm.matrixfile->answer, parm.group->answer);
+    A = open_files(parm.matrixfile->answer, parm.group->answer, 1);
    /* Spectral Matrix is stored in A now */
 
   /* Check matrix orthogonality 
@@ -115,17 +116,21 @@
    *          using spectral unmixing and modeling spectral mixtrues with 
    *          TM data. Photogrammetric Engineering & Remote Sensing, Vol.63, No6.
    */
-     
+    AT = open_files(parm.matrixfile->answer, parm.group->answer, 0);
     G_message("/* Check matrix orthogonality*/"); 
     for (i = 0; i < Ref.nfiles; i++) /* Ref.nfiles = matrixsize*/
     {
-     Avector = G_matvect_get_row(A, i);  /* go columnwise through matrix*/
+     Avector = G_matvect_get_column(AT, i);  /* go columnwise through matrix*/
+     G_message("Avector rows:%d cols:%d",Avector->rows,Avector->cols);
      for (j = 0; j < Ref.nfiles ; j++)
 	{
 	 if (j !=i)
 	    {
-	     b = G_matvect_get_row(A, j);      /* compare with next col in A */
+	     b = G_matvect_get_column(AT, j);      /* compare with next col in A */
+             G_message("b rows:%d cols:%d",b->rows,b->cols);
+             G_message("process spectangle %d %d %d",i,j, Ref.nfiles);
 	     spectangle = spectral_angle(Avector, b, RVEC);
+             G_message("processed spectangle");
 	     anglefield[i][j]= spectangle;
 	     G_vector_free(b);
 	    }

Modified: grass-addons/grass7/imagery/i.spec.sam/open.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/open.c	2015-12-02 09:49:26 UTC (rev 67002)
+++ grass-addons/grass7/imagery/i.spec.sam/open.c	2015-12-02 10:38:57 UTC (rev 67003)
@@ -13,7 +13,7 @@
 
 int G_matrix_read2(FILE * fp, mat_struct * out); /* Modified version of G_matrix_read(..). */
 
-mat_struct *open_files(char *matrixfile, char *img_grp)
+mat_struct *open_files(char *matrixfile, char *img_grp, int type)
 {
     char *name, *mapset;
     FILE *fp;
@@ -35,11 +35,14 @@
     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_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");
+    if(type==1)
+        A = G_matrix_transpose(&A_input);/*transposed for spec angle process*/
+    if(type==0){
+        A = G_matrix_copy(&A_input);/*for matrix orthogonality check*/
+        return A;
+    }
+/*    if(A->rows < A->cols)
+	G_fatal_error("Need m (rows) >= n (cols) to obtain least squares fit\n");*/
     /*Only for debug, so temporary disabled*/
     /*
     G_verbose_message("Your spectral matrix = ");
@@ -48,7 +51,7 @@
 	m_output(A);
     }
     */
-    matrixsize=A->rows; /*changed cols to rows*/
+    matrixsize=A->rows;
 
     G_message("/* open input files from group */");
 /* open input files from group */

Added: grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py	                        (rev 0)
+++ grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py	2015-12-02 10:38:57 UTC (rev 67003)
@@ -0,0 +1,65 @@
+#!/usr/bin/env
+import numpy as np
+
+#Spectral data to load as maps
+# Band: r g b i1 i2 i3
+# Spektren zeilenweise eingeben!
+# 1. Sagebrush 
+# 2. Saltbush
+# 3. Ground
+# 4. Dry Grass
+#row0:  8.87  13.14  11.71  35.85 
+#Matrix: 4 by 4
+#row0:  8.87  13.14  11.71  35.85
+#row1: 13.59  20.12  19.61  50.66
+#row2: 28.26  34.82  38.27  40.1
+#row3: 10.54  16.35  23.7   38.98
+
+#Define the spectral signatures for each land use class
+cls1=[8.87, 13.14, 11.71, 35.85]
+cls2=[13.59, 20.12, 19.61, 50.66]
+cls3=[28.26, 34.82, 38.27, 40.10]
+cls4=[10.54, 16.35, 23.70, 38.98]
+
+#Generate a disturbed class manually 
+cls1_wannabe=[9, 14, 12, 36]
+cls2_wannabe=[14, 21, 20, 51]
+cls3_wannabe=[29, 35, 39, 41]
+cls4_wannabe=[11, 17, 24, 39]
+
+
+#Define Norm of a Vector product 
+def prod(v1, v2):
+	return np.sum(np.multiply(v1,v2))
+
+#Define Euclidian Norm of a vector ("L2-norm")
+def l2(v):
+	return np.sqrt((np.multiply(v,v)).sum(axis=0))
+
+
+#Define Spectral Angle of a vector pair
+def specangle(v1,v2):
+	return np.arccos(prod(v1,v2)/(l2(v1)*l2(v2)))*180/np.pi
+
+
+#Run for all 4 classes
+print("cls1")
+print(cls1)
+print(cls1_wannabe)
+print (specangle(cls1,cls1_wannabe))
+
+print("cls2")
+print(cls2)
+print(cls2_wannabe)
+print (specangle(cls2,cls2_wannabe))
+
+print("cls3")
+print(cls3)
+print(cls3_wannabe)
+print (specangle(cls3,cls3_wannabe))
+
+print("cls4")
+print(cls4)
+print(cls4_wannabe)
+print (specangle(cls4,cls4_wannabe))
+



More information about the grass-commit mailing list