[GRASS-SVN] r31528 - in grass-addons/imagery: . i.spec.sam
i.spec.unmix i.spec.unmix/sample
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon May 26 10:37:42 EDT 2008
Author: neteler
Date: 2008-05-26 10:37:42 -0400 (Mon, 26 May 2008)
New Revision: 31528
Added:
grass-addons/imagery/i.spec.sam/
grass-addons/imagery/i.spec.sam/Makefile
grass-addons/imagery/i.spec.sam/description.html
grass-addons/imagery/i.spec.sam/err.h
grass-addons/imagery/i.spec.sam/global.h
grass-addons/imagery/i.spec.sam/hist.c
grass-addons/imagery/i.spec.sam/machine.h
grass-addons/imagery/i.spec.sam/main.c
grass-addons/imagery/i.spec.sam/matrix.h
grass-addons/imagery/i.spec.sam/matrix2.h
grass-addons/imagery/i.spec.sam/meminfo.h
grass-addons/imagery/i.spec.sam/open.c
grass-addons/imagery/i.spec.sam/spec_angle.c
grass-addons/imagery/i.spec.sam/spectrum.dat
grass-addons/imagery/i.spec.sam/spectrum2.dat
grass-addons/imagery/i.spec.sam/test.dat
grass-addons/imagery/i.spec.sam/view.sh
grass-addons/imagery/i.spec.unmix/
grass-addons/imagery/i.spec.unmix/.fig
grass-addons/imagery/i.spec.unmix/Makefile
grass-addons/imagery/i.spec.unmix/README
grass-addons/imagery/i.spec.unmix/description.html
grass-addons/imagery/i.spec.unmix/global.h
grass-addons/imagery/i.spec.unmix/hist.c
grass-addons/imagery/i.spec.unmix/histogram.c
grass-addons/imagery/i.spec.unmix/main.c
grass-addons/imagery/i.spec.unmix/main.c.meschach
grass-addons/imagery/i.spec.unmix/meminfo.h
grass-addons/imagery/i.spec.unmix/meschach_convert.txt
grass-addons/imagery/i.spec.unmix/meschach_grass.txt
grass-addons/imagery/i.spec.unmix/open.c
grass-addons/imagery/i.spec.unmix/sample/
grass-addons/imagery/i.spec.unmix/sample/bdg_farben.col
grass-addons/imagery/i.spec.unmix/sample/color.sh
grass-addons/imagery/i.spec.unmix/sample/diff_farben.col
grass-addons/imagery/i.spec.unmix/sample/diplom96.dat
grass-addons/imagery/i.spec.unmix/sample/diplom97.dat
grass-addons/imagery/i.spec.unmix/sample/noshadow.sh
grass-addons/imagery/i.spec.unmix/sample/run.sh
grass-addons/imagery/i.spec.unmix/sample/spectrum.dat
grass-addons/imagery/i.spec.unmix/sample/unmix_reclass.dat
grass-addons/imagery/i.spec.unmix/sample/unmix_reclass20.dat
grass-addons/imagery/i.spec.unmix/sample/unmix_reclass25.dat
grass-addons/imagery/i.spec.unmix/spec_angle.c
Log:
added spectral analysis tools from my master thesis (Brad Douglas should have updated the code, though)
Added: grass-addons/imagery/i.spec.sam/Makefile
===================================================================
--- grass-addons/imagery/i.spec.sam/Makefile (rev 0)
+++ grass-addons/imagery/i.spec.sam/Makefile 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.spec.sam
+
+EXTRALIB=meschach.a
+LIBES = $(IMAGERYLIB) $(GMATHLIB) $(GISLIB)
+DEPENDENCIES= $(IMAGERYDEP) $(GMATHDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/imagery/i.spec.sam/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/description.html
===================================================================
--- grass-addons/imagery/i.spec.sam/description.html (rev 0)
+++ grass-addons/imagery/i.spec.sam/description.html 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,26 @@
+<H2>DESCRIPTION</H2>
+
+<p><EM>i.spec.sam</EM> is used to perform Spectral Angle Mapping.
+
+<H2>EXAMPLES</H2>
+
+<div class="code"><pre>
+TODO
+</pre></div>
+
+<H2>REFERENCES</H2>
+
+Neteler, M. (1999): Spectral Mixture Analysis von Satellitendaten zur Bestimmung von Bodenbedeckungsgraden im Hinblick auf die Erosionsmodellierung. M.Sc. thesis, University of Hannover.
+
+Neteler, M., D. Grasso, I. Michelazzi, L. Miori, S. Merler, and C. Furlanello. New image processing tools for GRASS. In Proc. Free/Libre and Open Source Software for Geoinformatics: GIS-GRASS Users Conference 2004, Sept. 12-14, Bangkok, Thailand, 2004. http://gisws.media.osaka-cu.ac.jp/grass04/viewabstract.php?id=37 (PDF)
+
+Neteler, M., D. Grasso, I. Michelazzi, L. Miori, S. Merler, and C. Furlanello, 2005. An integrated toolbox for image registration, fusion and classification. International Journal of Geoinformatics. Special Issue on FOSS/GRASS 2004 and GIS-IDEAS 2004, 1(1), pp. 51-61, March 2005.
+
+<H2>SEE ALSO</H2>
+
+<H2>AUTHOR</H2>
+
+Markus Neteler, University of Hannover.
+
+<p>
+<i>Last changed: $Date: 2008/02/10 00:00:00 $</i>
Added: grass-addons/imagery/i.spec.sam/err.h
===================================================================
--- grass-addons/imagery/i.spec.sam/err.h (rev 0)
+++ grass-addons/imagery/i.spec.sam/err.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,183 @@
+
+/**************************************************************************
+**
+** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved.
+**
+** Meschach Library
+**
+** This Meschach Library is provided "as is" without any express
+** or implied warranty of any kind with respect to this software.
+** In particular the authors shall not be liable for any direct,
+** indirect, special, incidental or consequential damages arising
+** in any way from use of the software.
+**
+** Everyone is granted permission to copy, modify and redistribute this
+** Meschach Library, provided:
+** 1. All copies contain this copyright notice.
+** 2. All modified copies shall carry a notice stating who
+** made the last modification and the date of such modification.
+** 3. No charge is made for this software or works derived from it.
+** This clause shall not be construed as constraining other software
+** distributed on the same medium as this software, nor is a
+** distribution fee considered a charge.
+**
+***************************************************************************/
+
+
+/* err.h 28/09/1993 */
+
+/* RCS id: $Id: err.h,v 1.2 1995/01/30 14:48:05 des Exp $ */
+
+
+#ifndef ERRHEADER
+#define ERRHEADER
+
+
+#include <setjmp.h>
+#include "machine.h"
+
+/* Error recovery */
+
+extern jmp_buf restart;
+
+
+/* max. # of error lists */
+#define ERR_LIST_MAX_LEN 10
+
+/* main error functions */
+#ifndef ANSI_C
+extern int ev_err(); /* main error handler */
+extern int set_err_flag(); /* for different ways of handling
+ errors, returns old value */
+extern int count_errs(); /* to avoid "too many errors" */
+extern int err_list_attach(); /* for attaching a list of errors */
+extern int err_is_list_attached(); /* checking if a list is attached */
+extern int err_list_free(); /* freeing a list of errors */
+
+#else /* ANSI_C */
+
+extern int ev_err(char *,int,int,char *,int); /* main error handler */
+extern int set_err_flag(int flag); /* for different ways of handling
+ errors, returns old value */
+extern int count_errs(int true_false); /* to avoid "too many errors" */
+extern int err_list_attach(int list_num, int list_len,
+ char **err_ptr,int warn); /* for attaching a list of errors */
+extern int err_is_list_attached(int list_num); /* checking if a list
+ is attached */
+extern int err_list_free(int list_num); /* freeing a list of errors */
+
+#endif
+
+
+/* error(E_TYPE,"myfunc") raises error type E_TYPE for function my_func() */
+#define error(err_num,fn_name) ev_err(__FILE__,err_num,__LINE__,fn_name,0)
+
+/* warning(WARN_TYPE,"myfunc") raises warning type WARN_TYPE for
+ function my_func() */
+#define warning(err_num,fn_name) ev_err(__FILE__,err_num,__LINE__,fn_name,1)
+
+
+/* error flags */
+#define EF_EXIT 0 /* exit on error */
+#define EF_ABORT 1 /* abort (dump core) on error */
+#define EF_JUMP 2 /* jump on error */
+#define EF_SILENT 3 /* jump, but don't print message */
+#define ERREXIT() set_err_flag(EF_EXIT)
+#define ERRABORT() set_err_flag(EF_ABORT)
+/* don't print message */
+#define SILENTERR() if ( ! setjmp(restart) ) set_err_flag(EF_SILENT)
+/* return here on error */
+#define ON_ERROR() if ( ! setjmp(restart) ) set_err_flag(EF_JUMP)
+
+
+/* error types */
+#define E_UNKNOWN 0
+#define E_SIZES 1
+#define E_BOUNDS 2
+#define E_MEM 3
+#define E_SING 4
+#define E_POSDEF 5
+#define E_FORMAT 6
+#define E_INPUT 7
+#define E_NULL 8
+#define E_SQUARE 9
+#define E_RANGE 10
+#define E_INSITU2 11
+#define E_INSITU 12
+#define E_ITER 13
+#define E_CONV 14
+#define E_START 15
+#define E_SIGNAL 16
+#define E_INTERN 17
+#define E_EOF 18
+#define E_SHARED_VECS 19
+#define E_NEG 20
+#define E_OVERWRITE 21
+#define E_BREAKDOWN 22
+
+/* warning types */
+#define WARN_UNKNOWN 0
+#define WARN_WRONG_TYPE 1
+#define WARN_NO_MARK 2
+#define WARN_RES_LESS_0 3
+#define WARN_SHARED_VEC 4
+
+
+/* error catching macros */
+
+/* execute err_part if error errnum is raised while executing ok_part */
+#define catch(errnum,ok_part,err_part) \
+ { jmp_buf _save; int _err_num, _old_flag; \
+ _old_flag = set_err_flag(EF_SILENT); \
+ MEM_COPY(restart,_save,sizeof(jmp_buf)); \
+ if ( (_err_num=setjmp(restart)) == 0 ) \
+ { ok_part; \
+ set_err_flag(_old_flag); \
+ MEM_COPY(_save,restart,sizeof(jmp_buf)); } \
+ else if ( _err_num == errnum ) \
+ { set_err_flag(_old_flag); \
+ MEM_COPY(_save,restart,sizeof(jmp_buf)); \
+ err_part; } \
+ else { set_err_flag(_old_flag); \
+ MEM_COPY(_save,restart,sizeof(jmp_buf)); \
+ error(_err_num,"catch"); \
+ } \
+ }
+
+
+/* execute err_part if any error raised while executing ok_part */
+#define catchall(ok_part,err_part) \
+ { jmp_buf _save; int _err_num, _old_flag; \
+ _old_flag = set_err_flag(EF_SILENT); \
+ MEM_COPY(restart,_save,sizeof(jmp_buf)); \
+ if ( (_err_num=setjmp(restart)) == 0 ) \
+ { ok_part; \
+ set_err_flag(_old_flag); \
+ MEM_COPY(_save,restart,sizeof(jmp_buf)); } \
+ else \
+ { set_err_flag(_old_flag); \
+ MEM_COPY(_save,restart,sizeof(jmp_buf)); \
+ err_part; } \
+ }
+
+
+/* print message if error raised while executing ok_part,
+ then re-raise error to trace calls */
+#define tracecatch(ok_part,function) \
+ { jmp_buf _save; int _err_num, _old_flag; \
+ _old_flag = set_err_flag(EF_JUMP); \
+ MEM_COPY(restart,_save,sizeof(jmp_buf)); \
+ if ( (_err_num=setjmp(restart)) == 0 ) \
+ { ok_part; \
+ set_err_flag(_old_flag); \
+ MEM_COPY(_save,restart,sizeof(jmp_buf)); } \
+ else \
+ { set_err_flag(_old_flag); \
+ MEM_COPY(_save,restart,sizeof(jmp_buf)); \
+ error(_err_num,function); } \
+ }
+
+
+
+#endif /* ERRHEADER */
+
Property changes on: grass-addons/imagery/i.spec.sam/err.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/global.h
===================================================================
--- grass-addons/imagery/i.spec.sam/global.h (rev 0)
+++ grass-addons/imagery/i.spec.sam/global.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,35 @@
+#include "imagery.h"
+#include <math.h>
+#include "matrix.h"
+
+#ifndef GLOBAL
+#define GLOBAL extern
+#endif
+
+#define MAXFILES 255
+
+GLOBAL MAT *A;
+GLOBAL VEC *b, *Avector;
+GLOBAL int matrixsize;
+GLOBAL float curr_angle;
+
+GLOBAL char *group;
+GLOBAL struct Ref Ref;
+
+GLOBAL CELL **cell;
+GLOBAL int *cellfd;
+
+GLOBAL CELL **result_cell;
+GLOBAL int *resultfd;
+
+GLOBAL CELL **error_cell;
+GLOBAL int error_fd;
+
+GLOBAL char result_name[80];
+GLOBAL char *result_prefix, *matrixfile;
+
+GLOBAL struct
+ {
+ struct Flag *quiet;
+ } flag;
+
Property changes on: grass-addons/imagery/i.spec.sam/global.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/hist.c
===================================================================
--- grass-addons/imagery/i.spec.sam/hist.c (rev 0)
+++ grass-addons/imagery/i.spec.sam/hist.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,14 @@
+#include "gis.h"
+
+void make_history(name, group, matrixfile)
+ char *name, *group, *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);
+ }
+}
Property changes on: grass-addons/imagery/i.spec.sam/hist.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/machine.h
===================================================================
--- grass-addons/imagery/i.spec.sam/machine.h (rev 0)
+++ grass-addons/imagery/i.spec.sam/machine.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,193 @@
+/* machine.h. Generated automatically by configure. */
+/* Any machine specific stuff goes here */
+/* Add details necessary for your own installation here! */
+
+/* This is for use with "configure" -- if you are not using configure
+ then use machine.van for the "vanilla" version of machine.h */
+
+/* Note special macros: ANSI_C (ANSI C syntax)
+ SEGMENTED (segmented memory machine e.g. MS-DOS)
+ MALLOCDECL (declared if malloc() etc have
+ been declared) */
+
+/* #undef const */
+
+/* #undef MALLOCDECL */
+#define NOT_SEGMENTED 1
+/* #undef HAVE_COMPLEX_H */
+#define HAVE_MALLOC_H 1
+#define STDC_HEADERS 1
+#define HAVE_BCOPY 1
+#define HAVE_BZERO 1
+#define CHAR0ISDBL0 1
+/* #undef WORDS_BIGENDIAN */
+#define U_INT_DEF 1
+#define VARARGS 1
+
+
+/* for basic or larger versions */
+#define COMPLEX 1
+#define SPARSE 1
+
+/* for loop unrolling */
+/* #undef VUNROLL */
+/* #undef MUNROLL */
+
+/* for segmented memory */
+#ifndef NOT_SEGMENTED
+#define SEGMENTED
+#endif
+
+/* if the system has malloc.h */
+#ifdef HAVE_MALLOC_H
+#define MALLOCDECL 1
+#include <malloc.h>
+#endif
+
+/* any compiler should have this header */
+/* if not, change it */
+#include <stdio.h>
+
+
+/* Check for ANSI C memmove and memset */
+#ifdef STDC_HEADERS
+
+/* standard copy & zero functions */
+#define MEM_COPY(from,to,size) memmove((to),(from),(size))
+#define MEM_ZERO(where,size) memset((where),'\0',(size))
+
+#ifndef ANSI_C
+#define ANSI_C 1
+#endif
+
+#endif
+
+/* standard headers */
+#ifdef ANSI_C
+#include <stdlib.h>
+#include <stddef.h>
+#include <string.h>
+#include <float.h>
+#endif
+
+
+/* if have bcopy & bzero and no alternatives yet known, use them */
+#ifdef HAVE_BCOPY
+#ifndef MEM_COPY
+/* nonstandard copy function */
+#define MEM_COPY(from,to,size) bcopy((char *)(from),(char *)(to),(int)(size))
+#endif
+#endif
+
+#ifdef HAVE_BZERO
+#ifndef MEM_ZERO
+/* nonstandard zero function */
+#define MEM_ZERO(where,size) bzero((char *)(where),(int)(size))
+#endif
+#endif
+
+/* if the system has complex.h */
+#ifdef HAVE_COMPLEX_H
+#include <complex.h>
+#endif
+
+/* If prototypes are available & ANSI_C not yet defined, then define it,
+ but don't include any header files as the proper ANSI C headers
+ aren't here */
+#define HAVE_PROTOTYPES 1
+#ifdef HAVE_PROTOTYPES
+#ifndef ANSI_C
+#define ANSI_C 1
+#endif
+#endif
+
+/* floating point precision */
+
+/* you can choose single, double or long double (if available) precision */
+
+#define FLOAT 1
+#define DOUBLE 2
+#define LONG_DOUBLE 3
+
+/* #undef REAL_FLT */
+/* #undef REAL_DBL */
+
+/* if nothing is defined, choose double precision */
+#ifndef REAL_DBL
+#ifndef REAL_FLT
+#define REAL_DBL 1
+#endif
+#endif
+
+/* single precision */
+#ifdef REAL_FLT
+#define Real float
+#define LongReal float
+#define REAL FLOAT
+#define LONGREAL FLOAT
+#endif
+
+/* double precision */
+#ifdef REAL_DBL
+#define Real double
+#define LongReal double
+#define REAL DOUBLE
+#define LONGREAL DOUBLE
+#endif
+
+
+/* machine epsilon or unit roundoff error */
+/* This is correct on most IEEE Real precision systems */
+#ifdef DBL_EPSILON
+#if REAL == DOUBLE
+#define MACHEPS DBL_EPSILON
+#elif REAL == FLOAT
+#define MACHEPS FLT_EPSILON
+#elif REAL == LONGDOUBLE
+#define MACHEPS LDBL_EPSILON
+#endif
+#endif
+
+#define F_MACHEPS 1.19209e-07
+#define D_MACHEPS 2.22045e-16
+
+#ifndef MACHEPS
+#if REAL == DOUBLE
+#define MACHEPS D_MACHEPS
+#elif REAL == FLOAT
+#define MACHEPS F_MACHEPS
+#elif REAL == LONGDOUBLE
+#define MACHEPS D_MACHEPS
+#endif
+#endif
+
+/* #undef M_MACHEPS */
+
+/********************
+#ifdef DBL_EPSILON
+#define MACHEPS DBL_EPSILON
+#endif
+#ifdef M_MACHEPS
+#ifndef MACHEPS
+#define MACHEPS M_MACHEPS
+#endif
+#endif
+********************/
+
+#define M_MAX_INT 2147483647
+#ifdef M_MAX_INT
+#ifndef MAX_RAND
+#define MAX_RAND ((double)(M_MAX_INT))
+#endif
+#endif
+
+/* for non-ANSI systems */
+#ifndef HUGE_VAL
+#define HUGE_VAL HUGE
+#endif
+
+
+#ifdef ANSI_C
+extern int isatty(int);
+#endif
+
Property changes on: grass-addons/imagery/i.spec.sam/machine.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/main.c
===================================================================
--- grass-addons/imagery/i.spec.sam/main.c (rev 0)
+++ grass-addons/imagery/i.spec.sam/main.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,225 @@
+/* Spectral angle mapping
+ * (c) 1998 Markus Neteler, Hannover, Germany
+ * neteler at geog.uni-hannover.de
+ *
+ * V 0.1 - 26. Oct.1998
+ *
+ * * Based on Meschach Library (matrix operations)
+ * * Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+ *
+ * Cited references are from
+ * Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
+ * Proceedings of the centre for Mathematics and its Applicaions.
+ * The Australian National University. Vol. 32.
+ * ISBN 0 7315 1900 0
+ *
+ ********************************************************************/
+
+#define GLOBAL
+#include "global.h"
+#include <stdio.h>
+#include <strings.h>
+#include <math.h>
+#include "matrix.h"
+#include "matrix2.h"
+
+
+
+int main(argc,argv)
+char *argv[];
+{
+ int nrows, ncols;
+ int row, col;
+ int band;
+ int i, j, error=0;
+ VEC *svd_values;
+ char command[80];
+ float anglefield[255][255];
+ struct
+ {
+ struct Option *group, *matrixfile, *result;
+ } parm;
+
+ G_gisinit (argv[0]);
+
+ parm.group = G_define_option();
+ parm.group->key = "group";
+ parm.group->type = TYPE_STRING;
+ parm.group->required = YES;
+ parm.group->description = "Imagery group containing images to be analyzed with Spectral Mixture Analyis";
+
+ parm.matrixfile = G_define_option();
+ parm.matrixfile->key = "matrixfile";
+ parm.matrixfile->type = TYPE_STRING;
+ parm.matrixfile->required = YES;
+ parm.matrixfile->description = "Matrix file containing spectral signatures ";
+
+ parm.result = G_define_option();
+ parm.result->key = "result";
+ parm.result->type = TYPE_STRING;
+ parm.result->required = YES;
+ parm.result->description = "Raster map prefix to hold spectral angles";
+
+ flag.quiet = G_define_flag();
+ flag.quiet->key = 'q';
+ flag.quiet->description = "Run quietly";
+
+
+ if (G_parser(argc,argv))
+ exit(1);
+
+ result_prefix = parm.result->answer;
+ group = parm.group->answer;
+ matrixfile = parm.matrixfile->answer;
+
+
+/* here we go... */
+
+ open_files();
+ /* Spectral Matrix is stored in A now */
+
+ /* 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.
+ */
+
+ for (i = 0; i < Ref.nfiles; i++) /* Ref.nfiles = matrixsize*/
+ {
+ Avector = get_row(A, i, VNULL); /* go columnwise through matrix*/
+ for (j = 0; j < Ref.nfiles ; j++)
+ {
+ if (j !=i)
+ {
+ b = get_row(A, j, VNULL); /* compare with next col in A */
+ spectral_angle();
+ anglefield[i][j]= curr_angle;
+ V_FREE(b);
+ }
+ }
+ V_FREE(Avector);
+ }
+
+ /* print out the result */
+ fprintf(stderr,"Orthogonality check of Matrix A:\n");
+ for (i = 0; i < Ref.nfiles ; i++)
+ for (j = 0; j < Ref.nfiles ; j++)
+ {
+ if (j !=i)
+ fprintf(stderr," Angle between row %i and row %i: %g\n", (i+1), (j+1), anglefield[i][j]);
+ }
+ fprintf(stderr,"\n");
+
+ /* check it */
+ for (i = 0; i < Ref.nfiles ; i++)
+ for (j = 0; j < Ref.nfiles ; j++)
+ if (j !=i)
+ if (anglefield[i][j] < 10.0)
+ {
+ fprintf(stderr,"ERROR: Spectral entries row %i|%i in your matrix are linear dependent!\n",i,j);
+ error=1;
+ }
+
+ if (!error)
+ fprintf(stderr,"Spectral matrix is o.k. Proceeding...\n");
+
+/* check singular values of Matrix A
+ * Ref: Boardman, J.W. 1989: Inversion of imaging spectrometry data
+ * using singular value decomposition. IGARSS 1989: 12th Canadian
+ * symposium on Remote Sensing. Vol.4 pp.2069-2072
+ */
+ fprintf(stderr,"\n");
+ fprintf(stderr,"Singular values of Matrix A:");
+ svd_values = svd(A, MNULL, MNULL, VNULL);
+ v_output(svd_values);
+ fprintf(stderr,"\n");
+ if (error)
+ {
+ fprintf(stderr,"Exiting...\n");
+ exit(-1);
+ }
+
+ /* alright, start Spectral angle mapping */
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
+ if (!flag.quiet->answer)
+ {
+ fprintf (stderr, "Calculating for %i x %i = %i pixels:\n",nrows,ncols, (ncols * ncols));
+ fprintf (stderr, "%s ... ", G_program_name());
+ }
+ for (row = 0; row < nrows; row++) /* rows loop*/
+ {
+ if (!flag.quiet->answer)
+ G_percent(row, nrows, 2);
+ for (band = 0; band < Ref.nfiles; band++) /* get row for all bands*/
+ {
+ if (G_get_map_row (cellfd[band], cell[band], row) < 0)
+ exit(1);
+ }
+
+ for (col = 0; col < ncols; col++) /* cols loop, work pixelwise for all bands */
+ {
+ /* get pixel values of each band and store in b vector: */
+ b = v_get(A->m); /* dimension of vector = matrix size = Ref.nfiles*/
+ for (band = 0; band < Ref.nfiles; band++)
+ b->ve[band] = cell[band][col]; /* read input vector */
+
+ /* calculate spectral angle for current pixel
+ * between pixel spectrum and reference spectrum
+ * and write result in full degree */
+
+ for (i = 0; i < Ref.nfiles; i++) /* Ref.nfiles = matrixsize*/
+ {
+ Avector = get_row(A, i, VNULL); /* go row-wise through matrix*/
+ spectral_angle();
+ result_cell[i][col] = round (curr_angle);
+ V_FREE(Avector);
+ }
+
+ V_FREE(b);
+
+ } /* columns loop */
+
+ /* write the resulting rows: */
+ for (i = 0; i < Ref.nfiles; i++)
+ G_put_map_row (resultfd[i], result_cell[i], row);
+
+ } /* rows loop */
+
+ if (!flag.quiet->answer)
+ G_percent(row, nrows, 2);
+
+ /* close files */
+ for (i = 0; i < Ref.nfiles; i++)
+ {
+ G_close_cell (resultfd[i]);
+ G_unopen_cell(cellfd[i]);
+ /* make grey scale color table */
+ sprintf(result_name, "%s.%d", result_prefix, (i+1));
+ sprintf(command, "r.colors map=%s color=grey >/dev/null", result_name);
+ system(command);
+ /* write a color table */
+ }
+
+ M_FREE(A);
+ make_history(result_name, group, matrixfile);
+ exit(0);
+} /* main*/
+
+
+CELL round (x)
+ double x;
+ {
+ CELL n;
+
+ if (x >= 0.0)
+ n = x + .5;
+ else
+ {
+ n = -x + .5;
+ n = -n;
+ }
+ return n;
+ }
+
Property changes on: grass-addons/imagery/i.spec.sam/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/matrix.h
===================================================================
--- grass-addons/imagery/i.spec.sam/matrix.h (rev 0)
+++ grass-addons/imagery/i.spec.sam/matrix.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,665 @@
+
+/**************************************************************************
+**
+** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+**
+** Meschach Library
+**
+** This Meschach Library is provided "as is" without any express
+** or implied warranty of any kind with respect to this software.
+** In particular the authors shall not be liable for any direct,
+** indirect, special, incidental or consequential damages arising
+** in any way from use of the software.
+**
+** Everyone is granted permission to copy, modify and redistribute this
+** Meschach Library, provided:
+** 1. All copies contain this copyright notice.
+** 2. All modified copies shall carry a notice stating who
+** made the last modification and the date of such modification.
+** 3. No charge is made for this software or works derived from it.
+** This clause shall not be construed as constraining other software
+** distributed on the same medium as this software, nor is a
+** distribution fee considered a charge.
+**
+***************************************************************************/
+
+
+/*
+ Type definitions for general purpose maths package
+*/
+
+#ifndef MATRIXH
+
+/* RCS id: $Id: matrix.h,v 1.18 1994/04/16 00:33:37 des Exp $ */
+
+#define MATRIXH
+
+#include "machine.h"
+#include "err.h"
+#include "meminfo.h"
+
+/* unsigned integer type */
+#ifndef U_INT_DEF
+typedef unsigned int u_int;
+#define U_INT_DEF
+#endif
+
+/* vector definition */
+typedef struct {
+ u_int dim, max_dim;
+ Real *ve;
+ } VEC;
+
+/* matrix definition */
+typedef struct {
+ u_int m, n;
+ u_int max_m, max_n, max_size;
+ Real **me,*base; /* base is base of alloc'd mem */
+ } MAT;
+
+/* band matrix definition */
+typedef struct {
+ MAT *mat; /* matrix */
+ int lb,ub; /* lower and upper bandwidth */
+ } BAND;
+
+
+/* permutation definition */
+typedef struct {
+ u_int size, max_size, *pe;
+ } PERM;
+
+/* integer vector definition */
+typedef struct {
+ u_int dim, max_dim;
+ int *ive;
+ } IVEC;
+
+
+#ifndef MALLOCDECL
+#ifndef ANSI_C
+extern char *malloc(), *calloc(), *realloc();
+#else
+extern void *malloc(size_t),
+ *calloc(size_t,size_t),
+ *realloc(void *,size_t);
+#endif
+#endif
+
+#ifndef ANSI_C
+extern void m_version();
+#else
+void m_version( void );
+#endif
+
+#ifndef ANSI_C
+/* allocate one object of given type */
+#define NEW(type) ((type *)calloc(1,sizeof(type)))
+
+/* allocate num objects of given type */
+#define NEW_A(num,type) ((type *)calloc((unsigned)(num),sizeof(type)))
+
+ /* re-allocate arry to have num objects of the given type */
+#define RENEW(var,num,type) \
+ ((var)=(type *)((var) ? \
+ realloc((char *)(var),(unsigned)(num)*sizeof(type)) : \
+ calloc((unsigned)(num),sizeof(type))))
+
+#define MEMCOPY(from,to,n_items,type) \
+ MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type))
+
+#else
+/* allocate one object of given type */
+#define NEW(type) ((type *)calloc((size_t)1,(size_t)sizeof(type)))
+
+/* allocate num objects of given type */
+#define NEW_A(num,type) ((type *)calloc((size_t)(num),(size_t)sizeof(type)))
+
+ /* re-allocate arry to have num objects of the given type */
+#define RENEW(var,num,type) \
+ ((var)=(type *)((var) ? \
+ realloc((char *)(var),(size_t)((num)*sizeof(type))) : \
+ calloc((size_t)(num),(size_t)sizeof(type))))
+
+#define MEMCOPY(from,to,n_items,type) \
+ MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type))
+
+#endif
+
+/* type independent min and max operations */
+#ifndef max
+#define max(a,b) ((a) > (b) ? (a) : (b))
+#endif
+#ifndef min
+#define min(a,b) ((a) > (b) ? (b) : (a))
+#endif
+
+
+#undef TRUE
+#define TRUE 1
+#undef FALSE
+#define FALSE 0
+
+
+/* for input routines */
+#define MAXLINE 81
+
+
+/* Dynamic memory allocation */
+
+/* Should use M_FREE/V_FREE/PX_FREE in programs instead of m/v/px_free()
+ as this is considerably safer -- also provides a simple type check ! */
+
+#ifndef ANSI_C
+
+extern VEC *v_get(), *v_resize();
+extern MAT *m_get(), *m_resize();
+extern PERM *px_get(), *px_resize();
+extern IVEC *iv_get(), *iv_resize();
+extern int m_free(),v_free();
+extern int px_free();
+extern int iv_free();
+extern BAND *bd_get(), *bd_resize();
+extern int bd_free();
+
+#else
+
+/* get/resize vector to given dimension */
+extern VEC *v_get(int), *v_resize(VEC *,int);
+/* get/resize matrix to be m x n */
+extern MAT *m_get(int,int), *m_resize(MAT *,int,int);
+/* get/resize permutation to have the given size */
+extern PERM *px_get(int), *px_resize(PERM *,int);
+/* get/resize an integer vector to given dimension */
+extern IVEC *iv_get(int), *iv_resize(IVEC *,int);
+/* get/resize a band matrix to given dimension */
+extern BAND *bd_get(int,int,int), *bd_resize(BAND *,int,int,int);
+
+/* free (de-allocate) (band) matrices, vectors, permutations and
+ integer vectors */
+extern int iv_free(IVEC *);
+extern m_free(MAT *),v_free(VEC *),px_free(PERM *);
+extern int bd_free(BAND *);
+
+#endif
+
+
+/* MACROS */
+
+/* macros that also check types and sets pointers to NULL */
+#define M_FREE(mat) ( m_free(mat), (mat)=(MAT *)NULL )
+#define V_FREE(vec) ( v_free(vec), (vec)=(VEC *)NULL )
+#define PX_FREE(px) ( px_free(px), (px)=(PERM *)NULL )
+#define IV_FREE(iv) ( iv_free(iv), (iv)=(IVEC *)NULL )
+
+#define MAXDIM 2001
+
+
+/* Entry level access to data structures */
+#ifdef DEBUG
+
+/* returns x[i] */
+#define v_entry(x,i) (((i) < 0 || (i) >= (x)->dim) ? \
+ error(E_BOUNDS,"v_entry"), 0.0 : (x)->ve[i] )
+
+/* x[i] <- val */
+#define v_set_val(x,i,val) ((x)->ve[i] = ((i) < 0 || (i) >= (x)->dim) ? \
+ error(E_BOUNDS,"v_set_val"), 0.0 : (val))
+
+/* x[i] <- x[i] + val */
+#define v_add_val(x,i,val) ((x)->ve[i] += ((i) < 0 || (i) >= (x)->dim) ? \
+ error(E_BOUNDS,"v_add_val"), 0.0 : (val))
+
+/* x[i] <- x[i] - val */
+#define v_sub_val(x,i,val) ((x)->ve[i] -= ((i) < 0 || (i) >= (x)->dim) ? \
+ error(E_BOUNDS,"v_sub_val"), 0.0 : (val))
+
+/* returns A[i][j] */
+#define m_entry(A,i,j) (((i) < 0 || (i) >= (A)->m || \
+ (j) < 0 || (j) >= (A)->n) ? \
+ error(E_BOUNDS,"m_entry"), 0.0 : (A)->me[i][j] )
+
+/* A[i][j] <- val */
+#define m_set_val(A,i,j,val) ((A)->me[i][j] = ((i) < 0 || (i) >= (A)->m || \
+ (j) < 0 || (j) >= (A)->n) ? \
+ error(E_BOUNDS,"m_set_val"), 0.0 : (val) )
+
+/* A[i][j] <- A[i][j] + val */
+#define m_add_val(A,i,j,val) ((A)->me[i][j] += ((i) < 0 || (i) >= (A)->m || \
+ (j) < 0 || (j) >= (A)->n) ? \
+ error(E_BOUNDS,"m_add_val"), 0.0 : (val) )
+
+/* A[i][j] <- A[i][j] - val */
+#define m_sub_val(A,i,j,val) ((A)->me[i][j] -= ((i) < 0 || (i) >= (A)->m || \
+ (j) < 0 || (j) >= (A)->n) ? \
+ error(E_BOUNDS,"m_sub_val"), 0.0 : (val) )
+#else
+
+/* returns x[i] */
+#define v_entry(x,i) ((x)->ve[i])
+
+/* x[i] <- val */
+#define v_set_val(x,i,val) ((x)->ve[i] = (val))
+
+/* x[i] <- x[i] + val */
+#define v_add_val(x,i,val) ((x)->ve[i] += (val))
+
+ /* x[i] <- x[i] - val */
+#define v_sub_val(x,i,val) ((x)->ve[i] -= (val))
+
+/* returns A[i][j] */
+#define m_entry(A,i,j) ((A)->me[i][j])
+
+/* A[i][j] <- val */
+#define m_set_val(A,i,j,val) ((A)->me[i][j] = (val) )
+
+/* A[i][j] <- A[i][j] + val */
+#define m_add_val(A,i,j,val) ((A)->me[i][j] += (val) )
+
+/* A[i][j] <- A[i][j] - val */
+#define m_sub_val(A,i,j,val) ((A)->me[i][j] -= (val) )
+
+#endif
+
+
+/* I/O routines */
+#ifndef ANSI_C
+
+extern void v_foutput(),m_foutput(),px_foutput();
+extern void iv_foutput();
+extern VEC *v_finput();
+extern MAT *m_finput();
+extern PERM *px_finput();
+extern IVEC *iv_finput();
+extern int fy_or_n(), fin_int(), yn_dflt(), skipjunk();
+extern double fin_double();
+
+#else
+
+/* print x on file fp */
+void v_foutput(FILE *fp,VEC *x),
+ /* print A on file fp */
+ m_foutput(FILE *fp,MAT *A),
+ /* print px on file fp */
+ px_foutput(FILE *fp,PERM *px);
+/* print ix on file fp */
+void iv_foutput(FILE *fp,IVEC *ix);
+
+/* Note: if out is NULL, then returned object is newly allocated;
+ Also: if out is not NULL, then that size is assumed */
+
+/* read in vector from fp */
+VEC *v_finput(FILE *fp,VEC *out);
+/* read in matrix from fp */
+MAT *m_finput(FILE *fp,MAT *out);
+/* read in permutation from fp */
+PERM *px_finput(FILE *fp,PERM *out);
+/* read in int vector from fp */
+IVEC *iv_finput(FILE *fp,IVEC *out);
+
+/* fy_or_n -- yes-or-no to question in string s
+ -- question written to stderr, input from fp
+ -- if fp is NOT a tty then return y_n_dflt */
+int fy_or_n(FILE *fp,char *s);
+
+/* yn_dflt -- sets the value of y_n_dflt to val */
+int yn_dflt(int val);
+
+/* fin_int -- return integer read from file/stream fp
+ -- prompt s on stderr if fp is a tty
+ -- check that x lies between low and high: re-prompt if
+ fp is a tty, error exit otherwise
+ -- ignore check if low > high */
+int fin_int(FILE *fp,char *s,int low,int high);
+
+/* fin_double -- return double read from file/stream fp
+ -- prompt s on stderr if fp is a tty
+ -- check that x lies between low and high: re-prompt if
+ fp is a tty, error exit otherwise
+ -- ignore check if low > high */
+double fin_double(FILE *fp,char *s,double low,double high);
+
+/* it skips white spaces and strings of the form #....\n
+ Here .... is a comment string */
+int skipjunk(FILE *fp);
+
+#endif
+
+
+/* MACROS */
+
+/* macros to use stdout and stdin instead of explicit fp */
+#define v_output(vec) v_foutput(stdout,vec)
+#define v_input(vec) v_finput(stdin,vec)
+#define m_output(mat) m_foutput(stdout,mat)
+#define m_input(mat) m_finput(stdin,mat)
+#define px_output(px) px_foutput(stdout,px)
+#define px_input(px) px_finput(stdin,px)
+#define iv_output(iv) iv_foutput(stdout,iv)
+#define iv_input(iv) iv_finput(stdin,iv)
+
+/* general purpose input routine; skips comments # ... \n */
+#define finput(fp,prompt,fmt,var) \
+ ( ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) ), \
+ fscanf(fp,fmt,var) )
+#define input(prompt,fmt,var) finput(stdin,prompt,fmt,var)
+#define fprompter(fp,prompt) \
+ ( isatty(fileno(fp)) ? fprintf(stderr,prompt) : skipjunk(fp) )
+#define prompter(prompt) fprompter(stdin,prompt)
+#define y_or_n(s) fy_or_n(stdin,s)
+#define in_int(s,lo,hi) fin_int(stdin,s,lo,hi)
+#define in_double(s,lo,hi) fin_double(stdin,s,lo,hi)
+
+/* Copying routines */
+#ifndef ANSI_C
+extern MAT *_m_copy(), *m_move(), *vm_move();
+extern VEC *_v_copy(), *v_move(), *mv_move();
+extern PERM *px_copy();
+extern IVEC *iv_copy(), *iv_move();
+extern BAND *bd_copy();
+
+#else
+
+/* copy in to out starting at out[i0][j0] */
+extern MAT *_m_copy(MAT *in,MAT *out,u_int i0,u_int j0),
+ * m_move(MAT *in, int, int, int, int, MAT *out, int, int),
+ *vm_move(VEC *in, int, MAT *out, int, int, int, int);
+/* copy in to out starting at out[i0] */
+extern VEC *_v_copy(VEC *in,VEC *out,u_int i0),
+ * v_move(VEC *in, int, int, VEC *out, int),
+ *mv_move(MAT *in, int, int, int, int, VEC *out, int);
+extern PERM *px_copy(PERM *in,PERM *out);
+extern IVEC *iv_copy(IVEC *in,IVEC *out),
+ *iv_move(IVEC *in, int, int, IVEC *out, int);
+extern BAND *bd_copy(BAND *in,BAND *out);
+
+#endif
+
+
+/* MACROS */
+#define m_copy(in,out) _m_copy(in,out,0,0)
+#define v_copy(in,out) _v_copy(in,out,0)
+
+
+/* Initialisation routines -- to be zero, ones, random or identity */
+#ifndef ANSI_C
+extern VEC *v_zero(), *v_rand(), *v_ones();
+extern MAT *m_zero(), *m_ident(), *m_rand(), *m_ones();
+extern PERM *px_ident();
+extern IVEC *iv_zero();
+#else
+extern VEC *v_zero(VEC *), *v_rand(VEC *), *v_ones(VEC *);
+extern MAT *m_zero(MAT *), *m_ident(MAT *), *m_rand(MAT *),
+ *m_ones(MAT *);
+extern PERM *px_ident(PERM *);
+extern IVEC *iv_zero(IVEC *);
+#endif
+
+/* Basic vector operations */
+#ifndef ANSI_C
+extern VEC *sv_mlt(), *mv_mlt(), *vm_mlt(), *v_add(), *v_sub(),
+ *px_vec(), *pxinv_vec(), *v_mltadd(), *v_map(), *_v_map(),
+ *v_lincomb(), *v_linlist();
+extern double v_min(), v_max(), v_sum();
+extern VEC *v_star(), *v_slash(), *v_sort();
+extern double _in_prod(), __ip__();
+extern void __mltadd__(), __add__(), __sub__(),
+ __smlt__(), __zero__();
+#else
+
+extern VEC *sv_mlt(double,VEC *,VEC *), /* out <- s.x */
+ *mv_mlt(MAT *,VEC *,VEC *), /* out <- A.x */
+ *vm_mlt(MAT *,VEC *,VEC *), /* out^T <- x^T.A */
+ *v_add(VEC *,VEC *,VEC *), /* out <- x + y */
+ *v_sub(VEC *,VEC *,VEC *), /* out <- x - y */
+ *px_vec(PERM *,VEC *,VEC *), /* out <- P.x */
+ *pxinv_vec(PERM *,VEC *,VEC *), /* out <- P^{-1}.x */
+ *v_mltadd(VEC *,VEC *,double,VEC *), /* out <- x + s.y */
+#ifdef PROTOTYPES_IN_STRUCT
+ *v_map(double (*f)(double),VEC *,VEC *),
+ /* out[i] <- f(x[i]) */
+ *_v_map(double (*f)(void *,double),void *,VEC *,VEC *),
+#else
+ *v_map(double (*f)(),VEC *,VEC *), /* out[i] <- f(x[i]) */
+ *_v_map(double (*f)(),void *,VEC *,VEC *),
+#endif
+ *v_lincomb(int,VEC **,Real *,VEC *),
+ /* out <- sum_i s[i].x[i] */
+ *v_linlist(VEC *out,VEC *v1,double a1,...);
+ /* out <- s1.x1 + s2.x2 + ... */
+
+/* returns min_j x[j] (== x[i]) */
+extern double v_min(VEC *, int *),
+ /* returns max_j x[j] (== x[i]) */
+ v_max(VEC *, int *),
+ /* returns sum_i x[i] */
+ v_sum(VEC *);
+
+/* Hadamard product: out[i] <- x[i].y[i] */
+extern VEC *v_star(VEC *, VEC *, VEC *),
+ /* out[i] <- x[i] / y[i] */
+ *v_slash(VEC *, VEC *, VEC *),
+ /* sorts x, and sets order so that sorted x[i] = x[order[i]] */
+ *v_sort(VEC *, PERM *);
+
+/* returns inner product starting at component i0 */
+extern double _in_prod(VEC *x,VEC *y,u_int i0),
+ /* returns sum_{i=0}^{len-1} x[i].y[i] */
+ __ip__(Real *,Real *,int);
+
+/* see v_mltadd(), v_add(), v_sub() and v_zero() */
+extern void __mltadd__(Real *,Real *,double,int),
+ __add__(Real *,Real *,Real *,int),
+ __sub__(Real *,Real *,Real *,int),
+ __smlt__(Real *,double,Real *,int),
+ __zero__(Real *,int);
+
+#endif
+
+
+/* MACRO */
+/* usual way of computing the inner product */
+#define in_prod(a,b) _in_prod(a,b,0)
+
+/* Norms */
+/* scaled vector norms -- scale == NULL implies unscaled */
+#ifndef ANSI_C
+
+extern double _v_norm1(), _v_norm2(), _v_norm_inf(),
+ m_norm1(), m_norm_inf(), m_norm_frob();
+
+#else
+ /* returns sum_i |x[i]/scale[i]| */
+extern double _v_norm1(VEC *x,VEC *scale),
+ /* returns (scaled) Euclidean norm */
+ _v_norm2(VEC *x,VEC *scale),
+ /* returns max_i |x[i]/scale[i]| */
+ _v_norm_inf(VEC *x,VEC *scale);
+
+/* unscaled matrix norms */
+extern double m_norm1(MAT *A), m_norm_inf(MAT *A), m_norm_frob(MAT *A);
+
+#endif
+
+
+/* MACROS */
+/* unscaled vector norms */
+#define v_norm1(x) _v_norm1(x,VNULL)
+#define v_norm2(x) _v_norm2(x,VNULL)
+#define v_norm_inf(x) _v_norm_inf(x,VNULL)
+
+/* Basic matrix operations */
+#ifndef ANSI_C
+
+extern MAT *sm_mlt(), *m_mlt(), *mmtr_mlt(), *mtrm_mlt(), *m_add(), *m_sub(),
+ *sub_mat(), *m_transp(), *ms_mltadd();
+
+extern BAND *bd_transp();
+extern MAT *px_rows(), *px_cols(), *swap_rows(), *swap_cols(),
+ *_set_row(), *_set_col();
+extern VEC *get_row(), *get_col(), *sub_vec(),
+ *mv_mltadd(), *vm_mltadd();
+
+#else
+
+extern MAT *sm_mlt(double s,MAT *A,MAT *out), /* out <- s.A */
+ *m_mlt(MAT *A,MAT *B,MAT *out), /* out <- A.B */
+ *mmtr_mlt(MAT *A,MAT *B,MAT *out), /* out <- A.B^T */
+ *mtrm_mlt(MAT *A,MAT *B,MAT *out), /* out <- A^T.B */
+ *m_add(MAT *A,MAT *B,MAT *out), /* out <- A + B */
+ *m_sub(MAT *A,MAT *B,MAT *out), /* out <- A - B */
+ *sub_mat(MAT *A,u_int,u_int,u_int,u_int,MAT *out),
+ *m_transp(MAT *A,MAT *out), /* out <- A^T */
+ /* out <- A + s.B */
+ *ms_mltadd(MAT *A,MAT *B,double s,MAT *out);
+
+
+extern BAND *bd_transp(BAND *in, BAND *out); /* out <- A^T */
+extern MAT *px_rows(PERM *px,MAT *A,MAT *out), /* out <- P.A */
+ *px_cols(PERM *px,MAT *A,MAT *out), /* out <- A.P^T */
+ *swap_rows(MAT *,int,int,int,int),
+ *swap_cols(MAT *,int,int,int,int),
+ /* A[i][j] <- out[j], j >= j0 */
+ *_set_col(MAT *A,u_int i,VEC *out,u_int j0),
+ /* A[i][j] <- out[i], i >= i0 */
+ *_set_row(MAT *A,u_int j,VEC *out,u_int i0);
+
+extern VEC *get_row(MAT *,u_int,VEC *),
+ *get_col(MAT *,u_int,VEC *),
+ *sub_vec(VEC *,int,int,VEC *),
+ /* out <- x + s.A.y */
+ *mv_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out),
+ /* out^T <- x^T + s.y^T.A */
+ *vm_mltadd(VEC *x,VEC *y,MAT *A,double s,VEC *out);
+#endif
+
+
+/* MACROS */
+/* row i of A <- vec */
+#define set_row(mat,row,vec) _set_row(mat,row,vec,0)
+/* col j of A <- vec */
+#define set_col(mat,col,vec) _set_col(mat,col,vec,0)
+
+
+/* Basic permutation operations */
+#ifndef ANSI_C
+
+extern PERM *px_mlt(), *px_inv(), *px_transp();
+extern int px_sign();
+
+#else
+
+extern PERM *px_mlt(PERM *px1,PERM *px2,PERM *out), /* out <- px1.px2 */
+ *px_inv(PERM *px,PERM *out), /* out <- px^{-1} */
+ /* swap px[i] and px[j] */
+ *px_transp(PERM *px,u_int i,u_int j);
+
+ /* returns sign(px) = +1 if px product of even # transpositions
+ -1 if ps product of odd # transpositions */
+extern int px_sign(PERM *);
+
+#endif
+
+
+/* Basic integer vector operations */
+#ifndef ANSI_C
+
+extern IVEC *iv_add(), *iv_sub(), *iv_sort();
+
+#else
+
+extern IVEC *iv_add(IVEC *ix,IVEC *iy,IVEC *out), /* out <- ix + iy */
+ *iv_sub(IVEC *ix,IVEC *iy,IVEC *out), /* out <- ix - iy */
+ /* sorts ix & sets order so that sorted ix[i] = old ix[order[i]] */
+ *iv_sort(IVEC *ix, PERM *order);
+
+#endif
+
+
+/* miscellaneous functions */
+
+#ifndef ANSI_C
+
+extern double square(), cube(), mrand();
+extern void smrand(), mrandlist();
+extern void m_dump(), px_dump(), v_dump(), iv_dump();
+extern MAT *band2mat();
+extern BAND *mat2band();
+
+#else
+
+double square(double x), /* returns x^2 */
+ cube(double x), /* returns x^3 */
+ mrand(void); /* returns random # in [0,1) */
+
+void smrand(int seed), /* seeds mrand() */
+ mrandlist(Real *x, int len); /* generates len random numbers */
+
+void m_dump(FILE *fp,MAT *a), px_dump(FILE *,PERM *px),
+ v_dump(FILE *fp,VEC *x), iv_dump(FILE *fp, IVEC *ix);
+
+MAT *band2mat(BAND *bA, MAT *A);
+BAND *mat2band(MAT *A, int lb,int ub, BAND *bA);
+
+#endif
+
+
+/* miscellaneous constants */
+#define VNULL ((VEC *)NULL)
+#define MNULL ((MAT *)NULL)
+#define PNULL ((PERM *)NULL)
+#define IVNULL ((IVEC *)NULL)
+#define BDNULL ((BAND *)NULL)
+
+
+
+/* varying number of arguments */
+
+#ifdef ANSI_C
+#include <stdarg.h>
+
+/* prototypes */
+
+int v_get_vars(int dim,...);
+int iv_get_vars(int dim,...);
+int m_get_vars(int m,int n,...);
+int px_get_vars(int dim,...);
+
+int v_resize_vars(int new_dim,...);
+int iv_resize_vars(int new_dim,...);
+int m_resize_vars(int m,int n,...);
+int px_resize_vars(int new_dim,...);
+
+int v_free_vars(VEC **,...);
+int iv_free_vars(IVEC **,...);
+int px_free_vars(PERM **,...);
+int m_free_vars(MAT **,...);
+
+#elif VARARGS
+/* old varargs is used */
+
+#include <varargs.h>
+
+/* prototypes */
+
+int v_get_vars();
+int iv_get_vars();
+int m_get_vars();
+int px_get_vars();
+
+int v_resize_vars();
+int iv_resize_vars();
+int m_resize_vars();
+int px_resize_vars();
+
+int v_free_vars();
+int iv_free_vars();
+int px_free_vars();
+int m_free_vars();
+
+#endif
+
+
+#endif
+
+
Property changes on: grass-addons/imagery/i.spec.sam/matrix.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/matrix2.h
===================================================================
--- grass-addons/imagery/i.spec.sam/matrix2.h (rev 0)
+++ grass-addons/imagery/i.spec.sam/matrix2.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,229 @@
+
+/**************************************************************************
+**
+** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+**
+** Meschach Library
+**
+** This Meschach Library is provided "as is" without any express
+** or implied warranty of any kind with respect to this software.
+** In particular the authors shall not be liable for any direct,
+** indirect, special, incidental or consequential damages arising
+** in any way from use of the software.
+**
+** Everyone is granted permission to copy, modify and redistribute this
+** Meschach Library, provided:
+** 1. All copies contain this copyright notice.
+** 2. All modified copies shall carry a notice stating who
+** made the last modification and the date of such modification.
+** 3. No charge is made for this software or works derived from it.
+** This clause shall not be construed as constraining other software
+** distributed on the same medium as this software, nor is a
+** distribution fee considered a charge.
+**
+***************************************************************************/
+
+
+/*
+ Header file for ``matrix2.a'' library file
+*/
+
+
+#ifndef MATRIX2H
+#define MATRIX2H
+
+#include "matrix.h"
+
+/* Unless otherwise specified, factorisation routines overwrite the
+ matrix that is being factorised */
+
+#ifndef ANSI_C
+
+extern MAT *BKPfactor(), *CHfactor(), *LUfactor(), *QRfactor(),
+ *QRCPfactor(), *LDLfactor(), *Hfactor(), *MCHfactor(),
+ *m_inverse();
+extern double LUcondest(), QRcondest();
+extern MAT *makeQ(), *makeR(), *makeHQ(), *makeH();
+extern MAT *LDLupdate(), *QRupdate();
+
+extern VEC *BKPsolve(), *CHsolve(), *LUsolve(), *_Qsolve(), *QRsolve(),
+ *LDLsolve(), *Usolve(), *Lsolve(), *Dsolve(), *LTsolve(),
+ *UTsolve(), *LUTsolve(), *QRCPsolve();
+
+extern BAND *bdLUfactor(), *bdLDLfactor();
+extern VEC *bdLUsolve(), *bdLDLsolve();
+
+extern VEC *hhvec();
+extern VEC *hhtrvec();
+extern MAT *hhtrrows();
+extern MAT *hhtrcols();
+
+extern void givens();
+extern VEC *rot_vec(); /* in situ */
+extern MAT *rot_rows(); /* in situ */
+extern MAT *rot_cols(); /* in situ */
+
+
+/* eigenvalue routines */
+extern VEC *trieig(), *symmeig();
+extern MAT *schur();
+extern void schur_evals();
+extern MAT *schur_vecs();
+
+/* singular value decomposition */
+extern VEC *bisvd(), *svd();
+
+/* matrix powers and exponent */
+MAT *_m_pow();
+MAT *m_pow();
+MAT *m_exp(), *_m_exp();
+MAT *m_poly();
+
+/* FFT */
+void fft();
+void ifft();
+
+
+#else
+
+ /* forms Bunch-Kaufman-Parlett factorisation for
+ symmetric indefinite matrices */
+extern MAT *BKPfactor(MAT *A,PERM *pivot,PERM *blocks),
+ /* Cholesky factorisation of A
+ (symmetric, positive definite) */
+ *CHfactor(MAT *A),
+ /* LU factorisation of A (with partial pivoting) */
+ *LUfactor(MAT *A,PERM *pivot),
+ /* QR factorisation of A; need dim(diag) >= # rows of A */
+ *QRfactor(MAT *A,VEC *diag),
+ /* QR factorisation of A with column pivoting */
+ *QRCPfactor(MAT *A,VEC *diag,PERM *pivot),
+ /* L.D.L^T factorisation of A */
+ *LDLfactor(MAT *A),
+ /* Hessenberg factorisation of A -- for schur() */
+ *Hfactor(MAT *A,VEC *diag1,VEC *diag2),
+ /* modified Cholesky factorisation of A;
+ actually factors A+D, D diagonal with no
+ diagonal entry in the factor < sqrt(tol) */
+ *MCHfactor(MAT *A,double tol),
+ *m_inverse(MAT *A,MAT *out);
+
+ /* returns condition estimate for A after LUfactor() */
+extern double LUcondest(MAT *A,PERM *pivot),
+ /* returns condition estimate for Q after QRfactor() */
+ QRcondest(MAT *A);
+
+/* Note: The make..() and ..update() routines assume that the factorisation
+ has already been carried out */
+
+ /* Qout is the "Q" (orthongonal) matrix from QR factorisation */
+extern MAT *makeQ(MAT *A,VEC *diag,MAT *Qout),
+ /* Rout is the "R" (upper triangular) matrix
+ from QR factorisation */
+ *makeR(MAT *A,MAT *Rout),
+ /* Qout is orthogonal matrix in Hessenberg factorisation */
+ *makeHQ(MAT *A,VEC *diag1,VEC *diag2,MAT *Qout),
+ /* Hout is the Hessenberg matrix in Hessenberg factorisation */
+ *makeH(MAT *A,MAT *Hout);
+
+ /* updates L.D.L^T factorisation for A <- A + alpha.u.u^T */
+extern MAT *LDLupdate(MAT *A,VEC *u,double alpha),
+ /* updates QR factorisation for QR <- Q.(R+u.v^T)
+ Note: we need explicit Q & R matrices,
+ from makeQ() and makeR() */
+ *QRupdate(MAT *Q,MAT *R,VEC *u,VEC *v);
+
+/* Solve routines assume that the corresponding factorisation routine
+ has already been applied to the matrix along with auxiliary
+ objects (such as pivot permutations)
+
+ These solve the system A.x = b,
+ except for LUTsolve and QRTsolve which solve the transposed system
+ A^T.x. = b.
+ If x is NULL on entry, then it is created.
+*/
+
+extern VEC *BKPsolve(MAT *A,PERM *pivot,PERM *blocks,VEC *b,VEC *x),
+ *CHsolve(MAT *A,VEC *b,VEC *x),
+ *LDLsolve(MAT *A,VEC *b,VEC *x),
+ *LUsolve(MAT *A,PERM *pivot,VEC *b,VEC *x),
+ *_Qsolve(MAT *A,VEC *,VEC *,VEC *, VEC *),
+ *QRsolve(MAT *A,VEC *,VEC *b,VEC *x),
+ *QRTsolve(MAT *A,VEC *,VEC *b,VEC *x),
+
+
+ /* Triangular equations solve routines;
+ U for upper triangular, L for lower traingular, D for diagonal
+ if diag_val == 0.0 use that values in the matrix */
+
+ *Usolve(MAT *A,VEC *b,VEC *x,double diag_val),
+ *Lsolve(MAT *A,VEC *b,VEC *x,double diag_val),
+ *Dsolve(MAT *A,VEC *b,VEC *x),
+ *LTsolve(MAT *A,VEC *b,VEC *x,double diag_val),
+ *UTsolve(MAT *A,VEC *b,VEC *x,double diag_val),
+ *LUTsolve(MAT *A,PERM *,VEC *,VEC *),
+ *QRCPsolve(MAT *QR,VEC *diag,PERM *pivot,VEC *b,VEC *x);
+
+extern BAND *bdLUfactor(BAND *A,PERM *pivot),
+ *bdLDLfactor(BAND *A);
+extern VEC *bdLUsolve(BAND *A,PERM *pivot,VEC *b,VEC *x),
+ *bdLDLsolve(BAND *A,VEC *b,VEC *x);
+
+
+
+extern VEC *hhvec(VEC *,u_int,Real *,VEC *,Real *);
+extern VEC *hhtrvec(VEC *,double,u_int,VEC *,VEC *);
+extern MAT *hhtrrows(MAT *,u_int,u_int,VEC *,double);
+extern MAT *hhtrcols(MAT *,u_int,u_int,VEC *,double);
+
+extern void givens(double,double,Real *,Real *);
+extern VEC *rot_vec(VEC *,u_int,u_int,double,double,VEC *); /* in situ */
+extern MAT *rot_rows(MAT *,u_int,u_int,double,double,MAT *); /* in situ */
+extern MAT *rot_cols(MAT *,u_int,u_int,double,double,MAT *); /* in situ */
+
+
+/* eigenvalue routines */
+
+ /* compute eigenvalues of tridiagonal matrix
+ with diagonal entries a[i], super & sub diagonal entries
+ b[i]; eigenvectors stored in Q (if not NULL) */
+extern VEC *trieig(VEC *a,VEC *b,MAT *Q),
+ /* sets out to be vector of eigenvectors; eigenvectors
+ stored in Q (if not NULL). A is unchanged */
+ *symmeig(MAT *A,MAT *Q,VEC *out);
+
+ /* computes real Schur form = Q^T.A.Q */
+extern MAT *schur(MAT *A,MAT *Q);
+ /* computes real and imaginary parts of the eigenvalues
+ of A after schur() */
+extern void schur_evals(MAT *A,VEC *re_part,VEC *im_part);
+ /* computes real and imaginary parts of the eigenvectors
+ of A after schur() */
+extern MAT *schur_vecs(MAT *T,MAT *Q,MAT *X_re,MAT *X_im);
+
+
+/* singular value decomposition */
+
+ /* computes singular values of bi-diagonal matrix with
+ diagonal entries a[i] and superdiagonal entries b[i];
+ singular vectors stored in U and V (if not NULL) */
+VEC *bisvd(VEC *a,VEC *b,MAT *U,MAT *V),
+ /* sets out to be vector of singular values;
+ singular vectors stored in U and V */
+ *svd(MAT *A,MAT *U,MAT *V,VEC *out);
+
+/* matrix powers and exponent */
+MAT *_m_pow(MAT *,int,MAT *,MAT *);
+MAT *m_pow(MAT *,int, MAT *);
+MAT *m_exp(MAT *,double,MAT *);
+MAT *_m_exp(MAT *,double,MAT *,int *,int *);
+MAT *m_poly(MAT *,VEC *,MAT *);
+
+/* FFT */
+void fft(VEC *,VEC *);
+void ifft(VEC *,VEC *);
+
+#endif
+
+
+#endif
Property changes on: grass-addons/imagery/i.spec.sam/matrix2.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/meminfo.h
===================================================================
--- grass-addons/imagery/i.spec.sam/meminfo.h (rev 0)
+++ grass-addons/imagery/i.spec.sam/meminfo.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,155 @@
+
+/**************************************************************************
+**
+** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+**
+** Meschach Library
+**
+** This Meschach Library is provided "as is" without any express
+** or implied warranty of any kind with respect to this software.
+** In particular the authors shall not be liable for any direct,
+** indirect, special, incidental or consequential damages arising
+** in any way from use of the software.
+**
+** Everyone is granted permission to copy, modify and redistribute this
+** Meschach Library, provided:
+** 1. All copies contain this copyright notice.
+** 2. All modified copies shall carry a notice stating who
+** made the last modification and the date of such modification.
+** 3. No charge is made for this software or works derived from it.
+** This clause shall not be construed as constraining other software
+** distributed on the same medium as this software, nor is a
+** distribution fee considered a charge.
+**
+***************************************************************************/
+
+
+/* meminfo.h 26/08/93 */
+/* changed 11/12/93 */
+
+
+#ifndef MEM_INFOH
+#define MEM_INFOH
+
+
+
+/* for hash table in mem_stat.c */
+/* Note: the hash size should be a prime, or at very least odd */
+#define MEM_HASHSIZE 509
+#define MEM_HASHSIZE_FILE "meminfo.h"
+
+
+/* default: memory information is off */
+/* set it to 1 if you want it all the time */
+#define MEM_SWITCH_ON_DEF 0
+
+
+/* available standard types */
+#define TYPE_NULL (-1)
+#define TYPE_MAT 0
+#define TYPE_BAND 1
+#define TYPE_PERM 2
+#define TYPE_VEC 3
+#define TYPE_IVEC 4
+
+#ifdef SPARSE
+#define TYPE_ITER 5
+#define TYPE_SPROW 6
+#define TYPE_SPMAT 7
+#endif
+
+#ifdef COMPLEX
+#ifdef SPARSE
+#define TYPE_ZVEC 8
+#define TYPE_ZMAT 9
+#else
+#define TYPE_ZVEC 5
+#define TYPE_ZMAT 6
+#endif
+#endif
+
+/* structure for memory information */
+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_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);
+
+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_mark(int mark);
+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[]);
+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,...);
+
+#else
+int mem_info_is_on();
+int mem_info_on();
+
+long mem_info_bytes();
+int mem_info_numvar();
+void mem_info_file();
+
+void mem_bytes_list();
+void mem_numvar_list();
+
+int mem_stat_reg_list();
+int mem_stat_mark();
+int mem_stat_free_list();
+int mem_stat_show_mark();
+void mem_stat_dump();
+int mem_attach_list();
+int mem_free_vars();
+int mem_is_list_attached();
+void mem_dump_list();
+int mem_stat_reg_vars();
+
+#endif
+
+/* macros */
+
+#define mem_info() mem_info_file(stdout,0)
+
+#define mem_stat_reg(var,type) mem_stat_reg_list((void **)var,type,0)
+#define MEM_STAT_REG(var,type) mem_stat_reg_list((void **)&(var),type,0)
+#define mem_stat_free(mark) mem_stat_free_list(mark,0)
+
+#define mem_bytes(type,old_size,new_size) \
+ mem_bytes_list(type,old_size,new_size,0)
+
+#define mem_numvar(type,num) mem_numvar_list(type,num,0)
+
+
+/* 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 */
+} MEM_CONNECT;
+
+/* max number of lists of types */
+#define MEM_CONNECT_MAX_LISTS 5
+
+
+#endif
Property changes on: grass-addons/imagery/i.spec.sam/meminfo.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/open.c
===================================================================
--- grass-addons/imagery/i.spec.sam/open.c (rev 0)
+++ grass-addons/imagery/i.spec.sam/open.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,110 @@
+/* Spectral angle mapping
+ * (c) Oct/1998 Markus Neteler, Hannover
+ *
+ **************************************************************************
+ ** Matrix computations based on Meschach Library
+ ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+ **************************************************************************
+ *
+ * Cited references are from
+ Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
+ Proceedings of the centre for Mathematics and its Applicaions.
+ The Australian National University. Vol. 32.
+ ISBN 0 7315 1900 0
+*/
+
+#include "global.h"
+#include <stdio.h>
+#include "gis.h"
+#include <math.h>
+#include "matrix.h"
+
+int open_files()
+{
+ char *name, *mapset;
+ FILE *fp;
+ int i;
+
+/* read in matrix file with spectral library */
+
+ fp=fopen(matrixfile,"r");
+ if (fp == NULL)
+ {
+ fprintf(stderr,"ERROR: Matrixfile %s not found.\n",matrixfile);
+ exit(1);
+ }
+ A = m_finput(fp, MNULL);
+ fclose (fp);
+ if ( A->m < A->n )
+ {
+ fprintf(stderr, "Need m (rows) >= n (cols) to obtain least squares fit\n");
+ exit(1);
+ }
+ if (!flag.quiet->answer)
+ {
+ fprintf(stderr, "Your spectral matrix = ");
+ m_output(A);
+ }
+ matrixsize = A->n;
+
+/* open input files from group */
+ if (!I_find_group(group))
+ {
+ fprintf (stderr, "group=%s - not found\n", group);
+ exit(1);
+ }
+ I_get_group_ref(group, &Ref);
+ if (Ref.nfiles <= 1)
+ {
+ fprintf (stderr, "Group %s\n", group);
+ if (Ref.nfiles <= 0)
+ fprintf (stderr, "doesn't have any files\n");
+ else
+ fprintf (stderr, "only has 1 file\n");
+ fprintf (stderr, "The group must have at least 2 files\n");
+ exit(1);
+ }
+ /* Error check: input file number must be equal to matrix size */
+ if (Ref.nfiles != matrixsize)
+ {
+ fprintf (stderr, "Error: Number of %i input files in group <%s>\n", Ref.nfiles, group);
+ fprintf (stderr, " does not match matrix size (%i cols).\n", A->n);
+ exit(1);
+ }
+
+ /* 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();
+ name = Ref.file[i].name;
+ mapset = Ref.file[i].mapset;
+ if (!flag.quiet->answer)
+ fprintf (stderr,"Opening input file no. %i [%s]\n", (i+1), name);
+ if ((cellfd[i] = G_open_cell_old (name, mapset)) < 0)
+ {
+ fprintf (stderr, "Unable to proceed\n");
+ exit(1);
+ }
+ }
+
+ /* 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 < Ref.nfiles; i++)
+ {
+ sprintf(result_name, "%s.%d", result_prefix, (i+1));
+ if (!flag.quiet->answer)
+ fprintf (stderr,"Opening output file [%s]\n", result_name);
+ result_cell[i] = G_allocate_cell_buf();
+ if ((resultfd[i] = G_open_cell_new (result_name)) <0)
+ {
+ fprintf (stderr, "Unable to proceed\n");
+ exit(1) ;
+ }
+ }
+
+
+ return(matrixsize); /* give back number of output files (= Ref.nfiles) */
+}
Property changes on: grass-addons/imagery/i.spec.sam/open.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/spec_angle.c
===================================================================
--- grass-addons/imagery/i.spec.sam/spec_angle.c (rev 0)
+++ grass-addons/imagery/i.spec.sam/spec_angle.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,58 @@
+/* Spectral angle mapping
+ * (c) 1998 Markus Neteler, Hannover
+ *
+ * 26. Oct. 1998 - V. 0.1
+ *
+ ****************************************************************************
+ ** Based on Meschach Library
+ ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+ ****************************************************************************
+ *
+ * Cited references are from
+ * Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
+ * Proceedings of the centre for Mathematics and its Applicaions.
+ * The Australian National University. Vol. 32.
+ * ISBN 0 7315 1900 0
+ *****/
+
+#define GLOBAL
+#define MY_PI 3.141592653589793
+
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include "global.h"
+
+
+float spectral_angle() /* returns spectral angle*/
+{
+
+/* input MAT A, VEC Avector
+ * output cur_angle
+ *
+ * v_DN * v_reference
+ * cos alpha = ----------------------------
+ * ||v_DN|| * ||v_reference||
+ *
+ *
+ * b * Avector
+ * = -----------------------
+ * ||b|| * ||Avector||
+ */
+
+
+ VEC *vtmp1;
+ double norm1, norm2, norm3;
+
+ /* Measure spectral angle*/
+
+ vtmp1 = v_star(Avector, b, VNULL); /* multiply with b vector */
+ norm1 = v_norm1(vtmp1); /* calculate 1-norm */
+ norm2 = v_norm2(Avector); /* calculate 2-norm (Euclidean) */
+ norm3 = v_norm2(b); /* calculate 2-norm (Euclidean) */
+
+ V_FREE(vtmp1);
+
+ curr_angle = (acos(norm1/(norm2 * norm3)) * 180/MY_PI); /* Calculate angle */
+ /* return in degree*/
+}
\ No newline at end of file
Property changes on: grass-addons/imagery/i.spec.sam/spec_angle.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.sam/spectrum.dat
===================================================================
--- grass-addons/imagery/i.spec.sam/spectrum.dat (rev 0)
+++ grass-addons/imagery/i.spec.sam/spectrum.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,12 @@
+# Kanäle: r g b i1 i2 i3
+# Spektren zeilenweise eingeben!
+# 1. Sagebrush
+# 2. Saltbush
+# 3. Boden
+# 4. Trockengras
+#
+Matrix: 4 by 4
+row0: 8.87 13.14 11.71 35.85
+row1: 13.59 20.12 19.61 50.66
+row2: 28.26 34.82 38.27 40.1
+row3: 10.54 16.35 23.7 38.98
Added: grass-addons/imagery/i.spec.sam/spectrum2.dat
===================================================================
--- grass-addons/imagery/i.spec.sam/spectrum2.dat (rev 0)
+++ grass-addons/imagery/i.spec.sam/spectrum2.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,12 @@
+# Kanäle: r g b i1 i2 i3
+# Spektren zeilenweise eingeben!
+# 1. Wald
+# 2. helles Feld
+# 3. Boden
+# 4. Trockengras
+#
+Matrix: 4 by 4
+row0: 38 31 16 67
+row1: 41 28 17 118
+row2: 58 54 51 59
+row3: 44 44 90 107
Added: grass-addons/imagery/i.spec.sam/test.dat
===================================================================
--- grass-addons/imagery/i.spec.sam/test.dat (rev 0)
+++ grass-addons/imagery/i.spec.sam/test.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,9 @@
+# Kanäle: r g b i1 i2 i3
+# Spektren zeilenweise eingeben!
+#
+Matrix: 4 by 4
+row0: 54 83 44 13
+row1: 53 8 4 3
+row2: 48 45 34 11
+row3: 54 68 61 20
+
Added: grass-addons/imagery/i.spec.sam/view.sh
===================================================================
--- grass-addons/imagery/i.spec.sam/view.sh (rev 0)
+++ grass-addons/imagery/i.spec.sam/view.sh 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,4 @@
+d.mon x0; d.rast unmix2.1
+d.mon x1; d.rast unmix2.2
+d.mon x2; d.rast unmix2.3
+d.mon x3; d.rast unmix2.4
Property changes on: grass-addons/imagery/i.spec.sam/view.sh
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/.fig
===================================================================
--- grass-addons/imagery/i.spec.unmix/.fig (rev 0)
+++ grass-addons/imagery/i.spec.unmix/.fig 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,9 @@
+#FIG 3.2
+Landscape
+Center
+Metric
+A4
+100.00
+Single
+-2
+1200 2
Added: grass-addons/imagery/i.spec.unmix/Makefile
===================================================================
--- grass-addons/imagery/i.spec.unmix/Makefile (rev 0)
+++ grass-addons/imagery/i.spec.unmix/Makefile 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.spec.unmix
+
+LIBES = $(IMAGERYLIB) $(GMATHLIB) $(GISLIB)
+DEPENDENCIES= $(IMAGERYDEP) $(GMATHDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/imagery/i.spec.unmix/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/README
===================================================================
--- grass-addons/imagery/i.spec.unmix/README (rev 0)
+++ grass-addons/imagery/i.spec.unmix/README 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,6 @@
+horizontal vector x vertical vector = a number
+vertical vector x horizontal vector = matrix
+matrix x vertical vector = vertical vector
+horizontal vector x matrix = horizontal vector
+matrix A x matrix B = matrix C
+matrix B x matrix A = matrix D
Added: grass-addons/imagery/i.spec.unmix/description.html
===================================================================
--- grass-addons/imagery/i.spec.unmix/description.html (rev 0)
+++ grass-addons/imagery/i.spec.unmix/description.html 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,26 @@
+<H2>DESCRIPTION</H2>
+
+<p><EM>i.spec.unmix</EM> is used to Spectral Unmixing.
+
+<H2>EXAMPLES</H2>
+
+<div class="code"><pre>
+TODO
+</pre></div>
+
+<H2>REFERENCES</H2>
+
+Neteler, M. (1999): Spectral Mixture Analysis von Satellitendaten zur Bestimmung von Bodenbedeckungsgraden im Hinblick auf die Erosionsmodellierung. M.Sc. thesis, University of Hannover.
+
+Neteler, M., D. Grasso, I. Michelazzi, L. Miori, S. Merler, and C. Furlanello. New image processing tools for GRASS. In Proc. Free/Libre and Open Source Software for Geoinformatics: GIS-GRASS Users Conference 2004, Sept. 12-14, Bangkok, Thailand, 2004. http://gisws.media.osaka-cu.ac.jp/grass04/viewabstract.php?id=37 (PDF)
+
+Neteler, M., D. Grasso, I. Michelazzi, L. Miori, S. Merler, and C. Furlanello, 2005. An integrated toolbox for image registration, fusion and classification. International Journal of Geoinformatics. Special Issue on FOSS/GRASS 2004 and GIS-IDEAS 2004, 1(1), pp. 51-61, March 2005.
+
+<H2>SEE ALSO</H2>
+
+<H2>AUTHOR</H2>
+
+Markus Neteler, University of Hannover.
+
+<p>
+<i>Last changed: $Date: 2008/02/10 00:00:00 $</i>
Added: grass-addons/imagery/i.spec.unmix/global.h
===================================================================
--- grass-addons/imagery/i.spec.unmix/global.h (rev 0)
+++ grass-addons/imagery/i.spec.unmix/global.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,43 @@
+#include <math.h>
+#include <grass/imagery.h>
+#include <grass/la.h> /* LAPACK/BLAS */
+
+#ifndef GLOBAL
+#define GLOBAL extern
+#endif
+
+#define MAXFILES 255
+
+GLOBAL mat_struct *A, *A_tilde_trans;
+GLOBAL vec_struct *Avector1, *Avector2, *b, *fraction;
+GLOBAL int matrixsize;
+GLOBAL double svd_error;
+GLOBAL float curr_angle;
+
+GLOBAL char *group;
+GLOBAL struct Ref Ref;
+
+GLOBAL CELL **cell;
+GLOBAL int *cellfd;
+
+GLOBAL CELL **result_cell;
+GLOBAL int *resultfd;
+
+GLOBAL CELL *error_cell;
+GLOBAL int error_fd;
+
+GLOBAL CELL *iter_cell;
+GLOBAL int iter_fd;
+
+GLOBAL char result_name[80];
+GLOBAL char *result_prefix, *matrixfile, *error_name, *iter_name;
+
+GLOBAL struct
+ {
+ struct Flag *quiet;
+ } flag;
+
+GLOBAL struct
+ {
+ struct Flag *veryquiet;
+ } flag2;
Property changes on: grass-addons/imagery/i.spec.unmix/global.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/hist.c
===================================================================
--- grass-addons/imagery/i.spec.unmix/hist.c (rev 0)
+++ grass-addons/imagery/i.spec.unmix/hist.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,14 @@
+#include <grass/gis.h>
+
+void make_history(name, group, matrixfile)
+ char *name, *group, *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);
+ }
+}
Property changes on: grass-addons/imagery/i.spec.unmix/hist.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/histogram.c
===================================================================
--- grass-addons/imagery/i.spec.unmix/histogram.c (rev 0)
+++ grass-addons/imagery/i.spec.unmix/histogram.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,38 @@
+/* this code is from r.support */
+
+#include <grass/gis.h>
+
+int do_histogram (name, mapset)
+ char *name;
+ char *mapset;
+{
+ CELL *cell;
+ struct Cell_head cellhd;
+ struct Cell_stats statf;
+ int nrows, ncols, row;
+ int fd;
+
+ if (G_get_cellhd (name, mapset, &cellhd) < 0)
+ return 0;
+ G_set_window (&cellhd);
+ fd = G_open_cell_old (name, mapset);
+ if (fd < 0)
+ return 0;
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ cell = G_allocate_cell_buf();
+
+ 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_close_cell (fd);
+ free (cell);
+ if (row == nrows)
+ G_write_histogram_cs (name, &statf);
+ G_free_cell_stats (&statf);
+ return;
+}
Property changes on: grass-addons/imagery/i.spec.unmix/histogram.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/main.c
===================================================================
--- grass-addons/imagery/i.spec.unmix/main.c (rev 0)
+++ grass-addons/imagery/i.spec.unmix/main.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,410 @@
+/* Spectral mixture analysis of satellite/aerial images
+ * (c) 1999-2000 Markus Neteler, Hannover, Germany
+ *
+ * COPYRIGHT: (C) 2008 by the GRASS Development Team, Markus Neteler
+ * neteler cealp.it
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ * VERSION: Based on LAPACK/BLAS
+ *
+ * Module calculates
+ * -1
+ * x =A b
+ *
+ * A: matrix of reference spectra (endmembers)
+ * b: pixel vector from satellite image
+ * x: unknown fraction vector
+ *
+ * with two constraints:
+ * 1. SUM x_i = 1 (max. 100%) -> least square problem
+ * 2. 0 <= x_i <= 1 (no negative fractions) -> Steepest Descend of
+ * error surface
+ *
+ *
+ * IMPORTANT: Physical units of matrix file and image data set must fit
+ * otherwise the algorithm might run into endless loop!
+ *
+ * TODO: complete LAPACK/BLAS port. Check with Brad Douglas.
+ ********************************************************************/
+
+#define GLOBAL
+#include <stdio.h>
+#include <stdlib.h>
+#include <strings.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/la.h>
+#include "global.h"
+
+#define GAMMA 10 /* last row value in Matrix and last b vector element */
+ /* for constraint Sum xi = 1 (GAMMA=weight) */
+
+int open_files();
+void spectral_angle();
+
+double find_max(double x, double y)
+{
+ return((x*(x>y))+(y*(x<=y)));
+}
+
+
+CELL myround (x)
+ double x;
+ {
+ CELL n;
+
+ if (x >= 0.0)
+ n = x + .5;
+ else
+ {
+ n = -x + .5;
+ n = -n;
+ }
+ return n;
+ }
+
+
+int main(int argc, char *argv[])
+{
+ int nrows, ncols;
+ int row, col;
+ int band;
+ int i, j, k, iterations;
+ mat_struct *A_tilde;
+ vec_struct *b_gamma;
+ vec_struct *startvector, *A_times_startvector, *errorvector, *temp;
+ mat_struct *A_tilde_trans_mu;
+ int *index;
+ double max1, max2, max_total=0.0;
+ double change, mu, deviation;
+ float anglefield[255][255];
+ double error;
+ char command[80];
+ struct
+ {
+ struct Option *group, *matrixfile, *result, *error, *iter;
+ } parm;
+
+ G_gisinit (argv[0]);
+
+ parm.group = G_define_option();
+ parm.group->key = "group";
+ parm.group->type = TYPE_STRING;
+ parm.group->required = YES;
+ parm.group->description = "Imagery group to be analyzed with Spectral Mixture Analyis";
+
+ parm.matrixfile = G_define_option();
+ parm.matrixfile->key = "matrixfile";
+ parm.matrixfile->type = TYPE_STRING;
+ parm.matrixfile->required = YES;
+ parm.matrixfile->description = "Matrix file containing spectral signatures ";
+
+ parm.result = G_define_option();
+ parm.result->key = "result";
+ parm.result->type = TYPE_STRING;
+ parm.result->required = YES;
+ parm.result->description = "Raster map prefix to hold spectral unmixing results";
+
+ parm.error = G_define_option();
+ parm.error->key = "error";
+ parm.error->type = TYPE_STRING;
+ parm.error->required = YES;
+ parm.error->description = "Raster map to hold unmixing error";
+
+
+ parm.iter = G_define_option();
+ parm.iter->key = "iter";
+ parm.iter->type = TYPE_STRING;
+ parm.iter->required = YES;
+ parm.iter->description = "Raster map to hold number of iterations";
+
+ flag.quiet = G_define_flag();
+ flag.quiet->key = 'q';
+ flag.quiet->description = "Run quietly (but with percentage output)";
+
+ flag2.veryquiet = G_define_flag();
+ flag2.veryquiet->key = 's';
+ flag2.veryquiet->description = "Run silently (say nothing)";
+
+ if (G_parser(argc,argv))
+ exit(EXIT_FAILURE);
+
+ result_prefix = parm.result->answer;
+ error_name = parm.error->answer;
+ iter_name = parm.iter->answer;
+ group = parm.group->answer;
+ matrixfile = parm.matrixfile->answer;
+ if (flag2.veryquiet->answer)
+ flag.quiet->answer = 1;
+
+/* here we go... */
+
+ open_files(); /*A is created here */
+
+ /* ATTENTION: Internally we work here with col-oriented matrixfile,
+ * but the user has to enter the spectra row-wise for his/her's
+ * convenience... That means: Don't mix row- and col-orientation
+ * 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. */
+
+ for (i = 0; i < A->cols; i++) /* go columnwise through matrix*/
+ {
+ Avector1 = G_matvect_get_column(A, i);
+ max1 = G_vector_norm_maxval(Avector1, index); /* get the max. element of this vector */
+ for (j = 0; j < A->cols ; j++)
+ {
+ if (j !=i)
+ {
+ Avector2 = G_matvect_get_column(A, j); /* get next col in A */
+ max2 = G_vector_norm_maxval(Avector2, index); /* get the max. element of this vector */
+ max_total=(find_max(max1,max2),max_total); /* find max of matrix A */
+
+ spectral_angle(); /* check vector angle */
+ anglefield[i][j]= curr_angle; /* save angle in degree */
+ G_vector_free(Avector2);
+ }
+ }
+ G_vector_free(Avector1);
+ G_vector_free(Avector2);
+ }
+
+ /* print out the result */
+ if (!flag.quiet->answer)
+ {
+ fprintf(stderr,"Checking linear dependencies (orthogonality check) of Matrix A:\n");
+ for (i = 0; i < A->cols ; i++)
+ for (j = 0; j < A->cols ; j++)
+ {
+ if (j !=i)
+ fprintf(stderr," Angle between row %i and row %i: %g degree\n", \
+ (i+1), (j+1), anglefield[i][j]);
+ /* internally this is col and not row certainly */
+ }
+ fprintf(stderr,"\n");
+ }
+ /* check it */
+ error=0;
+ for (i = 0; i < A->cols ; i++)
+ for (j = 0; j < A->cols ; j++)
+ if (j !=i)
+ if (anglefield[i][j] < 8.0)
+ {
+ fprintf(stderr,"ERROR: Spectral entries row%i: and row%i: in your matrix \
+ are linear dependent!\n",i,j);
+ fprintf(stderr," You have to revise your reference spectra.\n");
+ error=1;
+ exit(EXIT_FAILURE);
+ }
+
+ if (!error)
+ if (!flag.quiet->answer)
+ fprintf(stderr,"Spectral matrix is o.k. Proceeding...\n\n");
+
+
+/* Begin calculations */
+ /* 1. contraint SUM xi = 1
+ * add last row "1" elements to Matrix A, store in A_tilde
+ * A_tilde is one row-dimension more than A */
+ A_tilde = m_get(A->rows+1, A->cols); /* memory allocation */
+
+ for (i=0; i < A->rows ; i++) /* copy rowwise */
+ for (j=0; j < A->cols; j++) /* copy colwise */
+ G_matrix_get_element(A_tilde, [i][j])= A->G_matrix_get_element[i][j];
+
+ /* fill last row with 1 elements */
+ for (j=0; j < A->cols ; j++)
+ G_matrix_get_element(A_tilde, [A->rows][j])= GAMMA;
+ /* now we have an overdetermined (non-square) system */
+
+
+/* We have a least square problem here: error minimization
+ * T -1 T
+ * unknown fraction = [A_tilde * A_tilde] * A_tilde * b
+ *
+ * A_tilde is the non-square matrix with first constraint in last row.
+ * b is pixel vector from satellite image
+ *
+ * Solve this by deriving above equation and searching the
+ * minimum of this error function in an iterative loop within
+ * both constraints.
+ */
+
+ /* calculate the transpose of A_tilde*/
+ A_tilde_trans = G_matrix_transpose(A_tilde);
+
+/* initialize some values */
+ /* step size must be small enough for covergence of iteration:
+ * mu=0.000001; step size for spectra in range of W/m^2/um
+ * mu=0.000000001; step size for spectra in range of mW/m^2/um
+ * mu=0.000001; step size for spectra in range of reflectance
+ ***/
+ /* check max_total for number of digits to configure mu size*/
+ mu=0.0001 * pow(10,-1*ceil(log10(max_total)));
+ startvector = v_get(A->cols); /* length: no. of spectra */
+ A_times_startvector = v_get(A_tilde->m); /* length: no. of bands */
+ errorvector = v_get(A_tilde->m); /* length: no. of bands */
+ temp = v_get(A_tilde->n); /* length: no. of spectra */
+ A_tilde_trans_mu = m_get(A_tilde->m,A_tilde->n);
+
+/* Now we can calculated the fractions pixelwise */
+ nrows = G_window_rows(); /* get geographical region */
+ ncols = G_window_cols();
+
+ if (!flag.quiet->answer)
+ {
+ fprintf (stderr, "Calculating for %i x %i pixels (%i bands) = %i pixelvectors:\n",\
+ nrows,ncols, Ref.nfiles, (ncols * ncols));
+ fprintf (stderr, "%s ... ", G_program_name());
+ }
+ for (row = 0; row < nrows; row++) /* rows loop in images */
+ {
+ if (!flag2.veryquiet->answer)
+ G_percent(row, nrows, 1);
+ for (band = 0; band < Ref.nfiles; band++) /* get one row for all bands*/
+ {
+ if (G_get_map_row (cellfd[band], cell[band], row) < 0)
+ exit(EXIT_FAILURE);
+ }
+
+ for (col = 0; col < ncols; col++)
+ /* cols loop, work pixelwise for all bands */
+ {
+
+ /* get pixel values of each band and store in b vector: */
+ b_gamma = v_get(A_tilde->m); /* length: no. of bands + 1 (GAMMA)*/
+ for (band = 0; band < Ref.nfiles; band++)
+ b_gamma->ve[band] = cell[band][col];
+ /* add GAMMA for 1. constraint as last element*/
+ b_gamma->ve[Ref.nfiles] = GAMMA;
+
+ /* calculate fraction vector for current pixel
+ * Result is stored in fractions vector
+ * with second constraint: Sum x_i = 1 */
+
+ change=1000; /* initialize */
+ deviation=1000;
+ iterations = 0;
+ for (k = 0; k < (A_tilde->n); k++) /* no. of spectra times */
+ startvector->ve[k] = (1./A_tilde->n);
+
+ /* get start vector and initialize it with equal fractions:
+ * using the neighbor pixelvector as startvector*/
+
+ /* solve with iterative solution: */
+ while( fabs(change) > 0.0001 )
+ {
+ /* go a small step into direction of negative gradient
+ * errorvector = A_tilde * startvector - b_gamma */
+ mv_mlt(A_tilde, startvector, A_times_startvector);
+ v_sub(A_times_startvector, b_gamma, errorvector);
+ sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
+ mv_mlt(A_tilde_trans_mu, errorvector, temp);
+ v_sub(startvector,temp,startvector); /* update startvector */
+
+ /* if one element gets negative, set it to zero */
+ for (k = 0; k < (A_tilde->n); k++) /* no. of spectra times */
+ if (startvector->ve[k] < 0)
+ startvector->ve[k] = 0;
+
+ /* Check the deviation */
+ change = deviation - G_vector_norm_euclid(errorvector);
+ deviation = G_vector_norm_euclid(errorvector);
+
+ /* debug output */
+ /*fprintf(stderr, "Change: %g - deviation: %g\n", \
+ *change, deviation); */
+
+ iterations++;
+ } /* while */
+
+ fraction=v_get(A->cols); /* length: no. of spectra */
+ error = deviation / G_vector_norm_euclid(b_gamma);
+ v_copy(startvector, fraction);
+
+ /*---------- end of second contraint -----------------------*/
+ /* store fractions in resulting rows of resulting files
+ * (number of bands = vector dimension) */
+
+ /* 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*/
+ error_cell[col] = (CELL) (100 * error);
+ iter_cell[col] = iterations;
+
+ G_vector_free(fraction);
+ G_vector_free(b);
+ } /* columns loop */
+
+ /* write the resulting rows into output files: */
+ for (i = 0; i < A->cols; i++) /* no. of spectra */
+ G_put_map_row (resultfd[i], result_cell[i]);
+ if (error_fd > 0)
+ G_put_map_row (error_fd, error_cell);
+ if (iter_fd > 0)
+ G_put_map_row (iter_fd, iter_cell);
+
+ } /* rows loop */
+
+ if (!flag2.veryquiet->answer)
+ G_percent(row, nrows, 2);
+
+ /* close files */
+ for (i = 0; i < Ref.nfiles; i++) /* no. of bands */
+ G_unopen_cell(cellfd[i]);
+
+ for (i = 0; i < A->cols; i++) /* no. of spectra */
+ {
+ G_close_cell (resultfd[i]);
+ /* make grey scale color table */
+ sprintf(result_name, "%s.%d", result_prefix, (i+1));
+ sprintf(command, "r.colors map=%s color=rules 2>&1> /dev/null <<EOF\n
+ 0 0 0 0 \n
+ 100 0 255 0\n
+ end\n
+ EOF", result_name);
+ system(command);
+ /* create histogram */
+ do_histogram(result_name, Ref.file[i].mapset);
+ }
+ if (error_fd > 0)
+ {
+ G_close_cell (error_fd);
+ sprintf(command, "r.colors map=%s color=gyr >/dev/null", error_name);
+ system(command);
+ }
+ if (iter_fd > 0)
+ {
+ G_close_cell (iter_fd);
+ /* sprintf(command, "r.colors map=%s color=gyr >/dev/null", iter_name);
+ system(command);*/
+ }
+
+ G_matrix_free(A);
+
+ make_history(result_name, group, matrixfile);
+ exit(EXIT_SUCCESS);
+} /* main*/
+
Property changes on: grass-addons/imagery/i.spec.unmix/main.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/main.c.meschach
===================================================================
--- grass-addons/imagery/i.spec.unmix/main.c.meschach (rev 0)
+++ grass-addons/imagery/i.spec.unmix/main.c.meschach 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,399 @@
+/* Spectral mixture analysis
+ * (c) 1999 Markus Neteler, Hannover, Germany
+ *
+ * V 1.2 - 20. Jan. 1999
+ * COPYRIGHT: (C) 2008 by the GRASS Development Team, Markus Neteler
+ * neteler cealp.it
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ * VERSION:
+ * * Based on Meschach Library (matrix operations)
+ * * Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+ *
+ * Reference:
+ * Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
+ * Proceedings of the centre for Mathematics and its Applications.
+ * The Australian National University. Vol. 32.
+ * ISBN 0 7315 1900 0
+ *
+ *
+ * Calculation of
+ * Ax = b
+ *
+ * A: matrix of reference spectra (endmembers)
+ * b: pixel vector from satellite image
+ * x: unknown fraction vector
+ *
+ * Two constraints:
+ * 1. SUM x_i = 1 (max. 100%) -> least square problem
+ * 2. 0 <= x_i <= 1 (no negative fractions) -> Steepest Descend of
+ * error surface
+ ********************************************************************/
+
+#define GLOBAL
+#include "global.h"
+#include <stdio.h>
+#include <strings.h>
+#include <math.h>
+#include "matrix.h"
+#include "matrix2.h"
+
+#define GAMMA 10 /* last row value in Matrix and last b vector element
+ * for constraint Sum xi = 1 (weight)*/
+
+int open_files();
+void spectral_angle();
+
+CELL myround (x)
+ double x;
+ {
+ CELL n;
+
+ if (x >= 0.0)
+ n = x + .5;
+ else
+ {
+ n = -x + .5;
+ n = -n;
+ }
+ return n;
+ }
+
+
+int main(argc,argv)
+char *argv[];
+{
+ int nrows, ncols;
+ int row, col;
+ int band;
+ int i, j, k, iterations;
+ MAT *A_tilde;
+ VEC *b_gamma;
+ VEC *startvector, *A_times_startvector, *errorvector, *temp;
+ MAT *A_tilde_trans_mu;
+ double change, mu, deviation;
+ float anglefield[255][255];
+ double error;
+ char command[80];
+ struct
+ {
+ struct Option *group, *matrixfile, *result, *error, *iter;
+ } parm;
+
+ G_gisinit (argv[0]);
+
+ parm.group = G_define_option();
+ parm.group->key = "group";
+ parm.group->type = TYPE_STRING;
+ parm.group->required = YES;
+ parm.group->description = "Imagery group to be analyzed with Spectral Mixture Analyis";
+
+ parm.matrixfile = G_define_option();
+ parm.matrixfile->key = "matrixfile";
+ parm.matrixfile->type = TYPE_STRING;
+ parm.matrixfile->required = YES;
+ parm.matrixfile->description = "Matrix file containing spectral signatures ";
+
+ parm.result = G_define_option();
+ parm.result->key = "result";
+ parm.result->type = TYPE_STRING;
+ parm.result->required = YES;
+ parm.result->description = "Raster map prefix to hold spectral unmixing results";
+
+ parm.error = G_define_option();
+ parm.error->key = "error";
+ parm.error->type = TYPE_STRING;
+ parm.error->required = YES;
+ parm.error->description = "Raster map to hold unmixing error";
+
+
+ parm.iter = G_define_option();
+ parm.iter->key = "iter";
+ parm.iter->type = TYPE_STRING;
+ parm.iter->required = YES;
+ parm.iter->description = "Raster map to hold number of iterations";
+
+ flag.quiet = G_define_flag();
+ flag.quiet->key = 'q';
+ flag.quiet->description = "Run quietly (but with percentage output)";
+
+ flag2.veryquiet = G_define_flag();
+ flag2.veryquiet->key = 's';
+ flag2.veryquiet->description = "Run silently (say nothing)";
+
+ if (G_parser(argc,argv))
+ exit(1);
+
+ result_prefix = parm.result->answer;
+ error_name = parm.error->answer;
+ iter_name = parm.iter->answer;
+ group = parm.group->answer;
+ matrixfile = parm.matrixfile->answer;
+ if (flag2.veryquiet->answer)
+ flag.quiet->answer = 1;
+
+/* here we go... */
+
+ open_files();
+ /* Spectral Matrix is stored in A now
+ * Generally: n: cols .... for matrix A
+ * m: rows
+ * |
+ * |
+ **/
+
+ /* 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.
+ *
+ * We work here internally with col-oriented matrixfile, but the user
+ * has to enter the spectra row-wise for his/her's convenience...
+ * That means: Don't mix rows and cols in the source code and the modules
+ * messages output!
+ */
+
+ for (i = 0; i < A->n; i++) /* go columnwise through matrix*/
+ {
+ Avector1 = get_col(A, i, VNULL);
+ for (j = 0; j < A->n ; j++)
+ {
+ if (j !=i)
+ {
+ Avector2 = get_col(A, j, VNULL); /* compare with next col in A */
+ spectral_angle();
+ anglefield[i][j]= curr_angle; /* in degree */
+ V_FREE(Avector2);
+ }
+ }
+ V_FREE(Avector1);
+ V_FREE(Avector2);
+ }
+
+ /* print out the result */
+ if (!flag.quiet->answer)
+ {
+ fprintf(stderr,"Checking linear dependencies (orthogonality check) of Matrix A:\n");
+ for (i = 0; i < A->n ; i++)
+ for (j = 0; j < A->n ; j++)
+ {
+ if (j !=i)
+ fprintf(stderr," Angle between row %i and row %i: %g degree\n", \
+ (i+1), (j+1), anglefield[i][j]);
+ /* internally this is col and not row certainly */
+ }
+ fprintf(stderr,"\n");
+ }
+ /* check it */
+ error=0;
+ for (i = 0; i < A->n ; i++)
+ for (j = 0; j < A->n ; j++)
+ if (j !=i)
+ if (anglefield[i][j] < 8.0)
+ {
+ fprintf(stderr,"ERROR: Spectral entries row%i: and row%i: in your matrix \
+ are linear dependent!\n",i,j);
+ fprintf(stderr," You have to revise your reference spectra.\n");
+ error=1;
+ exit(1);
+ }
+
+ if (!error)
+ if (!flag.quiet->answer)
+ fprintf(stderr,"Spectral matrix is o.k. Proceeding...\n\n");
+
+
+/* Begin calculations */
+ /* 1. contraint SUM xi = 1
+ * add last row "1" elements to Matrix A, store in A_tilde
+ * A_tilde is one row-dimension more than A */
+ A_tilde = m_get( (A->m+1), A->n);
+
+ for (i=0; i < A->m ; i++) /* copy rowwise */
+ for (j=0; j < A->n; j++) /* copy colwise */
+ A_tilde->me[i][j]= A->me[i][j];
+
+ /* fill last row with 1 elements */
+ for (j=0; j < A->n ; j++)
+ A_tilde->me[A->m][j]= GAMMA;
+ /* 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
+ */
+
+ /* calculate the transpose of A_tilde*/
+ A_tilde_trans = m_transp(A_tilde, MNULL);
+
+ /* calculate the matrix product of A_tilde_trans and A_tilde
+ * T
+ * MM = A_tilde * A_tilde */
+/* MM = m_mlt(A_tilde_trans, A_tilde, MNULL);*/
+
+ /* Calculate the inverse of matrix M with SVD */
+/* calc_inverse(MM);*/
+ /* The inverted matrix of [A_tilde_trans*A_tilde] is
+ * stored in MM_inverse now */
+
+ /* initialize some values */
+ /* mu=0.000001; */ /* step size for spectra in range of W/m^2/um */
+ mu=0.000000001; /* step size for spectra in range of 0.01W/m^2/um */
+ /* step size must be small enough for covergence of iteration*/
+
+ startvector = v_get(A->n); /* length: no. of spectra */
+ A_times_startvector = v_get(A_tilde->m); /* length: no. of bands */
+ errorvector = v_get(A_tilde->m); /* length: no. of bands */
+ temp = v_get(A_tilde->n); /* length: no. of spectra */
+ A_tilde_trans_mu = m_get(A_tilde->m,A_tilde->n);
+
+/* Now we can calculated the fractions pixelwise */
+ nrows = G_window_rows(); /* get geographical region */
+ ncols = G_window_cols();
+
+ if (!flag.quiet->answer)
+ {
+ fprintf (stderr, "Calculating for %i x %i pixels (%i bands) = %i pixelvectors:\n",\
+ nrows,ncols, Ref.nfiles, (ncols * ncols));
+ fprintf (stderr, "%s ... ", G_program_name());
+ }
+ for (row = 0; row < nrows; row++) /* rows loop in images */
+ {
+ if (!flag2.veryquiet->answer)
+ G_percent(row, nrows, 1);
+ for (band = 0; band < Ref.nfiles; band++) /* get one row for all bands*/
+ {
+ if (G_get_map_row (cellfd[band], cell[band], row) < 0)
+ exit(1);
+ }
+
+ for (col = 0; col < ncols; col++)
+ /* cols loop, work pixelwise for all bands */
+ {
+
+ /* get pixel values of each band and store in b vector: */
+ b_gamma = v_get(A_tilde->m); /* length: no. of bands + 1 (GAMMA)*/
+ for (band = 0; band < Ref.nfiles; band++)
+ b_gamma->ve[band] = cell[band][col];
+ /* add GAMMA for 1. constraint as last element*/
+ b_gamma->ve[Ref.nfiles] = GAMMA;
+
+ /* calculate fraction vector for current pixel
+ * Result is stored in fractions vector
+ * with second constraint: Sum x_i = 1 */
+
+ change=1000; /* initialize */
+ deviation=1000;
+ iterations = 0;
+ for (k = 0; k < (A_tilde->n); k++) /* no. of spectra times */
+ startvector->ve[k] = (1./A_tilde->n);
+
+ /* get start vector and initialize it with equal fractions:
+ * using the neighbor pixelvector as startvector*/
+
+ /* solve with iterative solution: */
+ while( fabs(change) > 0.0001 )
+ {
+ /* go a small step into direction of negative gradient
+ * errorvector = A_tilde * startvector - b_gamma */
+ mv_mlt(A_tilde, startvector, A_times_startvector);
+ v_sub(A_times_startvector, b_gamma, errorvector);
+ sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
+ mv_mlt(A_tilde_trans_mu, errorvector, temp);
+ v_sub(startvector,temp,startvector); /* update startvector */
+
+ /* if one element gets negative, set it to zero */
+ for (k = 0; k < (A_tilde->n); k++) /* no. of spectra times */
+ if (startvector->ve[k] < 0)
+ startvector->ve[k] = 0;
+
+ /* Check the deviation */
+ change = deviation - v_norm2(errorvector);
+ deviation = v_norm2(errorvector);
+
+ /* debug output */
+ /* fprintf(stderr, "Change: %g - deviation: %g\n", \
+ change, deviation); */
+
+ iterations++;
+ } /* while */
+
+ fraction=v_get(A->n); /* length: no. of spectra */
+ error = deviation / v_norm2(b_gamma);
+ v_copy(startvector, fraction);
+
+ /*---------- end of second contraint -----------------------*/
+ /* store fractions in resulting rows of resulting files
+ * (number of bands = vector dimension) */
+
+ /* write result in full percent */
+ for (i = 0; i < A->n; 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;
+
+ V_FREE(fraction);
+ V_FREE(b);
+ } /* columns loop */
+
+ /* write the resulting rows into output files: */
+ for (i = 0; i < A->n; i++) /* no. of spectra */
+ G_put_map_row (resultfd[i], result_cell[i], row);
+ if (error_fd > 0)
+ G_put_map_row (error_fd, error_cell, row);
+ if (iter_fd > 0)
+ G_put_map_row (iter_fd, iter_cell, row);
+
+ } /* rows loop */
+
+ if (!flag2.veryquiet->answer)
+ G_percent(row, nrows, 2);
+
+ /* close files */
+ for (i = 0; i < Ref.nfiles; i++) /* no. of bands */
+ G_unopen_cell(cellfd[i]);
+
+ for (i = 0; i < A->n; i++) /* no. of spectra */
+ {
+ G_close_cell (resultfd[i]);
+ /* make grey scale color table */
+ sprintf(result_name, "%s.%d", result_prefix, (i+1));
+ sprintf(command, "r.colors map=%s color=rules 2>&1> /dev/null <<EOF\n
+ 0 0 0 0 \n
+ 100 0 255 0\n
+ end\n
+ EOF", result_name);
+ system(command);
+ /* create histogram */
+ do_histogram(result_name, Ref.file[i].mapset);
+ }
+ if (error_fd > 0)
+ {
+ G_close_cell (error_fd);
+ sprintf(command, "r.colors map=%s color=gyr >/dev/null", error_name);
+ system(command);
+ }
+ if (iter_fd > 0)
+ {
+ G_close_cell (iter_fd);
+ /* sprintf(command, "r.colors map=%s color=gyr >/dev/null", iter_name);
+ system(command);*/
+ }
+
+ M_FREE(A);
+
+ make_history(result_name, group, matrixfile);
+ exit(0);
+} /* main*/
+
Added: grass-addons/imagery/i.spec.unmix/meminfo.h
===================================================================
--- grass-addons/imagery/i.spec.unmix/meminfo.h (rev 0)
+++ grass-addons/imagery/i.spec.unmix/meminfo.h 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,155 @@
+
+/**************************************************************************
+**
+** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+**
+** Meschach Library
+**
+** This Meschach Library is provided "as is" without any express
+** or implied warranty of any kind with respect to this software.
+** In particular the authors shall not be liable for any direct,
+** indirect, special, incidental or consequential damages arising
+** in any way from use of the software.
+**
+** Everyone is granted permission to copy, modify and redistribute this
+** Meschach Library, provided:
+** 1. All copies contain this copyright notice.
+** 2. All modified copies shall carry a notice stating who
+** made the last modification and the date of such modification.
+** 3. No charge is made for this software or works derived from it.
+** This clause shall not be construed as constraining other software
+** distributed on the same medium as this software, nor is a
+** distribution fee considered a charge.
+**
+***************************************************************************/
+
+
+/* meminfo.h 26/08/93 */
+/* changed 11/12/93 */
+
+
+#ifndef MEM_INFOH
+#define MEM_INFOH
+
+
+
+/* for hash table in mem_stat.c */
+/* Note: the hash size should be a prime, or at very least odd */
+#define MEM_HASHSIZE 509
+#define MEM_HASHSIZE_FILE "meminfo.h"
+
+
+/* default: memory information is off */
+/* set it to 1 if you want it all the time */
+#define MEM_SWITCH_ON_DEF 0
+
+
+/* available standard types */
+#define TYPE_NULL (-1)
+#define TYPE_MAT 0
+#define TYPE_BAND 1
+#define TYPE_PERM 2
+#define TYPE_VEC 3
+#define TYPE_IVEC 4
+
+#ifdef SPARSE
+#define TYPE_ITER 5
+#define TYPE_SPROW 6
+#define TYPE_SPMAT 7
+#endif
+
+#ifdef COMPLEX
+#ifdef SPARSE
+#define TYPE_ZVEC 8
+#define TYPE_ZMAT 9
+#else
+#define TYPE_ZVEC 5
+#define TYPE_ZMAT 6
+#endif
+#endif
+
+/* structure for memory information */
+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_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);
+
+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_mark(int mark);
+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[]);
+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,...);
+
+#else
+int mem_info_is_on();
+int mem_info_on();
+
+long mem_info_bytes();
+int mem_info_numvar();
+void mem_info_file();
+
+void mem_bytes_list();
+void mem_numvar_list();
+
+int mem_stat_reg_list();
+int mem_stat_mark();
+int mem_stat_free_list();
+int mem_stat_show_mark();
+void mem_stat_dump();
+int mem_attach_list();
+int mem_free_vars();
+int mem_is_list_attached();
+void mem_dump_list();
+int mem_stat_reg_vars();
+
+#endif
+
+/* macros */
+
+#define mem_info() mem_info_file(stdout,0)
+
+#define mem_stat_reg(var,type) mem_stat_reg_list((void **)var,type,0)
+#define MEM_STAT_REG(var,type) mem_stat_reg_list((void **)&(var),type,0)
+#define mem_stat_free(mark) mem_stat_free_list(mark,0)
+
+#define mem_bytes(type,old_size,new_size) \
+ mem_bytes_list(type,old_size,new_size,0)
+
+#define mem_numvar(type,num) mem_numvar_list(type,num,0)
+
+
+/* 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 */
+} MEM_CONNECT;
+
+/* max number of lists of types */
+#define MEM_CONNECT_MAX_LISTS 5
+
+
+#endif
Property changes on: grass-addons/imagery/i.spec.unmix/meminfo.h
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/meschach_convert.txt
===================================================================
--- grass-addons/imagery/i.spec.unmix/meschach_convert.txt (rev 0)
+++ grass-addons/imagery/i.spec.unmix/meschach_convert.txt 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,5 @@
+For instance, to copy a vector from
+x to y it is sufficient to write y = v_copy(x,VNULL). The VNULL is the
+NULL vector, and usually tells the routine called to create a vector for
+output.
+
Property changes on: grass-addons/imagery/i.spec.unmix/meschach_convert.txt
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/meschach_grass.txt
===================================================================
--- grass-addons/imagery/i.spec.unmix/meschach_grass.txt (rev 0)
+++ grass-addons/imagery/i.spec.unmix/meschach_grass.txt 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,57 @@
+$Id: meschach_grass.txt 11 2000-12-05 21:00:14Z neteler $
+
+missing:
+- G_matrix_norm_maxval()
+- G_vector_free()
+- find max of matrix A (find_max)
+- Subtract two vectors (v_sub)
+- Scalar-matrix multiplication (sm_mlt)
+- stdin read-function "m_input"
+- stdin read in vector from stdin (v_input)
+
+VNULL: null vector
+
+* G\_vector\_free() can be done by G\_matrix\_free().
+* Ax calculations can be done with G\_matrix\_multiply()
+
+/* get matrix element shortcut*/
+#define ME(x, y, z) G_matrix_get_element( (x), (y), (z) )
+
+-------------------------------------------------------
+MESCHACH vs. GRASS implementation
+
+Vector:
+ V_FREE(Avector1): G_vector_free(Avector1)
+ v_norm2: G_vector_norm_euclid
+ v_copy: G_vector_copy
+ v_sub:
+ v_get(A->cols): G_vector_init(?,A->cols, CVEC)
+ ve: G_vector_get_element()
+ v_max(Avector1, index): G_vector_norm_maxval(Avector1, index) ?
+
+Matrix-Vector:
+ get_col(A, i, VNULL): G_matvect_get_column(A, i)
+
+Matrix:
+ M_FREE: G_matrix_free
+ me: G_matrix_get_element A_tilde->me[A->m][j]
+ m_get((A->m+1), A->n) G_matrix_init(A->rows+1, A->cols)
+ A->m: A->rows (matrix dimensions)
+ A->n: A->cols
+
+ m_transp(A_tilde, MNULL): G_matrix_transpose(A_tilde)
+ find_max:
+ m_copy: G_matrix_copy
+ mv_mlt: G_matrix_product
+ ?: G_matrix_print
+ ?: G_matrix_inverse
+ ?: G_matrix_LU_solve
+
+sm_mlt:
+
+m_input: Do you mean something that would
+initialise a matrix from an array like x[i][j]?
+
+v_input:Probably get params:
+[rows,cols] and G_matrix_init(), then set values with
+G_set_matrix_element. See my remarks above re. G_vector_init()
Property changes on: grass-addons/imagery/i.spec.unmix/meschach_grass.txt
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/open.c
===================================================================
--- grass-addons/imagery/i.spec.unmix/open.c (rev 0)
+++ grass-addons/imagery/i.spec.unmix/open.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,155 @@
+/* Spextral unmixing with Singular Value Decomposition */
+/* (c) 15. Jan. 1999 Markus Neteler, Hannover*/
+
+/**************************************************************************
+ ** Matrix computations based on Meschach Library
+ ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+ **************************************************************************/
+
+/* Cited references are from
+ Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
+ Proceedings of the centre for Mathematics and its Applicaions.
+ The Australian National University. Vol. 32.
+ ISBN 0 7315 1900 0
+*/
+
+#include "global.h"
+#include <stdio.h>
+#include <math.h>
+#include <grass/gis.h>
+#include "matrix.h"
+#include "matrix2.h"
+
+int open_files()
+{
+ char *name, *mapset;
+ FILE *fp;
+ int i;
+ MAT *A_input;
+
+/* Read in matrix file with spectral library.
+ Input matrix must contain spectra row-wise (for user's convenience)!
+ Transposed here to col-wise orientation (for modules/mathematical
+ convenience).
+ */
+
+ fp=fopen(matrixfile,"r");
+ if (fp == NULL)
+ {
+ fprintf(stderr,"ERROR: Matrixfile %s not found.\n",matrixfile);
+ exit(EXIT_FAILURE);
+ }
+ A_input = m_finput(fp, MNULL);
+ fclose (fp);
+
+ if (!flag.quiet->answer)
+ {
+ fprintf(stderr, "Your spectral matrix = ");
+ m_output(A_input);
+ }
+
+/* transpose input matrix from row orientation to col orientation.
+ * Don't mix rows and cols in the source code and the modules
+ * messages output! */
+
+ A=m_get(A_input->m, A_input->n);
+ m_transp(A_input, A);
+ M_FREE(A_input);
+
+ if ( A->m < A->n )
+ {
+ fprintf(stderr, "ERROR: Need number of cols >= rows to perform least squares fitting.\n");
+ exit(EXIT_FAILURE);
+ }
+ matrixsize = A->m; /* number of rows must be equivalent to no. of bands */
+
+/* open input files from group */
+ if (!I_find_group(group))
+ {
+ fprintf (stderr, "group=%s - not found\n", group);
+ exit(EXIT_FAILURE);
+ }
+ I_get_group_ref(group, &Ref);
+ if (Ref.nfiles <= 1)
+ {
+ fprintf (stderr, "ERROR: Group %s\n", group);
+ if (Ref.nfiles <= 0)
+ fprintf (stderr, "doesn't have any files\n");
+ else
+ fprintf (stderr, "only has 1 file\n");
+ fprintf (stderr, "The group must have at least 2 files\n");
+ exit(EXIT_FAILURE);
+ }
+ /* Error check: input file number must be equal to matrix size */
+ if (Ref.nfiles != matrixsize)
+ {
+ fprintf (stderr, "ERROR: Number of %i input files in group <%s>\n", Ref.nfiles, group);
+ fprintf (stderr, " does not match no. of spectra in matrix \
+ (contains only %i cols).\n", A->m);
+ exit(1);
+ }
+
+ /* 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();
+ name = Ref.file[i].name;
+ mapset = Ref.file[i].mapset;
+ if (!flag.quiet->answer)
+ fprintf (stderr,"Opening input file no. %i [%s]\n", (i+1), name);
+ if ((cellfd[i] = G_open_cell_old (name, mapset)) < 0)
+ {
+ fprintf (stderr, "Unable to proceed\n");
+ exit(1);
+ }
+ }
+
+/* open files for results*/
+ result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+ for (i=0; i < (A->n); i++) /* no. of spectra */
+ {
+ sprintf(result_name, "%s.%d", result_prefix, (i+1));
+ if (!flag.quiet->answer)
+ fprintf (stderr,"Opening output file [%s]\n", result_name);
+ result_cell[i] = G_allocate_cell_buf();
+ if ((resultfd[i] = G_open_cell_new (result_name)) <0)
+ {
+ fprintf (stderr, "GRASS-Database internal error: Unable to proceed\n");
+ exit(1) ;
+ }
+ }
+
+
+/* open file containing SMA error*/
+
+ error_cell = (CELL *) G_malloc (sizeof (CELL *));
+ if (error_name)
+ {
+ error_fd = G_open_cell_new (error_name);
+ if (!flag.quiet->answer)
+ fprintf (stderr,"Opening error file [%s]\n", error_name);
+ if (error_fd < 0)
+ fprintf (stderr, "Unable to create error layer [%s]", error_name);
+ else
+ error_cell = G_allocate_cell_buf();
+ }
+
+/* open file containing number of iterations */
+
+ iter_cell = (CELL *) G_malloc (sizeof (CELL *));
+ if (iter_name)
+ {
+ iter_fd = G_open_cell_new (iter_name);
+ if (!flag.quiet->answer)
+ fprintf (stderr,"Opening iteration file [%s]\n", iter_name);
+ if (iter_fd < 0)
+ fprintf (stderr, "Unable to create iterations layer [%s]", iter_name);
+ else
+ iter_cell = G_allocate_cell_buf();
+ }
+
+ return(matrixsize); /* give back number of output files (= Ref.nfiles) */
+}
Property changes on: grass-addons/imagery/i.spec.unmix/open.c
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/sample/bdg_farben.col
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/bdg_farben.col (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/bdg_farben.col 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,12 @@
+0 0 0 0
+10% 0 51 0
+20% 0 100 0
+30% 0 170 0
+40% 0 170 0
+50% 140 200 0
+60% 160 220 0
+70% 50 230 50
+80% 30 255 80
+90% 150 255 150
+100% 0 255 0
+end
Added: grass-addons/imagery/i.spec.unmix/sample/color.sh
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/color.sh (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/color.sh 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,5 @@
+r.colors map=$1 col=rules 2>&1 > /dev/null << EOF
+0 0 0 0
+100 0 255 0
+end
+EOF
Property changes on: grass-addons/imagery/i.spec.unmix/sample/color.sh
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/sample/diff_farben.col
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/diff_farben.col (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/diff_farben.col 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,12 @@
+0 0 0 0
+1% 0 0 0
+20% 0 100 0
+30% 0 170 0
+40% 100 200 0
+50% 250 250 0
+60% 240 220 0
+70% 255 220 0
+80% 255 180 0
+90% 255 130 0
+100% 255 0 0
+end
\ No newline at end of file
Added: grass-addons/imagery/i.spec.unmix/sample/diplom96.dat
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/diplom96.dat (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/diplom96.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,12 @@
+# Kanäle: r g b i1 i2 i3
+# Spektren zeilenweise eingeben!
+#
+# 0: Vegetation
+# 1: unbedeckter Boden
+# 2: ? im (Norden) Wald bei Bodenburg
+#
+# Zeilen, Spalten
+Matrix: 3 by 6
+row0: 0 4 0 83 20 0
+row1: 13 24 26 51 53 48
+row2: 5 5 11 41 32 21
Added: grass-addons/imagery/i.spec.unmix/sample/diplom97.dat
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/diplom97.dat (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/diplom97.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,22 @@
+# Kanäle: r g b i1 i2 i3
+# Spektren zeilenweise eingeben!
+#
+#3570637.5 5766432.58928571
+#3571837.5 5764210.76785714
+#3570962.5 5763711.48214286
+#3573662.5 5763287.08928571
+#3573112.5 5766058.125
+#
+# 0: Wasser ?
+# 1: Vegetation volles Wachstum
+# 2: Vegetation auflaufend
+# 3: unbedeckter Boden
+# 4: Boden/Vegetation?
+#
+# Zeilen, Spalten
+Matrix: 5 by 6
+row0: 0 0 0 14 7 0
+row1: 1 4 2 81 21 7
+row2: 1 5 4 52 35 8
+row3: 14 20 27 54 56 49
+row4: 2 3 5 33 40 22
Added: grass-addons/imagery/i.spec.unmix/sample/noshadow.sh
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/noshadow.sh (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/noshadow.sh 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,8 @@
+ r.mapcalc unmix1_noshade = "(unmix.1/(unmix.1-unmix.3)) + unmix.1"
+ r.mapcalc unmix2_noshade = "(unmix.2/(unmix.2-unmix.3)) + unmix.2"
+cat unmix_reclass25.dat |r.reclass i=unmix1_noshade o=unmix1_noshade.recl
+cat unmix_reclass25.dat |r.reclass i=unmix2_noshade o=unmix2_noshade.recl
+cat unmix_reclass25.dat |r.reclass i=bodenbedeckungsgrad95 o=bdgrad95.recl
+r.kappa -mwz classification=unmix1_noshade.recl reference=bdgrad95.recl output=kappa.1
+r.kappa -mwz classification=unmix2_noshade.recl reference=bdgrad95.recl output=kappa.2
+textedit kappa.1&textedit kappa.2&
Property changes on: grass-addons/imagery/i.spec.unmix/sample/noshadow.sh
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/sample/run.sh
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/run.sh (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/run.sh 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,17 @@
+i.spec.unmix gr=96corr err=error97 iter=iter97 res=unmix97 mat=diplom97.dat
+cat unmix_reclass25.dat |r.reclass i=unmix97.1 o=unmix971.recl
+cat unmix_reclass25.dat |r.reclass i=unmix97.2 o=unmix972.recl
+cat unmix_reclass25.dat |r.reclass i=bdg_wrede o=bdg_wrede.recl
+
+cat unmix_reclass25.dat |r.reclass i=unmix97.3 o=unmix973.recl
+cat unmix_reclass25.dat |r.reclass i=unmix97.4 o=unmix974.recl
+cat unmix_reclass25.dat |r.reclass i=unmix97.5 o=unmix975.recl
+#cat unmix_reclass25.dat |r.reclass i=unmix.6 o=unmix6.recl
+
+r.kappa -mwz classification=unmix971.recl reference=bdg_wrede.recl output=kappa.1
+r.kappa -mwz classification=unmix972.recl reference=bdg_wrede.recl output=kappa.2
+#r.kappa -mwz classification=unmix3.recl reference=bodenbedeckungsgrad95 output=kappa.3
+#r.kappa -mwz classification=unmix4.recl reference=bodenbedeckungsgrad95 output=kappa.4
+#r.kappa -mwz classification=unmix5.recl reference=bodenbedeckungsgrad95 output=kappa.5
+#r.kappa -mwz classification=unmix6.recl reference=bodenbedeckungsgrad95 output=kappa.6
+textedit kappa.1&textedit kappa.2&
Property changes on: grass-addons/imagery/i.spec.unmix/sample/run.sh
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:eol-style
+ native
Added: grass-addons/imagery/i.spec.unmix/sample/spectrum.dat
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/spectrum.dat (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/spectrum.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,12 @@
+# Kanäle: r g b i1 i2 i3
+# Spektren zeilenweise eingeben!
+# 1. Sagebrush
+# 2. Saltbush
+# 3. Boden
+# 4. Trockengras
+#
+Matrix: 4 by 6
+row0: 8.87 13.14 11.71 35.85 28.26 10.54
+row1: 13.59 20.12 19.61 70.66 34.82 16.35
+row2: 28.26 34.82 38.27 40.1 38.27 23.7
+row3: 10.54 16.35 23.7 38.98 40.1 38.98
Added: grass-addons/imagery/i.spec.unmix/sample/unmix_reclass.dat
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/unmix_reclass.dat (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/unmix_reclass.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,10 @@
+1 thru 10 = 10
+11 thru 20 = 20
+21 thru 30 = 30
+31 thru 40 = 40
+41 thru 50 = 50
+51 thru 60 = 60
+61 thru 70 = 70
+71 thru 80 = 80
+81 thru 90 = 90
+91 thru 100 = 100
Added: grass-addons/imagery/i.spec.unmix/sample/unmix_reclass20.dat
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/unmix_reclass20.dat (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/unmix_reclass20.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,5 @@
+0 thru 20 = 10
+21 thru 40 = 30
+41 thru 60 = 50
+61 thru 80 = 70
+81 thru 100 = 90
Added: grass-addons/imagery/i.spec.unmix/sample/unmix_reclass25.dat
===================================================================
--- grass-addons/imagery/i.spec.unmix/sample/unmix_reclass25.dat (rev 0)
+++ grass-addons/imagery/i.spec.unmix/sample/unmix_reclass25.dat 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,4 @@
+0 thru 25 = 13
+26 thru 50 = 38
+51 thru 75 = 63
+76 thru 100 = 88
Added: grass-addons/imagery/i.spec.unmix/spec_angle.c
===================================================================
--- grass-addons/imagery/i.spec.unmix/spec_angle.c (rev 0)
+++ grass-addons/imagery/i.spec.unmix/spec_angle.c 2008-05-26 14:37:42 UTC (rev 31528)
@@ -0,0 +1,64 @@
+/* Spectral angle mapping
+ * (c) 1999 Markus Neteler, Hannover
+ *
+ * 25. Nov. 1998 - V. 0.2
+ *
+ ****************************************************************************
+ ** Based on Meschach Library
+ ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
+ ****************************************************************************
+ *
+ * Cited references are from
+ * Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
+ * Proceedings of the centre for Mathematics and its Applicaions.
+ * The Australian National University. Vol. 32.
+ * ISBN 0 7315 1900 0
+ *****/
+
+#define GLOBAL
+#define MY_PI 3.141592653589793
+
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include "global.h"
+
+
+void spectral_angle() /* returns spectral angle globally*/
+{
+
+/* input MAT A, VEC Avector1, Avector2
+ * output cur_angle
+ *
+ * v_DN * v_reference
+ * cos alpha = ----------------------------
+ * ||v_DN|| * ||v_reference||
+ *
+ *
+ * Avector1 * Avector2
+ * = ---------------------------
+ * ||Avector1|| * ||Avector2||
+ *
+ *
+ * Ref.: van der Meer, F. 1997: Mineral mapping and Landsat Thematic Mapper
+ * image Classification using spectral unmixing.
+ * Geocarto International, Vol.12, no.3 (Sept.). pp. 27-40
+ */
+
+
+ VEC *vtmp1;
+ double norm1, norm2, norm3;
+
+
+/* Measure spectral angle*/
+
+ vtmp1 = v_star(Avector1, Avector2, VNULL); /* multiply one A column with second */
+ norm1 = v_norm1(vtmp1); /* calculate 1-norm */
+ norm2 = v_norm2(Avector1); /* calculate 2-norm (Euclidean) */
+ norm3 = v_norm2(Avector2); /* calculate 2-norm (Euclidean) */
+
+ V_FREE(vtmp1);
+
+ curr_angle = (acos(norm1/(norm2 * norm3)) * 180/MY_PI); /* Calculate angle */
+ /* return in degree globally*/
+}
Property changes on: grass-addons/imagery/i.spec.unmix/spec_angle.c
___________________________________________________________________
Name: svn:eol-style
+ native
More information about the grass-commit
mailing list