[GRASS-SVN] r54197 - in grass-addons: grass6/imagery/i.spec.unmix grass7/imagery/i.spec.unmix grass7/imagery/i.spec.unmix/unused
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Dec 5 03:06:36 PST 2012
Author: rashadkm
Date: 2012-12-05 03:06:36 -0800 (Wed, 05 Dec 2012)
New Revision: 54197
Modified:
grass-addons/grass6/imagery/i.spec.unmix/la_extra.c
grass-addons/grass6/imagery/i.spec.unmix/main.c
grass-addons/grass6/imagery/i.spec.unmix/open.c
grass-addons/grass7/imagery/i.spec.unmix/histogram.c
grass-addons/grass7/imagery/i.spec.unmix/la_extra.c
grass-addons/grass7/imagery/i.spec.unmix/main.c
grass-addons/grass7/imagery/i.spec.unmix/open.c
grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c
grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h
Log:
updated C-comment style and error check for libLAPACK and libLAS
Modified: grass-addons/grass6/imagery/i.spec.unmix/la_extra.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/la_extra.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass6/imagery/i.spec.unmix/la_extra.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -42,7 +42,7 @@
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);
}
@@ -70,8 +70,7 @@
sprintf(numbuf, "%14.6f", element);
strcat(buf, numbuf);
- //if( j < mt->cols - 1 )
- //strcat(buf, ", ");
+
}
G_message("%s", buf);
}
@@ -79,24 +78,6 @@
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");
-
- */
}
@@ -115,13 +96,11 @@
- //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;
}
}
@@ -158,7 +137,7 @@
m = matrix->rows;
n = matrix->cols;
for (i = 0; i < m; i++) {
- //__smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
+ /* __smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n); */
/**************************************************/
for (j = 0; j < n; j++) {
@@ -166,9 +145,6 @@
}
}
- //G_matrix_print(matrix,"matrix");
- //G_matrix_print(matrix,"out");
-
/**************************************************/
return (out);
}
@@ -231,11 +207,10 @@
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);
@@ -264,12 +239,12 @@
/* register Real sum; */
- //G_vec_print(out, "xx");
+ /*
- // if ( A==(MAT *)NULL || b==(VEC *)NULL )
- // error(E_NULL,"mv_mlt");
+ if ( A==(MAT *)NULL || b==(VEC *)NULL )
+ error(E_NULL,"mv_mlt");
+ */
- //G_matrix_print(A, "A(mlt)");
if (A->cols != b->dim)
@@ -278,7 +253,7 @@
if (b == out)
G_fatal_error("mv_mlt2(error)");
- //error(E_INSITU,"mv_mlt");
+
if (!out) {
G_fatal_error("mv_mltsss3(error)");
exit(1);
@@ -289,39 +264,33 @@
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;
}
@@ -383,7 +352,7 @@
}
vector->dim = vector->max_dim = size;
- //G_vec_print(vector, "vector");
+
return vector;
}
@@ -391,9 +360,9 @@
VEC *G_vec_get2(int size, VEC * vector)
{
- //VEC *vector;
+
vector = (VEC *) G_malloc(sizeof(VEC));
@@ -408,7 +377,7 @@
}
vector->dim = vector->max_dim = size;
- //G_vec_print(vector, "vector");
+
return vector;
}
@@ -420,8 +389,8 @@
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];
@@ -448,14 +417,8 @@
- //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;
@@ -488,7 +451,7 @@
return NULL;
}
-#if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS) //&& defined(HAVE_G2C_H)
+#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;
@@ -539,7 +502,6 @@
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) {
@@ -557,6 +519,6 @@
G_matrix_set_element(out, i, j, val);
}
}
- //G_matrix_print(out);
+
return 0;
}
Modified: grass-addons/grass6/imagery/i.spec.unmix/main.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/main.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass6/imagery/i.spec.unmix/main.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -19,6 +19,16 @@
#define GLOBAL
#include <grass/config.h>
+
+
+#ifndef HAVE_LIBBLAS
+#error GRASS is not configured with libLAS
+#endif
+
+#ifndef HAVE_LIBLAPACK
+#error GRASS is not configured with LAPACK
+#endif
+
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
Modified: grass-addons/grass6/imagery/i.spec.unmix/open.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/open.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass6/imagery/i.spec.unmix/open.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -27,54 +27,39 @@
mat_struct A_input;
- // mat_struct A_input1;
- /*
+ /* 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);
- //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);
- //*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->cols;
- // open input files from group
+
if (!I_find_group (img_grp))
G_fatal_error (_("Unable to find imagery group %s."), img_grp);
@@ -89,14 +74,14 @@
"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 only %i cols)."),
Ref.nfiles, img_grp, A->cols);
- // 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++)
@@ -109,11 +94,10 @@
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));
- for (i = 0; i < A->cols; i++) // no. of spectra
+ for (i = 0; i < A->cols; i++)
{
sprintf (result_name, "%s.%d", result_prefix, (i + 1));
G_message (_("Opening output file [%s]"), result_name);
@@ -123,7 +107,6 @@
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)
{
@@ -135,7 +118,7 @@
error_cell = G_allocate_cell_buf ();
}
- // open file containing number of iterations
+
iter_cell = (CELL *) G_malloc (sizeof(CELL *));
if (iter_name)
{
@@ -147,7 +130,7 @@
iter_cell = G_allocate_cell_buf ();
}
- //G_matrix_print(A);
+
return matrixsize;
*/
@@ -173,9 +156,7 @@
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
@@ -185,18 +166,16 @@
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
@@ -205,9 +184,7 @@
* 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);
@@ -217,17 +194,14 @@
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
+ /* number of rows must be equivalent to no. of bands */
matrixsize = A->rows;
- // open input files from group
+ /* open input files from group */
if (!I_find_group(img_grp))
G_fatal_error(_("Unable to find imagery group %s."), img_grp);
@@ -243,13 +217,13 @@
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)
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
+ /* 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++) {
@@ -267,11 +241,11 @@
- // 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->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));
@@ -282,7 +256,7 @@
G_fatal_error(_("GRASS-DB internal error: Unable to proceed."));
}
}
- // open file containing SMA error
+ /* open file containing SMA error */
error_cell = (CELL *) G_malloc(sizeof(CELL *));
if (error_name) {
G_message(_("Opening error file [%s]"), error_name);
@@ -293,7 +267,7 @@
error_cell = G_allocate_cell_buf();
}
- // open file containing number of iterations
+ /* open file containing number of iterations */
iter_cell = (CELL *) G_malloc(sizeof(CELL *));
if (iter_name) {
G_message(_("Opening iteration file [%s]"), iter_name);
Modified: grass-addons/grass7/imagery/i.spec.unmix/histogram.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/histogram.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass7/imagery/i.spec.unmix/histogram.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -12,8 +12,9 @@
struct Cell_head region;
Rast_get_cellhd(name, mapset, &cellhd);
- // if (Rast_get_cellhd (name, mapset, &cellhd) < 0)
- // return 0;
+ /* if (Rast_get_cellhd (name, mapset, &cellhd) < 0)
+ return 0;
+ */
G_set_window(&cellhd);
fd = Rast_open_old(name, mapset);
@@ -27,14 +28,15 @@
nrows = region.rows;
ncols = region.cols;
- // nrows = G_window_rows ();
- // ncols = G_window_cols ();
+ /* nrows = G_window_rows ();
+ ncols = G_window_cols ();
+ */
cell = Rast_allocate_c_buf();
Rast_init_cell_stats(&statf);
for (row = 0; row < nrows; row++) {
Rast_get_c_row_nomask(fd, cell, row);
- // break;
+ /* break; */
Rast_update_cell_stats(cell, ncols, &statf);
}
Modified: grass-addons/grass7/imagery/i.spec.unmix/la_extra.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/la_extra.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass7/imagery/i.spec.unmix/la_extra.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -40,7 +40,7 @@
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);
}
@@ -68,8 +68,7 @@
sprintf(numbuf, "%14.6f", element);
strcat(buf, numbuf);
- //if( j < mt->cols - 1 )
- //strcat(buf, ", ");
+
}
G_message("%s", buf);
}
@@ -77,24 +76,7 @@
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");
-
- */
}
@@ -113,13 +95,11 @@
- //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;
}
}
@@ -156,7 +136,7 @@
m = matrix->rows;
n = matrix->cols;
for (i = 0; i < m; i++) {
- //__smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
+ /* __smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n); */
/**************************************************/
for (j = 0; j < n; j++) {
@@ -164,9 +144,6 @@
}
}
- //G_matrix_print(matrix,"matrix");
- //G_matrix_print(matrix,"out");
-
/**************************************************/
return (out);
}
@@ -229,9 +206,7 @@
if (out == NULL)
out = G_vec_get(vec1->dim);
- //out = v_resize(out,A->m);
- //G_vec_print(vec1, "vec1");
if (out->dim != vec1->dim)
@@ -263,12 +238,6 @@
/* register Real sum; */
-
- // 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)");
@@ -276,7 +245,7 @@
if (b == out)
G_fatal_error("mv_mlt2(error)");
- //error(E_INSITU,"mv_mlt");
+
if (!out) {
G_fatal_error("mv_mltsss3(error)");
exit(1);
@@ -287,14 +256,8 @@
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;
@@ -306,15 +269,12 @@
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);
-
}
}
@@ -389,9 +349,6 @@
VEC *G_vec_get2(int size, VEC * vector)
{
- //VEC *vector;
-
-
vector = (VEC *) G_malloc(sizeof(VEC));
@@ -418,7 +375,6 @@
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];
@@ -445,15 +401,6 @@
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;
@@ -486,7 +433,7 @@
return NULL;
}
-#if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS) //&& defined(HAVE_G2C_H)
+#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;
@@ -537,8 +484,6 @@
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"));
@@ -555,6 +500,6 @@
G_matrix_set_element(out, i, j, val);
}
}
- //G_matrix_print(out);
+
return 0;
}
Modified: grass-addons/grass7/imagery/i.spec.unmix/main.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/main.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass7/imagery/i.spec.unmix/main.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -15,10 +15,26 @@
* comes with GRASS for details.
*
*****************************************************************************/
+
+
+
+
+
#define GLOBAL
#include <grass/config.h>
+
+
+#ifndef HAVE_LIBBLAS
+#error GRASS is not configured with libLAS
+#endif
+
+#ifndef HAVE_LIBLAPACK
+#error GRASS is not configured with LAPACK
+#endif
+
+
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
@@ -44,8 +60,10 @@
}
+
int main(int argc, char *argv[])
{
+
char result_name[80];
int nrows, ncols;
int row;
@@ -114,8 +132,7 @@
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
@@ -150,41 +167,36 @@
Avector1 = G_matvect_get_column2(A, i);
- // G_matrix_print(Avector1);
- //G_warning("Avector1: %d", Avector1->rows);
-
- // get the max. element of this vector
+ /* 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
+ /* get next col in A */
Avector2 = G_matvect_get_column(A, j);
- //G_matrix_print(Avector2);
- // get the max. element of this vector
+
+ /* 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
+ /* find max of matrix A */
+ /* max_total = (find_max (max1, max2), max_total); */
+
+ /* save angle in degree */
anglefield[i][j] = spectral_angle(Avector1, Avector2);
- //G_warning("anglefield[i][j]: %lf", anglefield[i][j]);
+
}
}
@@ -197,7 +209,7 @@
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
+ /* 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]);
@@ -213,24 +225,21 @@
if (!error)
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
+ /* 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
+ */
-
- // memory allocation
+ /* 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"));
- //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,
@@ -238,49 +247,48 @@
- // fill last row with 1 elements
+ /* 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.
+ /* 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);
+ /* calculate the transpose of 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
- //
+ /* 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
+ /* 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);
@@ -290,25 +298,25 @@
- 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 */
- // length: no. of bands
+ /* length: no. of bands */
+
if (A_times_startvector == NULL)
G_fatal_error(_("Unable to allocate memory for vector"));
- // length: no. of bands
+ /* length: no. of bands */
if (errorvector == NULL)
G_fatal_error(_("Unable to allocate memory for vector"));
- // length: no. of spectra
+ /* length: no. of spectra */
if (temp == NULL)
G_fatal_error(_("Unable to allocate memory for vector"));
@@ -318,12 +326,11 @@
if (A_tilde_trans_mu == NULL)
G_fatal_error(_("Unable to allocate memory for matrix"));
- // Now we can calculated the fractions pixelwise
- //nrows = G_region_rows (); // get geographical region
- //ncols = G_window_cols ();
+ /* Now we can calculated the fractions pixelwise */
- G_get_window(®ion);
+ G_get_window(®ion); /* get geographical region */
+
nrows = region.rows;
ncols = region.cols;
@@ -337,20 +344,11 @@
G_percent(row, nrows, 1);
- // get one row for all bands
+ /* get one row for all bands */
for (band = 0; band < Ref.nfiles; band++)
Rast_get_c_row(cellfd[band], cell[band], row);
-
- // if (Rast_get_c_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 (Rast_get_c_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++) {
@@ -360,51 +358,42 @@
int iterations = 0;
- // get pixel values of each band and store in b vector:
- // length: no. of bands + 1 (GAMMA)
+ /* 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"));
-
- //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
+ /* 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
+ /* calculate fraction vector for current pixel
+ Result is stored in fractions vector
+ with second constraint: Sum x_i = 1
+ */
+ /* get start vector and initialize it with equal fractions:
+ using the neighbor pixelvector as startvector
+ */
- // get start vector and initialize it with equal fractions:
- // using the neighbor pixelvector as startvector
-
-
- // solve with iterative solution:
+ /* solve with iterative solution: */
while (fabs(change) > 0.0001) {
@@ -420,18 +409,15 @@
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 //
+ 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 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 //
+ /* Check the deviation */
double norm2 = v_norm2(errorvector);
change = deviation - norm2;
@@ -439,131 +425,47 @@
iterations++;
- // if(fabs (change) > 0.0001)
- //G_message("change=%lf, deviation=%lf",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 //
+ /* 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 //
+ /* 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;
- //V_FREE(fraction);
- //V_FREE(b);
+ /****V_FREE(fraction);
+ V_FREE(b);
+ *****/
- // }
+ } /* end cols loop */
- //---------- 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
+ /* write the resulting rows into output files: */
+ for (i = 0; i < A->cols; i++) /* no. of spectra */
Rast_put_c_row(resultfd[i], result_cell[i]);
if (error_fd > 0)
@@ -572,29 +474,29 @@
if (iter_fd > 0)
Rast_put_c_row(iter_fd, iter_cell);
- } // rows loop
+ } /* rows loop */
G_percent(row, nrows, 2);
- // close files
- for (i = 0; i < Ref.nfiles; i++) // no. of bands
+ /* close files */
+ for (i = 0; i < Ref.nfiles; i++) /* no. of bands */
Rast_unopen(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];
Rast_close(resultfd[i]);
- // make grey scale color table
+ /* 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);
+ /* G_message(command); */
+ /* G_system (command); */
- // create histogram
+ /* create histogram */
do_histogram(result_name, Ref.file[i].mapset);
}
@@ -602,8 +504,8 @@
char command[80];
Rast_close(error_fd);
- //sprintf (command, "r.colors map=%s color=gyr >/dev/null", parm.error->answer);
- //G_system (command);
+ /* sprintf (command, "r.colors map=%s color=gyr >/dev/null", parm.error->answer);*/
+ /* G_system (command); */
}
if (iter_fd > 0)
@@ -613,7 +515,5 @@
make_history(result_name, parm.group->answer, parm.matrixfile->answer);
-/*********************
-*********************/
exit(EXIT_SUCCESS);
}
Modified: grass-addons/grass7/imagery/i.spec.unmix/open.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/open.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass7/imagery/i.spec.unmix/open.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -20,61 +20,36 @@
char *img_grp,
char *iter_name, char *error_name, mat_struct * A)
{
+ /*
+
char result_name[80];
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);
- //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);
-
-
- //*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->cols;
- // open input files from group
if (!I_find_group (img_grp))
G_fatal_error (_("Unable to find imagery group %s."), img_grp);
@@ -89,14 +64,12 @@
"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 only %i cols)."),
Ref.nfiles, img_grp, A->cols);
- // 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++)
@@ -109,11 +82,10 @@
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));
- for (i = 0; i < A->cols; i++) // no. of spectra
+ for (i = 0; i < A->cols; i++)
{
sprintf (result_name, "%s.%d", result_prefix, (i + 1));
G_message (_("Opening output file [%s]"), result_name);
@@ -123,7 +95,6 @@
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)
{
@@ -135,7 +106,6 @@
error_cell = Rast_allocate_c_buf();
}
- // open file containing number of iterations
iter_cell = (CELL *) G_malloc (sizeof(CELL *));
if (iter_name)
{
@@ -147,20 +117,11 @@
iter_cell = Rast_allocate_c_buf();
}
- //G_matrix_print(A);
-
return matrixsize;
*/
return 0;
}
-
-
-
-
-
-
-
mat_struct *open_files2(char *matrixfile,
char *img_grp,
char *result_prefix,
@@ -173,61 +134,44 @@
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). */
+ * 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! */
+ * 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
+ /* number of rows must be equivalent to no. of bands */
matrixsize = A->rows;
- // open input files from group
+ /* open input files from group */
if (!I_find_group(img_grp))
G_fatal_error(_("Unable to find imagery group %s."), img_grp);
@@ -243,13 +187,13 @@
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)
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
+ /* 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++) {
@@ -267,11 +211,11 @@
- // 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->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));
@@ -282,7 +226,7 @@
G_fatal_error(_("GRASS-DB internal error: Unable to proceed."));
}
}
- // open file containing SMA error
+ /* open file containing SMA error */
error_cell = (CELL *) G_malloc(sizeof(CELL *));
if (error_name) {
G_message(_("Opening error file [%s]"), error_name);
@@ -293,7 +237,7 @@
error_cell = Rast_allocate_c_buf();
}
- // open file containing number of iterations
+ /* open file containing number of iterations */
iter_cell = (CELL *) G_malloc(sizeof(CELL *));
if (iter_name) {
G_message(_("Opening iteration file [%s]"), iter_name);
Modified: grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c 2012-12-05 11:06:36 UTC (rev 54197)
@@ -37,7 +37,7 @@
/* Measure spectral angle */
/* multiply one A column with second */
- // vtmp1 = G_vector_init (0, 0, RVEC);
+
vtmp1 = G_vector_product(Avector1, Avector2);
norm1 = G_vector_norm1(vtmp1); /* calculate 1-norm */
norm2 = G_vector_norm_euclid(Avector1); /* calculate 2-norm (Euclidean) */
Modified: grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h 2012-12-05 10:42:08 UTC (rev 54196)
+++ grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h 2012-12-05 11:06:36 UTC (rev 54197)
@@ -168,7 +168,7 @@
MAT *m_poly(MAT *,VEC *,MAT *);
/* FFT */
-//void fft(VEC *,VEC *);
-//void ifft(VEC *,VEC *);
+/*void fft(VEC *,VEC *);*/
+/*void ifft(VEC *,VEC *);*/
#endif /* __MATRIX2_H__ */
More information about the grass-commit
mailing list