[GRASS-SVN] r54312 - grass-addons/grass7/imagery/i.spec.unmix
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Dec 16 02:44:22 PST 2012
Author: rashadkm
Date: 2012-12-16 02:44:22 -0800 (Sun, 16 Dec 2012)
New Revision: 54312
Removed:
grass-addons/grass7/imagery/i.spec.unmix/open.h
Modified:
grass-addons/grass7/imagery/i.spec.unmix/global.h
grass-addons/grass7/imagery/i.spec.unmix/la_extra.c
grass-addons/grass7/imagery/i.spec.unmix/la_extra.h
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
Log:
cleaning up gmath specific code; la_extra ready to be merged to lib/gmath/la.c
Modified: grass-addons/grass7/imagery/i.spec.unmix/global.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/global.h 2012-12-16 10:33:35 UTC (rev 54311)
+++ grass-addons/grass7/imagery/i.spec.unmix/global.h 2012-12-16 10:44:22 UTC (rev 54312)
@@ -29,6 +29,6 @@
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 mat_struct *open_files(char *matrixfile, char *img_grp, char *result_prefix, char *iter_name, char *error_name);
#endif /* __GLOBAL_H__ */
Modified: grass-addons/grass7/imagery/i.spec.unmix/la_extra.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/la_extra.c 2012-12-16 10:33:35 UTC (rev 54311)
+++ grass-addons/grass7/imagery/i.spec.unmix/la_extra.c 2012-12-16 10:44:22 UTC (rev 54312)
@@ -1,4 +1,4 @@
-/* TODO: move to GMATHLIB? */
+/* TODO: move to GMATHLIB */
#include <stdio.h> /* needed here for ifdef/else */
#include <stdlib.h>
@@ -13,436 +13,153 @@
#include "la_extra.h"
-vec_struct *G_matvect_get_column2(mat_struct * mt, int col)
-{
- int i; /* loop */
- vec_struct *vc1;
-
- if (col < 0 || col >= mt->cols) {
- G_warning(_("Specified matrix column index is outside range"));
- return NULL;
- }
-
- if (!mt->is_init) {
- G_warning(_("Matrix is not initialised"));
- return NULL;
- }
-
- if ((vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) == NULL) {
- G_warning(_("Could not allocate space for vector structure"));
- return NULL;
- }
-
- for (i = 0; i < mt->rows; i++) {
- double dd = G_matrix_get_element(mt, i, col);
-
-
- G_matrix_set_element((mat_struct *) vc1, i, 0, dd);
-
- }
-
- return vc1;
-}
-
-
-void G_matrix_print2(mat_struct * mt, const char *name)
-{
- int i, j;
-
-
- if (mt != NULL) {
- G_message("start matrix(%s)", name);
- G_message("Size: %d x %d", mt->rows, mt->cols);
-
- for (i = 0; i < mt->rows; i++) {
- char buf[2048], numbuf[640];
-
- sprintf(buf, "row%d: ", i);
- for (j = 0; j < mt->cols; j++) {
-
- double element = G_matrix_get_element(mt, i, j);
-
- sprintf(numbuf, "%14.6f", element);
- strcat(buf, numbuf);
-
- }
- G_message("%s", buf);
- }
-
- G_message("end matrix(%s)", name);
- }
-
-
-}
-
-
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++)
+ {
+ for (j = 0; j < cols; j++)
+ matrix->vals[index++] = in->vals[i + j * cols];
+ }
- mat_struct *matrix;
+ 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;
- matrix = G_matrix_init(rows, cols, rows);
-
-
- int i, j, p, index = 0;
-
- for (i = 0; i < rows; i++) {
-
-
-
-
- for (j = 0; j < cols; j++) {
-
- matrix->vals[index++] = in->vals[i + j * cols];
-
- }
-
- }
-
- 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;
+ return (matrix);
}
-
-mat_struct *sm_mlt(double scalar, mat_struct * matrix, mat_struct * out)
+mat_struct *G_matrix_scalar_mul(double scalar, mat_struct * matrix, mat_struct * out)
{
- int m, n, i, j;
+ int m, n, i, j;
+ int index = 0;
- int index = 0;
+ if (matrix == NULL)
+ {
+ G_warning (_("Input matrix is uninitialized"));
+ return NULL;
+ }
- if (matrix == NULL)
- G_fatal_error("sm_mlt1(error)");
+ if (out == NULL)
+ out = G_matrix_init(matrix->rows, matrix->cols, matrix->rows);
- 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 (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++)
+ {
+ for (j = 0; j < n; j++)
+ {
+ doublereal value = scalar * G_matrix_get_element(matrix, i, j);
+ G_matrix_set_element (out, i,j, value);
+ }
+ }
- 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];
- }
- }
-
- /**************************************************/
- return (out);
+ return (out);
}
-VEC *G_vec_copy(VEC * in)
+vec_struct* G_mat_vector_product(mat_struct * A, vec_struct * b,vec_struct *out)
{
- int i;
- VEC *out;
+ unsigned int i, m, n, j;
+ register doublereal sum;
- if (!in)
- G_fatal_error("v_copy(error1)");
+ if (A->cols != b->ldim) {
+ G_warning (_("Input matrix and vector have differing dimensions"));
+ return NULL;
+ }
+ if (!out) {
+ G_warning (_("Output vector is uninitialized"));
+ return NULL;
+ }
+ if (out->ldim != A->rows) {
+ G_warning (_("Output vector has incorrect dimension"));
+ return NULL;
+ }
- int dim = in->dim;
- out = G_vec_get(dim);
+ m = A->rows;
+ n = A->cols;
-
- for (i = 0; i < dim; i++) {
- out->ve[i] = in->ve[i];
- }
-
- return (out);
-}
-
-
-double v_norm2(VEC * x)
-{
- int i, dim;
- double s, sum;
-
- if (!x)
- G_fatal_error("v_norm2(error1)");
-
-
- dim = x->dim;
-
+ for (i = 0; i < m; i++) {
sum = 0.0;
-
- for (i = 0; i < dim; i++) {
- sum += x->ve[i] * x->ve[i];
- }
-
-
-
- return sqrt(sum);
+ int width = A->rows;
+ for (j = 0; j < n; j++) {
+ sum += A->vals[i + j * width] * b->vals[j];
+ out->vals[i] = sum;
+ }
+ }
+ return (out);
}
-VEC *v_sub(VEC * vec1, VEC * vec2, VEC * out)
+vec_struct *
+G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out)
{
- /* u_int i, dim; */
- /* Real *out_ve, *vec1_ve, *vec2_ve; */
-
- if (!vec1 || !vec2)
- G_fatal_error("v_sub1(error)");
-
- if (vec1->dim != vec2->dim)
- G_fatal_error("v_sub2(error)");
-
- if (out == NULL)
- out = G_vec_get(vec1->dim);
-
- //G_vec_print(vec2, "vec2");
-
- if (out->dim != vec1->dim)
- out = G_vec_resize(out, vec1->dim);
-
- int i;
-
- for (i = 0; i < vec1->dim; i++)
- out->ve[i] = vec1->ve[i] - vec2->ve[i];
-
- /************************************************************
- dim = vec1->dim;
- out_ve = out->ve; vec1_ve = vec1->ve; vec2_ve = vec2->ve;
- for ( i=0; i<dim; i++ )
- out->ve[i] = vec1->ve[i]-vec2->ve[i];
- (*out_ve++) = (*vec1_ve++) - (*vec2_ve++);
- ************************************************************/
-
- return (out);
-}
-
-
-
-VEC *mv_mlt(mat_struct * A, VEC * b, VEC * out)
-{
- unsigned int i, m, n, j;
- double **A_v, *b_v /*, *A_row */ ;
-
- /* register Real sum; */
-
-
- if (A->cols != b->dim)
-
- G_fatal_error("mv_mlt1(error)");
-
-
- if (b == out)
- G_fatal_error("mv_mlt2(error)");
-
- if (!out) {
- G_fatal_error("mv_mltsss3(error)");
- out = G_vec_get2(A->rows, out);
- }
- if (out->dim != A->rows) {
- G_fatal_error("mv_mlt3(error)");
- out = G_vec_resize(out, A->rows);
- }
-
-
- m = A->rows;
- n = A->cols;
- A_v = A->vals;
- b_v = b->ve;
-
-
-
- for (i = 0; i < m; i++) {
- double sum = 0.0;
- int width = A->rows;
-
- for (j = 0; j < n; j++) {
-
- sum += A->vals[i + j * width] * b->ve[j];
- out->ve[i] = sum;
-
-
- }
- }
-
-
- return out;
-}
-
-
-VEC *G_vec_resize(VEC * in, int size)
-{
-
- VEC *vector;
-
-
- vector = (VEC *) G_malloc(sizeof(VEC));
-
-
- vector->ve = (double *)G_malloc(size * sizeof(double));
- int i, j;
-
- G_message(":%d", in->dim);
- for (i = 0; i < in->dim; i++) {
-
- vector->ve[i] = in->ve[i];
- G_message("ss:%lf", in->ve[i]);
-
- }
-
- if (size > in->dim)
- for (j = i; j < size; j++)
- vector->ve[j] = 0.0;
-
- vector->dim = vector->max_dim = size;
-
- return vector;
-}
-
-
-
-
-
-
-
-VEC *G_vec_get(int size)
-{
-
- VEC *vector;
-
-
- vector = (VEC *) G_malloc(sizeof(VEC));
-
-
- vector->ve = (double *)G_malloc(size * sizeof(double));
- int i;
-
- for (i = 0; i < size; i++) {
-
- vector->ve[i] = 0.0;
-
-
- }
-
- vector->dim = vector->max_dim = size;
- //G_vec_print(vector, "vector");
- return vector;
-}
-
-
-VEC *G_vec_get2(int size, VEC * vector)
-{
-
- vector = (VEC *) G_malloc(sizeof(VEC));
-
-
- vector->ve = (double *)G_malloc(size * sizeof(double));
- int i;
-
- for (i = 0; i < size; i++) {
-
- vector->ve[i] = 0.0;
-
-
- }
-
- vector->dim = vector->max_dim = size;
- //G_vec_print(vector, "vector");
- return vector;
-}
-
-
-void G_vec_print(VEC * vector, const char *name)
-{
- int i;
-
-
- if (vector != NULL) {
- G_message("start vector(%s)", name);
-
- for (i = 0; i < vector->dim; i++) {
- char buf[2048], numbuf[640];
-
- sprintf(buf, "%lf ", vector->ve[i]);
-
- G_message("%s", buf);
- }
-
- G_message("end vector(%s)", name);
- }
-
-}
-
-
-
-
-vec_struct *G_vector_product(vec_struct * v1, vec_struct * v2)
-{
int idx1, idx2, idx0;
int i;
-
- vec_struct *out = G_vector_init(v1->rows, v1->ldim, CVEC);
-
-
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)
+ 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));
+ 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));
}
- 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;
@@ -450,60 +167,3 @@
-
-int G_matrix_read2(FILE * fp, mat_struct * out)
-{
- char buff[100];
- int rows, cols;
- int i, j, row;
- double val;
-
- /* skip comments */
- for (;;) {
- if (!G_getl(buff, sizeof(buff), fp))
- return -1;
- if (buff[0] != '#')
- break;
- }
-
- if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) {
- G_warning(_("Input format error1"));
- return -1;
- }
-
-
- G_matrix_set(out, rows, cols, rows);
-
-
- 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);
- }
- }
-
- return 0;
-}
-
-int vec_free(VEC *vec)
-{
- if ( !vec || (vec->dim) < 0 )
- /* don't trust it */
- return -1;
-
- //G_free(vec->ve);
-
- G_free(vec);
-
- return 0;
-}
Modified: grass-addons/grass7/imagery/i.spec.unmix/la_extra.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/la_extra.h 2012-12-16 10:33:35 UTC (rev 54311)
+++ grass-addons/grass7/imagery/i.spec.unmix/la_extra.h 2012-12-16 10:44:22 UTC (rev 54312)
@@ -1,3 +1,5 @@
+/* TODO: move to GMATHLIB */
+
#include <grass/config.h>
#include <stdio.h>
#include <grass/blas.h>
@@ -2,24 +4,6 @@
#include <grass/lapack.h>
-typedef struct VEC_
-{
- int dim, max_dim;
- double *ve;
-} VEC;
-int G_matrix_read2(FILE * fp, mat_struct * out);
-void G_matrix_print2(mat_struct * mt, const char *name);
-vec_struct *G_matvect_get_column2(mat_struct * mt, int col);
-vec_struct *G_vector_product(vec_struct *, vec_struct *);
+mat_struct *G_matrix_scalar_mul(double scalar, mat_struct * matrix, mat_struct * out);
+vec_struct* G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out);
+vec_struct* G_mat_vector_product(mat_struct * A, vec_struct * b,vec_struct *out);
-VEC *G_vec_get(int size);
-VEC *G_vec_get2(int size, VEC * vector);
-void G_vec_print(VEC * vector, const char *name);
-VEC *G_vec_resize(VEC * in, int size);
-
-VEC *v_sub(VEC * vec1, VEC * vec2, VEC * out);
-VEC *mv_mlt(mat_struct * A, VEC * b, VEC * out);
-mat_struct *sm_mlt(double scalar, mat_struct * matrix, mat_struct * out);
-mat_struct *G_matrix_resize(mat_struct * in, int rows, int cols);
-double v_norm2(VEC * x);
-VEC *G_vec_copy(VEC * in);
-int vec_free(VEC *vec);
Modified: grass-addons/grass7/imagery/i.spec.unmix/main.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/main.c 2012-12-16 10:33:35 UTC (rev 54311)
+++ grass-addons/grass7/imagery/i.spec.unmix/main.c 2012-12-16 10:44:22 UTC (rev 54312)
@@ -45,14 +45,13 @@
#include <grass/gmath.h>
#include <grass/glocale.h>
#include "global.h"
-#include "open.h"
#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)
{
@@ -66,9 +65,13 @@
int nrows, ncols;
int row;
int i, j, k;
- VEC *b, *b_gamma;
+ vec_struct *b_gamma;
struct Cell_head region;
- VEC *startvector, *A_times_startvector, *errorvector, *temp;
+
+ // VEC *startvector, *A_times_startvector, *errorvector, *temp;
+
+ vec_struct *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;
@@ -126,10 +129,10 @@
/* 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);
+ A = open_files(parm.matrixfile->answer,
+ parm.group->answer,
+ parm.result->answer,
+ parm.iter->answer, parm.error->answer);
/* ATTENTION: Internally we work here with col-oriented matrixfile,
@@ -163,7 +166,7 @@
vec_struct *Avector1, *Avector2;
double max1, max2;
- Avector1 = G_matvect_get_column2(A, i);
+ Avector1 = G_matvect_get_column(A, i);
/* get the max. element of this vector */
@@ -181,24 +184,24 @@
if (max2 > max1)
temp = max2;
-
if (temp > max_total)
max_total = temp;
/* find max of matrix A */
- /* max_total = (find_max (max1, max2), max_total); */
+ /* max_total = (find_max (max1, max2), max_total); */
/* save angle in degree */
anglefield[i][j] = spectral_angle(Avector1, Avector2);
-
-
}
}
G_vector_free(Avector1);
G_vector_free(Avector2);
+
}
- G_message(_("Checking linear dependencies (orthogonality check) of Matrix A..."));
+ G_message("%s",
+ _
+ ("Checking linear dependencies (orthogonality check) of Matrix A..."));
for (i = 0; i < A->cols; i++)
for (j = 0; j < A->cols; j++)
@@ -236,9 +239,7 @@
G_matrix_get_element(A, i, j));
/* fill last row with 1 elements */
-
for (j = 0; j < A_tilde->cols; j++) {
-
G_matrix_set_element(A_tilde, i, j, GAMMA);
}
@@ -272,18 +273,20 @@
mu = 0.0001 * pow(10, -1 * ceil(log10(max_total)));
/* TODO: Missing? startvector = G_vector_init (0, 0, RVEC); */
- startvector = G_vec_get2(A->cols, startvector);
+ startvector = G_vector_init(A->cols, A->cols, CVEC);
+
if (startvector == NULL)
G_fatal_error(_("Unable to allocate memory for vector"));
/* TODO: Missing? A_times_startvector = G_vector_init (0, 0, RVEC); */
- A_times_startvector = G_vec_get2(A_tilde->rows, A_times_startvector); /* length: no. of bands */
+ A_times_startvector = G_vector_init(A_tilde->rows, A_tilde->rows, CVEC); /* length: no. of bands */
/* TODO: Missing? errorvector = G_vector_init (0, 0, RVEC); */
- errorvector = G_vec_get2(A_tilde->rows, errorvector); /* length: no. of bands */
+ errorvector = G_vector_init(A_tilde->rows, A_tilde->rows, CVEC); /* length: no. of bands */
/* TODO: Missing? temp = G_vector_init (0, 0, RVEC); */
- temp = G_vec_get2(A_tilde->cols, temp); /* length: no. of spectra */
+ temp = G_vector_init(A_tilde->cols, A_tilde->cols, CVEC); /* length: no. of spectra */
+
/* length: no. of bands */
if (A_times_startvector == NULL)
G_fatal_error(_("Unable to allocate memory for vector"));
@@ -307,7 +310,8 @@
nrows = region.rows;
ncols = region.cols;
- G_message(_("Calculating for %i x %i pixels (%i bands) = %i pixelvectors."),
+ G_message(_
+ ("Calculating for %i x %i pixels (%i bands) = %i pixelvectors."),
nrows, ncols, Ref.nfiles, (ncols * ncols));
for (row = 0; row < nrows; row++) {
@@ -328,22 +332,23 @@
/* 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, A_tilde->rows, CVEC);
if (b_gamma == NULL)
G_fatal_error(_("Unable to allocate memory for matrix"));
- for (band = 0; band < Ref.nfiles; band++) {
- b_gamma->ve[band] = cell[band][col];
- }
+ for (band = 0; band < Ref.nfiles; band++)
+ G_matrix_set_element(b_gamma, band, 0, cell[band][col]);
/* add GAMMA for 1. constraint as last element */
- b_gamma->ve[Ref.nfiles] = GAMMA;
+ G_matrix_set_element(b_gamma, Ref.nfiles, 0, 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));
/* calculate fraction vector for current pixel
Result is stored in fractions vector
@@ -357,53 +362,64 @@
/* solve with iterative solution: */
while (fabs(change) > 0.0001) {
- A_times_startvector =
- mv_mlt(A_tilde, startvector, A_times_startvector);
+ G_mat_vector_product(A_tilde, startvector,
+ A_times_startvector);
- errorvector =
- v_sub(A_times_startvector, b_gamma, errorvector);
+ G_vector_sub(A_times_startvector, b_gamma, errorvector);
+
A_tilde_trans_mu =
- sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
+ G_matrix_scalar_mul(mu, A_tilde_trans, A_tilde_trans_mu);
- temp = mv_mlt(A_tilde_trans_mu, errorvector, temp);
- startvector = v_sub(startvector, temp, startvector); /* update startvector */
+ G_mat_vector_product(A_tilde_trans_mu, errorvector, temp);
+ G_vector_sub(startvector, temp, startvector); /* update startvector */
- /* if one element gets negative, set it to zero */
- for (k = 0; k < (A_tilde->cols); k++) /* no. of spectra times */
- if (startvector->ve[k] < 0)
- startvector->ve[k] = 0;
+ for (k = 0; k < A_tilde->cols; k++)
+ /* no. of spectra times */
+ /* if one element gets negative, set it to zero */
+ if ((G_matrix_get_element(startvector, k, 0) < 0))
+ G_matrix_set_element(startvector, k, 0, 0);
+
/* Check the deviation */
- double norm2 = v_norm2(errorvector);
+ double norm = G_vector_norm_euclid(errorvector);
- change = deviation - norm2;
- deviation = norm2;
+ change = deviation - norm;
+ deviation = norm;
+
iterations++;
- /* G_message("change=%lf, deviation=%lf",change, 0.0001); */
-
+ /* G_message("change=%lf, norm2=%lf",change, norm); */
}
- VEC *fraction;
+ vec_struct *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);
+ fraction = G_vector_init(A->cols, A->cols, CVEC); /* length: no. of spectra */
+ error = deviation / G_vector_norm_euclid(b_gamma);
+ fraction = G_vector_copy(startvector, NO_COMPACT);
+
+
/* write result in full percent */
for (i = 0; i < A->cols; i++) /* no. of spectra */
- result_cell[i][col] = (CELL) (100 * fraction->ve[i]);
+ result_cell[i][col] =
+ (CELL) (100 * G_matrix_get_element(fraction, i, 0));
/* save error and iterations */
error_cell[col] = (CELL) (100 * error);
iter_cell[col] = iterations;
- vec_free(fraction);
+
+ G_vector_free(fraction);
+
+
} /* end cols loop */
+
+
/* 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]);
@@ -414,6 +430,7 @@
if (iter_fd > 0)
Rast_put_c_row(iter_fd, iter_cell);
+
} /* rows loop */
G_percent(row, nrows, 2);
@@ -425,11 +442,9 @@
char command[1080];
Rast_close(resultfd[i]);
-
- vec_free(startvector);
- vec_free(A_times_startvector);
+
/* 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"
@@ -454,13 +469,17 @@
Rast_close(iter_fd);
G_matrix_free(A);
-
- vec_free(errorvector);
- vec_free(temp);
-
- vec_free(b_gamma);
+ G_vector_free(errorvector);
+ G_vector_free(temp);
+
+ G_vector_free(b_gamma);
+
+ G_vector_free(startvector);
+
+ G_vector_free(A_times_startvector);
+
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-16 10:33:35 UTC (rev 54311)
+++ grass-addons/grass7/imagery/i.spec.unmix/open.c 2012-12-16 10:44:22 UTC (rev 54312)
@@ -16,113 +16,10 @@
#include "global.h"
#include "la_extra.h"
-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";
- FILE *fp;
- int i, matrixsize;
- mat_struct A_input;
+int G_matrix_read2(FILE * fp, mat_struct * out); /* Modified version of G_matrix_read(..). */
- if ((fp = fopen (matrixfile, "r")) == NULL)
- G_fatal_error (_("Matrix file %s not found."), matrixfile);
-
-
-
- if ((G_matrix_read (fp, &A_input) < 0))
- G_fatal_error (_("Unable to read matrix file %s."), matrixfile);
- fclose (fp);
-
- #if 0
- G_message(_("Your spectral matrix = %d"), m_output(A_input));
- #endif
-
- G_matrix_clone2(&A_input, A);
-
- A = G_matrix_transpose (&A_input);
-
- if ((A->rows) < (A->cols))
- G_fatal_error (_("Need number of cols >= rows to perform least squares fitting."));
-
- matrixsize = A->cols;
-
- 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);
- }
-
- 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);
-
- 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] = Rast_allocate_c_buf();
-
- G_message (_("Opening input file no. %i [%s]"), (i + 1), Ref.file[i].name);
-
- if ((cellfd[i] = Rast_open_old (Ref.file[i].name, Ref.file[i].mapset)) < 0)
- G_fatal_error (_("Unable to open <%s>"), Ref.file[i].name);
- }
-
- result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
- resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
-
- for (i = 0; i < A->cols; i++)
- {
- sprintf (result_name, "%s.%d", result_prefix, (i + 1));
- G_message (_("Opening output file [%s]"), result_name);
-
- result_cell[i] = Rast_allocate_c_buf();
- if ((resultfd[i] = Rast_open_c_new (result_name)) < 0)
- G_fatal_error (_("GRASS-DB internal error: Unable to proceed."));
- }
-
- error_cell = (CELL *) G_malloc (sizeof(CELL *));
- if (error_name)
- {
- G_message (_("Opening error file [%s]"), error_name);
-
- if ((error_fd = Rast_open_c_new (error_name)) < 0)
- G_fatal_error (_("Unable to create error layer [%s]"), error_name);
- else
- error_cell = Rast_allocate_c_buf();
- }
-
- iter_cell = (CELL *) G_malloc (sizeof(CELL *));
- if (iter_name)
- {
- G_message (_("Opening iteration file [%s]"), iter_name);
-
- if ((iter_fd = Rast_open_c_new (iter_name)) < 0)
- G_fatal_error (_("Unable to create iterations layer [%s]"), iter_name);
- else
- iter_cell = Rast_allocate_c_buf();
- }
-
- return matrixsize;
- */
- return 0;
-}
-
-mat_struct *open_files2(char *matrixfile,
+mat_struct *open_files(char *matrixfile,
char *img_grp,
char *result_prefix,
char *iter_name, char *error_name)
@@ -251,3 +148,47 @@
/* give back number of output files (= Ref.nfiles) */
return A;
}
+
+int G_matrix_read2(FILE * fp, mat_struct * out)
+{
+ char buff[100];
+ int rows, cols;
+ int i, j, row;
+ double val;
+
+ /* skip comments */
+ for (;;) {
+ if (!G_getl(buff, sizeof(buff), fp))
+ return -1;
+ if (buff[0] != '#')
+ break;
+ }
+
+ if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) {
+ G_warning(_("Input format error1"));
+ return -1;
+ }
+
+
+ G_matrix_set(out, rows, cols, rows);
+
+
+ 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);
+ }
+ }
+
+ return 0;
+}
Deleted: grass-addons/grass7/imagery/i.spec.unmix/open.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/open.h 2012-12-16 10:33:35 UTC (rev 54311)
+++ grass-addons/grass7/imagery/i.spec.unmix/open.h 2012-12-16 10:44:22 UTC (rev 54312)
@@ -1 +0,0 @@
-mat_struct *open_files2(char *, char *, char *, char *, char *);
Modified: grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c 2012-12-16 10:33:35 UTC (rev 54311)
+++ grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c 2012-12-16 10:44:22 UTC (rev 54312)
@@ -9,8 +9,8 @@
#include <math.h>
#include <grass/gmath.h>
#include "global.h"
-#include "la_extra.h"
+
/* input mat_struct A, vec_struct Avector1, Avector2
* output cur_angle
*
@@ -37,8 +37,8 @@
/* Measure spectral angle */
/* multiply one A column with second */
-
- vtmp1 = G_vector_product(Avector1, Avector2);
+ vtmp1 = G_vector_init (Avector1->rows, Avector1->rows, CVEC);
+ vtmp1 = G_vector_product(Avector1, Avector2,vtmp1);
norm1 = G_vector_norm1(vtmp1); /* calculate 1-norm */
norm2 = G_vector_norm_euclid(Avector1); /* calculate 2-norm (Euclidean) */
norm3 = G_vector_norm_euclid(Avector2); /* calculate 2-norm (Euclidean) */
More information about the grass-commit
mailing list