[GRASS-SVN] r67490 - grass-addons/grass7/imagery/i.spec.sam
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 5 13:27:32 PST 2016
Author: ychemin
Date: 2016-01-05 13:27:32 -0800 (Tue, 05 Jan 2016)
New Revision: 67490
Modified:
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
grass-addons/grass7/imagery/i.spec.sam/test1.py
grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py
Log:
Getting similar results as with python test
Modified: grass-addons/grass7/imagery/i.spec.sam/main.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/main.c 2016-01-05 19:21:16 UTC (rev 67489)
+++ grass-addons/grass7/imagery/i.spec.sam/main.c 2016-01-05 21:27:32 UTC (rev 67490)
@@ -57,7 +57,7 @@
char result_name[80];
char *result_prefix;
-mat_struct *open_files(char * matrixfile, char *img_grp, int type);
+mat_struct *open_files(char * matrixfile, char *img_grp);
float spectral_angle(vec_struct * Avector1, vec_struct * Avector2, int vtype);
DCELL myround(double x);
@@ -108,26 +108,34 @@
G_message("%s",result_prefix);
/*Creating A, the spectral signature matrix here*/
- A = open_files(parm.matrixfile->answer, parm.group->answer, 1);
- /* Spectral Matrix is stored in A now */
-
+ A = open_files(parm.matrixfile->answer, parm.group->answer);
+ /* Spectral Matrix is stored in A now */
+ G_message("Your incoming spectral signature matrix");
+ for (i=0; i<A->rows; i++)
+ {
+ for (j=0; j<A->cols; j++)
+ {
+ G_message("%f ", A->vals[i*A->rows+j]);
+ }
+ G_message("\n");
+ }
/* Check matrix orthogonality
* Ref: Youngsinn Sohn, Roger M. McCoy 1997: Mapping desert shrub rangeland
* 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);
+ AT = open_files(parm.matrixfile->answer, parm.group->answer);
G_message("/* Check matrix orthogonality*/");
for (i = 0; i < Ref.nfiles; i++) /* Ref.nfiles = matrixsize*/
{
Avector = G_matvect_get_column(AT, i); /* go columnwise through matrix*/
- G_message("Avector rows:%d cols:%d",Avector->rows,Avector->cols);
+ G_message("Avector rows:%d cols:%d, vals %f %f %f %f",Avector->rows,Avector->cols, Avector->vals[0], Avector->vals[1],Avector->vals[2],Avector->vals[3]);
for (j = 0; j < Ref.nfiles ; j++)
{
if (j !=i)
{
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("b rows:%d cols:%d, vals %f %f %f %f",b->rows,b->cols,b->vals[0],b->vals[1],b->vals[2],b->vals[3]);
G_message("process spectangle %d %d %d",i,j, Ref.nfiles);
spectangle = spectral_angle(Avector, b, RVEC);
G_message("processed spectangle");
@@ -203,9 +211,10 @@
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*/
+ /* Yann: b is a spectral signature extracted for a given pixel */
+ /*b = v_get(A->m); //m=rows; dim of vector=matrix size=Ref.nfiles*/
b = G_vector_init(A->cols,A->cols,CVEC);
- for (band = 0; band < Ref.nfiles-1; band++)
+ for (band = 0; band < Ref.nfiles; band++)
b->vals[band] = cell[band][col]; /* read input vector */
/* calculate spectral angle for current pixel
@@ -215,6 +224,7 @@
{
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]);
+ G_verbose_message("b: %f %f %f %f",b->vals[0],b->vals[1],b->vals[2],b->vals[3]);
spectangle = spectral_angle(Avector, b, CVEC);
result_cell[i][col] = myround (spectangle);
G_vector_free(Avector);
@@ -225,7 +235,6 @@
/* write the resulting rows: */
for (i = 0; i < Ref.nfiles; i++)
- for (i = 0; i < Ref.nfiles; i++)
Rast_put_d_row (resultfd[i], result_cell[i]);
} /* rows loop */
Modified: grass-addons/grass7/imagery/i.spec.sam/open.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/open.c 2016-01-05 19:21:16 UTC (rev 67489)
+++ grass-addons/grass7/imagery/i.spec.sam/open.c 2016-01-05 21:27:32 UTC (rev 67490)
@@ -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, int type)
+mat_struct *open_files(char *matrixfile, char *img_grp)
{
char *name, *mapset;
FILE *fp;
@@ -35,22 +35,23 @@
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"));
- 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;
- }
+ A = G_matrix_transpose(&A_input);/*transposed for spec angle process*/
+ /*if(type==1){
+ 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 = ");
- if (G_verbose() > G_verbose_std())
+ G_message("Your spectral matrix");
+ int j;
+ for (i=0; i<A->rows; i++)
{
- m_output(A);
+ for (j=0; j<A->cols; j++)
+ {
+ G_message("%f ", A->vals[i*A->rows+j]);
+ }
+ G_message("\n");
}
- */
matrixsize=A->rows;
G_message("/* open input files from group */");
Modified: grass-addons/grass7/imagery/i.spec.sam/spec_angle.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/spec_angle.c 2016-01-05 19:21:16 UTC (rev 67489)
+++ grass-addons/grass7/imagery/i.spec.sam/spec_angle.c 2016-01-05 21:27:32 UTC (rev 67490)
@@ -7,6 +7,7 @@
*****/
#include <stdio.h>
+#include <stdlib.h>
#include <math.h>
#include <grass/gmath.h>
#include "global.h"
@@ -47,13 +48,25 @@
/*vtmp1 = G_vector_init (Avector1->cols, Avector1->cols, CVEC);*/
vtmp1 = G_vector_copy(Avector1, CVEC);
}
+ int i;
+ G_verbose_message("spec_angle.c: Avector1");
+ for (i=0; i<Avector1->rows; i++)
+ G_verbose_message("%f ", Avector1->vals[i]);
+ G_verbose_message("spec_angle.c: Avector2");
+ for (i=0; i<Avector1->rows; i++)
+ G_verbose_message("%f ", Avector2->vals[i]);
+ /*Yann: Manually checked vtmp1 output: OK*/
vtmp1 = G_vector_product(Avector1, Avector2, vtmp1);
+ G_verbose_message("spec_angle.c: vtmp1");
+ for (i=0; i<vtmp1->rows; i++)
+ G_verbose_message("%f ", vtmp1->vals[i]);
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) */
+ norm2 = G_vector_norm_euclid(Avector1);/* calculate 2-norm (Euclidean) */
+ norm3 = G_vector_norm_euclid(Avector2);/* calculate 2-norm (Euclidean) */
G_vector_free(vtmp1);
-
+ G_verbose_message("norms: %f, %f, %f",norm1,norm2,norm3);
+ G_verbose_message("specangle: %f", acos(norm1 / (norm2 * norm3)) * 180 / M_PI);
/* Calculate angle and return in degree globally */
return (acos(norm1 / (norm2 * norm3)) * 180 / M_PI);
}
Modified: grass-addons/grass7/imagery/i.spec.sam/test1.py
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/test1.py 2016-01-05 19:21:16 UTC (rev 67489)
+++ grass-addons/grass7/imagery/i.spec.sam/test1.py 2016-01-05 21:27:32 UTC (rev 67490)
@@ -1,5 +1,9 @@
#!/usr/bin/env python
+#This creates 4 spectral maps
+#Each map has 4 horizontal sections
+#Each horizontal section is a test spectral class
+
# -*- coding: utf-8 -*-
from __future__ import (nested_scopes, generators, division, absolute_import,
with_statement, print_function, unicode_literals)
Modified: grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py 2016-01-05 19:21:16 UTC (rev 67489)
+++ grass-addons/grass7/imagery/i.spec.sam/test1_spec_angle.py 2016-01-05 21:27:32 UTC (rev 67490)
@@ -37,7 +37,7 @@
return np.sqrt((np.multiply(v,v)).sum(axis=0))
-#Define Spectral Angle of a vector pair
+#Define Spectral Angle (degrees) of a vector pair
def specangle(v1,v2):
return np.arccos(prod(v1,v2)/(l2(v1)*l2(v2)))*180/np.pi
@@ -46,21 +46,25 @@
print("cls1")
print(cls1)
print(cls1_wannabe)
+print(prod(cls1,cls1_wannabe))
print (specangle(cls1,cls1_wannabe))
print("cls2")
print(cls2)
print(cls2_wannabe)
+print(prod(cls2,cls2_wannabe))
print (specangle(cls2,cls2_wannabe))
print("cls3")
print(cls3)
print(cls3_wannabe)
+print(prod(cls3,cls3_wannabe))
print (specangle(cls3,cls3_wannabe))
print("cls4")
print(cls4)
print(cls4_wannabe)
+print(prod(cls4,cls4_wannabe))
print (specangle(cls4,cls4_wannabe))
More information about the grass-commit
mailing list