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

[svn_grass at osgeo.org]

Author: neteler
Date: 2015-03-03 01:50:48 -0800 (Tue, 03 Mar 2015)
New Revision: 64788

Modified:
   grass-addons/grass7/imagery/i.spec.sam/global.h
   grass-addons/grass7/imagery/i.spec.sam/spec_angle.c
Log:
i.spec.sam Addon: partial update with code from i.spec.unmix

Modified: grass-addons/grass7/imagery/i.spec.sam/global.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/global.h	2015-03-03 09:46:34 UTC (rev 64787)
+++ grass-addons/grass7/imagery/i.spec.sam/global.h	2015-03-03 09:50:48 UTC (rev 64788)
@@ -1,7 +1,15 @@
-#include <math.h>
-#include <grass/imagery.h>
-#include "matrix.h"
+#ifndef __GLOBAL_H__
+#define __GLOBAL_H__
 
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/la.h>
+
+#ifndef GLOBAL
+#define GLOBAL extern
+#endif
+
 #define MAXFILES 255
 
 extern MAT *A;

Modified: grass-addons/grass7/imagery/i.spec.sam/spec_angle.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.sam/spec_angle.c	2015-03-03 09:46:34 UTC (rev 64787)
+++ grass-addons/grass7/imagery/i.spec.sam/spec_angle.c	2015-03-03 09:50:48 UTC (rev 64788)
@@ -1,33 +1,18 @@
 /* Spectral angle mapping
- * (c) 1998 Markus Neteler, Hannover
+ * (c) 1999 Markus Neteler, Hannover, GPL >= 2.0
  *
- * 26. Oct. 1998 - V. 0.1
+ * 25. Nov. 1998 - V. 0.2
  *
- ****************************************************************************
- ** Based on Meschach Library
- ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- ****************************************************************************
- *
- * Cited references are from
- *    Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
- *       Proceedings of the centre for Mathematics and its Applicaions.
- *       The Australian National University. Vol. 32.
- *       ISBN 0 7315 1900 0
+ * Updated to LAPACK/BLAS by Mohammed Rashad, 2012
  *****/
 
-
-#define MY_PI 3.141592653589793
-
 #include <stdio.h>
 #include <math.h>
-#include <meschach/matrix.h>
+#include <grass/gmath.h>
 #include "global.h"
 
 
-float spectral_angle() /* returns spectral angle*/
-{
-
-/* input MAT A, VEC Avector
+/* input mat_struct A, vec_struct Avector1, Avector2
  * output cur_angle
  *
  *                   v_DN * v_reference
@@ -35,24 +20,32 @@
  *               ||v_DN|| * ||v_reference||
  * 
  *
- *                 b * Avector
- *            = -----------------------
- *              ||b|| * ||Avector||
+ *                  Avector1 * Avector2
+ *            = ---------------------------
+ *              ||Avector1|| * ||Avector2||
+ *
+ * 
+ * Ref.: van der Meer, F. 1997: Mineral mapping and Landsat Thematic Mapper 
+ *                             image Classification using spectral unmixing.
+ *       Geocarto International, Vol.12, no.3 (Sept.). pp. 27-40
  */
 
+float spectral_angle(vec_struct * Avector1, vec_struct * Avector2)
+{
+    vec_struct *vtmp1;
+    double norm1, norm2, norm3;
 
-  VEC *vtmp1;
-  double norm1, norm2, norm3;
-        
-  /* Measure spectral angle*/
+    /* Measure spectral angle */
 
-   vtmp1 = v_star(Avector, b, VNULL); /* multiply with b vector */
-   norm1 = v_norm1(vtmp1);       /* calculate 1-norm */
-   norm2 = v_norm2(Avector);     /* calculate 2-norm (Euclidean) */
-   norm3 = v_norm2(b);           /* calculate 2-norm (Euclidean) */
+    /* multiply one A column with second */
+    vtmp1 = G_vector_init (Avector1->rows, Avector1->rows, 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) */
 
-   V_FREE(vtmp1);
-      
-   curr_angle = (acos(norm1/(norm2 * norm3)) * 180/MY_PI);  /* Calculate angle */
-        /* return in degree*/
+    G_vector_free(vtmp1);
+
+    /* Calculate angle and return in degree globally */
+    return (acos(norm1 / (norm2 * norm3)) * 180 / M_PI);
 }