[GRASS-SVN] r54183 - in grass-addons: grass6/imagery/i.spec.unmix grass7/imagery grass7/imagery/i.spec.unmix grass7/imagery/i.spec.unmix/test grass7/imagery/i.spec.unmix/unused
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Dec 4 10:05:01 PST 2012
Author: rashadkm
Date: 2012-12-04 10:05:00 -0800 (Tue, 04 Dec 2012)
New Revision: 54183
Added:
grass-addons/grass7/imagery/i.spec.unmix/
grass-addons/grass7/imagery/i.spec.unmix/.fig
grass-addons/grass7/imagery/i.spec.unmix/Makefile
grass-addons/grass7/imagery/i.spec.unmix/README
grass-addons/grass7/imagery/i.spec.unmix/bdg_farben.col
grass-addons/grass7/imagery/i.spec.unmix/color.sh
grass-addons/grass7/imagery/i.spec.unmix/diff_farben.col
grass-addons/grass7/imagery/i.spec.unmix/diplom96.dat
grass-addons/grass7/imagery/i.spec.unmix/diplom97.dat
grass-addons/grass7/imagery/i.spec.unmix/docu.lg
grass-addons/grass7/imagery/i.spec.unmix/global.h
grass-addons/grass7/imagery/i.spec.unmix/hist.c
grass-addons/grass7/imagery/i.spec.unmix/histogram.c
grass-addons/grass7/imagery/i.spec.unmix/i.spec.unmix.html
grass-addons/grass7/imagery/i.spec.unmix/kappa.1
grass-addons/grass7/imagery/i.spec.unmix/kappa.2
grass-addons/grass7/imagery/i.spec.unmix/kappa.3
grass-addons/grass7/imagery/i.spec.unmix/kappa.4
grass-addons/grass7/imagery/i.spec.unmix/la_extra.c
grass-addons/grass7/imagery/i.spec.unmix/la_extra.h
grass-addons/grass7/imagery/i.spec.unmix/main.c
grass-addons/grass7/imagery/i.spec.unmix/makedocu.sh
grass-addons/grass7/imagery/i.spec.unmix/meschach_grass.txt
grass-addons/grass7/imagery/i.spec.unmix/noshadow.sh
grass-addons/grass7/imagery/i.spec.unmix/open.c
grass-addons/grass7/imagery/i.spec.unmix/region.unmix.txt
grass-addons/grass7/imagery/i.spec.unmix/run.sh
grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c
grass-addons/grass7/imagery/i.spec.unmix/spectrum.dat
grass-addons/grass7/imagery/i.spec.unmix/test/
grass-addons/grass7/imagery/i.spec.unmix/test/Gmakefile
grass-addons/grass7/imagery/i.spec.unmix/test/gmath_test.c
grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass.dat
grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass20.dat
grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass25.dat
grass-addons/grass7/imagery/i.spec.unmix/unused/
grass-addons/grass7/imagery/i.spec.unmix/unused/main.c.old.s
grass-addons/grass7/imagery/i.spec.unmix/unused/matrix.h
grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h
grass-addons/grass7/imagery/i.spec.unmix/unused/meminfo.h
grass-addons/grass7/imagery/i.spec.unmix/unused/meschach_convert.txt
grass-addons/grass7/imagery/i.spec.unmix/unused/svd_calc.c
Modified:
grass-addons/grass6/imagery/i.spec.unmix/main.c
Log:
ported i.spec.unmix to grass7
Modified: grass-addons/grass6/imagery/i.spec.unmix/main.c
===================================================================
--- grass-addons/grass6/imagery/i.spec.unmix/main.c 2012-12-04 17:49:25 UTC (rev 54182)
+++ grass-addons/grass6/imagery/i.spec.unmix/main.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -4,7 +4,7 @@
* MODULE: i.spec.unmix
*
* AUTHOR(S): Markus Neteler <neteler itc.it>
- * Mohammed Rashad <rashadkm gmail.com>
+ * Mohammed Rashad <rashadkm gmail.com>(updates)
*
* PURPOSE: Spectral mixture analysis of satellite/aerial images
*
Added: grass-addons/grass7/imagery/i.spec.unmix/.fig
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/.fig (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/.fig 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,9 @@
+#FIG 3.2
+Landscape
+Center
+Metric
+A4
+100.00
+Single
+-2
+1200 2
Added: grass-addons/grass7/imagery/i.spec.unmix/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/Makefile 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.spec.unmix
+
+
+LIBES = $(GISLIB) $(GMATHLIB) $(IMAGERYLIB) $(RASTERLIB)
+DEPENDENCIES= $(GISDEP) $(GMATHDEP) $(IMAGERYDEP) $(RASTERDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/imagery/i.spec.unmix/README
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/README (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/README 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/bdg_farben.col
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/bdg_farben.col (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/bdg_farben.col 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/color.sh
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/color.sh (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/color.sh 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/color.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/imagery/i.spec.unmix/diff_farben.col
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/diff_farben.col (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/diff_farben.col 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/diplom96.dat
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/diplom96.dat (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/diplom96.dat 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/diplom97.dat
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/diplom97.dat (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/diplom97.dat 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/docu.lg
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/docu.lg (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/docu.lg 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,17 @@
+\subsection*{global.h}
+%< global.h
+
+\subsection*{open.c}
+%< open.c
+
+\subsection*{hist.c}
+%< hist.c
+
+\subsection*{histogram.c}
+%< histogram.c
+
+\subsection*{spec\underline{~}angle.c}
+%< spec_angle.c
+
+\subsection*{main.c}
+%< main.c
Added: grass-addons/grass7/imagery/i.spec.unmix/global.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/global.h (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/global.h 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,34 @@
+#ifndef __GLOBAL_H__
+#define __GLOBAL_H__
+
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/la.h>
+
+#ifndef GLOBAL
+#define GLOBAL extern
+#endif
+
+
+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 float spectral_angle(vec_struct *, vec_struct *);
+GLOBAL int do_histogram(char *, char *);
+GLOBAL void make_history(char *, char *, char *);
+GLOBAL int open_files(char *, char *, char *, char *, mat_struct *A);
+
+#endif /* __GLOBAL_H__ */
Added: grass-addons/grass7/imagery/i.spec.unmix/hist.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/hist.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/hist.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,16 @@
+#include <stdio.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+
+
+void make_history (char *name, char *group, char *matrixfile)
+{
+ struct History hist;
+
+ if (Rast_read_history (name, G_mapset(), &hist) >= 0)
+ {
+ sprintf (hist.fields[1], "Group: %s", group);
+ sprintf (hist.fields[2], "Matrix file: %s", matrixfile);
+ Rast_write_history (name, &hist);
+ }
+}
Added: grass-addons/grass7/imagery/i.spec.unmix/histogram.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/histogram.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/histogram.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,51 @@
+#include <grass/gis.h>
+
+#include <grass/raster.h>
+
+int do_histogram (char *name, char *mapset)
+{
+ CELL *cell;
+ struct Cell_head cellhd;
+ struct Cell_stats statf;
+ int nrows, ncols, row;
+ int fd;
+ struct Cell_head region;
+
+Rast_get_cellhd (name, mapset, &cellhd);
+// if (Rast_get_cellhd (name, mapset, &cellhd) < 0)
+ // return 0;
+
+ G_set_window (&cellhd);
+ fd = Rast_open_old (name, mapset);
+ if (fd < 0)
+ return 0;
+
+
+
+ G_get_window(®ion);
+
+ nrows = region.rows;
+ ncols = region.cols;
+
+// nrows = G_window_rows ();
+ // ncols = G_window_cols ();
+ cell = Rast_allocate_c_buf ();
+
+ Rast_init_cell_stats (&statf);
+ for (row = 0; row < nrows; row++)
+ {
+ Rast_get_c_row_nomask (fd, cell, row);
+ // break;
+
+ Rast_update_cell_stats (cell, ncols, &statf);
+ }
+ Rast_close (fd);
+ G_free (cell);
+
+ if (row == nrows)
+ Rast_write_histogram_cs (name, &statf);
+
+ Rast_free_cell_stats (&statf);
+
+ return 0;
+}
Added: grass-addons/grass7/imagery/i.spec.unmix/i.spec.unmix.html
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/i.spec.unmix.html (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/i.spec.unmix.html 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,80 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS manual: i.spec.unmix</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>i.spec.unmix</b></em> - Perfroms Spectral mixture analysis of satellite/aerial images
+<h2>KEYWORDS</h2>
+imagery, spectral unmixing
+<h2>SYNOPSIS</h2>
+<b>i.spec.unmix</b><br>
+<b>i.spec.unmix help</b><br>
+<b>i.spec.unmix</b> <b>group</b>=<em>name</em> <b>matrix</b>=<em>string</em> [<b>result</b>=<em>string</em>] <b>error</b>=<em>name</em> <b>iter</b>=<em>name</em> [--<b>overwrite</b>] [--<b>verbose</b>] [--<b>quiet</b>]
+
+<h3>Flags:</h3>
+<DL>
+<DT><b>--overwrite</b></DT>
+<DD>Allow output files to overwrite existing files</DD>
+<DT><b>--verbose</b></DT>
+<DD>Verbose module output</DD>
+<DT><b>--quiet</b></DT>
+<DD>Quiet module output</DD>
+</DL>
+
+<h3>Parameters:</h3>
+<DL>
+<DT><b>group</b>=<em>name</em></DT>
+<DD>Name of input imagery group</DD>
+
+<DT><b>matrix</b>=<em>string</em></DT>
+<DD>Open Matrix file</DD>
+<DD>Matrix file containing spectral signatures</DD>
+
+<DT><b>result</b>=<em>string</em></DT>
+<DD>Name of raster map prefix to hold spectral unmixing results</DD>
+
+<DT><b>error</b>=<em>name</em></DT>
+<DD>Name of raster map to hold unmixing error</DD>
+
+<DT><b>iter</b>=<em>name</em></DT>
+<DD>Raster map to hold number of iterations</DD>
+
+</DL>
+<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: 2011-11-09 03:26:45 +0530 (Wed, 09 Nov 2011) $</i>
+<HR>
+<P><a href="index.html">Main index</a> - <a href="imagery.html">imagery index</a> - <a href="full_index.html">Full index</a></P>
+<P>© 2003-2011 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>
Added: grass-addons/grass7/imagery/i.spec.unmix/kappa.1
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/kappa.1 (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/kappa.1 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,40 @@
+ ACCURACY ASSESSMENT
+LOCATION: geosum Fri Jan 29 17:07:20 1999
+MASK: wrede_maske in geosum, categories 1
+MAPS: MAP1 = Reclass of bdg_wrede in geosum (bdg_wrede.recl in geosum)
+ MAP2 = Reclass of unmix97.1 in geosum (unmix971.recl in geosum)
+
+Error Matrix
+Panel #1 of 1
+ MAP1
+ cat# 13 38 63 88 Row Sum
+ M 13 59856 104 0 0 59960
+ A 38 4758 0 0 0 4758
+ P 63 76177 81545 0 0 157722
+ 2 88 53421 6876 0 0 60297
+Col Sum 194212 88525 0 0 282737
+
+
+Cats % Commission % Ommission Estimated Kappa
+13 99.826551 30.819929 0.994460
+38 0.000000 0.000000 -0.455816
+63 NA NA NA
+88 NA NA NA
+
+Kappa Kappa Variance
+0.071564 0.000000
+
+Obs Correct Total Obs % Observed Correct
+59856 282737 21.170204
+
+MAP1 Category Description
+13: (no description)
+38: (no description)
+63: (no description)
+88: (no description)
+
+MAP2 Category Description
+13: (no description)
+38: (no description)
+63: (no description)
+88: (no description)
Added: grass-addons/grass7/imagery/i.spec.unmix/kappa.2
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/kappa.2 (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/kappa.2 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,40 @@
+ ACCURACY ASSESSMENT
+LOCATION: geosum Fri Jan 29 17:06:25 1999
+MASK: wrede_maske in geosum, categories 1
+MAPS: MAP1 = Reclass of bdg_wrede in geosum (bdg_wrede.recl in geosum)
+ MAP2 = Reclass of unmix97.2 in geosum (unmix972.recl in geosum)
+
+Error Matrix
+Panel #1 of 1
+ MAP1
+ cat# 13 38 63 88 Row Sum
+ M 13 53797 5965 198 0 59960
+ A 38 3318 1145 295 0 4758
+ P 63 5922 39483 112317 0 157722
+ 2 88 13404 24993 17702 4198 60297
+Col Sum 76441 71586 130512 4198 282737
+
+
+Cats % Commission % Ommission Estimated Kappa
+13 89.721481 70.377154 0.859129
+38 24.064733 1.599475 -0.016794
+63 71.212006 86.058753 0.465303
+88 6.962204 100.000000 0.055600
+
+Kappa Kappa Variance
+0.419272 0.000017
+
+Obs Correct Total Obs % Observed Correct
+171457 282737 60.641869
+
+MAP1 Category Description
+13: (no description)
+38: (no description)
+63: (no description)
+88: (no description)
+
+MAP2 Category Description
+13: (no description)
+38: (no description)
+63: (no description)
+88: (no description)
Added: grass-addons/grass7/imagery/i.spec.unmix/kappa.3
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/kappa.3 (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/kappa.3 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,89 @@
+ ACCURACY ASSESSMENT
+LOCATION: geosum Thu Jan 28 11:43:43 1999
+MASK: feldmaske in geosum, categories 1
+MAPS: MAP1 = Reclass of rfn95 in geosum (bodenbedeckungsgrad95 in geosum)
+ MAP2 = Reclass of unmix.3 in geosum (unmix3.recl in geosum)
+
+Error Matrix
+Panel #1 of 2
+ MAP1
+ cat# -1 0 2 5 10 20 30 40 50
+ M -1 0 0 0 0 0 0 570 5755 4
+ A 0 0 0 0 0 0 1 197 3297 1
+ P 2 0 0 0 0 0 0 0 226 0
+ 2 5 0 0 0 0 2 6 58 3892 0
+ 10 0 0 0 0 0 0 30 184 0
+ 20 0 0 0 0 0 0 0 0 0
+ 30 0 0 0 0 0 0 0 0 0
+ 40 0 0 0 0 0 0 0 0 0
+ 50 0 0 0 0 0 0 33 314 0
+ 70 0 0 0 0 0 4 964 8291 5
+ 90 0 0 0 0 0 0 44 118 0
+ 100 0 0 0 0 2 35 588 2632 0
+Col Sum 0 0 0 0 4 46 2484 24709 10
+
+Panel #2 of 2
+ MAP1
+ cat# 70 90 100 Row Sum
+ M -1 0 0 0 6329
+ A 0 0 0 0 3496
+ P 2 0 0 0 226
+ 2 5 0 0 0 3958
+ 10 0 0 0 214
+ 20 0 0 0 0
+ 30 0 0 0 0
+ 40 0 0 0 0
+ 50 0 0 0 347
+ 70 0 0 0 9264
+ 90 0 0 0 162
+ 100 0 0 0 3257
+Col Sum 0 0 0 27253
+
+
+Cats % Commission % Ommission Estimated Kappa
+-1 0.000000 NaN -999.000000
+0 NA NA NA
+2 NA NA NA
+5 NA NA NA
+10 0.000000 0.000000 -0.000147
+20 NA NA NA
+30 NA NA NA
+40 NA NA NA
+50 0.000000 0.000000 -0.000367
+70 NA NA NA
+90 NA NA NA
+100 NA NA NA
+
+Kappa Kappa Variance
+-0.000006 0.000000
+
+Obs Correct Total Obs % Observed Correct
+0 27253 0.000000
+
+MAP1 Category Description
+-1: no data
+0: % Bodenbedeckungsgrad
+2: % Bodenbedeckungsgrad
+5: % Bodenbedeckungsgrad
+10: % Bodenbedeckungsgrad
+20: (no description)
+30: (no description)
+40: (no description)
+50: % Bodenbedeckungsgrad
+70: % Bodenbedeckungsgrad
+90: % Bodenbedeckungsgrad
+100: % Bodenbedeckungsgrad
+
+MAP2 Category Description
+-1: (no description)
+0: (no description)
+2: (no description)
+5: (no description)
+10: (no description)
+20: (no description)
+30: (no description)
+40: (no description)
+50: (no description)
+70: (no description)
+90: (no description)
+100: (no description)
Added: grass-addons/grass7/imagery/i.spec.unmix/kappa.4
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/kappa.4 (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/kappa.4 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,40 @@
+ ACCURACY ASSESSMENT
+LOCATION: geosum Fri Jan 29 17:06:18 1999
+MASK: wrede_maske in geosum, categories 1
+MAPS: MAP1 = Reclass of bdg_wrede in geosum (bdg_wrede.recl in geosum)
+ MAP2 = Reclass of unmix97.4 in geosum (unmix974.recl in geosum)
+
+Error Matrix
+Panel #1 of 1
+ MAP1
+ cat# 13 38 63 88 Row Sum
+ M 13 7292 32193 20124 351 59960
+ A 38 1224 1951 1534 49 4758
+ P 63 153501 4052 169 0 157722
+ 2 88 51002 8270 1025 0 60297
+Col Sum 213019 46466 22852 400 282737
+
+
+Cats % Commission % Ommission Estimated Kappa
+13 12.161441 3.423169 -2.562238
+38 41.004624 4.198769 0.294024
+63 0.107151 0.739541 -0.086765
+88 0.000000 0.000000 -0.001417
+
+Kappa Kappa Variance
+-0.220489 -0.000001
+
+Obs Correct Total Obs % Observed Correct
+9412 282737 3.328889
+
+MAP1 Category Description
+13: (no description)
+38: (no description)
+63: (no description)
+88: (no description)
+
+MAP2 Category Description
+13: (no description)
+38: (no description)
+63: (no description)
+88: (no description)
Added: grass-addons/grass7/imagery/i.spec.unmix/la_extra.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/la_extra.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/la_extra.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,575 @@
+
+
+#include <stdio.h> /* needed here for ifdef/else */
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+
+/********
+ ******** only compile this LAPACK/BLAS wrapper file if g2c.h is present!
+ ********/
+#include <grass/config.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/la.h>
+#include "la_extra.h"
+
+
+vec_struct *
+G_matvect_get_column2(mat_struct *mt, int col)
+{
+ int i; /* loop */
+ vec_struct *vc1;
+
+ if(col < 0 || col >= mt->cols) {
+ G_warning(_("Specified matrix column index is outside range"));
+ return NULL;
+ }
+
+ if(!mt->is_init) {
+ G_warning(_("Matrix is not initialised"));
+ return NULL;
+ }
+
+ if( (vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) == NULL ) {
+ G_warning(_("Could not allocate space for vector structure"));
+ return NULL;
+ }
+
+ for ( i = 0; i < mt->rows; i++ )
+ {
+ double dd = G_matrix_get_element(mt, i, col);
+ //G_warning("element at row=%d, col=%d is %lf",i,col,dd);
+ G_matrix_set_element( (mat_struct *)vc1, i, 0, dd);
+
+ }
+
+ return vc1;
+}
+
+
+void
+G_matrix_print2 (mat_struct *mt, const char* name)
+{
+ int i, j;
+
+
+ if(mt!=NULL)
+ {
+ G_message ("start matrix(%s)", name);
+ G_message ("Size: %d x %d", mt->rows, mt->cols);
+
+ for( i = 0; i < mt->rows; i++ )
+ {
+ char buf[2048], numbuf[640];
+ sprintf( buf, "row%d: ", i);
+ for( j = 0; j < mt->cols; j++ )
+ {
+
+ double element = G_matrix_get_element(mt, i, j);
+ sprintf( numbuf, "%14.6f",element);
+ strcat(buf, numbuf);
+ //if( j < mt->cols - 1 )
+ //strcat(buf, ", ");
+ }
+ G_message ("%s", buf);
+ }
+
+ G_message ("end matrix(%s)",name);
+ }
+
+/*
+ for( i = 0; i < mt->rows; i++ ) {
+ strcpy(buf, "");
+
+ for( j = 0; j < mt->cols; j++ ) {
+
+ sprintf( numbuf, "%14.6f", G_matrix_get_element(mt, i, j) );
+ strcat(buf, numbuf);
+ if( j < mt->cols - 1 )
+ strcat(buf, ", ");
+ }
+
+ G_message ("%s", buf);
+ }
+
+ fprintf (stderr, "\n");
+
+*/
+}
+
+
+mat_struct *G_matrix_resize(mat_struct * in, int rows, int cols)
+{
+
+ mat_struct *matrix;
+
+
+ matrix = G_matrix_init(rows,cols,rows);
+
+
+ int i,j,p, index = 0;
+ for ( i=0; i<rows; i++ )
+ {
+
+
+
+ //A_row = A_v[i]; b_v = b->ve;
+ for ( j=0; j<cols; j++ )
+ {
+
+ matrix->vals[index++] = in->vals[i + j * cols];
+
+ // sum += in->vals[i + j * cols] *b->ve[j];
+ //out->ve[i] = sum;
+ }
+
+ }
+
+ int old_size = in->rows * in->cols;
+
+ int new_size = rows * cols;
+ if(new_size > old_size)
+ for(p=old_size;p<new_size;p++)
+ matrix->vals[p] = 0.0;
+
+
+ return matrix;
+}
+
+
+
+mat_struct *sm_mlt(double scalar,mat_struct *matrix, mat_struct *out)
+
+{
+ int m,n,i,j;
+
+ int index = 0;
+
+ if ( matrix==NULL )
+ G_fatal_error("sm_mlt1(error)");
+
+ if ( out==NULL)
+ out = G_matrix_init(matrix->rows, matrix->cols,matrix->rows);
+
+ if (out->rows != matrix->rows || out->cols != matrix->cols )
+ out = G_matrix_resize(out,matrix->rows,matrix->cols);
+
+m = matrix->rows;
+n = matrix->cols;
+ for ( i=0; i<m; i++ )
+ {
+ //__smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
+ /**************************************************/
+ for ( j=0; j<n; j++ )
+ {
+ out->vals[index++] = scalar*matrix->vals[i+ j * m];
+ }}
+
+ //G_matrix_print(matrix,"matrix");
+ //G_matrix_print(matrix,"out");
+ /**************************************************/
+ return (out);
+}
+
+VEC *G_vec_copy(VEC *in)
+{
+ int i;
+ VEC *out;
+ if ( !in )
+ G_fatal_error("v_copy(error1)");
+
+ int dim = in->dim;
+
+ out = G_vec_get(dim);
+
+
+ for( i = 0; i < dim; i++ )
+ {
+ out->ve[i]= in->ve[i];
+ }
+
+ return (out);
+}
+
+
+double v_norm2(VEC *x)
+{
+ int i, dim;
+ double s, sum;
+
+ if ( !x )
+ G_fatal_error("v_norm2(error1)");
+
+
+ dim = x->dim;
+
+ sum = 0.0;
+
+ for ( i = 0; i < dim; i++ )
+ {
+ sum += x->ve[i] * x->ve[i];
+ }
+
+
+
+ return sqrt(sum);
+}
+
+
+
+VEC *v_sub(VEC *vec1,VEC *vec2,VEC *out)
+{
+ /* u_int i, dim; */
+ /* Real *out_ve, *vec1_ve, *vec2_ve; */
+
+ if ( !vec1 || !vec2 )
+ G_fatal_error("v_sub1(error)");
+
+ if ( vec1->dim != vec2->dim )
+ G_fatal_error("v_sub2(error)");
+
+ if ( out == NULL)
+ out = G_vec_get(vec1->dim);
+ //out = v_resize(out,A->m);
+
+ //G_vec_print(vec1, "vec1");
+
+
+ if ( out->dim != vec1->dim )
+ out = G_vec_resize(out,vec1->dim);
+
+ int i;
+ for ( i = 0; i < vec1->dim; i++ )
+ out->ve[i] = vec1->ve[i] - vec2->ve[i];
+
+ /************************************************************
+ dim = vec1->dim;
+ out_ve = out->ve; vec1_ve = vec1->ve; vec2_ve = vec2->ve;
+ for ( i=0; i<dim; i++ )
+ out->ve[i] = vec1->ve[i]-vec2->ve[i];
+ (*out_ve++) = (*vec1_ve++) - (*vec2_ve++);
+ ************************************************************/
+
+ return (out);
+}
+
+
+
+VEC *mv_mlt(mat_struct *A, VEC *b, VEC *out)
+{
+ unsigned int i, m, n,j;
+ double **A_v, *b_v /*, *A_row */;
+ /* register Real sum; */
+
+
+
+// if ( A==(MAT *)NULL || b==(VEC *)NULL )
+ // error(E_NULL,"mv_mlt");
+
+//G_matrix_print(A, "A(mlt)");
+
+ if ( A->cols != b->dim )
+
+ G_fatal_error("mv_mlt1(error)");
+
+
+ if ( b == out )
+ G_fatal_error("mv_mlt2(error)");
+ //error(E_INSITU,"mv_mlt");
+ if(!out)
+ {
+ G_fatal_error("mv_mltsss3(error)");
+ exit(1);
+ out = G_vec_get2(A->rows, out);
+ }
+ if ( out->dim != A->rows)
+ {
+ G_fatal_error("mv_mlt3(error)");
+ exit(1);
+ out = G_vec_resize(out,A->rows);
+ }
+ //out = G_vec_get(A->rows);
+ //out = v_resize(out,A->m);
+
+
+
+ //if ( out->dim != A->rows )
+
+
+ m = A->rows; n = A->cols;
+ A_v = A->vals; b_v = b->ve;
+
+
+
+ for ( i=0; i<m; i++ )
+ {
+ double sum=0.0;
+int width = A->rows;
+
+ //A_row = A_v[i]; b_v = b->ve;
+ for ( j=0; j<n; j++ )
+ {
+
+ sum += A->vals[i + j * width] *b->ve[j];
+ out->ve[i] = sum;
+
+
+//G_message("sum: %lf of j=%d", b->ve[j],j);
+
+}
+ }
+
+
+ return out;
+}
+
+
+
+
+
+VEC *G_vec_resize(VEC * in, int size)
+{
+
+ VEC *vector;
+
+
+ vector = (VEC *)G_malloc( sizeof(VEC) );
+
+
+ vector->ve = (double *)G_malloc(size * sizeof(double) );
+ int i,j;
+
+ G_message(":%d",in->dim);
+for( i = 0; i < in->dim; i++ )
+ {
+
+ vector->ve[i]= in->ve[i];
+ G_message("ss:%lf",in->ve[i]);
+
+ }
+
+ if(size > in->dim)
+ for(j=i;j<size;j++)
+ vector->ve[j]= 0.0;
+
+ vector->dim = vector->max_dim = size;
+
+ return vector;
+}
+
+
+
+
+
+
+
+VEC *G_vec_get(int size)
+{
+
+ VEC *vector;
+
+
+ vector = (VEC *)G_malloc( sizeof(VEC) );
+
+
+ vector->ve = (double *)G_malloc(size * sizeof(double) );
+ int i;
+for( i = 0; i < size; i++ )
+ {
+
+ vector->ve[i]= 0.0;
+
+
+ }
+
+ vector->dim = vector->max_dim = size;
+
+ return vector;
+}
+
+
+VEC *G_vec_get2(int size, VEC *vector)
+{
+
+ //VEC *vector;
+
+
+ vector = (VEC *)G_malloc( sizeof(VEC) );
+
+
+ vector->ve = (double *)G_malloc(size * sizeof(double) );
+ int i;
+for( i = 0; i < size; i++ )
+ {
+
+ vector->ve[i]= 0.0;
+
+
+ }
+
+ vector->dim = vector->max_dim = size;
+
+ return vector;
+}
+
+
+void G_vec_print (VEC *vector, const char* name)
+{
+ int i;
+
+
+ if(vector!=NULL)
+ {
+ G_message ("start vector(%s)", name);
+ //G_message ("Size: %d x %d", vector->dim);
+
+ for( i = 0; i < vector->dim; i++ )
+ {
+ char buf[2048], numbuf[640];
+ sprintf( buf, "%lf ", vector->ve[i]);
+
+ G_message ("%s", buf);
+ }
+
+ G_message ("end vector(%s)",name);
+ }
+
+}
+
+
+
+
+vec_struct *
+G_vector_product (vec_struct *v1, vec_struct *v2)
+{
+ int idx1, idx2, idx0;
+ int i;
+
+
+ vec_struct *out = G_vector_init (v1->rows, v1->ldim, CVEC);
+
+
+
+ //G_warning("Avector1->rows: %d",Avector1->rows);
+ //G_warning("Avector1->cols: %d",Avector1->cols);
+
+ //G_warning("vtmp1->rows1: %d",vtmp1->rows);
+ //G_warning("vtmp1->cols1: %d",out->cols);
+
+ //vtmp1.type = Avector1->type;
+
+ if (!out->is_init) {
+ G_warning (_("Output vector is uninitialized"));
+ return NULL;
+ }
+
+ if (v1->type != v2->type) {
+ G_warning (_("Vectors are not of the same type"));
+ return NULL;
+ }
+
+ if (v1->type != out->type) {
+ G_warning (_("Output vector is of incorrect type"));
+ return NULL;
+ }
+
+ if (v1->type == MATRIX_) {
+ G_warning (_("Matrices not allowed"));
+ return NULL;
+ }
+
+ if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
+ (v1->type == COLVEC_ && v1->rows != v2->rows))
+ {
+ G_warning (_("Vectors have differing dimensions"));
+ return NULL;
+ }
+
+ if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
+ (v1->type == COLVEC_ && v1->rows != out->rows))
+ {
+ G_warning (_("Output vector has incorrect dimension"));
+ return NULL;
+ }
+
+#if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS) //&& defined(HAVE_G2C_H)
+ f77_dhad (v1->cols, 1.0, v1->vals, 1, v2->vals, 1, 0.0, out->vals, 1.0);
+#else
+ idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
+ idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
+ idx0 = (out->v_indx > 0) ? out->v_indx : 0;
+
+ if (v1->type == ROWVEC_) {
+ for (i = 0; i < v1->cols; i++)
+ G_matrix_set_element(out, idx0, i,
+ G_matrix_get_element(v1, idx1, i) *
+ G_matrix_get_element(v2, idx2, i));
+ } else {
+ for (i = 0; i < v1->rows; i++)
+ G_matrix_set_element(out, i, idx0,
+ G_matrix_get_element(v1, i, idx1) *
+ G_matrix_get_element(v2, i, idx2));
+ }
+#endif
+
+ return out;
+}
+
+
+
+
+int
+G_matrix_read2(FILE *fp, mat_struct *out)
+{
+ char buff[100];
+ int rows, cols;
+ int i, j, row;
+ double val;
+
+ /* skip comments */
+ for (;;) {
+ if (!G_getl(buff, sizeof(buff), fp))
+ return -1;
+ if (buff[0] != '#')
+ break;
+ }
+
+ if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) {
+ G_warning(_("Input format error1"));
+ return -1;
+ }
+
+
+ G_matrix_set(out, rows, cols, rows);
+
+
+ //G_warning("row: %d",row);
+
+ for (i = 0; i < rows; i++)
+ {
+ if (fscanf(fp, "row%d:", &row) != 1)
+ {
+ G_warning(_("Input format error"));
+ return -1;
+ }
+
+ for (j = 0; j < cols; j++)
+ {
+ if (fscanf(fp, "%lf:", &val) != 1)
+ {
+ G_warning(_("Input format error"));
+ return -1;
+ }
+
+ fgetc(fp);
+ G_matrix_set_element(out, i, j, val);
+ }
+ }
+//G_matrix_print(out);
+ return 0;
+}
+
+
Added: grass-addons/grass7/imagery/i.spec.unmix/la_extra.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/la_extra.h (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/la_extra.h 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,31 @@
+#include <grass/config.h>
+#include <stdio.h>
+#include <grass/blas.h>
+#include <grass/lapack.h>
+typedef struct VEC_ {
+ int dim, max_dim;
+ double *ve;
+ } VEC;
+
+int
+G_matrix_read2(FILE *fp, mat_struct *out);
+void
+G_matrix_print2 (mat_struct *mt, const char* name);
+vec_struct *
+G_matvect_get_column2(mat_struct *mt, int col);
+
+
+
+VEC *G_vec_get(int size);
+VEC* G_vec_get2(int size, VEC *vector);
+void G_vec_print (VEC *vector, const char* name);
+VEC *G_vec_resize(VEC * in, int size);
+
+
+VEC *v_sub(VEC *vec1,VEC *vec2,VEC *out);
+VEC *mv_mlt(mat_struct *A, VEC *b, VEC *out);
+mat_struct *sm_mlt(double scalar,mat_struct *matrix, mat_struct *out);
+mat_struct *G_matrix_resize(mat_struct * in, int rows, int cols);
+double v_norm2(VEC *x);
+VEC *G_vec_copy(VEC *in);
+
Added: grass-addons/grass7/imagery/i.spec.unmix/main.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/main.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/main.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,618 @@
+/****************************************************************************
+ *
+ * MODULE: i.spec.unmix
+ *
+ * AUTHOR(S): Markus Neteler <neteler itc.it>
+ * Mohammed Rashad <rashadkm gmail.com>(updates)
+ *
+ * PURPOSE: Spectral mixture analysis of satellite/aerial images
+ *
+ * COPYRIGHT: (C) 2006-2012 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General
+ * Public License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
+ *
+ *****************************************************************************/
+
+#define GLOBAL
+
+#include <grass/config.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <strings.h>
+#include <math.h>
+#include <grass/imagery.h>
+#include <grass/gmath.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+#include <grass/gis.h>
+
+#include "la_extra.h"
+
+
+#define GAMMA 10 /* last row value in Matrix and last b vector element
+ * for constraint Sum xi = 1 (GAMMA=weight) */
+
+
+static double find_max (double x, double y)
+{
+ return ((x * (x > y)) + (y * (x <= y)));
+}
+
+
+int main (int argc, char *argv[])
+{
+ char result_name[80];
+ int nrows, ncols;
+ int row;
+ int i, j, k;
+ VEC *b, *b_gamma;
+ struct Cell_head region;
+ VEC *startvector, *A_times_startvector, *errorvector, *temp;
+ mat_struct *A, *A_tilde, *A_tilde_trans_mu, *A_tilde_trans;
+
+ mat_struct B, B_tilde, B_tilde_trans_mu;
+
+ struct GModule *module;
+
+ double max_total = 0.0;
+ double mu;
+ float anglefield[255][255];
+ double error = 0.0;
+ struct {
+ struct Option *group, *matrixfile, *result, *error, *iter;
+ } parm;
+
+ /* initialize GIS engine */
+ G_gisinit (argv[0]);
+ module = G_define_module();
+
+ module->keywords = _("imagery, spectral unmixing");
+ module->description =_("Perfroms Spectral mixture analysis of satellite/aerial images");
+
+ parm.group = G_define_standard_option (G_OPT_I_GROUP);
+
+
+ parm.matrixfile = G_define_standard_option(G_OPT_F_INPUT);
+ parm.matrixfile->key = "matrix";
+ parm.matrixfile->type = TYPE_STRING;
+ parm.matrixfile->required = YES;
+ parm.matrixfile->label = _("Open Matrix file");
+ parm.matrixfile->description = _("Matrix file containing spectral signatures");
+
+
+
+
+ parm.result = G_define_option ();
+ parm.result->key = "result";
+ parm.result->description = _("Name of raster map prefix to hold spectral unmixing results");
+ parm.result->guisection = _("Required");
+
+ parm.error = G_define_standard_option (G_OPT_R_OUTPUT);
+ parm.error->key = "error";
+ parm.error->description = _("Name of raster map to hold unmixing error");
+
+ parm.iter = G_define_standard_option (G_OPT_R_OUTPUT);
+ parm.iter->key = "iter";
+ parm.iter->description = _("Raster map to hold number of iterations");
+
+ if (G_parser (argc, argv))
+
+ exit (EXIT_FAILURE);
+
+
+ /* here we go... A is created here */
+ A = open_files2 (parm.matrixfile->answer,
+ parm.group->answer,
+ parm.result->answer,
+ parm.iter->answer,
+ parm.error->answer);
+ //G_warning("printing A******");
+ //G_matrix_print2(A,"A");
+
+ /* 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. */
+
+ /* go columnwise through matrix */
+
+
+
+
+ for (i = 0; i < A->cols; i++)
+ {
+ vec_struct *Avector1, *Avector2;
+ double max1, max2;
+
+ Avector1 = G_matvect_get_column2 (A, i);
+
+ // G_matrix_print(Avector1);
+
+//G_warning("Avector1: %d", Avector1->rows);
+
+ // get the max. element of this vector
+ max1 = G_vector_norm_maxval (Avector1, 1);
+
+double temp = max1;
+
+ for (j = 0; j < A->cols; j++)
+ {
+ if (j != i)
+ {
+ // get next col in A
+ Avector2 = G_matvect_get_column (A, j);
+//G_matrix_print(Avector2);
+ // get the max. element of this vector
+ max2 = G_vector_norm_maxval (Avector2, 1);
+
+//G_warning("max2: %lf", max2);
+
+
+
+if(max2 > max1)
+temp = max2;
+
+//G_warning("temp: %lf", temp);
+
+if(temp > max_total)
+max_total = temp;
+ // find max of matrix A
+// max_total = (find_max (max1, max2), max_total);
+//G_warning("max_total: %lf", max_total);
+ // save angle in degree
+ anglefield[i][j] = spectral_angle (Avector1, Avector2);
+ //G_warning("anglefield[i][j]: %lf", anglefield[i][j]);
+
+ }
+ }
+
+ G_vector_free (Avector1);
+ G_vector_free (Avector2);
+ }
+
+ G_message (_("Checking linear dependencies (orthogonality check) of Matrix A..."));
+
+ for (i = 0; i < A->cols; i++)
+ for (j = 0; j < A->cols; j++)
+ if (j != i)
+ // internally this is col and not row certainly
+ G_message (_("Angle between row %i and row %i: %g degree"), (i + 1), (j + 1), anglefield[i][j]);
+
+ for (i = 0; i < A->cols ; i++)
+ for (j = 0; j < A->cols ; j++)
+ if (j != i)
+ if (anglefield[i][j] < 8.0)
+ G_fatal_error (_("Spectral entries row %i: and row %i: in "
+ "your matrix are linear dependent!\nYou "
+ "have to revise your reference spectra."),
+ i, j);
+
+ if (!error)
+ G_message (_("Spectral matrix is o.k. Proceeding..."));
+
+ // Begin calculations
+ // 1. contraint SUM xi = 1
+ // add last row "1" elements to Matrix A, store in A_tilde
+ // A_tilde is one row-dimension more than A
+
+
+
+
+
+ // memory allocation
+ A_tilde = G_matrix_init (A->rows + 1, A->cols, A->rows+1);
+ if (A_tilde == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
+
+
+ //G_message("rrrr%d", A_tilde->rows);
+
+
+ for (i = 0; i < A->rows ; i++)
+ for (j = 0; j < A->cols; j++)
+ G_matrix_set_element (A_tilde, i, j, G_matrix_get_element (A, i, j));
+
+
+
+ // fill last row with 1 elements
+
+ for (j = 0; j < A_tilde->cols; j++)
+ {
+ //G_message("Row: %d, Col:%d",i,j);
+ G_matrix_set_element (A_tilde, i, j, GAMMA);
+ }
+//G_matrix_print2(A_tilde, "A_tilde");
+
+
+ // now we have an overdetermined (non-square) system
+
+ // We have a least square problem here: error minimization
+ // T -1 T
+ // unknown fraction = [A_tilde * A_tilde] * A_tilde * b
+ //
+ // A_tilde is the non-square matrix with first constraint in last row.
+ // b is pixel vector from satellite image
+ //
+ // Solve this by deriving above equation and searching the
+ // minimum of this error function in an iterative loop within
+ // both constraints.
+
+
+ // calculate the transpose of A_tilde
+//G_matrix_print(A_tilde);
+
+ A_tilde_trans = G_matrix_transpose (A_tilde);
+//G_matrix_print2(A_tilde_trans, "A_tilde_trans");
+
+
+ // initialize some values
+
+ // step size must be small enough for covergence of iteration:
+ // mu = 0.000001; step size for spectra in range of W/m^2/um
+ // mu = 0.000000001; step size for spectra in range of mW/m^2/um
+ // mu = 0.000001; step size for spectra in range of reflectance
+ //
+
+ // check max_total for number of digits to configure mu size
+ mu = 0.0001 * pow (10, -1 * ceil (log10 (max_total)));
+
+
+ //G_message("mu = %lf", mu);
+startvector = G_vec_get2(A->cols,startvector);
+
+
+
+ if (startvector == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+
+
+ A_times_startvector = G_vec_get2(A_tilde->rows,A_times_startvector); // length: no. of bands //
+ errorvector = G_vec_get2(A_tilde->rows,errorvector); // length: no. of bands //
+ temp = G_vec_get2(A_tilde->cols,temp); // length: no. of spectra //
+ // A_tilde_trans_mu = m_get(A_tilde->m,A_tilde->n);
+
+
+
+
+ // length: no. of bands
+
+ if (A_times_startvector == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+ // length: no. of bands
+
+ if (errorvector == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+ // length: no. of spectra
+
+ if (temp == NULL)
+ G_fatal_error (_("Unable to allocate memory for vector"));
+
+ A_tilde_trans_mu = G_matrix_init (A_tilde->rows, A_tilde->cols, A_tilde->rows);
+ if (A_tilde_trans_mu == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
+
+ // Now we can calculated the fractions pixelwise
+ //nrows = G_region_rows (); // get geographical region
+ //ncols = G_window_cols ();
+
+ G_get_window(®ion);
+
+ nrows = region.rows;
+ ncols = region.cols;
+
+ G_message (_("Calculating for %i x %i pixels (%i bands) = %i pixelvectors."),
+ nrows, ncols, Ref.nfiles, (ncols * ncols));
+
+
+
+ for (row = 0; row < nrows; row++)
+ {
+ int col, band;
+
+ G_percent (row, nrows, 1);
+
+ // get one row for all bands
+ for (band = 0; band < Ref.nfiles; band++)
+ Rast_get_c_row (cellfd[band], cell[band], row);
+
+
+// if (Rast_get_c_row (cellfd[band], cell[band], row) < 0)
+ // G_fatal_error (_("Unable to get map row [%d]"), row);
+
+ // for (band = 0; band < Ref.nfiles; band++)
+ // {
+ // if (Rast_get_c_row (cellfd[band], cell[band], row) < 0)
+ // G_fatal_error (_("Unable to get map row [%d]"), row);
+
+ //G_message("row: %d, nrows: %d", row, nrows);
+ for (col = 0; col < ncols; col++)
+ {
+
+
+
+ double change = 1000;
+ double deviation = 1000;
+ int iterations = 0;
+
+
+ // get pixel values of each band and store in b vector:
+ // length: no. of bands + 1 (GAMMA)
+
+ b_gamma = G_vec_get2(A_tilde->rows,b_gamma);
+
+// b_gamma = G_vector_init (A_tilde->rows, 1, RVEC);
+ if (b_gamma == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
+
+
+ //G_message("%d", A_tilde->rows);
+ for (band = 0; band < Ref.nfiles; band++)
+ {
+ b_gamma->ve[band] = cell[band][col];
+ //G_message("band: %d col: %d", band,col);
+ //G_matrix_set_element (b_gamma, 0, band, cell[band][col]);
+ }
+
+
+ // add GAMMA for 1. constraint as last element
+ b_gamma->ve[Ref.nfiles] = GAMMA;
+/// G_matrix_set_element (b_gamma, 0, Ref.nfiles, GAMMA);
+
+
+//G_matrix_print(b_gamma,"b_gamma");
+
+ for (k = 0; k < A_tilde->cols; k++)
+ startvector->ve[k] = (1.0 / A_tilde->cols);
+ // G_matrix_set_element (startvector, k,0, (1.0 / A_tilde->cols));
+//G_matrix_print(startvector,"startvector1");
+
+
+//G_matrix_print(b_gamma, "b_gama");
+ // calculate fraction vector for current pixel
+ // Result is stored in fractions vector
+ // with second constraint: Sum x_i = 1
+
+
+
+
+
+ // get start vector and initialize it with equal fractions:
+ // using the neighbor pixelvector as startvector
+
+
+ // solve with iterative solution:
+ while (fabs (change) > 0.0001)
+ {
+
+
+
+ A_times_startvector = mv_mlt(A_tilde, startvector, A_times_startvector);
+
+
+ errorvector = v_sub(A_times_startvector, b_gamma, errorvector);
+
+ A_tilde_trans_mu = sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
+
+
+//G_matrix_print(A_tilde_trans_mu,"A_tilde_trans_mu");
+
+
+ temp = mv_mlt(A_tilde_trans_mu, errorvector, temp);
+ startvector = v_sub(startvector,temp,startvector); // update startvector //
+
+ // if one element gets negative, set it to zero //
+ for (k = 0; k < (A_tilde->cols); k++) // no. of spectra times //
+ if (startvector->ve[k] < 0)
+ startvector->ve[k] = 0;
+
+ // Check the deviation //
+ double norm2 = v_norm2(errorvector);
+ change = deviation - norm2;
+ deviation = norm2;
+
+ iterations++;
+
+ // if(fabs (change) > 0.0001)
+ //G_message("change=%lf, deviation=%lf",change, 0.0001);
+
+ /********************/
+
+
+ // go a small step into direction of negative gradient
+ // errorvector = A_tilde * startvector - b_gamma
+ //
+// mv_mlt(A_tilde, startvector, A_times_startvector);
+
+/*
+
+//G_matrix_print(A_times_startvector,"A_times_startvector");
+//G_matrix_print(b_gamma, "b_gamma");
+//G_matrix_print(errorvector, "errorvector");
+ G_vector_sub (A_times_startvector, b_gamma, errorvector);
+
+
+// sm_mlt(mu, A_tilde_trans, A_tilde_trans_mu);
+// mv_mlt(A_tilde_trans_mu, errorvector, temp);
+ G_vector_sub (startvector, temp, startvector);
+
+ // if one element gets negative, set it to zero
+ for (k = 0; k < A_tilde->cols; k++) // no. of spectra times
+ {
+// if (startvector->ve[k] < 0)
+// startvector->ve[k] = 0;
+
+//G_message("get_element: %lf", G_matrix_get_element (startvector, startvector->cols, k));
+
+ if ((G_matrix_get_element (startvector, startvector->cols, k) < 0))
+ {
+ G_matrix_set_element (startvector, startvector->cols, k, 0);
+ //G_message("A_tilde->cols: %d", A_tilde->cols); //4
+ }
+ }
+ // Check the deviation
+ double norm_euclid = G_vector_norm_euclid (errorvector);
+
+ // G_message("norm_euclid : %lf", norm_euclid );
+ change = deviation - norm_euclid;
+ deviation = G_vector_norm_euclid (errorvector);
+
+ // G_message ("Change: %g - deviation: %g", change, deviation);
+
+ G_debug (5, "Change: %g - deviation: %g",
+ change, deviation);
+*/
+
+
+ }
+
+ // if(fabs (change) > 0.0001)
+
+
+
+
+ VEC *fraction;
+//G_message("fcol %d and A->cols %d", startvector->dim, A->cols);
+ fraction=G_vec_get(A->cols); // length: no. of spectra //
+ error = deviation / v_norm2(b_gamma);
+ fraction = G_vec_copy (startvector);
+
+
+
+ // write result in full percent //
+ for (i = 0; i < A->cols; i++) // no. of spectra //
+ result_cell[i][col] = (CELL)(100 * fraction->ve[i]);
+
+
+
+ // save error and iterations//
+ error_cell[col] = (CELL) (100 * error);
+ iter_cell[col] = iterations;
+
+ //V_FREE(fraction);
+ //V_FREE(b);
+
+ // }
+
+
+
+
+ //---------- end of second contraint -----------------------
+ // store fractions in resulting rows of resulting files
+ // (number of bands = vector dimension)
+
+ // write result in full percent
+ //G_matrix_print(fraction,"fraction"); //same as startvector
+
+ //G_message ("fraction->rows: %d",fraction->rows);
+ //G_message ("i=%d, col=%d",i,col);
+
+/*
+
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ double dd = G_matrix_get_element (fraction, i, 0);
+ dd = 100*dd;
+ //G_message ("i=%d, col=%d",i,col);
+ result_cell[i][col] = (CELL) dd;
+ //result_cell[i][col] = (CELL)(100 * G_matrix_get_element (fraction, fraction->rows-1, i));
+ }
+
+ // save error and iterations
+ error_cell[col] = (CELL) (100 * error);
+ iter_cell[col] = iterations;
+
+ G_vector_free (fraction);
+ // G_vector_free (b);
+
+*/
+ } //end cols loop
+
+ //G_message("finished %d of %d", row,nrows);
+ // }
+ // }
+
+
+
+ // write the resulting rows into output files:
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ Rast_put_c_row (resultfd[i], result_cell[i]);
+
+ if (error_fd > 0)
+ Rast_put_c_row (error_fd, error_cell);
+
+ if (iter_fd > 0)
+ Rast_put_c_row (iter_fd, iter_cell);
+
+ } // rows loop
+ G_percent (row, nrows, 2);
+
+ // close files
+ for (i = 0; i < Ref.nfiles; i++) // no. of bands
+ Rast_unopen (cellfd[i]);
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ char command[1080];
+
+ Rast_close (resultfd[i]);
+
+ // make grey scale color table
+ sprintf (result_name, "%s.%d", parm.result->answer, (i+1));
+ sprintf (command, "r.colors map=%s color=rules <<EOF\n"
+ "0 0 0 0 \n"
+ "201 0 255 0\n"
+ "end\n"
+ "EOF", result_name);
+
+
+ //G_message(command);
+ // G_system (command);
+
+ // create histogram
+ do_histogram (result_name, Ref.file[i].mapset);
+ }
+
+ if (error_fd > 0)
+ {
+ char command[80];
+
+ Rast_close (error_fd);
+ //sprintf (command, "r.colors map=%s color=gyr >/dev/null", parm.error->answer);
+ //G_system (command);
+ }
+
+ if (iter_fd > 0)
+ Rast_close (iter_fd);
+
+ G_matrix_free (A);
+
+ make_history (result_name, parm.group->answer, parm.matrixfile->answer);
+
+/*********************
+*********************/
+ exit (EXIT_SUCCESS);
+}
Added: grass-addons/grass7/imagery/i.spec.unmix/makedocu.sh
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/makedocu.sh (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/makedocu.sh 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,4 @@
+lgrind -e docu.lg > docu_ispecunmix.tex
+#latex docu_ispecunmix.tex
+rm -rf OBJ.linux statistics.txt *.aux *.log *.idx
+mv docu_ispecunmix.tex ../..
Property changes on: grass-addons/grass7/imagery/i.spec.unmix/makedocu.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/imagery/i.spec.unmix/meschach_grass.txt
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/meschach_grass.txt (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/meschach_grass.txt 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,54 @@
+$Id: meschach_grass.txt 11 2000-12-05 21:00:14Z neteler $
+
+missing:
+- 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
+
+* 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: G_vector_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, i, 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()
Added: grass-addons/grass7/imagery/i.spec.unmix/noshadow.sh
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/noshadow.sh (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/noshadow.sh 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/noshadow.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/imagery/i.spec.unmix/open.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/open.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/open.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,316 @@
+/* Spextral unmixing with Singular Value Decomposition */
+/* (c) 15. Jan. 1999 Markus Neteler, Hannover*/
+
+/* Cited references are from
+ Steward, D.E, Leyk, Z. 1994: Meschach: Matrix computations in C.
+ Proceedings of the centre for Mathematics and its Applicaions.
+ The Australian National University. Vol. 32.
+ ISBN 0 7315 1900 0
+*/
+
+#include <stdio.h>
+#include <math.h>
+#include <grass/imagery.h>
+#include <grass/gmath.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+
+int open_files (char *matrixfile,
+ char *img_grp,
+ char *iter_name,
+ char *error_name,
+ mat_struct *A)
+{
+ char result_name[80];
+ char *result_prefix="out";
+ FILE *fp;
+ int i, matrixsize;
+ mat_struct A_input;
+
+
+// mat_struct A_input1;
+/*
+
+ if ((fp = fopen (matrixfile, "r")) == NULL)
+ G_fatal_error (_("Matrix file %s not found."), matrixfile);
+
+
+ //G_matrix_init2( &A_input,3,3,3);
+ // G_warning("matrix read start");
+// G_warning("cols=%d",A_input->cols);
+
+ if ((G_matrix_read (fp, &A_input) < 0))
+ G_fatal_error (_("Unable to read matrix file %s."), matrixfile);
+ fclose (fp);
+
+ //G_matrix_print(&A_input);
+
+
+ // G_warning("matrix read done");
+#if 0
+ G_message(_("Your spectral matrix = %d"), m_output(A_input));
+#endif
+
+
+
+// A = m_get(A_input->rows, A_input->cols);
+
+ G_matrix_clone2(&A_input, A);
+
+
+
+ //*A = *G_matrix_init (A_input.rows, A_input.cols, A_input.rows);
+// if (A == NULL)
+ // G_fatal_error (_("Unable to allocate memory for matrix"));
+
+ A = G_matrix_transpose (&A_input);
+
+ // G_matrix_free (&A_input);
+
+// G_matrix_print(A);
+
+ if ((A->rows) < (A->cols))
+ G_fatal_error (_("Need number of cols >= rows to perform least squares fitting."));
+
+ // number of rows must be equivalent to no. of bands
+ matrixsize = A->cols;
+
+ // open input files from group
+ if (!I_find_group (img_grp))
+ G_fatal_error (_("Unable to find imagery group %s."), img_grp);
+
+ I_get_group_ref (img_grp, &Ref);
+ if (Ref.nfiles <= 1)
+ {
+ if (Ref.nfiles <= 0)
+ G_fatal_error (_("Group %s does not have any rasters. "
+ "The group must have at least 2 rasters."), img_grp);
+ else
+ G_fatal_error (_("Group %s only has 1 raster. "
+ "The group must have at least 2 rasters."), img_grp);
+ }
+
+ // Error check: input file number must be equal to matrix size
+ if (Ref.nfiles != matrixsize)
+ G_fatal_error (_("Number of input files (%i) in group <%s> "
+ "does not match number of spectra in matrix. "
+ "(contains only %i cols)."),
+ Ref.nfiles, img_grp, A->cols);
+
+ // get memory for input files
+ cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ cellfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+ for (i = 0; i < Ref.nfiles; i++)
+ {
+ cell[i] = Rast_allocate_c_buf();
+
+ G_message (_("Opening input file no. %i [%s]"), (i + 1), Ref.file[i].name);
+
+ if ((cellfd[i] = Rast_open_old (Ref.file[i].name, Ref.file[i].mapset)) < 0)
+ G_fatal_error (_("Unable to open <%s>"), Ref.file[i].name);
+ }
+
+ // open files for results
+ result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ sprintf (result_name, "%s.%d", result_prefix, (i + 1));
+ G_message (_("Opening output file [%s]"), result_name);
+
+ result_cell[i] = Rast_allocate_c_buf();
+ if ((resultfd[i] = Rast_open_c_new (result_name)) < 0)
+ G_fatal_error (_("GRASS-DB internal error: Unable to proceed."));
+ }
+
+ // open file containing SMA error
+ error_cell = (CELL *) G_malloc (sizeof(CELL *));
+ if (error_name)
+ {
+ G_message (_("Opening error file [%s]"), error_name);
+
+ if ((error_fd = Rast_open_c_new (error_name)) < 0)
+ G_fatal_error (_("Unable to create error layer [%s]"), error_name);
+ else
+ error_cell = Rast_allocate_c_buf();
+ }
+
+ // open file containing number of iterations
+ iter_cell = (CELL *) G_malloc (sizeof(CELL *));
+ if (iter_name)
+ {
+ G_message (_("Opening iteration file [%s]"), iter_name);
+
+ if ((iter_fd = Rast_open_c_new (iter_name)) < 0)
+ G_fatal_error (_("Unable to create iterations layer [%s]"), iter_name);
+ else
+ iter_cell = Rast_allocate_c_buf();
+ }
+
+ //G_matrix_print(A);
+
+ return matrixsize;
+ */
+ return 0;
+}
+
+
+
+
+
+
+
+
+mat_struct* open_files2 (char *matrixfile,
+ char *img_grp,
+ char *result_prefix,
+ char *iter_name,
+ char *error_name)
+{
+ char result_name[80];
+
+ FILE *fp;
+ int i, matrixsize;
+ mat_struct A_input, *A;
+
+
+// mat_struct A_input1;
+
+
+ /* Read in matrix file with spectral library.
+ * Input matrix must contain spectra row-wise (for user's convenience)!
+ * Transposed here to col-wise orientation (for modules/mathematical
+ * convenience). */
+
+ if ((fp = fopen (matrixfile, "r")) == NULL)
+ G_fatal_error (_("Matrix file %s not found."), matrixfile);
+
+
+ //G_matrix_init2( &A_input,3,3,3);
+ // G_warning("matrix read start");
+// G_warning("cols=%d",A_input->cols);
+ /* Read data and close file */
+ if ((G_matrix_read2 (fp, &A_input) < 0))
+ G_fatal_error (_("Unable to read matrix file %s."), matrixfile);
+ fclose (fp);
+
+ //G_matrix_print2(&A_input, "A_input");
+
+
+ // G_warning("matrix read done");
+#if 0
+ G_message(_("Your spectral matrix = %d"), m_output(A_input));
+#endif
+
+ /* transpose input matrix from row orientation to col orientation.
+ * Don't mix rows and cols in the source code and the modules
+ * messages output! */
+
+// A = m_get(A_input->rows, A_input->cols);
+
+ //G_matrix_clone2(&A_input, A);
+
+
+
+ A = G_matrix_init (A_input.rows, A_input.cols, A_input.rows);
+ if (A == NULL)
+ G_fatal_error (_("Unable to allocate memory for matrix"));
+
+ A = G_matrix_transpose (&A_input);
+
+ // G_matrix_free (&A_input);
+
+// G_matrix_print(A);
+
+ if ((A->rows) < (A->cols))
+ G_fatal_error (_("Need number of cols >= rows to perform least squares fitting."));
+
+ // number of rows must be equivalent to no. of bands
+ matrixsize = A->rows;
+
+ // open input files from group
+ if (!I_find_group (img_grp))
+ G_fatal_error (_("Unable to find imagery group %s."), img_grp);
+
+ I_get_group_ref (img_grp, &Ref);
+ if (Ref.nfiles <= 1)
+ {
+ if (Ref.nfiles <= 0)
+ G_fatal_error (_("Group %s does not have any rasters. "
+ "The group must have at least 2 rasters."), img_grp);
+ else
+ G_fatal_error (_("Group %s only has 1 raster. "
+ "The group must have at least 2 rasters."), img_grp);
+ }
+
+ // Error check: input file number must be equal to matrix size
+ if (Ref.nfiles != matrixsize)
+ G_fatal_error (_("Number of input files (%i) in group <%s> "
+ "does not match number of spectra in matrix. "
+ "(contains %i cols)."),
+ Ref.nfiles, img_grp, A->rows);
+
+ // get memory for input files
+ cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ cellfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+ for (i = 0; i < Ref.nfiles; i++)
+ {
+ cell[i] = Rast_allocate_c_buf();
+
+ G_message (_("Opening input file no. %i [%s]"), (i + 1), Ref.file[i].name);
+
+ if ((cellfd[i] = Rast_open_old (Ref.file[i].name, Ref.file[i].mapset)) < 0)
+ G_fatal_error (_("Unable to open <%s>"), Ref.file[i].name);
+
+
+ }
+
+
+
+ // open files for results
+ result_cell = (CELL **) G_malloc (Ref.nfiles * sizeof (CELL *));
+ resultfd = (int *) G_malloc (Ref.nfiles * sizeof (int));
+
+ for (i = 0; i < A->cols; i++) // no. of spectra
+ {
+ if (result_prefix)
+ {
+ sprintf (result_name, "%s.%d", result_prefix, (i + 1));
+ G_message (_("Opening output file [%s]"), result_name);
+
+ result_cell[i] = Rast_allocate_c_buf();
+ if ((resultfd[i] = Rast_open_c_new (result_name)) < 0)
+ G_fatal_error (_("GRASS-DB internal error: Unable to proceed."));
+ }
+ }
+ // open file containing SMA error
+ error_cell = (CELL *) G_malloc (sizeof(CELL *));
+ if (error_name)
+ {
+ G_message (_("Opening error file [%s]"), error_name);
+
+ if ((error_fd = Rast_open_c_new (error_name)) < 0)
+ G_fatal_error (_("Unable to create error layer [%s]"), error_name);
+ else
+ error_cell = Rast_allocate_c_buf();
+ }
+
+ // open file containing number of iterations
+ iter_cell = (CELL *) G_malloc (sizeof(CELL *));
+ if (iter_name)
+ {
+ G_message (_("Opening iteration file [%s]"), iter_name);
+
+ if ((iter_fd = Rast_open_c_new (iter_name)) < 0)
+ G_fatal_error (_("Unable to create iterations layer [%s]"), iter_name);
+ else
+ iter_cell = Rast_allocate_c_buf();
+ }
+
+
+
+ /* give back number of output files (= Ref.nfiles) */
+ return A;
+}
Added: grass-addons/grass7/imagery/i.spec.unmix/region.unmix.txt
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/region.unmix.txt (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/region.unmix.txt 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,10 @@
+projection: 99 (Transverse Mercator)
+zone: 0
+north: 5767525
+south: 5763000
+east: 3574025
+west: 3569500
+nsres: 25
+ewres: 25
+rows: 181
+cols: 181
Added: grass-addons/grass7/imagery/i.spec.unmix/run.sh
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/run.sh (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/run.sh 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/run.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/spec_angle.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,50 @@
+/* Spectral angle mapping
+ * (c) 1999 Markus Neteler, Hannover
+ *
+ * 25. Nov. 1998 - V. 0.2
+ *
+ *****/
+
+#include <stdio.h>
+#include <math.h>
+#include <grass/gmath.h>
+#include "global.h"
+
+
+/* input mat_struct A, vec_struct 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
+ */
+
+float spectral_angle (vec_struct *Avector1, vec_struct *Avector2)
+{
+ vec_struct *vtmp1;
+ double norm1, norm2, norm3;
+
+ /* Measure spectral angle */
+
+ /* multiply one A column with second */
+// vtmp1 = G_vector_init (0, 0, RVEC);
+ vtmp1 = G_vector_product (Avector1, Avector2, vtmp1);
+ norm1 = G_vector_norm1 (vtmp1); /* calculate 1-norm */
+ norm2 = G_vector_norm_euclid (Avector1); /* calculate 2-norm (Euclidean) */
+ norm3 = G_vector_norm_euclid (Avector2); /* calculate 2-norm (Euclidean) */
+
+ G_vector_free (vtmp1);
+
+ /* Calculate angle and return in degree globally */
+ return (acos (norm1 / (norm2 * norm3)) * 180 / M_PI);
+}
Added: grass-addons/grass7/imagery/i.spec.unmix/spectrum.dat
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/spectrum.dat (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/spectrum.dat 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/test/Gmakefile
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/test/Gmakefile (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/test/Gmakefile 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,18 @@
+
+OFILES = \
+ gmath_test.o
+
+
+PROGS = gmath_test
+G2C = -lg2c
+LA = -lblas -llapack
+
+all: $(PROGS)
+
+gmath_test: gmath_test.o
+ $(CC) $(LDFLAGS) -o $@ $< $(LA) $(G2C) \
+ $(GMATHLIB) $(GISLIB) $(MATHLIB) $(XDRLIB)
+
+
+$(GISLIB): #
+$(GMATHLIB): #
Added: grass-addons/grass7/imagery/i.spec.unmix/test/gmath_test.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/test/gmath_test.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/test/gmath_test.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,99 @@
+/* Test the facility for adding, subtracting, multiplying matrices */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "gis.h"
+#include "la.h"
+
+int main(int argc, char *argv[]) {
+
+ int i, j;
+ mat_struct *m1, *m2, *m_sum, *m_sub, *m_scale, *m3, *m4;
+ vec_struct *v1;
+
+ double testmat1[5][3] = { 1.0, 3.4, 8.1,
+ 1.5, 1.5, 2.3,
+ 2.0, 2.2, 1.7,
+ 2.5, 6.3, 6.1,
+ 3.0, 1.6, 5.0
+ };
+
+ double testmat2[5][3] = { 7.3, 0.5, 2.6,
+ 6.9, 1.2, 2.8,
+ 5.5, 1.9, 5.1,
+ 9.3, 7.7, 7.0,
+ 0.5, 3.1, 3.8
+ };
+
+ double testvec1[5] = { 1 2 3 4 5 }
+
+ double testc = 4.3;
+
+
+ /* Initialise the matrix structures */
+
+ m1 = G_matrix_init(5, 3, 5);
+ m2 = G_matrix_init(5, 3, 5);
+
+ for( i = 0; i < 5; i++ )
+ for( j = 0; j < 3; j++ ) {
+
+ G_matrix_set_element(m1, i, j, testmat1[i][j]);
+ G_matrix_set_element(m2, i, j, testmat2[i][j]);
+ }
+
+ /* Add the matrices */
+
+ m_sum = G_matrix_add(m1, m2);
+
+ /* Subtract */
+
+ m_sub = G_matrix_subtract(m1, m2);
+
+ /* Scale the matrix by a given scalar value */
+
+ m_scale = G_matrix_scale(m1, testc);
+
+ /* Multiply two matrices */
+
+ m3 = G_matrix_transpose(m1);
+
+ m4 = G_matrix_product(m3, m2);
+
+
+ /* Print out the results */
+
+ printf("*** TEST OF MATRIX WRAPPER FUNCTIONS ***\n\n");
+
+ printf("1. Simple matrix manipulations\n\n");
+
+ printf(" Matrix 1:\n");
+ G_matrix_print(m1);
+
+ printf(" Matrix 2:\n");
+ G_matrix_print(m2);
+
+ printf(" Sum (m1 + m2):\n");
+ G_matrix_print(m_sum);
+
+ printf(" Difference (m1 - m2):\n");
+ G_matrix_print(m_sub);
+
+ printf(" Scale (c x m1):\n");
+ G_matrix_print(m_scale);
+
+ printf(" Multiply transpose of M1 by M2 (m1~ x m2):\n");
+ G_matrix_print(m4);
+
+ G_matrix_free(m1);
+ G_matrix_free(m2);
+ G_matrix_free(m_sum);
+ G_matrix_free(m_sub);
+ G_matrix_free(m_scale);
+ G_matrix_free(m3);
+ G_matrix_free(m4);
+
+ return 0;
+
+}
Added: grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass.dat
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass.dat (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass.dat 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/unmix_reclass20.dat
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass20.dat (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass20.dat 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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/grass7/imagery/i.spec.unmix/unmix_reclass25.dat
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass25.dat (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unmix_reclass25.dat 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,4 @@
+0 thru 25 = 13
+26 thru 50 = 38
+51 thru 75 = 63
+76 thru 100 = 88
Added: grass-addons/grass7/imagery/i.spec.unmix/unused/main.c.old.s
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unused/main.c.old.s (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unused/main.c.old.s 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,390 @@
+/* Spectral mixture analysis
+ * (c) 1999 Markus Neteler, Hannover, Germany
+ * neteler@@geog.uni-hannover.de
+ *
+ * V 1.2 - 20. Jan. 1999
+ *
+ * * 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)*/
+
+CELL round (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/grass7/imagery/i.spec.unmix/unused/matrix.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unused/matrix.h (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unused/matrix.h 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,474 @@
+/**************************************************************************
+**
+** 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.
+**
+***************************************************************************/
+
+
+#ifndef __MATRIX_H__
+#define __MATRIX_H__
+
+
+#include <grass/config.h>
+#include <grass/gis.h>
+
+/* vector definition */
+typedef struct {
+ unsigned int dim, max_dim;
+ double *ve;
+} VEC;
+
+/* matrix definition */
+typedef struct {
+ unsigned int m, n;
+ unsigned int max_m, max_n, max_size;
+ double **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 {
+ unsigned int size, max_size, *pe;
+} PERM;
+
+/* integer vector definition */
+typedef struct {
+ unsigned int dim, max_dim;
+ int *ive;
+} IVEC;
+
+
+void m_version( void );
+
+/* allocate one object of given type */
+#define NEW(type) ((type *)G_calloc((size_t)1, sizeof(type)))
+
+/* allocate num objects of given type */
+#define NEW_A(num,type) ((type *)G_calloc((size_t)(num), sizeof(type)))
+
+ /* re-allocate arry to have num objects of the given type */
+#define RENEW(var,num,type) \
+ ((var)=(type *)((var) ? \
+ G_realloc((char *)(var), (num)*sizeof(type))) : \
+ G_calloc((size_t)(num), sizeof(type))))
+
+#define MEMCOPY(from,to,n_items,type) \
+ MEM_COPY((char *)(from), (char *)(to), (unsigned)(n_items)*sizeof(type))
+
+
+/* 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
+
+
+/* for input routines */
+#define MAXLINE 81
+
+
+#if 0
+/* 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 ! */
+
+/* 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 int m_free(MAT *);
+extern int v_free(VEC *);
+extern int px_free(PERM *);
+extern int bd_free(BAND *);
+
+
+/* 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
+
+#endif
+
+/* 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 */
+/* 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);
+
+
+/* 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 */
+/* copy in to out starting at out[i0][j0] */
+extern MAT *_m_copy(MAT *in,MAT *out,unsigned int i0,unsigned 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,unsigned 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);
+
+
+/* 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 */
+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 *);
+
+/* Basic vector operations */
+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 **,double *,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,unsigned int i0),
+ /* returns sum_{i=0}^{len-1} x[i].y[i] */
+ __ip__(double *,double *,int);
+
+/* see v_mltadd(), v_add(), v_sub() and v_zero() */
+extern void __mltadd__(double *,double *, double, int),
+ __add__(double *, double *, double *, int),
+ __sub__(double *, double *, double *, int),
+ __smlt__(double *, double, double *, int),
+ __zero__(double *,int);
+
+
+/* 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 */
+ /* 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);
+
+
+/* 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 */
+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,unsigned int,unsigned int,unsigned int,unsigned 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,unsigned int i,VEC *out,unsigned int j0),
+ /* A[i][j] <- out[i], i >= i0 */
+ *_set_row(MAT *A,unsigned int j,VEC *out,unsigned int i0);
+
+extern VEC *get_row(MAT *,unsigned int,VEC *),
+ *get_col(MAT *,unsigned 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);
+
+
+/* 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 */
+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,unsigned int i,unsigned int j);
+
+ /* returns sign(px) = +1 if px product of even # transpositions
+ -1 if ps product of odd # transpositions */
+extern int px_sign(PERM *);
+
+
+/* Basic integer vector operations */
+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);
+
+
+/* miscellaneous functions */
+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(double *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);
+
+
+/* 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 */
+#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 **,...);
+
+
+#endif /* __MATRIX_H__ */
+
Added: grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unused/matrix2.h 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,174 @@
+/**************************************************************************
+**
+** 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 __MATRIX2_H__
+#define __MATRIX2_H__
+
+
+#include "matrix.h"
+
+/* Unless otherwise specified, factorisation routines overwrite the
+ matrix that is being factorised */
+
+ /* 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 *,unsigned int, double *,VEC *, double *);
+extern VEC *hhtrvec(VEC *,double,unsigned int,VEC *,VEC *);
+extern MAT *hhtrrows(MAT *,unsigned int,unsigned int,VEC *,double);
+extern MAT *hhtrcols(MAT *,unsigned int,unsigned int,VEC *,double);
+
+extern void givens(double,double, double *, double *);
+extern VEC *rot_vec(VEC *,unsigned int,unsigned int,double,double,VEC *); /* in situ */
+extern MAT *rot_rows(MAT *,unsigned int,unsigned int,double,double,MAT *); /* in situ */
+extern MAT *rot_cols(MAT *,unsigned int,unsigned 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 /* __MATRIX2_H__ */
Added: grass-addons/grass7/imagery/i.spec.unmix/unused/meminfo.h
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unused/meminfo.h (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unused/meminfo.h 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,124 @@
+/**************************************************************************
+**
+** 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.
+**
+***************************************************************************/
+
+
+#ifndef __MEM_INFO_H__
+#define __MEM_INFO_H__
+
+
+/* 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;
+
+
+
+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,...);
+
+
+/* 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 /* __MEM_INFO_H__ */
Added: grass-addons/grass7/imagery/i.spec.unmix/unused/meschach_convert.txt
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unused/meschach_convert.txt (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unused/meschach_convert.txt 2012-12-04 18:05:00 UTC (rev 54183)
@@ -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.
+
Added: grass-addons/grass7/imagery/i.spec.unmix/unused/svd_calc.c
===================================================================
--- grass-addons/grass7/imagery/i.spec.unmix/unused/svd_calc.c (rev 0)
+++ grass-addons/grass7/imagery/i.spec.unmix/unused/svd_calc.c 2012-12-04 18:05:00 UTC (rev 54183)
@@ -0,0 +1,145 @@
+/* Matrix inversion with Singular Value Decomposition
+ * (c) 1998 Markus Neteler, Hannover
+ *
+ * 20. Nov. 1998 - V 1.0
+ *
+ ****************************************************************************
+ ** 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 abs(x) ((x)<0?-(x):(x))
+
+#define GLOBAL
+#define errmesg(mesg) fprintf(stderr,"Error: %s error: line %d\n",mesg,__LINE__)
+
+#include <stdio.h>
+#include <math.h>
+#include "matrix.h"
+#include "matrix2.h"
+#include "global.h"
+
+
+void calc_inverse(M) /* returns inverted matrix MM_inverse, which is globally defined*/
+MAT *M;
+{
+ MAT *U, *U_trans, *V, *Temp1, *Temp2, *MM_inverse_pre, *Svdvals_mat_inv, *C, *D;
+ VEC *svdvals;
+ int i, j, m_dim, n_dim;
+
+ m_dim = M->m;
+ n_dim = M->n;
+ U= m_get(m_dim,m_dim); /* create empty matrices */
+ V= m_get(n_dim,n_dim);
+ svdvals = svd(M,U,V,VNULL);
+
+/* 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
+ */
+
+ /* debug output */
+ if (!flag.quiet->answer)
+ {
+ fprintf(stderr, "SVD-Vector of M: ");
+ v_output(svdvals);
+ }
+
+/* Now we have U, V, svdvals */
+
+/* convert vector svdvals to diagonal matrix 1/W: */
+ Svdvals_mat_inv = m_get(m_dim,n_dim);
+ for (i=0; i < m_dim; i++)
+ Svdvals_mat_inv->me[i][i] = (double)(1/svdvals->ve[i]);
+
+/* calculate transpose of U */
+ U_trans = m_transp(U, MNULL);
+
+ Temp1 = m_mlt(Svdvals_mat_inv, V, MNULL);
+ Temp2 = m_mlt(U_trans, Temp1, MNULL); /* why that? no idea */
+ MM_inverse_pre = m_transp(Temp2, MNULL);
+ MM_inverse = m_get(m_dim,n_dim);
+
+/* set values nearly zero to to zero */
+
+ for (i=0; i < MM_inverse_pre->m ; i++) /* rowwise */
+ for (j=0; j < MM_inverse_pre->n; j++) /* colwise */
+ {
+ if (abs(MM_inverse_pre->me[i][j]) > 0.00000000000001)
+ MM_inverse->me[i][j]= MM_inverse_pre->me[i][j];
+ else
+ MM_inverse->me[i][j]= 0;
+ }
+
+ M_FREE(Svdvals_mat_inv);
+ M_FREE(U_trans);
+ M_FREE(Temp1);
+ M_FREE(Temp2);
+ M_FREE(MM_inverse_pre);
+
+ /***************************************************************
+ * Error check for SVD
+ * check reconstruction of M - taken from torture.c (Meschach)
+ ****/
+
+ C = m_get(M->m,M->n);
+ D = m_get(M->m,M->n);
+ for ( i = 0; i < min(M->m,M->n); i++ )
+ m_set_val(D,i,i,v_entry(svdvals,i));
+ mtrm_mlt(U,D,C);
+ m_mlt(C,V,D);
+ m_sub(M,D,D);
+ if ( m_norm1(D) >= MACHEPS*m_norm1(V)*m_norm_inf(U)*m_norm1(M) )
+ if (!flag.quiet->answer)
+ fprintf(stderr, "ERROR: SVD reconstruction error = %g [allowed MACHEPS = %g]\n",\
+ m_norm1(D), MACHEPS);
+
+ /* check orthogonality of U and V */
+ M_FREE(D);
+ D = m_resize(D,U->n,U->n);
+ mtrm_mlt(U,U,D);
+ for ( i = 0; i < D->m; i++ )
+ m_set_val(D,i,i,m_entry(D,i,i)-1.0);
+ if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(U)*5 )
+ if (!flag.quiet->answer)
+ fprintf(stderr, "ERROR: SVD orthogonality error (V) = %g [allowed MACHEPS = %g\n",\
+ m_norm1(D), MACHEPS);
+
+ M_FREE(D);
+ D = m_resize(D,V->n,V->n);
+ mtrm_mlt(V,V,D);
+ for ( i = 0; i < D->m; i++ )
+ m_set_val(D,i,i,m_entry(D,i,i)-1.0);
+ if ( m_norm1(D) >= MACHEPS*m_norm1(V)*m_norm_inf(V)*5 )
+ if (!flag.quiet->answer)
+ fprintf(stderr, "ERROR: SVD orthogonality error (U) = %g [allowed MACHEPS = %g\n",\
+ m_norm1(D), MACHEPS);
+
+ for ( i = 0; i < svdvals->dim; i++ )
+ if ( v_entry(svdvals,i) < 0 || (i < svdvals->dim-1 &&
+ v_entry(svdvals,i+1) > v_entry(svdvals,i)) )
+ break;
+
+ /* check next only, when svdvals->dim high enough: */
+ if (( i < svdvals->dim ) && (svdvals->dim > 2))
+ if (!flag.quiet->answer)
+ fprintf(stderr, "ERROR: SVD sorting error\n");
+
+
+ /*free the memory*/
+ M_FREE(V);
+ M_FREE(U);
+ V_FREE(svdvals);
+ M_FREE(C);
+ M_FREE(D);
+}
More information about the grass-commit
mailing list