[GRASS-SVN] r54176 - grass-addons/grass6/imagery/i.spec.unmix
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Dec 4 02:17:41 PST 2012
Author: rashadkm
Date: 2012-12-04 02:17:41 -0800 (Tue, 04 Dec 2012)
New Revision: 54176
Added:
grass-addons/grass6/imagery/i.spec.unmix/la_extra.c
grass-addons/grass6/imagery/i.spec.unmix/la_extra.h
Modified:
grass-addons/grass6/imagery/i.spec.unmix/Makefile
grass-addons/grass6/imagery/i.spec.unmix/global.h
grass-addons/grass6/imagery/i.spec.unmix/hist.c
grass-addons/grass6/imagery/i.spec.unmix/histogram.c
grass-addons/grass6/imagery/i.spec.unmix/main.c
grass-addons/grass6/imagery/i.spec.unmix/open.c
grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c
Log:
ported i.spec.unmix to grass6.5
Modified: grass-addons/grass6/imagery/i.spec.unmix/Makefile
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/Makefile 2012-12-04 08:11:17 UTC (rev 54175)
+++ grass-addons/grass6/imagery/i.spec.unmix/Makefile 2012-12-04 10:17:41 UTC (rev 54176)
@@ -2,8 +2,8 @@
PGM = i.spec.unmix
-LIBES = $(IMAGERYLIB) $(GMATHLIB) $(GISLIB)
-DEPENDENCIES= $(IMAGERYDEP) $(GMATHDEP) $(GISDEP)
+LIBES = $(IMAGERYLIB) $(GMATHLIB) $(GISLIB) $(CURSES)
+DEPENDENCIES= $(IMAGERYDEP) $(GMATHDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass-addons/grass6/imagery/i.spec.unmix/global.h
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/global.h 2012-12-04 08:11:17 UTC (rev 54175)
+++ grass-addons/grass6/imagery/i.spec.unmix/global.h 2012-12-04 10:17:41 UTC (rev 54176)
@@ -1,20 +1,16 @@
-#include <math.h>
-#include <grass/imagery.h>
-#include <grass/la.h> /* LAPACK/BLAS */
+#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
-GLOBAL mat_struct *A, *A_tilde_trans;
-GLOBAL vec_struct *Avector1, *Avector2, *b, *fraction;
-GLOBAL int matrixsize;
-GLOBAL double svd_error;
-GLOBAL float curr_angle;
-
-GLOBAL char *group;
GLOBAL struct Ref Ref;
GLOBAL CELL **cell;
@@ -29,15 +25,10 @@
GLOBAL CELL *iter_cell;
GLOBAL int iter_fd;
-GLOBAL char result_name[80];
-GLOBAL char *result_prefix, *matrixfile, *error_name, *iter_name;
-GLOBAL struct
- {
- struct Flag *quiet;
- } flag;
-
-GLOBAL struct
- {
- struct Flag *veryquiet;
- } flag2;
+GLOBAL float spectral_angle(vec_struct *, vec_struct *);
+GLOBAL int do_histogram(char *, char *);
+GLOBAL void make_history(char *, char *, char *);
+GLOBAL int open_files(char *, char *, char *, char *, mat_struct *A);
+
+#endif /* __GLOBAL_H__ */
Modified: grass-addons/grass6/imagery/i.spec.unmix/hist.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/hist.c 2012-12-04 08:11:17 UTC (rev 54175)
+++ grass-addons/grass6/imagery/i.spec.unmix/hist.c 2012-12-04 10:17:41 UTC (rev 54176)
@@ -1,11 +1,12 @@
+#include <stdio.h>
#include <grass/gis.h>
-void make_history(name, group, matrixfile)
- char *name, *group, *matrixfile;
+
+void make_history (char *name, char *group, char *matrixfile)
{
struct History hist;
- if(G_read_history (name, G_mapset(), &hist) >= 0)
+ if (G_read_history (name, G_mapset(), &hist) >= 0)
{
sprintf (hist.datsrc_1, "Group: %s", group);
sprintf (hist.datsrc_2, "Matrix file: %s", matrixfile);
Modified: grass-addons/grass6/imagery/i.spec.unmix/histogram.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/histogram.c 2012-12-04 08:11:17 UTC (rev 54175)
+++ grass-addons/grass6/imagery/i.spec.unmix/histogram.c 2012-12-04 10:17:41 UTC (rev 54176)
@@ -1,10 +1,7 @@
-/* this code is from r.support */
-
#include <grass/gis.h>
-int do_histogram (name, mapset)
- char *name;
- char *mapset;
+
+int do_histogram (char *name, char *mapset)
{
CELL *cell;
struct Cell_head cellhd;
@@ -13,26 +10,32 @@
int fd;
if (G_get_cellhd (name, mapset, &cellhd) < 0)
- return 0;
+ return 0;
+
G_set_window (&cellhd);
fd = G_open_cell_old (name, mapset);
if (fd < 0)
- return 0;
- nrows = G_window_rows();
- ncols = G_window_cols();
- cell = G_allocate_cell_buf();
+ return 0;
+ nrows = G_window_rows ();
+ ncols = G_window_cols ();
+ cell = G_allocate_cell_buf ();
+
G_init_cell_stats (&statf);
for (row = 0; row < nrows; row++)
{
- if (G_get_map_row_nomask (fd, cell, row) < 0)
- break;
- G_update_cell_stats (cell, ncols, &statf);
+ if (G_get_map_row_nomask (fd, cell, row) < 0)
+ break;
+
+ G_update_cell_stats (cell, ncols, &statf);
}
G_close_cell (fd);
- free (cell);
+ G_free (cell);
+
if (row == nrows)
- G_write_histogram_cs (name, &statf);
+ G_write_histogram_cs (name, &statf);
+
G_free_cell_stats (&statf);
- return;
+
+ return 0;
}
Added: grass-addons/grass6/imagery/i.spec.unmix/la_extra.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/la_extra.c (rev 0)
+++ grass-addons/grass6/imagery/i.spec.unmix/la_extra.c 2012-12-04 10:17:41 UTC (rev 54176)
@@ -0,0 +1,577 @@
+
+
+#include <stdio.h> /* needed here for ifdef/else */
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+
+/********
+ ******** only compile this LAPACK/BLAS wrapper file if g2c.h is present!
+ ********/
+#include <grass/config.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/la.h>
+#include "la_extra.h"
+
+
+
+
+vec_struct *
+G_matvect_get_column2(mat_struct *mt, int col)
+{
+ int i; /* loop */
+ vec_struct *vc1;
+
+ if(col < 0 || col >= mt->cols) {
+ G_warning(_("Specified matrix column index is outside range"));
+ return NULL;
+ }
+
+ if(!mt->is_init) {
+ G_warning(_("Matrix is not initialised"));
+ return NULL;
+ }
+
+ if( (vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) == NULL ) {
+ G_warning(_("Could not allocate space for vector structure"));
+ return NULL;
+ }
+
+ for ( i = 0; i < mt->rows; i++ )
+ {
+ double dd = G_matrix_get_element(mt, i, col);
+ //G_warning("element at row=%d, col=%d is %lf",i,col,dd);
+ G_matrix_set_element( (mat_struct *)vc1, i, 0, dd);
+
+ }
+
+ return vc1;
+}
+
+
+void
+G_matrix_print2 (mat_struct *mt, const char* name)
+{
+ int i, j;
+
+
+ if(mt!=NULL)
+ {
+ G_message ("start matrix(%s)", name);
+ G_message ("Size: %d x %d", mt->rows, mt->cols);
+
+ for( i = 0; i < mt->rows; i++ )
+ {
+ char buf[2048], numbuf[640];
+ sprintf( buf, "row%d: ", i);
+ for( j = 0; j < mt->cols; j++ )
+ {
+
+ double element = G_matrix_get_element(mt, i, j);
+ sprintf( numbuf, "%14.6f",element);
+ strcat(buf, numbuf);
+ //if( j < mt->cols - 1 )
+ //strcat(buf, ", ");
+ }
+ G_message ("%s", buf);
+ }
+
+ G_message ("end matrix(%s)",name);
+ }
+
+/*
+ for( i = 0; i < mt->rows; i++ ) {
+ strcpy(buf, "");
+
+ for( j = 0; j < mt->cols; j++ ) {
+
+ sprintf( numbuf, "%14.6f", G_matrix_get_element(mt, i, j) );
+ strcat(buf, numbuf);
+ if( j < mt->cols - 1 )
+ strcat(buf, ", ");
+ }
+
+ G_message ("%s", buf);
+ }
+
+ fprintf (stderr, "\n");
+
+*/
+}
+
+
+mat_struct *G_matrix_resize(mat_struct * in, int rows, int cols)
+{
+
+ mat_struct *matrix;
+
+
+ matrix = G_matrix_init(rows,cols,rows);
+
+
+ int i,j,p, index = 0;
+ for ( i=0; i<rows; i++ )
+ {
+
+
+
+ //A_row = A_v[i]; b_v = b->ve;
+ for ( j=0; j<cols; j++ )
+ {
+
+ matrix->vals[index++] = in->vals[i + j * cols];
+
+ // sum += in->vals[i + j * cols] *b->ve[j];
+ //out->ve[i] = sum;
+ }
+
+ }
+
+ int old_size = in->rows * in->cols;
+
+ int new_size = rows * cols;
+ if(new_size > old_size)
+ for(p=old_size;p<new_size;p++)
+ matrix->vals[p] = 0.0;
+
+
+ return matrix;
+}
+
+
+
+mat_struct *sm_mlt(double scalar,mat_struct *matrix, mat_struct *out)
+
+{
+ int m,n,i,j;
+
+ int index = 0;
+
+ if ( matrix==NULL )
+ G_fatal_error("sm_mlt1(error)");
+
+ if ( out==NULL)
+ out = G_matrix_init(matrix->rows, matrix->cols,matrix->rows);
+
+ if (out->rows != matrix->rows || out->cols != matrix->cols )
+ out = G_matrix_resize(out,matrix->rows,matrix->cols);
+
+m = matrix->rows;
+n = matrix->cols;
+ for ( i=0; i<m; i++ )
+ {
+ //__smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
+ /**************************************************/
+ for ( j=0; j<n; j++ )
+ {
+ out->vals[index++] = scalar*matrix->vals[i+ j * m];
+ }}
+
+ //G_matrix_print(matrix,"matrix");
+ //G_matrix_print(matrix,"out");
+ /**************************************************/
+ return (out);
+}
+
+VEC *G_vec_copy(VEC *in)
+{
+ int i;
+ VEC *out;
+ if ( !in )
+ G_fatal_error("v_copy(error1)");
+
+ int dim = in->dim;
+
+ out = G_vec_get(dim);
+
+
+ for( i = 0; i < dim; i++ )
+ {
+ out->ve[i]= in->ve[i];
+ }
+
+ return (out);
+}
+
+
+double v_norm2(VEC *x)
+{
+ int i, dim;
+ double s, sum;
+
+ if ( !x )
+ G_fatal_error("v_norm2(error1)");
+
+
+ dim = x->dim;
+
+ sum = 0.0;
+
+ for ( i = 0; i < dim; i++ )
+ {
+ sum += x->ve[i] * x->ve[i];
+ }
+
+
+
+ return sqrt(sum);
+}
+
+
+
+VEC *v_sub(VEC *vec1,VEC *vec2,VEC *out)
+{
+ /* u_int i, dim; */
+ /* Real *out_ve, *vec1_ve, *vec2_ve; */
+
+ if ( !vec1 || !vec2 )
+ G_fatal_error("v_sub1(error)");
+
+ if ( vec1->dim != vec2->dim )
+ G_fatal_error("v_sub2(error)");
+
+ if ( out == NULL)
+ out = G_vec_get(vec1->dim);
+ //out = v_resize(out,A->m);
+
+ //G_vec_print(vec1, "vec1");
+ //G_vec_print(vec2, "vec2");
+
+ if ( out->dim != vec1->dim )
+ out = G_vec_resize(out,vec1->dim);
+
+ int i;
+ for ( i = 0; i < vec1->dim; i++ )
+ out->ve[i] = vec1->ve[i] - vec2->ve[i];
+
+ /************************************************************
+ dim = vec1->dim;
+ out_ve = out->ve; vec1_ve = vec1->ve; vec2_ve = vec2->ve;
+ for ( i=0; i<dim; i++ )
+ out->ve[i] = vec1->ve[i]-vec2->ve[i];
+ (*out_ve++) = (*vec1_ve++) - (*vec2_ve++);
+ ************************************************************/
+
+ return (out);
+}
+
+
+
+VEC *mv_mlt(mat_struct *A, VEC *b, VEC *out)
+{
+ unsigned int i, m, n,j;
+ double **A_v, *b_v /*, *A_row */;
+ /* register Real sum; */
+
+ //G_vec_print(out, "xx");
+
+// if ( A==(MAT *)NULL || b==(VEC *)NULL )
+ // error(E_NULL,"mv_mlt");
+
+//G_matrix_print(A, "A(mlt)");
+
+ if ( A->cols != b->dim )
+
+ G_fatal_error("mv_mlt1(error)");
+
+
+ if ( b == out )
+ G_fatal_error("mv_mlt2(error)");
+ //error(E_INSITU,"mv_mlt");
+ if(!out)
+ {
+ G_fatal_error("mv_mltsss3(error)");
+ exit(1);
+ out = G_vec_get2(A->rows, out);
+ }
+ if ( out->dim != A->rows)
+ {
+ G_fatal_error("mv_mlt3(error)");
+ exit(1);
+ out = G_vec_resize(out,A->rows);
+ }
+ //out = G_vec_get(A->rows);
+ //out = v_resize(out,A->m);
+
+
+
+ //if ( out->dim != A->rows )
+
+
+ m = A->rows; n = A->cols;
+ A_v = A->vals; b_v = b->ve;
+
+ //G_vec_print(out, "xx");
+
+ for ( i=0; i<m; i++ )
+ {
+ double sum=0.0;
+int width = A->rows;
+
+ //A_row = A_v[i]; b_v = b->ve;
+ for ( j=0; j<n; j++ )
+ {
+
+ sum += A->vals[i + j * width] *b->ve[j];
+ out->ve[i] = sum;
+
+
+//G_message("sum: %lf of j=%d", b->ve[j],j);
+
+}
+ }
+
+ //G_vec_print(out, "xx");
+
+ return out;
+}
+
+
+
+
+VEC *G_vec_resize(VEC * in, int size)
+{
+
+ VEC *vector;
+
+
+ vector = (VEC *)G_malloc( sizeof(VEC) );
+
+
+ vector->ve = (double *)G_malloc(size * sizeof(double) );
+ int i,j;
+
+ G_message(":%d",in->dim);
+for( i = 0; i < in->dim; i++ )
+ {
+
+ vector->ve[i]= in->ve[i];
+ G_message("ss:%lf",in->ve[i]);
+
+ }
+
+ if(size > in->dim)
+ for(j=i;j<size;j++)
+ vector->ve[j]= 0.0;
+
+ vector->dim = vector->max_dim = size;
+
+ return vector;
+}
+
+
+
+
+
+
+
+VEC *G_vec_get(int size)
+{
+
+ VEC *vector;
+
+
+ vector = (VEC *)G_malloc( sizeof(VEC) );
+
+
+ vector->ve = (double *)G_malloc(size * sizeof(double) );
+ int i;
+for( i = 0; i < size; i++ )
+ {
+
+ vector->ve[i]= 0.0;
+
+
+ }
+
+ vector->dim = vector->max_dim = size;
+ //G_vec_print(vector, "vector");
+ return vector;
+}
+
+
+VEC *G_vec_get2(int size, VEC *vector)
+{
+
+ //VEC *vector;
+
+
+ vector = (VEC *)G_malloc( sizeof(VEC) );
+
+
+ vector->ve = (double *)G_malloc(size * sizeof(double) );
+ int i;
+for( i = 0; i < size; i++ )
+ {
+
+ vector->ve[i]= 0.0;
+
+
+ }
+
+ vector->dim = vector->max_dim = size;
+ //G_vec_print(vector, "vector");
+ return vector;
+}
+
+
+void G_vec_print (VEC *vector, const char* name)
+{
+ int i;
+
+
+ if(vector!=NULL)
+ {
+ G_message ("start vector(%s)", name);
+ //G_message ("Size: %d x %d", vector->dim);
+
+ for( i = 0; i < vector->dim; i++ )
+ {
+ char buf[2048], numbuf[640];
+ sprintf( buf, "%lf ", vector->ve[i]);
+
+ G_message ("%s", buf);
+ }
+
+ G_message ("end vector(%s)",name);
+ }
+
+}
+
+
+
+
+vec_struct *
+G_vector_product (vec_struct *v1, vec_struct *v2)
+{
+ int idx1, idx2, idx0;
+ int i;
+
+
+ vec_struct *out = G_vector_init (v1->rows, v1->ldim, CVEC);
+
+
+
+ //G_warning("Avector1->rows: %d",Avector1->rows);
+ //G_warning("Avector1->cols: %d",Avector1->cols);
+
+ //G_warning("vtmp1->rows1: %d",vtmp1->rows);
+ //G_warning("vtmp1->cols1: %d",out->cols);
+
+ //vtmp1.type = Avector1->type;
+
+ if (!out->is_init) {
+ G_warning (_("Output vector is uninitialized"));
+ return NULL;
+ }
+
+ if (v1->type != v2->type) {
+ G_warning (_("Vectors are not of the same type"));
+ return NULL;
+ }
+
+ if (v1->type != out->type) {
+ G_warning (_("Output vector is of incorrect type"));
+ return NULL;
+ }
+
+ if (v1->type == MATRIX_) {
+ G_warning (_("Matrices not allowed"));
+ return NULL;
+ }
+
+ if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
+ (v1->type == COLVEC_ && v1->rows != v2->rows))
+ {
+ G_warning (_("Vectors have differing dimensions"));
+ return NULL;
+ }
+
+ if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
+ (v1->type == COLVEC_ && v1->rows != out->rows))
+ {
+ G_warning (_("Output vector has incorrect dimension"));
+ return NULL;
+ }
+
+#if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS) //&& defined(HAVE_G2C_H)
+ f77_dhad (v1->cols, 1.0, v1->vals, 1, v2->vals, 1, 0.0, out->vals, 1.0);
+#else
+ idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
+ idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
+ idx0 = (out->v_indx > 0) ? out->v_indx : 0;
+
+ if (v1->type == ROWVEC_) {
+ for (i = 0; i < v1->cols; i++)
+ G_matrix_set_element(out, idx0, i,
+ G_matrix_get_element(v1, idx1, i) *
+ G_matrix_get_element(v2, idx2, i));
+ } else {
+ for (i = 0; i < v1->rows; i++)
+ G_matrix_set_element(out, i, idx0,
+ G_matrix_get_element(v1, i, idx1) *
+ G_matrix_get_element(v2, i, idx2));
+ }
+#endif
+
+ return out;
+}
+
+
+
+
+int
+G_matrix_read2(FILE *fp, mat_struct *out)
+{
+ char buff[100];
+ int rows, cols;
+ int i, j, row;
+ double val;
+
+ /* skip comments */
+ for (;;) {
+ if (!G_getl(buff, sizeof(buff), fp))
+ return -1;
+ if (buff[0] != '#')
+ break;
+ }
+
+ if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) {
+ G_warning(_("Input format error1"));
+ return -1;
+ }
+
+
+ G_matrix_set(out, rows, cols, rows);
+
+
+ //G_warning("row: %d",row);
+
+ for (i = 0; i < rows; i++)
+ {
+ if (fscanf(fp, "row%d:", &row) != 1)
+ {
+ G_warning(_("Input format error"));
+ return -1;
+ }
+
+ for (j = 0; j < cols; j++)
+ {
+ if (fscanf(fp, "%lf:", &val) != 1)
+ {
+ G_warning(_("Input format error"));
+ return -1;
+ }
+
+ fgetc(fp);
+ G_matrix_set_element(out, i, j, val);
+ }
+ }
+//G_matrix_print(out);
+ return 0;
+}
+
+
Added: grass-addons/grass6/imagery/i.spec.unmix/la_extra.h
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/la_extra.h (rev 0)
+++ grass-addons/grass6/imagery/i.spec.unmix/la_extra.h 2012-12-04 10:17:41 UTC (rev 54176)
@@ -0,0 +1,31 @@
+#include <grass/config.h>
+#include <stdio.h>
+#include <grass/blas.h>
+#include <grass/lapack.h>
+typedef struct VEC_ {
+ int dim, max_dim;
+ double *ve;
+ } VEC;
+
+int
+G_matrix_read2(FILE *fp, mat_struct *out);
+void
+G_matrix_print2 (mat_struct *mt, const char* name);
+vec_struct *
+G_matvect_get_column2(mat_struct *mt, int col);
+
+
+
+
+VEC *G_vec_get(int size);
+VEC* G_vec_get2(int size, VEC *vector);
+void G_vec_print (VEC *vector, const char* name);
+VEC *G_vec_resize(VEC * in, int size);
+
+VEC *v_sub(VEC *vec1,VEC *vec2,VEC *out);
+VEC *mv_mlt(mat_struct *A, VEC *b, VEC *out);
+mat_struct *sm_mlt(double scalar,mat_struct *matrix, mat_struct *out);
+mat_struct *G_matrix_resize(mat_struct * in, int rows, int cols);
+double v_norm2(VEC *x);
+VEC *G_vec_copy(VEC *in);
+
Modified: grass-addons/grass6/imagery/i.spec.unmix/main.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/main.c 2012-12-04 08:11:17 UTC (rev 54175)
+++ grass-addons/grass6/imagery/i.spec.unmix/main.c 2012-12-04 10:17:41 UTC (rev 54176)
@@ -1,149 +1,117 @@
-/* Spectral mixture analysis of satellite/aerial images
- * (c) 1999-2000 Markus Neteler, Hannover, Germany
+/****************************************************************************
*
- * COPYRIGHT: (C) 2008 by the GRASS Development Team, Markus Neteler
- * neteler cealp.it
+ * MODULE: i.spec.unmix
*
- * This program is free software under the GNU General Public
- * License (>=v2). Read the file COPYING that comes with GRASS
- * for details.
+ * AUTHOR(S): Markus Neteler <neteler itc.it>
+ * Mohammed Rashad <rashadkm gmail.com>
*
- * VERSION: Based on LAPACK/BLAS
+ * PURPOSE: Spectral mixture analysis of satellite/aerial images
*
- * Module calculates
- * -1
- * x =A b
+ * COPYRIGHT: (C) 2006-2012 by the GRASS Development Team
*
- * A: matrix of reference spectra (endmembers)
- * b: pixel vector from satellite image
- * x: unknown fraction vector
+ * This program is free software under the GNU General
+ * Public License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
*
- * with two constraints:
- * 1. SUM x_i = 1 (max. 100%) -> least square problem
- * 2. 0 <= x_i <= 1 (no negative fractions) -> Steepest Descend of
- * error surface
- *
- *
- * IMPORTANT: Physical units of matrix file and image data set must fit
- * otherwise the algorithm might run into endless loop!
- *
- * TODO: complete LAPACK/BLAS port. Check with Brad Douglas.
- ********************************************************************/
+ *****************************************************************************/
#define GLOBAL
+
+#include <grass/config.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <math.h>
-#include <grass/gis.h>
-#include <grass/la.h>
+#include <grass/imagery.h>
+#include <grass/gmath.h>
+#include <grass/glocale.h>
#include "global.h"
-#define GAMMA 10 /* last row value in Matrix and last b vector element */
- /* for constraint Sum xi = 1 (GAMMA=weight) */
+#include "la_extra.h"
-int open_files();
-void spectral_angle();
-double find_max(double x, double y)
+#define GAMMA 10 /* last row value in Matrix and last b vector element
+ * for constraint Sum xi = 1 (GAMMA=weight) */
+
+
+static double find_max (double x, double y)
{
- return((x*(x>y))+(y*(x<=y)));
+ return ((x * (x > y)) + (y * (x <= y)));
}
-CELL myround (x)
- double x;
- {
- CELL n;
-
- if (x >= 0.0)
- n = x + .5;
- else
- {
- n = -x + .5;
- n = -n;
- }
- return n;
- }
-
-
-int main(int argc, char *argv[])
+int main (int argc, char *argv[])
{
+ char result_name[80];
int nrows, ncols;
- int row, col;
- int band;
- int i, j, k, iterations;
- mat_struct *A_tilde;
- vec_struct *b_gamma;
- vec_struct *startvector, *A_times_startvector, *errorvector, *temp;
- mat_struct *A_tilde_trans_mu;
- int *index;
- double max1, max2, max_total=0.0;
- double change, mu, deviation;
+ int row;
+ int i, j, k;
+ VEC *b, *b_gamma;
+ VEC *startvector, *A_times_startvector, *errorvector, *temp;
+ mat_struct *A, *A_tilde, *A_tilde_trans_mu, *A_tilde_trans;
+
+ mat_struct B, B_tilde, B_tilde_trans_mu;
+
+ struct GModule *module;
+
+ double max_total = 0.0;
+ double mu;
float anglefield[255][255];
- double error;
- char command[80];
- struct
- {
+ double error = 0.0;
+ struct {
struct Option *group, *matrixfile, *result, *error, *iter;
} parm;
+ /* initialize GIS engine */
G_gisinit (argv[0]);
+
+ module = G_define_module();
+
+ module->keywords = _("imagery, spectral unmixing");
+ module->description =
+ _("Perfroms Spectral mixture analysis of satellite/aerial images");
- parm.group = G_define_option();
- parm.group->key = "group";
- parm.group->type = TYPE_STRING;
- parm.group->required = YES;
- parm.group->description = "Imagery group to be analyzed with Spectral Mixture Analyis";
+ parm.group = G_define_standard_option (G_OPT_I_GROUP);
parm.matrixfile = G_define_option();
- parm.matrixfile->key = "matrixfile";
- parm.matrixfile->type = TYPE_STRING;
+ parm.matrixfile->key = "matrix";
+ parm.matrixfile->type = TYPE_STRING;
parm.matrixfile->required = YES;
- parm.matrixfile->description = "Matrix file containing spectral signatures ";
+ parm.matrixfile->gisprompt = "old_file,file";
+ parm.matrixfile->label = _("Open Matrix file");
+ parm.matrixfile->description = _("Matrix file containing spectral signatures");
+
+
- parm.result = G_define_option();
- parm.result->key = "result";
- parm.result->type = TYPE_STRING;
- parm.result->required = YES;
- parm.result->description = "Raster map prefix to hold spectral unmixing results";
+
- parm.error = G_define_option();
- parm.error->key = "error";
- parm.error->type = TYPE_STRING;
- parm.error->required = YES;
- parm.error->description = "Raster map to hold unmixing error";
+ parm.result = G_define_option ();
+ parm.result->key = "result";
+ parm.result->description = _("Name of raster map prefix to hold spectral unmixing results");
+ parm.result->guisection = _("Required");
-
- parm.iter = G_define_option();
- parm.iter->key = "iter";
- parm.iter->type = TYPE_STRING;
- parm.iter->required = YES;
- parm.iter->description = "Raster map to hold number of iterations";
+ parm.error = G_define_standard_option (G_OPT_R_OUTPUT);
+ parm.error->key = "error";
+ parm.error->description = _("Name of raster map to hold unmixing error");
- flag.quiet = G_define_flag();
- flag.quiet->key = 'q';
- flag.quiet->description = "Run quietly (but with percentage output)";
+ parm.iter = G_define_standard_option (G_OPT_R_OUTPUT);
+ parm.iter->key = "iter";
+ parm.iter->description = _("Raster map to hold number of iterations");
- flag2.veryquiet = G_define_flag();
- flag2.veryquiet->key = 's';
- flag2.veryquiet->description = "Run silently (say nothing)";
+ if (G_parser (argc, argv))
+
+ exit (EXIT_FAILURE);
+
- if (G_parser(argc,argv))
- exit(EXIT_FAILURE);
+ /* here we go... A is created here */
+ A = open_files2 (parm.matrixfile->answer,
+ parm.group->answer,
+ parm.result->answer,
+ parm.iter->answer,
+ parm.error->answer);
+ //G_warning("printing A******");
+ //G_matrix_print2(A,"A");
- result_prefix = parm.result->answer;
- error_name = parm.error->answer;
- iter_name = parm.iter->answer;
- group = parm.group->answer;
- matrixfile = parm.matrixfile->answer;
- if (flag2.veryquiet->answer)
- flag.quiet->answer = 1;
-
-/* here we go... */
-
- open_files(); /*A is created here */
-
/* ATTENTION: Internally we work here with col-oriented matrixfile,
* but the user has to enter the spectra row-wise for his/her's
* convenience... That means: Don't mix row- and col-orientation
@@ -165,246 +133,480 @@
*
* 2. Beside checking matrix orthogonality we find out the maximum
* entry of the matrix for configuring stepsize mu later. */
-
- for (i = 0; i < A->cols; i++) /* go columnwise through matrix*/
+
+ /* go columnwise through matrix */
+
+
+
+
+ for (i = 0; i < A->cols; i++)
+ {
+ vec_struct *Avector1, *Avector2;
+ double max1, max2;
+
+ Avector1 = G_matvect_get_column2 (A, i);
+
+ // G_matrix_print(Avector1);
+
+//G_warning("Avector1: %d", Avector1->rows);
+
+ // get the max. element of this vector
+ max1 = G_vector_norm_maxval (Avector1, 1);
+
+double temp = max1;
+
+ for (j = 0; j < A->cols; j++)
{
- Avector1 = G_matvect_get_column(A, i);
- max1 = G_vector_norm_maxval(Avector1, index); /* get the max. element of this vector */
- for (j = 0; j < A->cols ; j++)
- {
- if (j !=i)
+ if (j != i)
{
- Avector2 = G_matvect_get_column(A, j); /* get next col in A */
- max2 = G_vector_norm_maxval(Avector2, index); /* get the max. element of this vector */
- max_total=(find_max(max1,max2),max_total); /* find max of matrix A */
+ // get next col in A
+ Avector2 = G_matvect_get_column (A, j);
+//G_matrix_print(Avector2);
+ // get the max. element of this vector
+ max2 = G_vector_norm_maxval (Avector2, 1);
+
+//G_warning("max2: %lf", max2);
- spectral_angle(); /* check vector angle */
- anglefield[i][j]= curr_angle; /* save angle in degree */
- G_vector_free(Avector2);
- }
- }
- G_vector_free(Avector1);
- G_vector_free(Avector2);
+
+
+if(max2 > max1)
+temp = max2;
+
+//G_warning("temp: %lf", temp);
+
+if(temp > max_total)
+max_total = temp;
+ // find max of matrix A
+// max_total = (find_max (max1, max2), max_total);
+//G_warning("max_total: %lf", max_total);
+ // save angle in degree
+ anglefield[i][j] = spectral_angle (Avector1, Avector2);
+ //G_warning("anglefield[i][j]: %lf", anglefield[i][j]);
+
+ }
}
- /* print out the result */
- if (!flag.quiet->answer)
- {
- fprintf(stderr,"Checking linear dependencies (orthogonality check) of Matrix A:\n");
- for (i = 0; i < A->cols ; i++)
- for (j = 0; j < A->cols ; j++)
- {
- if (j !=i)
- fprintf(stderr," Angle between row %i and row %i: %g degree\n", \
- (i+1), (j+1), anglefield[i][j]);
- /* internally this is col and not row certainly */
- }
- fprintf(stderr,"\n");
- }
- /* check it */
- error=0;
- for (i = 0; i < A->cols ; i++)
- for (j = 0; j < A->cols ; j++)
- if (j !=i)
- if (anglefield[i][j] < 8.0)
- {
- fprintf(stderr,"ERROR: Spectral entries row%i: and row%i: in your matrix \
- are linear dependent!\n",i,j);
- fprintf(stderr," You have to revise your reference spectra.\n");
- error=1;
- exit(EXIT_FAILURE);
- }
+ G_vector_free (Avector1);
+ G_vector_free (Avector2);
+ }
- if (!error)
- if (!flag.quiet->answer)
- fprintf(stderr,"Spectral matrix is o.k. Proceeding...\n\n");
+ G_message (_("Checking linear dependencies (orthogonality check) of Matrix A..."));
+ for (i = 0; i < A->cols; i++)
+ for (j = 0; j < A->cols; j++)
+ if (j != i)
+ // internally this is col and not row certainly
+ G_message (_("Angle between row %i and row %i: %g degree"), (i + 1), (j + 1), anglefield[i][j]);
-/* Begin calculations */
- /* 1. contraint SUM xi = 1
- * add last row "1" elements to Matrix A, store in A_tilde
- * A_tilde is one row-dimension more than A */
- A_tilde = m_get(A->rows+1, A->cols); /* memory allocation */
+ for (i = 0; i < A->cols ; i++)
+ for (j = 0; j < A->cols ; j++)
+ if (j != i)
+ if (anglefield[i][j] < 8.0)
+ G_fatal_error (_("Spectral entries row %i: and row %i: in "
+ "your matrix are linear dependent!\nYou "
+ "have to revise your reference spectra."),
+ i, j);
- for (i=0; i < A->rows ; i++) /* copy rowwise */
- for (j=0; j < A->cols; j++) /* copy colwise */
- G_matrix_get_element(A_tilde, i, j)= A->G_matrix_get_element[i][j];
+ if (!error)
+ G_message (_("Spectral matrix is o.k. Proceeding..."));
- /* fill last row with 1 elements */
- for (j=0; j < A->cols ; j++)
- G_matrix_get_element(A_tilde, [A->rows][j])= GAMMA;
- /* now we have an overdetermined (non-square) system */
+ // Begin calculations
+ // 1. contraint SUM xi = 1
+ // add last row "1" elements to Matrix A, store in A_tilde
+ // A_tilde is one row-dimension more than A
+
+
-/* We have a least square problem here: error minimization
- * T -1 T
- * unknown fraction = [A_tilde * A_tilde] * A_tilde * b
- *
- * A_tilde is the non-square matrix with first constraint in last row.
- * b is pixel vector from satellite image
- *
- * Solve this by deriving above equation and searching the
- * minimum of this error function in an iterative loop within
- * both constraints.
- */
- /* calculate the transpose of A_tilde*/
- A_tilde_trans = G_matrix_transpose(A_tilde);
+ // memory allocation
+ A_tilde = G_matrix_init (A->rows + 1, A->cols, A->rows+1);
+ if (A_tilde == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
-/* initialize some values */
- /* step size must be small enough for covergence of iteration:
- * mu=0.000001; step size for spectra in range of W/m^2/um
- * mu=0.000000001; step size for spectra in range of mW/m^2/um
- * mu=0.000001; step size for spectra in range of reflectance
- ***/
- /* check max_total for number of digits to configure mu size*/
- mu=0.0001 * pow(10,-1*ceil(log10(max_total)));
- startvector = v_get(A->cols); /* length: no. of spectra */
- A_times_startvector = v_get(A_tilde->m); /* length: no. of bands */
- errorvector = v_get(A_tilde->m); /* length: no. of bands */
- temp = v_get(A_tilde->n); /* length: no. of spectra */
- A_tilde_trans_mu = m_get(A_tilde->m,A_tilde->n);
-/* Now we can calculated the fractions pixelwise */
- nrows = G_window_rows(); /* get geographical region */
- ncols = G_window_cols();
+ //G_message("rrrr%d", A_tilde->rows);
+
+
+ for (i = 0; i < A->rows ; i++)
+ for (j = 0; j < A->cols; j++)
+ G_matrix_set_element (A_tilde, i, j, G_matrix_get_element (A, i, j));
+
+
+
+ // fill last row with 1 elements
+
+ for (j = 0; j < A_tilde->cols; j++)
+ {
+ //G_message("Row: %d, Col:%d",i,j);
+ G_matrix_set_element (A_tilde, i, j, GAMMA);
+ }
+//G_matrix_print2(A_tilde, "A_tilde");
+
+
+ // now we have an overdetermined (non-square) system
+
+ // We have a least square problem here: error minimization
+ // T -1 T
+ // unknown fraction = [A_tilde * A_tilde] * A_tilde * b
+ //
+ // A_tilde is the non-square matrix with first constraint in last row.
+ // b is pixel vector from satellite image
+ //
+ // Solve this by deriving above equation and searching the
+ // minimum of this error function in an iterative loop within
+ // both constraints.
+
+
+ // calculate the transpose of A_tilde
+//G_matrix_print(A_tilde);
- if (!flag.quiet->answer)
- {
- fprintf (stderr, "Calculating for %i x %i pixels (%i bands) = %i pixelvectors:\n",\
- nrows,ncols, Ref.nfiles, (ncols * ncols));
- fprintf (stderr, "%s ... ", G_program_name());
- }
- for (row = 0; row < nrows; row++) /* rows loop in images */
+ A_tilde_trans = G_matrix_transpose (A_tilde);
+//G_matrix_print2(A_tilde_trans, "A_tilde_trans");
+
+
+ // initialize some values
+
+ // step size must be small enough for covergence of iteration:
+ // mu = 0.000001; step size for spectra in range of W/m^2/um
+ // mu = 0.000000001; step size for spectra in range of mW/m^2/um
+ // mu = 0.000001; step size for spectra in range of reflectance
+ //
+
+ // check max_total for number of digits to configure mu size
+ mu = 0.0001 * pow (10, -1 * ceil (log10 (max_total)));
+
+
+ //G_message("mu = %lf", mu);
+startvector = G_vec_get2(A->cols,startvector);
+
+
+
+ if (startvector == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+
+
+ A_times_startvector = G_vec_get2(A_tilde->rows,A_times_startvector); // length: no. of bands //
+ errorvector = G_vec_get2(A_tilde->rows,errorvector); // length: no. of bands //
+ temp = G_vec_get2(A_tilde->cols,temp); // length: no. of spectra //
+ // A_tilde_trans_mu = m_get(A_tilde->m,A_tilde->n);
+
+
+
+
+ // length: no. of bands
+
+ if (A_times_startvector == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+ // length: no. of bands
+
+ if (errorvector == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+ // length: no. of spectra
+
+ if (temp == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+ A_tilde_trans_mu = G_matrix_init (A_tilde->rows, A_tilde->cols, A_tilde->rows);
+ if (A_tilde_trans_mu == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
+
+ // Now we can calculated the fractions pixelwise
+ nrows = G_window_rows (); // get geographical region
+ ncols = G_window_cols ();
+
+ G_message (_("Calculating for %i x %i pixels (%i bands) = %i pixelvectors."),
+ nrows, ncols, Ref.nfiles, (ncols * ncols));
+
+
+ //G_vec_print(A_times_startvector, "A_times_startvectorq");
+
+ for (row = 0; row < nrows; row++)
{
- if (!flag2.veryquiet->answer)
- G_percent(row, nrows, 1);
- for (band = 0; band < Ref.nfiles; band++) /* get one row for all bands*/
- {
- if (G_get_map_row (cellfd[band], cell[band], row) < 0)
- exit(EXIT_FAILURE);
- }
+ int col, band;
- for (col = 0; col < ncols; col++)
- /* cols loop, work pixelwise for all bands */
- {
+ G_percent (row, nrows, 1);
- /* get pixel values of each band and store in b vector: */
- b_gamma = v_get(A_tilde->m); /* length: no. of bands + 1 (GAMMA)*/
- for (band = 0; band < Ref.nfiles; band++)
- b_gamma->ve[band] = cell[band][col];
- /* add GAMMA for 1. constraint as last element*/
- b_gamma->ve[Ref.nfiles] = GAMMA;
+ // get one row for all bands
+ for (band = 0; band < Ref.nfiles; band++)
+ if (G_get_map_row (cellfd[band], cell[band], row) < 0)
+ G_fatal_error (_("Unable to get map row [%d]"), row);
+
+ // for (band = 0; band < Ref.nfiles; band++)
+ // {
+ // if (G_get_map_row (cellfd[band], cell[band], row) < 0)
+ // G_fatal_error (_("Unable to get map row [%d]"), row);
+
+ //G_message("row: %d, nrows: %d", row, nrows);
+ for (col = 0; col < ncols; col++)
+ {
+
+
+
+ double change = 1000;
+ double deviation = 1000;
+ int iterations = 0;
- /* calculate fraction vector for current pixel
- * Result is stored in fractions vector
- * with second constraint: Sum x_i = 1 */
- change=1000; /* initialize */
- deviation=1000;
- iterations = 0;
- for (k = 0; k < (A_tilde->n); k++) /* no. of spectra times */
- startvector->ve[k] = (1./A_tilde->n);
+ // get pixel values of each band and store in b vector:
+ // length: no. of bands + 1 (GAMMA)
+
+ b_gamma = G_vec_get2(A_tilde->rows,b_gamma);
+
+// b_gamma = G_vector_init (A_tilde->rows, 1, RVEC);
+ if (b_gamma == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
- /* get start vector and initialize it with equal fractions:
- * using the neighbor pixelvector as startvector*/
- /* solve with iterative solution: */
- while( fabs(change) > 0.0001 )
- {
- /* go a small step into direction of negative gradient
- * errorvector = A_tilde * startvector - b_gamma */
- mv_mlt(A_tilde, startvector, A_times_startvector);
- v_sub(A_times_startvector, b_gamma, errorvector);
- sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
- mv_mlt(A_tilde_trans_mu, errorvector, temp);
- v_sub(startvector,temp,startvector); /* update startvector */
+ //G_message("%d", A_tilde->rows);
+ for (band = 0; band < Ref.nfiles; band++)
+ {
+ b_gamma->ve[band] = cell[band][col];
+ //G_message("band: %d col: %d", band,col);
+ //G_matrix_set_element (b_gamma, 0, band, cell[band][col]);
+ }
+
- /* if one element gets negative, set it to zero */
- for (k = 0; k < (A_tilde->n); k++) /* no. of spectra times */
+ // add GAMMA for 1. constraint as last element
+ b_gamma->ve[Ref.nfiles] = GAMMA;
+/// G_matrix_set_element (b_gamma, 0, Ref.nfiles, GAMMA);
+
+
+//G_matrix_print(b_gamma,"b_gamma");
+
+ for (k = 0; k < A_tilde->cols; k++)
+ startvector->ve[k] = (1.0 / A_tilde->cols);
+ // G_matrix_set_element (startvector, k,0, (1.0 / A_tilde->cols));
+//G_matrix_print(startvector,"startvector1");
+
+
+//G_matrix_print(b_gamma, "b_gama");
+ // calculate fraction vector for current pixel
+ // Result is stored in fractions vector
+ // with second constraint: Sum x_i = 1
+
+
+
+
+//G_vec_print(startvector, "startvector");
+ // get start vector and initialize it with equal fractions:
+ // using the neighbor pixelvector as startvector
+
+
+ // solve with iterative solution:
+ while (fabs (change) > 0.0001)
+ {
+
+
+
+ A_times_startvector = mv_mlt(A_tilde, startvector, A_times_startvector);
+
+//G_vec_print(A_times_startvector, "xx");
+ errorvector = v_sub(A_times_startvector, b_gamma, errorvector);
+ //G_vec_print(A_times_startvector, "errorvector");
+
+ A_tilde_trans_mu = sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
+
+
+//G_matrix_print(A_tilde_trans_mu,"A_tilde_trans_mu");
+
+
+ temp = mv_mlt(A_tilde_trans_mu, errorvector, temp);
+ startvector = v_sub(startvector,temp,startvector); // update startvector //
+
+ // if one element gets negative, set it to zero //
+ for (k = 0; k < (A_tilde->cols); k++) // no. of spectra times //
if (startvector->ve[k] < 0)
startvector->ve[k] = 0;
- /* Check the deviation */
- change = deviation - G_vector_norm_euclid(errorvector);
- deviation = G_vector_norm_euclid(errorvector);
+ // Check the deviation //
+ double norm2 = v_norm2(errorvector);
+ change = deviation - norm2;
+ deviation = norm2;
+
+ iterations++;
+
+ // if(fabs (change) > 0.0001)
+ //G_message("change=%lf, deviation=%lf",change, 0.0001);
+
+ /********************/
+
+
+ // go a small step into direction of negative gradient
+ // errorvector = A_tilde * startvector - b_gamma
+ //
+// mv_mlt(A_tilde, startvector, A_times_startvector);
- /* debug output */
- /*fprintf(stderr, "Change: %g - deviation: %g\n", \
- *change, deviation); */
+/*
- iterations++;
- } /* while */
+//G_matrix_print(A_times_startvector,"A_times_startvector");
+//G_matrix_print(b_gamma, "b_gamma");
+//G_matrix_print(errorvector, "errorvector");
+ G_vector_sub (A_times_startvector, b_gamma, errorvector);
+
+
+// sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
+// mv_mlt(A_tilde_trans_mu, errorvector, temp);
+ G_vector_sub (startvector, temp, startvector);
- fraction=v_get(A->cols); /* length: no. of spectra */
- error = deviation / G_vector_norm_euclid(b_gamma);
- v_copy(startvector, fraction);
+ // if one element gets negative, set it to zero
+ for (k = 0; k < A_tilde->cols; k++) // no. of spectra times
+ {
+// if (startvector->ve[k] < 0)
+// startvector->ve[k] = 0;
- /*---------- end of second contraint -----------------------*/
- /* store fractions in resulting rows of resulting files
- * (number of bands = vector dimension) */
+//G_message("get_element: %lf", G_matrix_get_element (startvector, startvector->cols, k));
- /* write result in full percent */
- for (i = 0; i < A->cols; i++) /* no. of spectra */
+ if ((G_matrix_get_element (startvector, startvector->cols, k) < 0))
+ {
+ G_matrix_set_element (startvector, startvector->cols, k, 0);
+ //G_message("A_tilde->cols: %d", A_tilde->cols); //4
+ }
+ }
+ // Check the deviation
+ double norm_euclid = G_vector_norm_euclid (errorvector);
+
+ // G_message("norm_euclid : %lf", norm_euclid );
+ change = deviation - norm_euclid;
+ deviation = G_vector_norm_euclid (errorvector);
+
+ // G_message ("Change: %g - deviation: %g", change, deviation);
+
+ G_debug (5, "Change: %g - deviation: %g",
+ change, deviation);
+*/
+
+
+ }
+
+ // if(fabs (change) > 0.0001)
+
+
+
+
+ VEC *fraction;
+//G_message("fcol %d and A->cols %d", startvector->dim, A->cols);
+ fraction=G_vec_get(A->cols); // length: no. of spectra //
+ error = deviation / v_norm2(b_gamma);
+ fraction = G_vec_copy (startvector);
+
+
+
+ // write result in full percent //
+ for (i = 0; i < A->cols; i++) // no. of spectra //
result_cell[i][col] = (CELL)(100 * fraction->ve[i]);
- /* save error and iterations*/
+ // save error and iterations//
error_cell[col] = (CELL) (100 * error);
iter_cell[col] = iterations;
- G_vector_free(fraction);
- G_vector_free(b);
- } /* columns loop */
+ //V_FREE(fraction);
+ //V_FREE(b);
+
+ // }
+
- /* write the resulting rows into output files: */
- for (i = 0; i < A->cols; i++) /* no. of spectra */
- G_put_map_row (resultfd[i], result_cell[i]);
- if (error_fd > 0)
- G_put_map_row (error_fd, error_cell);
- if (iter_fd > 0)
- G_put_map_row (iter_fd, iter_cell);
- } /* rows loop */
- if (!flag2.veryquiet->answer)
- G_percent(row, nrows, 2);
-
- /* close files */
- for (i = 0; i < Ref.nfiles; i++) /* no. of bands */
- G_unopen_cell(cellfd[i]);
-
- for (i = 0; i < A->cols; i++) /* no. of spectra */
- {
- G_close_cell (resultfd[i]);
- /* make grey scale color table */
- sprintf(result_name, "%s.%d", result_prefix, (i+1));
- sprintf(command, "r.colors map=%s color=rules 2>&1> /dev/null <<EOF\n
- 0 0 0 0 \n
- 100 0 255 0\n
- end\n
- EOF", result_name);
- system(command);
- /* create histogram */
- do_histogram(result_name, Ref.file[i].mapset);
- }
+ //---------- end of second contraint -----------------------
+ // store fractions in resulting rows of resulting files
+ // (number of bands = vector dimension)
+
+ // write result in full percent
+ //G_matrix_print(fraction,"fraction"); //same as startvector
+
+ //G_message ("fraction->rows: %d",fraction->rows);
+ //G_message ("i=%d, col=%d",i,col);
+
+/*
+
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ double dd = G_matrix_get_element (fraction, i, 0);
+ dd = 100*dd;
+ //G_message ("i=%d, col=%d",i,col);
+ result_cell[i][col] = (CELL) dd;
+ //result_cell[i][col] = (CELL)(100 * G_matrix_get_element (fraction, fraction->rows-1, i));
+ }
+
+ // save error and iterations
+ error_cell[col] = (CELL) (100 * error);
+ iter_cell[col] = iterations;
+
+ G_vector_free (fraction);
+ // G_vector_free (b);
+
+*/
+ } //end cols loop
+
+ //G_message("finished %d of %d", row,nrows);
+ // }
+ // }
+
+
+
+ // write the resulting rows into output files:
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ G_put_map_row (resultfd[i], result_cell[i]);
+
+ if (error_fd > 0)
+ G_put_map_row (error_fd, error_cell);
+
+ if (iter_fd > 0)
+ G_put_map_row (iter_fd, iter_cell);
+
+ } // rows loop
+ G_percent (row, nrows, 2);
+
+ // close files
+ for (i = 0; i < Ref.nfiles; i++) // no. of bands
+ G_unopen_cell (cellfd[i]);
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ char command[1080];
+
+ G_close_cell (resultfd[i]);
+
+ // make grey scale color table
+ sprintf (result_name, "%s.%d", parm.result->answer, (i+1));
+ sprintf (command, "r.colors map=%s color=rules <<EOF\n"
+ "0 0 0 0 \n"
+ "201 0 255 0\n"
+ "end\n"
+ "EOF", result_name);
+
+
+ //G_message(command);
+ //G_system (command);
+
+ // create histogram
+ do_histogram (result_name, Ref.file[i].mapset);
+ }
+
if (error_fd > 0)
- {
- G_close_cell (error_fd);
- sprintf(command, "r.colors map=%s color=gyr >/dev/null", error_name);
- system(command);
- }
+ {
+ char command[80];
+
+ G_close_cell (error_fd);
+ sprintf (command, "r.colors map=%s color=gyr >/dev/null", parm.error->answer);
+ //G_system (command);
+ }
+
if (iter_fd > 0)
- {
- G_close_cell (iter_fd);
- /* sprintf(command, "r.colors map=%s color=gyr >/dev/null", iter_name);
- system(command);*/
- }
+ G_close_cell (iter_fd);
- G_matrix_free(A);
+ G_matrix_free (A);
- make_history(result_name, group, matrixfile);
- exit(EXIT_SUCCESS);
-} /* main*/
-
+ make_history (result_name, parm.group->answer, parm.matrixfile->answer);
+
+/*********************
+*********************/
+ exit (EXIT_SUCCESS);
+}
Modified: grass-addons/grass6/imagery/i.spec.unmix/open.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/open.c 2012-12-04 08:11:17 UTC (rev 54175)
+++ grass-addons/grass6/imagery/i.spec.unmix/open.c 2012-12-04 10:17:41 UTC (rev 54176)
@@ -1,11 +1,6 @@
/* Spextral unmixing with Singular Value Decomposition */
/* (c) 15. Jan. 1999 Markus Neteler, Hannover*/
-/**************************************************************************
- ** Matrix computations 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.
@@ -13,143 +8,309 @@
ISBN 0 7315 1900 0
*/
-#include "global.h"
#include <stdio.h>
#include <math.h>
-#include <grass/gis.h>
-#include "matrix.h"
-#include "matrix2.h"
+#include <grass/imagery.h>
+#include <grass/gmath.h>
+#include <grass/glocale.h>
+#include "global.h"
-int open_files()
+
+int open_files (char *matrixfile,
+ char *img_grp,
+ char *iter_name,
+ char *error_name,
+ mat_struct *A)
{
- char *name, *mapset;
+ char result_name[80];
+ char *result_prefix="out";
FILE *fp;
- int i;
- MAT *A_input;
+ int i, matrixsize;
+ mat_struct A_input;
+
-/* Read in matrix file with spectral library.
- Input matrix must contain spectra row-wise (for user's convenience)!
- Transposed here to col-wise orientation (for modules/mathematical
- convenience).
- */
+// mat_struct A_input1;
+/*
- fp=fopen(matrixfile,"r");
- if (fp == NULL)
- {
- fprintf(stderr,"ERROR: Matrixfile %s not found.\n",matrixfile);
- exit(EXIT_FAILURE);
- }
- A_input = m_finput(fp, MNULL);
+ if ((fp = fopen (matrixfile, "r")) == NULL)
+ G_fatal_error (_("Matrix file %s not found."), matrixfile);
+
+
+ //G_matrix_init2( &A_input,3,3,3);
+ // G_warning("matrix read start");
+// G_warning("cols=%d",A_input->cols);
+
+ if ((G_matrix_read (fp, &A_input) < 0))
+ G_fatal_error (_("Unable to read matrix file %s."), matrixfile);
fclose (fp);
- if (!flag.quiet->answer)
- {
- fprintf(stderr, "Your spectral matrix = ");
- m_output(A_input);
- }
+ //G_matrix_print(&A_input);
+
+
+ // G_warning("matrix read done");
+#if 0
+ G_message(_("Your spectral matrix = %d"), m_output(A_input));
+#endif
+
+
+
+// A = m_get(A_input->rows, A_input->cols);
+
+ G_matrix_clone2(&A_input, A);
+
+
-/* transpose input matrix from row orientation to col orientation.
- * Don't mix rows and cols in the source code and the modules
- * messages output! */
-
- A=m_get(A_input->m, A_input->n);
- m_transp(A_input, A);
- M_FREE(A_input);
+ //*A = *G_matrix_init (A_input.rows, A_input.cols, A_input.rows);
+// if (A == NULL)
+ // G_fatal_error (_("Unable to allocate memory for matrix"));
+
+ A = G_matrix_transpose (&A_input);
- if ( A->m < A->n )
- {
- fprintf(stderr, "ERROR: Need number of cols >= rows to perform least squares fitting.\n");
- exit(EXIT_FAILURE);
- }
- matrixsize = A->m; /* number of rows must be equivalent to no. of bands */
+ // G_matrix_free (&A_input);
-/* open input files from group */
- if (!I_find_group(group))
- {
- fprintf (stderr, "group=%s - not found\n", group);
- exit(EXIT_FAILURE);
- }
- I_get_group_ref(group, &Ref);
+// G_matrix_print(A);
+
+ if ((A->rows) < (A->cols))
+ G_fatal_error (_("Need number of cols >= rows to perform least squares fitting."));
+
+ // number of rows must be equivalent to no. of bands
+ matrixsize = A->cols;
+
+ // open input files from group
+ if (!I_find_group (img_grp))
+ G_fatal_error (_("Unable to find imagery group %s."), img_grp);
+
+ I_get_group_ref (img_grp, &Ref);
if (Ref.nfiles <= 1)
{
- fprintf (stderr, "ERROR: Group %s\n", group);
if (Ref.nfiles <= 0)
- fprintf (stderr, "doesn't have any files\n");
+ G_fatal_error (_("Group %s does not have any rasters. "
+ "The group must have at least 2 rasters."), img_grp);
else
- fprintf (stderr, "only has 1 file\n");
- fprintf (stderr, "The group must have at least 2 files\n");
- exit(EXIT_FAILURE);
+ G_fatal_error (_("Group %s only has 1 raster. "
+ "The group must have at least 2 rasters."), img_grp);
}
- /* Error check: input file number must be equal to matrix size */
+
+ // Error check: input file number must be equal to matrix size
if (Ref.nfiles != matrixsize)
- {
- fprintf (stderr, "ERROR: Number of %i input files in group <%s>\n", Ref.nfiles, group);
- fprintf (stderr, " does not match no. of spectra in matrix \
- (contains only %i cols).\n", A->m);
- exit(1);
- }
+ G_fatal_error (_("Number of input files (%i) in group <%s> "
+ "does not match number of spectra in matrix. "
+ "(contains only %i cols)."),
+ Ref.nfiles, img_grp, A->cols);
- /* get memory for input files */
+ // get memory for input files
cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
cellfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
- for (i=0; i < Ref.nfiles; i++)
+ for (i = 0; i < Ref.nfiles; i++)
{
- cell[i] = G_allocate_cell_buf();
- name = Ref.file[i].name;
- mapset = Ref.file[i].mapset;
- if (!flag.quiet->answer)
- fprintf (stderr,"Opening input file no. %i [%s]\n", (i+1), name);
- if ((cellfd[i] = G_open_cell_old (name, mapset)) < 0)
- {
- fprintf (stderr, "Unable to proceed\n");
- exit(1);
- }
+ cell[i] = G_allocate_cell_buf ();
+
+ G_message (_("Opening input file no. %i [%s]"), (i + 1), Ref.file[i].name);
+
+ if ((cellfd[i] = G_open_cell_old (Ref.file[i].name, Ref.file[i].mapset)) < 0)
+ G_fatal_error (_("Unable to open <%s>"), Ref.file[i].name);
}
-/* open files for results*/
+ // open files for results
result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
- for (i=0; i < (A->n); i++) /* no. of spectra */
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
{
- sprintf(result_name, "%s.%d", result_prefix, (i+1));
- if (!flag.quiet->answer)
- fprintf (stderr,"Opening output file [%s]\n", result_name);
- result_cell[i] = G_allocate_cell_buf();
- if ((resultfd[i] = G_open_cell_new (result_name)) <0)
- {
- fprintf (stderr, "GRASS-Database internal error: Unable to proceed\n");
- exit(1) ;
- }
+ sprintf (result_name, "%s.%d", result_prefix, (i + 1));
+ G_message (_("Opening output file [%s]"), result_name);
+
+ result_cell[i] = G_allocate_cell_buf ();
+ if ((resultfd[i] = G_open_cell_new (result_name)) < 0)
+ G_fatal_error (_("GRASS-DB internal error: Unable to proceed."));
}
+ // open file containing SMA error
+ error_cell = (CELL *) G_malloc (sizeof(CELL *));
+ if (error_name)
+ {
+ G_message (_("Opening error file [%s]"), error_name);
-/* open file containing SMA error*/
+ if ((error_fd = G_open_cell_new (error_name)) < 0)
+ G_fatal_error (_("Unable to create error layer [%s]"), error_name);
+ else
+ error_cell = G_allocate_cell_buf ();
+ }
- error_cell = (CELL *) G_malloc (sizeof (CELL *));
- if (error_name)
+ // open file containing number of iterations
+ iter_cell = (CELL *) G_malloc (sizeof(CELL *));
+ if (iter_name)
{
- error_fd = G_open_cell_new (error_name);
- if (!flag.quiet->answer)
- fprintf (stderr,"Opening error file [%s]\n", error_name);
- if (error_fd < 0)
- fprintf (stderr, "Unable to create error layer [%s]", error_name);
+ G_message (_("Opening iteration file [%s]"), iter_name);
+
+ if ((iter_fd = G_open_cell_new (iter_name)) < 0)
+ G_fatal_error (_("Unable to create iterations layer [%s]"), iter_name);
else
- error_cell = G_allocate_cell_buf();
+ iter_cell = G_allocate_cell_buf ();
}
+
+ //G_matrix_print(A);
+
+ return matrixsize;
+ */
+ return 0;
+}
+
+
+
+
+
+
+
+
+mat_struct* open_files2 (char *matrixfile,
+ char *img_grp,
+ char *result_prefix,
+ char *iter_name,
+ char *error_name)
+{
+ char result_name[80];
+
+ FILE *fp;
+ int i, matrixsize;
+ mat_struct A_input, *A;
+
+
+// mat_struct A_input1;
+
+
+ /* Read in matrix file with spectral library.
+ * Input matrix must contain spectra row-wise (for user's convenience)!
+ * Transposed here to col-wise orientation (for modules/mathematical
+ * convenience). */
+
+ if ((fp = fopen (matrixfile, "r")) == NULL)
+ G_fatal_error (_("Matrix file %s not found."), matrixfile);
+
+
+ //G_matrix_init2( &A_input,3,3,3);
+ // G_warning("matrix read start");
+// G_warning("cols=%d",A_input->cols);
+ /* Read data and close file */
+ if ((G_matrix_read2 (fp, &A_input) < 0))
+ G_fatal_error (_("Unable to read matrix file %s."), matrixfile);
+ fclose (fp);
+
+ //G_matrix_print2(&A_input, "A_input");
+
+
+ // G_warning("matrix read done");
+#if 0
+ G_message(_("Your spectral matrix = %d"), m_output(A_input));
+#endif
+
+ /* transpose input matrix from row orientation to col orientation.
+ * Don't mix rows and cols in the source code and the modules
+ * messages output! */
+
+// A = m_get(A_input->rows, A_input->cols);
+
+ //G_matrix_clone2(&A_input, A);
+
+
+
+ A = G_matrix_init (A_input.rows, A_input.cols, A_input.rows);
+ if (A == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
+
+ A = G_matrix_transpose (&A_input);
+
+ // G_matrix_free (&A_input);
+
+// G_matrix_print(A);
+
+ if ((A->rows) < (A->cols))
+ G_fatal_error (_("Need number of cols >= rows to perform least squares fitting."));
+
+ // number of rows must be equivalent to no. of bands
+ matrixsize = A->rows;
+
+ // open input files from group
+ if (!I_find_group (img_grp))
+ G_fatal_error (_("Unable to find imagery group %s."), img_grp);
+
+ I_get_group_ref (img_grp, &Ref);
+ if (Ref.nfiles <= 1)
+ {
+ if (Ref.nfiles <= 0)
+ G_fatal_error (_("Group %s does not have any rasters. "
+ "The group must have at least 2 rasters."), img_grp);
+ else
+ G_fatal_error (_("Group %s only has 1 raster. "
+ "The group must have at least 2 rasters."), img_grp);
+ }
+
+ // Error check: input file number must be equal to matrix size
+ if (Ref.nfiles != matrixsize)
+ G_fatal_error (_("Number of input files (%i) in group <%s> "
+ "does not match number of spectra in matrix. "
+ "(contains %i cols)."),
+ Ref.nfiles, img_grp, A->rows);
+
+ // get memory for input files
+ cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ cellfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+ for (i = 0; i < Ref.nfiles; i++)
+ {
+ cell[i] = G_allocate_cell_buf ();
+
+ G_message (_("Opening input file no. %i [%s]"), (i + 1), Ref.file[i].name);
+
+ if ((cellfd[i] = G_open_cell_old (Ref.file[i].name, Ref.file[i].mapset)) < 0)
+ G_fatal_error (_("Unable to open <%s>"), Ref.file[i].name);
+
+
+ }
+
-/* open file containing number of iterations */
- iter_cell = (CELL *) G_malloc (sizeof (CELL *));
+ // open files for results
+ result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ if (result_prefix)
+ {
+ sprintf (result_name, "%s.%d", result_prefix, (i + 1));
+ G_message (_("Opening output file [%s]"), result_name);
+
+ result_cell[i] = G_allocate_cell_buf ();
+ if ((resultfd[i] = G_open_cell_new (result_name)) < 0)
+ G_fatal_error (_("GRASS-DB internal error: Unable to proceed."));
+ }
+ }
+ // open file containing SMA error
+ error_cell = (CELL *) G_malloc (sizeof(CELL *));
+ if (error_name)
+ {
+ G_message (_("Opening error file [%s]"), error_name);
+
+ if ((error_fd = G_open_cell_new (error_name)) < 0)
+ G_fatal_error (_("Unable to create error layer [%s]"), error_name);
+ else
+ error_cell = G_allocate_cell_buf ();
+ }
+
+ // open file containing number of iterations
+ iter_cell = (CELL *) G_malloc (sizeof(CELL *));
if (iter_name)
{
- iter_fd = G_open_cell_new (iter_name);
- if (!flag.quiet->answer)
- fprintf (stderr,"Opening iteration file [%s]\n", iter_name);
- if (iter_fd < 0)
- fprintf (stderr, "Unable to create iterations layer [%s]", iter_name);
+ G_message (_("Opening iteration file [%s]"), iter_name);
+
+ if ((iter_fd = G_open_cell_new (iter_name)) < 0)
+ G_fatal_error (_("Unable to create iterations layer [%s]"), iter_name);
else
- iter_cell = G_allocate_cell_buf();
+ iter_cell = G_allocate_cell_buf ();
}
- return(matrixsize); /* give back number of output files (= Ref.nfiles) */
+
+
+ /* give back number of output files (= Ref.nfiles) */
+ return A;
}
Modified: grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c 2012-12-04 08:11:17 UTC (rev 54175)
+++ grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c 2012-12-04 10:17:41 UTC (rev 54176)
@@ -3,31 +3,15 @@
*
* 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
*****/
-#define GLOBAL
-#define MY_PI 3.141592653589793
-
#include <stdio.h>
#include <math.h>
-#include "matrix.h"
+#include <grass/gmath.h>
#include "global.h"
-void spectral_angle() /* returns spectral angle globally*/
-{
-
-/* input MAT A, VEC Avector1, Avector2
+/* input mat_struct A, vec_struct Avector1, Avector2
* output cur_angle
*
* v_DN * v_reference
@@ -45,21 +29,22 @@
* 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 */
+ /* multiply one A column with second */
+// vtmp1 = G_vector_init (0, 0, RVEC);
+ 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) */
-/* Measure spectral angle*/
+ G_vector_free (vtmp1);
- vtmp1 = v_star(Avector1, Avector2, VNULL); /* multiply one A column with second */
- norm1 = v_norm1(vtmp1); /* calculate 1-norm */
- norm2 = v_norm2(Avector1); /* calculate 2-norm (Euclidean) */
- norm3 = v_norm2(Avector2); /* calculate 2-norm (Euclidean) */
-
- V_FREE(vtmp1);
-
- curr_angle = (acos(norm1/(norm2 * norm3)) * 180/MY_PI); /* Calculate angle */
- /* return in degree globally*/
+ /* Calculate angle and return in degree globally */
+ return (acos (norm1 / (norm2 * norm3)) * 180 / M_PI);
}
-
More information about the grass-commit
mailing list