[GRASS-SVN] r54180 - grass-addons/grass6/imagery/i.spec.unmix
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Dec 4 09:37:01 PST 2012
Author: neteler
Date: 2012-12-04 09:37:01 -0800 (Tue, 04 Dec 2012)
New Revision: 54180
Modified:
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/la_extra.c
grass-addons/grass6/imagery/i.spec.unmix/la_extra.h
grass-addons/grass6/imagery/i.spec.unmix/main.c
grass-addons/grass6/imagery/i.spec.unmix/meminfo.h
grass-addons/grass6/imagery/i.spec.unmix/open.c
grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c
Log:
proper indentation with tools/grass_indent.sh
Modified: grass-addons/grass6/imagery/i.spec.unmix/global.h
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/global.h 2012-12-04 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/global.h 2012-12-04 17:37:01 UTC (rev 54180)
@@ -20,15 +20,15 @@
GLOBAL int *resultfd;
GLOBAL CELL *error_cell;
-GLOBAL int error_fd;
+GLOBAL int error_fd;
GLOBAL CELL *iter_cell;
-GLOBAL int iter_fd;
+GLOBAL int iter_fd;
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);
+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 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/hist.c 2012-12-04 17:37:01 UTC (rev 54180)
@@ -2,14 +2,13 @@
#include <grass/gis.h>
-void make_history (char *name, char *group, char *matrixfile)
+void make_history(char *name, char *group, char *matrixfile)
{
struct History hist;
- if (G_read_history (name, G_mapset(), &hist) >= 0)
- {
- sprintf (hist.datsrc_1, "Group: %s", group);
- sprintf (hist.datsrc_2, "Matrix file: %s", matrixfile);
- G_write_history (name, &hist);
+ if (G_read_history(name, G_mapset(), &hist) >= 0) {
+ sprintf(hist.datsrc_1, "Group: %s", group);
+ sprintf(hist.datsrc_2, "Matrix file: %s", matrixfile);
+ G_write_history(name, &hist);
}
}
Modified: grass-addons/grass6/imagery/i.spec.unmix/histogram.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/histogram.c 2012-12-04 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/histogram.c 2012-12-04 17:37:01 UTC (rev 54180)
@@ -1,7 +1,7 @@
#include <grass/gis.h>
-int do_histogram (char *name, char *mapset)
+int do_histogram(char *name, char *mapset)
{
CELL *cell;
struct Cell_head cellhd;
@@ -9,33 +9,32 @@
int nrows, ncols, row;
int fd;
- if (G_get_cellhd (name, mapset, &cellhd) < 0)
- return 0;
+ if (G_get_cellhd(name, mapset, &cellhd) < 0)
+ return 0;
- G_set_window (&cellhd);
- fd = G_open_cell_old (name, mapset);
+ G_set_window(&cellhd);
+ fd = G_open_cell_old(name, mapset);
if (fd < 0)
- return 0;
+ return 0;
- nrows = G_window_rows ();
- ncols = G_window_cols ();
- cell = G_allocate_cell_buf ();
+ 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_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);
+ G_update_cell_stats(cell, ncols, &statf);
}
- G_close_cell (fd);
- G_free (cell);
+ G_close_cell(fd);
+ G_free(cell);
if (row == nrows)
- G_write_histogram_cs (name, &statf);
+ G_write_histogram_cs(name, &statf);
- G_free_cell_stats (&statf);
+ G_free_cell_stats(&statf);
return 0;
}
Modified: grass-addons/grass6/imagery/i.spec.unmix/la_extra.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/la_extra.c 2012-12-04 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/la_extra.c 2012-12-04 17:37:01 UTC (rev 54180)
@@ -1,6 +1,6 @@
-#include <stdio.h> /* needed here for ifdef/else */
+#include <stdio.h> /* needed here for ifdef/else */
#include <stdlib.h>
#include <string.h>
#include <math.h>
@@ -19,232 +19,229 @@
-vec_struct *
-G_matvect_get_column2(mat_struct *mt, int col)
+vec_struct *G_matvect_get_column2(mat_struct * mt, int col)
{
- int i; /* loop */
- vec_struct *vc1;
+ int i; /* loop */
+ vec_struct *vc1;
- if(col < 0 || col >= mt->cols) {
- G_warning(_("Specified matrix column index is outside range"));
- return NULL;
- }
+ 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 (!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;
- }
+ 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);
-
- }
+ for (i = 0; i < mt->rows; i++) {
+ double dd = G_matrix_get_element(mt, i, col);
- return vc1;
+ //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)
+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);
+ int i, j;
- 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);
- }
+ 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++ ) {
- strcpy(buf, "");
+ for (i = 0; i < mt->rows; i++) {
+ char buf[2048], numbuf[640];
- for( j = 0; j < mt->cols; j++ ) {
+ sprintf(buf, "row%d: ", i);
+ 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, ", ");
+ 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);
}
- G_message ("%s", buf);
- }
+ /*
+ for( i = 0; i < mt->rows; i++ ) {
+ strcpy(buf, "");
- fprintf (stderr, "\n");
-
-*/
+ 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++ )
- {
-
+ mat_struct *matrix;
- //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;
-}
+ matrix = G_matrix_init(rows, cols, rows);
+ int i, j, p, index = 0;
-mat_struct *sm_mlt(double scalar,mat_struct *matrix, mat_struct *out)
+ 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;
+ int m, n, i, j;
- if ( matrix==NULL )
- G_fatal_error("sm_mlt1(error)");
+ int index = 0;
- 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);
+ if (matrix == NULL)
+ G_fatal_error("sm_mlt1(error)");
-m = matrix->rows;
-n = matrix->cols;
- for ( i=0; i<m; i++ )
- {
- //__smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
+ 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");
+ 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);
+ return (out);
}
-VEC *G_vec_copy(VEC *in)
+VEC *G_vec_copy(VEC * in)
{
- int i;
- VEC *out;
- if ( !in )
- G_fatal_error("v_copy(error1)");
-
- int dim = in->dim;
+ int i;
+ VEC *out;
- out = G_vec_get(dim);
-
-
- for( i = 0; i < dim; i++ )
- {
- out->ve[i]= in->ve[i];
- }
+ if (!in)
+ G_fatal_error("v_copy(error1)");
- return (out);
+ 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)
+double v_norm2(VEC * x)
{
- int i, dim;
- double s, sum;
+ int i, dim;
+ double s, sum;
- if ( !x )
- G_fatal_error("v_norm2(error1)");
-
-
- dim = x->dim;
+ if (!x)
+ G_fatal_error("v_norm2(error1)");
- sum = 0.0;
- for ( i = 0; i < dim; i++ )
- {
- sum += x->ve[i] * x->ve[i];
- }
-
-
-
- return sqrt(sum);
+ 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)
+VEC *v_sub(VEC * vec1, VEC * vec2, VEC * out)
{
- /* u_int i, dim; */
- /* Real *out_ve, *vec1_ve, *vec2_ve; */
+ /* u_int i, dim; */
+ /* Real *out_ve, *vec1_ve, *vec2_ve; */
- if ( !vec1 || !vec2 )
+ 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++ )
+ 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];
/************************************************************
@@ -255,78 +252,77 @@
(*out_ve++) = (*vec1_ve++) - (*vec2_ve++);
************************************************************/
- return (out);
+ return (out);
}
-VEC *mv_mlt(mat_struct *A, VEC *b, VEC *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");
+ unsigned int i, m, n, j;
+ double **A_v, *b_v /*, *A_row */ ;
-// 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 )
+ /* 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 )
-
+ //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);
- 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);
-}
+ //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;
+ }
+
+ //G_vec_print(out, "xx");
+
+ return out;
}
@@ -334,32 +330,31 @@
VEC *G_vec_resize(VEC * in, int size)
{
-
- VEC *vector;
-
- vector = (VEC *)G_malloc( sizeof(VEC) );
+ VEC *vector;
-
- 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]);
-
+
+ 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;
-
+ if (size > in->dim)
+ for (j = i; j < size; j++)
+ vector->ve[j] = 0.0;
+
vector->dim = vector->max_dim = size;
- return vector;
+ return vector;
}
@@ -370,151 +365,148 @@
VEC *G_vec_get(int size)
{
-
- VEC *vector;
-
- vector = (VEC *)G_malloc( sizeof(VEC) );
+ VEC *vector;
-
- vector->ve = (double *)G_malloc(size * sizeof(double) );
- int i;
-for( i = 0; i < size; i++ )
- {
-
- vector->ve[i]= 0.0;
+ 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;
+ //G_vec_print(vector, "vector");
+ return vector;
}
-VEC *G_vec_get2(int size, VEC *vector)
+VEC *G_vec_get2(int size, VEC * vector)
{
-
- //VEC *vector;
-
- vector = (VEC *)G_malloc( sizeof(VEC) );
+ //VEC *vector;
-
- vector->ve = (double *)G_malloc(size * sizeof(double) );
- int i;
-for( i = 0; i < size; i++ )
- {
-
- vector->ve[i]= 0.0;
+ 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;
+ //G_vec_print(vector, "vector");
+ return vector;
}
-void G_vec_print (VEC *vector, const char* name)
+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);
+ int i;
- for( i = 0; i < vector->dim; i++ )
- {
- char buf[2048], numbuf[640];
- sprintf( buf, "%lf ", vector->ve[i]);
- G_message ("%s", buf);
+ 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);
}
-
- G_message ("end vector(%s)",name);
- }
}
-vec_struct *
-G_vector_product (vec_struct *v1, vec_struct *v2)
+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);
-
+
+
+ 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;
+ G_warning(_("Output vector is uninitialized"));
+ return NULL;
}
if (v1->type != v2->type) {
- G_warning (_("Vectors are not of the same type"));
- return NULL;
+ 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;
+ G_warning(_("Output vector is of incorrect type"));
+ return NULL;
}
if (v1->type == MATRIX_) {
- G_warning (_("Matrices not allowed"));
- return NULL;
+ 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;
+ (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;
+ (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);
+#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));
+ 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;
@@ -523,55 +515,48 @@
-int
-G_matrix_read2(FILE *fp, mat_struct *out)
+int G_matrix_read2(FILE * fp, mat_struct * out)
{
- char buff[100];
- int rows, cols;
- int i, j, row;
- double val;
+ 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;
- }
+ /* 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;
- }
-
+ 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);
+ G_matrix_set(out, rows, cols, rows);
- 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;
- }
+ //G_warning("row: %d",row);
- fgetc(fp);
- G_matrix_set_element(out, i, j, val);
- }
- }
-//G_matrix_print(out);
- return 0;
-}
+ 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;
+}
Modified: grass-addons/grass6/imagery/i.spec.unmix/la_extra.h
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/la_extra.h 2012-12-04 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/la_extra.h 2012-12-04 17:37:01 UTC (rev 54180)
@@ -2,30 +2,27 @@
#include <stdio.h>
#include <grass/blas.h>
#include <grass/lapack.h>
-typedef struct VEC_ {
- int dim, max_dim;
- double *ve;
- } VEC;
+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);
+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_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);
+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);
-
+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 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/main.c 2012-12-04 17:37:01 UTC (rev 54180)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: i.spec.unmix
@@ -14,7 +15,7 @@
* comes with GRASS for details.
*
*****************************************************************************/
-
+
#define GLOBAL
#include <grass/config.h>
@@ -30,17 +31,17 @@
#include "la_extra.h"
-#define GAMMA 10 /* last row value in Matrix and last b vector element
- * for constraint Sum xi = 1 (GAMMA=weight) */
+#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)
+static double find_max(double x, double y)
{
return ((x * (x > y)) + (y * (x <= y)));
}
-int main (int argc, char *argv[])
+int main(int argc, char *argv[])
{
char result_name[80];
int nrows, ncols;
@@ -49,248 +50,248 @@
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;
-
+
+ struct GModule *module;
+
double max_total = 0.0;
double mu;
float anglefield[255][255];
double error = 0.0;
- struct {
+ 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");
+ G_gisinit(argv[0]);
- parm.group = G_define_standard_option (G_OPT_I_GROUP);
+ module = G_define_module();
+ module->keywords = _("imagery, spectral unmixing");
+ module->description =
+ _("Perfroms Spectral mixture analysis of satellite/aerial images");
+
+ parm.group = G_define_standard_option(G_OPT_I_GROUP);
+
parm.matrixfile = G_define_option();
- parm.matrixfile->key = "matrix";
- parm.matrixfile->type = TYPE_STRING;
+ parm.matrixfile->key = "matrix";
+ parm.matrixfile->type = TYPE_STRING;
parm.matrixfile->required = YES;
parm.matrixfile->gisprompt = "old_file,file";
- parm.matrixfile->label = _("Open Matrix file");
- parm.matrixfile->description = _("Matrix file containing spectral signatures");
-
-
+ 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->description = _("Name of raster map prefix to hold spectral unmixing results");
+
+
+
+ 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.error = G_define_standard_option (G_OPT_R_OUTPUT);
- parm.error->key = "error";
+ 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");
- parm.iter = G_define_standard_option (G_OPT_R_OUTPUT);
- parm.iter->key = "iter";
+ parm.iter = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.iter->key = "iter";
parm.iter->description = _("Raster map to hold number of iterations");
- 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******");
+ 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");
- /* 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
- * in source code and modules messages output!
- *
- * Spectral Matrix is stored in A now (diagonally flipped to input
- * file) Generally: n: cols .... for matrix A
- * m: rows
- * |
- * |
- **/
+ /* 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
+ * in source code and modules messages output!
+ *
+ * Spectral Matrix is stored in A now (diagonally flipped to input
+ * file) Generally: n: cols .... for matrix A
+ * m: rows
+ * |
+ * |
+ **/
- /* 1. 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.
- *
- *
- * 2. Beside checking matrix orthogonality we find out the maximum
- * entry of the matrix for configuring stepsize mu later. */
+ /* 1. 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.
+ *
+ *
+ * 2. Beside checking matrix orthogonality we find out the maximum
+ * entry of the matrix for configuring stepsize mu later. */
/* 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++)
- {
- if (j != i)
- {
- // 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);
+ 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);
-if(max2 > max1)
-temp = max2;
+ //G_warning("Avector1: %d", Avector1->rows);
-//G_warning("temp: %lf", temp);
+ // get the max. element of this vector
+ max1 = G_vector_norm_maxval(Avector1, 1);
-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]);
-
- }
+ double temp = max1;
+
+ for (j = 0; j < A->cols; j++) {
+ if (j != i) {
+ // 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);
+
+
+
+ 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]);
+
+ }
+ }
+
+ G_vector_free(Avector1);
+ G_vector_free(Avector2);
}
- G_vector_free (Avector1);
- G_vector_free (Avector2);
- }
+ G_message(_("Checking linear dependencies (orthogonality check) of Matrix A..."));
- 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]);
- 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]);
+ 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->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);
-
if (!error)
- G_message (_("Spectral matrix is o.k. Proceeding..."));
+ G_message(_("Spectral matrix is o.k. Proceeding..."));
// 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
-
-
+ // add last row "1" elements to Matrix A, store in A_tilde
+ // A_tilde is one row-dimension more than A
+
+
// memory allocation
- A_tilde = G_matrix_init (A->rows + 1, A->cols, A->rows+1);
+ 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"));
+ G_fatal_error(_("Unable to allocate memory for matrix"));
- //G_message("rrrr%d", A_tilde->rows);
+ //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));
+ 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);
+ 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");
+ //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.
-
+ // 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);
-
- A_tilde_trans = G_matrix_transpose (A_tilde);
-//G_matrix_print2(A_tilde_trans, "A_tilde_trans");
+ //G_matrix_print(A_tilde);
+ 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
- //
+ // 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)));
-
-
+ mu = 0.0001 * pow(10, -1 * ceil(log10(max_total)));
+
+
//G_message("mu = %lf", mu);
-startvector = G_vec_get2(A->cols,startvector);
+ startvector = G_vec_get2(A->cols, startvector);
-
+
if (startvector == NULL)
- G_fatal_error (_("Unable to allocate memory for vector"));
+ 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);
+ 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);
@@ -298,315 +299,314 @@
// length: no. of bands
if (A_times_startvector == NULL)
- G_fatal_error (_("Unable to allocate memory for vector"));
+ 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"));
+ 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"));
+ 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);
+ 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"));
+ 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 ();
+ 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");
+ G_message(_("Calculating for %i x %i pixels (%i bands) = %i pixelvectors."),
+ nrows, ncols, Ref.nfiles, (ncols * ncols));
- for (row = 0; row < nrows; row++)
- {
- int col, band;
- G_percent (row, nrows, 1);
+ //G_vec_print(A_times_startvector, "A_times_startvectorq");
- // 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;
+ for (row = 0; row < nrows; row++) {
+ int col, band;
+ G_percent(row, nrows, 1);
- // 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 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("%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]);
- }
-
+ //G_message("row: %d, nrows: %d", row, nrows);
+ for (col = 0; col < ncols; col++) {
- // 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
+ double change = 1000;
+ double deviation = 1000;
+ int iterations = 0;
+ // 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);
-//G_vec_print(startvector, "startvector");
+ // b_gamma = G_vector_init (A_tilde->rows, 1, RVEC);
+ if (b_gamma == NULL)
+ G_fatal_error(_("Unable to allocate memory for matrix"));
+
+
+ //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]);
+ }
+
+
+ // 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
-
+ // 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 //
+ while (fabs(change) > 0.0001) {
- // 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 //
- 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);
-/*
+ A_times_startvector =
+ mv_mlt(A_tilde, startvector, A_times_startvector);
-//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);
+ //G_vec_print(A_times_startvector, "xx");
+ errorvector =
+ v_sub(A_times_startvector, b_gamma, errorvector);
+ //G_vec_print(A_times_startvector, "errorvector");
- // 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;
+ A_tilde_trans_mu =
+ sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
-//G_message("get_element: %lf", G_matrix_get_element (startvector, startvector->cols, k));
- 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_matrix_print(A_tilde_trans_mu,"A_tilde_trans_mu");
- G_debug (5, "Change: %g - deviation: %g",
- change, deviation);
-*/
-
-
- }
-
- // if(fabs (change) > 0.0001)
-
+ temp = mv_mlt(A_tilde_trans_mu, errorvector, temp);
+ startvector = v_sub(startvector, temp, startvector); // update startvector //
-
- 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);
+ // 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 //
+ 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);
+
+ /*
+
+ //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);
+
+ // 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;
+
+ //G_message("get_element: %lf", G_matrix_get_element (startvector, startvector->cols, k));
+
+ 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]);
-
+ for (i = 0; i < A->cols; i++) // no. of spectra //
+ result_cell[i][col] = (CELL) (100 * fraction->ve[i]);
+
// save error and iterations//
- error_cell[col] = (CELL) (100 * error);
- iter_cell[col] = iterations;
+ error_cell[col] = (CELL) (100 * error);
+ iter_cell[col] = iterations;
- //V_FREE(fraction);
- //V_FREE(b);
-
- // }
-
+ //V_FREE(fraction);
+ //V_FREE(b);
+ // }
- //---------- end of second contraint -----------------------
- // store fractions in resulting rows of resulting files
- // (number of bands = vector dimension)
+
+ //---------- 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 ("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);
- // }
- // }
-
+ 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));
+ }
- // 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]);
+ // save error and iterations
+ error_cell[col] = (CELL) (100 * error);
+ iter_cell[col] = iterations;
- if (error_fd > 0)
- G_put_map_row (error_fd, error_cell);
+ G_vector_free (fraction);
+ // G_vector_free (b);
- if (iter_fd > 0)
- G_put_map_row (iter_fd, iter_cell);
-
- } // rows loop
- G_percent (row, nrows, 2);
+ */
+ } //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 < Ref.nfiles; i++) // no. of bands
+ G_unopen_cell(cellfd[i]);
- for (i = 0; i < A->cols; i++) // no. of spectra
+ for (i = 0; i < A->cols; i++) // no. of spectra
{
- char command[1080];
+ char command[1080];
- G_close_cell (resultfd[i]);
+ 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);
+ // 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);
- // create histogram
- do_histogram (result_name, Ref.file[i].mapset);
+
+ //G_message(command);
+ //G_system (command);
+
+ // create histogram
+ do_histogram(result_name, Ref.file[i].mapset);
}
- if (error_fd > 0)
- {
- char command[80];
+ if (error_fd > 0) {
+ char command[80];
- G_close_cell (error_fd);
- sprintf (command, "r.colors map=%s color=gyr >/dev/null", parm.error->answer);
- //G_system (command);
+ 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);
+ G_close_cell(iter_fd);
- G_matrix_free (A);
+ G_matrix_free(A);
- make_history (result_name, parm.group->answer, parm.matrixfile->answer);
-
+ make_history(result_name, parm.group->answer, parm.matrixfile->answer);
+
/*********************
*********************/
- exit (EXIT_SUCCESS);
+ exit(EXIT_SUCCESS);
}
Modified: grass-addons/grass6/imagery/i.spec.unmix/meminfo.h
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/meminfo.h 2012-12-04 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/meminfo.h 2012-12-04 17:37:01 UTC (rev 54180)
@@ -69,37 +69,37 @@
#endif
/* structure for memory information */
-typedef struct {
- long bytes; /* # of allocated bytes for each type (summary) */
- int numvar; /* # of allocated variables for each type */
+typedef struct
+{
+ long bytes; /* # of allocated bytes for each type (summary) */
+ int numvar; /* # of allocated variables for each type */
} MEM_ARRAY;
#ifdef ANSI_C
-int mem_info_is_on(void);
+int mem_info_is_on(void);
int mem_info_on(int sw);
-long mem_info_bytes(int type,int list);
-int mem_info_numvar(int type,int list);
-void mem_info_file(FILE * fp,int list);
+long mem_info_bytes(int type, int list);
+int mem_info_numvar(int type, int list);
+void mem_info_file(FILE * fp, int list);
-void mem_bytes_list(int type,int old_size,int new_size,
- int list);
+void mem_bytes_list(int type, int old_size, int new_size, int list);
void mem_numvar_list(int type, int num, int list);
-int mem_stat_reg_list(void **var,int type,int list);
+int mem_stat_reg_list(void **var, int type, int list);
int mem_stat_mark(int mark);
-int mem_stat_free_list(int mark,int list);
+int mem_stat_free_list(int mark, int list);
int mem_stat_show_mark(void);
-void mem_stat_dump(FILE *fp,int list);
-int mem_attach_list(int list,int ntypes,char *type_names[],
- int (*free_funcs[])(), MEM_ARRAY info_sum[]);
+void mem_stat_dump(FILE * fp, int list);
+int mem_attach_list(int list, int ntypes, char *type_names[],
+ int (*free_funcs[])(), MEM_ARRAY info_sum[]);
int mem_free_vars(int list);
int mem_is_list_attached(int list);
-void mem_dump_list(FILE *fp,int list);
-int mem_stat_reg_vars(int list,int type,...);
+void mem_dump_list(FILE * fp, int list);
+int mem_stat_reg_vars(int list, int type, ...);
#else
int mem_info_is_on();
@@ -123,7 +123,7 @@
void mem_dump_list();
int mem_stat_reg_vars();
-#endif
+#endif
/* macros */
@@ -141,11 +141,12 @@
/* internal type */
-typedef struct {
- char **type_names; /* array of names of types (strings) */
- int (**free_funcs)(); /* array of functions for releasing types */
- unsigned ntypes; /* max number of types */
- MEM_ARRAY *info_sum; /* local array for keeping track of memory */
+typedef struct
+{
+ char **type_names; /* array of names of types (strings) */
+ int (**free_funcs) (); /* array of functions for releasing types */
+ unsigned ntypes; /* max number of types */
+ MEM_ARRAY *info_sum; /* local array for keeping track of memory */
} MEM_CONNECT;
/* max number of lists of types */
Modified: grass-addons/grass6/imagery/i.spec.unmix/open.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/open.c 2012-12-04 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/open.c 2012-12-04 17:37:01 UTC (rev 54180)
@@ -1,12 +1,12 @@
/* Spextral unmixing with Singular Value Decomposition */
-/* (c) 15. Jan. 1999 Markus Neteler, Hannover*/
+/* (c) 15. Jan. 1999 Markus Neteler, Hannover */
/* 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
-*/
+ 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
+ */
#include <stdio.h>
#include <math.h>
@@ -16,143 +16,141 @@
#include "global.h"
-int open_files (char *matrixfile,
- char *img_grp,
- char *iter_name,
- char *error_name,
- mat_struct *A)
+int open_files(char *matrixfile,
+ char *img_grp,
+ char *iter_name, char *error_name, mat_struct * A)
{
char result_name[80];
- char *result_prefix="out";
+ char *result_prefix = "out";
FILE *fp;
int i, matrixsize;
mat_struct A_input;
-
-
-// mat_struct A_input1;
-/*
- 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);
+ // mat_struct A_input1;
+ /*
- //G_matrix_print(&A_input);
+ if ((fp = fopen (matrixfile, "r")) == NULL)
+ G_fatal_error (_("Matrix file %s not found."), matrixfile);
- // G_warning("matrix read done");
-#if 0
- G_message(_("Your spectral matrix = %d"), m_output(A_input));
-#endif
+ //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);
+ //G_matrix_print(&A_input);
-// A = m_get(A_input->rows, A_input->cols);
- G_matrix_clone2(&A_input, A);
-
+ // G_warning("matrix read done");
+ #if 0
+ G_message(_("Your spectral matrix = %d"), m_output(A_input));
+ #endif
-
- //*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);
+ // A = m_get(A_input->rows, A_input->cols);
- if ((A->rows) < (A->cols))
- G_fatal_error (_("Need number of cols >= rows to perform least squares fitting."));
+ G_matrix_clone2(&A_input, A);
- // 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)
- {
- 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);
- }
+ //*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"));
- // 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 only %i cols)."),
- Ref.nfiles, img_grp, A->cols);
+ A = G_matrix_transpose (&A_input);
- // 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_matrix_free (&A_input);
- G_message (_("Opening input file no. %i [%s]"), (i + 1), Ref.file[i].name);
+ // G_matrix_print(A);
- 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);
- }
+ if ((A->rows) < (A->cols))
+ G_fatal_error (_("Need number of cols >= rows to perform least squares fitting."));
- // open files for results
- result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
- resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+ // number of rows must be equivalent to no. of bands
+ matrixsize = A->cols;
- for (i = 0; i < A->cols; i++) // no. of spectra
- {
- sprintf (result_name, "%s.%d", result_prefix, (i + 1));
- G_message (_("Opening output file [%s]"), result_name);
+ // open input files from group
+ if (!I_find_group (img_grp))
+ G_fatal_error (_("Unable to find imagery group %s."), img_grp);
- 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."));
- }
+ 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);
+ }
- // open file containing SMA error
- error_cell = (CELL *) G_malloc (sizeof(CELL *));
- if (error_name)
- {
- G_message (_("Opening error file [%s]"), error_name);
+ // 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 only %i cols)."),
+ Ref.nfiles, img_grp, A->cols);
- 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 ();
- }
+ // 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 ();
- // open file containing number of iterations
- iter_cell = (CELL *) G_malloc (sizeof(CELL *));
- if (iter_name)
- {
- G_message (_("Opening iteration file [%s]"), iter_name);
+ G_message (_("Opening input file no. %i [%s]"), (i + 1), Ref.file[i].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 ();
- }
+ 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);
+ }
- //G_matrix_print(A);
+ // open files for results
+ result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
- return matrixsize;
- */
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ 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)
+ {
+ 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 ();
+ }
+
+ //G_matrix_print(A);
+
+ return matrixsize;
+ */
return 0;
}
@@ -163,43 +161,42 @@
-mat_struct* open_files2 (char *matrixfile,
- char *img_grp,
- char *result_prefix,
- char *iter_name,
- char *error_name)
+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;
-
-
+
+
+ // 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);
+ 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);
+ 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");
+ // G_warning("matrix read done");
#if 0
G_message(_("Your spectral matrix = %d"), m_output(A_input));
#endif
@@ -208,105 +205,104 @@
* Don't mix rows and cols in the source code and the modules
* messages output! */
-// A = m_get(A_input->rows, A_input->cols);
+ // 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);
+
+
+ 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"));
+ G_fatal_error(_("Unable to allocate memory for matrix"));
- A = G_matrix_transpose (&A_input);
-
- // G_matrix_free (&A_input);
+ A = G_matrix_transpose(&A_input);
-// G_matrix_print(A);
+ // 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."));
+ 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);
+ 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)
- {
+ 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);
+ 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);
+ 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);
+ 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 ();
+ 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);
+ 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);
+ 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
- result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
- resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+ 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
+ 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);
+ 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."));
- }
+ 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);
+ 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 ();
+ 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)
- {
- G_message (_("Opening iteration file [%s]"), iter_name);
+ iter_cell = (CELL *) G_malloc(sizeof(CELL *));
+ if (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);
+ 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();
}
Modified: grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c 2012-12-04 17:35:01 UTC (rev 54179)
+++ grass-addons/grass6/imagery/i.spec.unmix/spec_angle.c 2012-12-04 17:37:01 UTC (rev 54180)
@@ -29,7 +29,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)
{
vec_struct *vtmp1;
double norm1, norm2, norm3;
@@ -37,14 +37,14 @@
/* 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) */
+ // 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) */
- G_vector_free (vtmp1);
+ G_vector_free(vtmp1);
/* Calculate angle and return in degree globally */
- return (acos (norm1 / (norm2 * norm3)) * 180 / M_PI);
+ return (acos(norm1 / (norm2 * norm3)) * 180 / M_PI);
}
More information about the grass-commit
mailing list