[GRASS-SVN] r51067 - in grass-addons: grass6/imagery/gipe
grass7/imagery grass7/imagery/i.evapo.zk grass7/imagery/i.lmf
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Mar 15 06:16:18 EDT 2012
Author: ychemin
Date: 2012-03-15 03:16:18 -0700 (Thu, 15 Mar 2012)
New Revision: 51067
Added:
grass-addons/grass7/imagery/i.evapo.zk/
grass-addons/grass7/imagery/i.lmf/
grass-addons/grass7/imagery/i.lmf/Makefile
grass-addons/grass7/imagery/i.lmf/fitting.c
grass-addons/grass7/imagery/i.lmf/i.lmf.html
grass-addons/grass7/imagery/i.lmf/invert_matrix.c
grass-addons/grass7/imagery/i.lmf/lmf.c
grass-addons/grass7/imagery/i.lmf/main.c
grass-addons/grass7/imagery/i.lmf/make_matrix.c
grass-addons/grass7/imagery/i.lmf/maxmin.c
grass-addons/grass7/imagery/i.lmf/minmax.c
Removed:
grass-addons/grass6/imagery/gipe/i.evapo.zk/
grass-addons/grass6/imagery/gipe/i.lmf/
Modified:
grass-addons/grass6/imagery/gipe/readme.gipe
grass-addons/grass7/imagery/i.evapo.zk/main.c
Log:
more maintenance on gipe, upgrade some modules to G7 addons
Modified: grass-addons/grass6/imagery/gipe/readme.gipe
===================================================================
--- grass-addons/grass6/imagery/gipe/readme.gipe 2012-03-15 09:48:13 UTC (rev 51066)
+++ grass-addons/grass6/imagery/gipe/readme.gipe 2012-03-15 10:16:18 UTC (rev 51067)
@@ -29,6 +29,7 @@
20120315: Created i.vi.mpi in addons/grass7/imagery/, tested with "ndvi"
: Moved i.evapo.potrad/senay to grass7, removed m.gem and i.vi.mpi from grass6.
: Moved i.eb.h_SEBAL95 to grass7
+ : Moved i.lmf to grass7
Imagery modules
---------------
@@ -36,7 +37,7 @@
**i.atcorr: (ported from GRASS 5.x, now in Main SVN) Atmospheric correction
**i.biomass: (Now in Main SVN) Calculates biomass growth of a plant/crop
**i.emissivity: (now in main SVN) Calculates the generic emissivity from NDVI
-i.lmf: Calculates NDVI temporal splinning after Sawada et al. (2001)
+**i.lmf: (Grass7 version in addons) Calculates NDVI temporal splinning after Sawada et al. (2001)
**i.vi: (now in Main SVN) Calculates any of 13 types of vegetation indices
**i.vi.mpi: (Grass7 version in addons) Educational version in MPICH for HPC
i.vi.grid: Educational version in Ninf/G for grid computing
Modified: grass-addons/grass7/imagery/i.evapo.zk/main.c
===================================================================
--- grass-addons/grass6/imagery/gipe/i.evapo.zk/main.c 2012-03-14 23:30:31 UTC (rev 51063)
+++ grass-addons/grass7/imagery/i.evapo.zk/main.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -4,9 +4,9 @@
* AUTHOR: Yann Chemin yann.chemin at gmail.com
*
* PURPOSE: To estimate the daily evapotranspiration by means
-* of Prestley and Taylor method (1972).
+* of Zhang and Kimberley.
*
-* COPYRIGHT: (C) 2007-2011 by the GRASS Development Team
+* COPYRIGHT: (C) 2007-2012 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* Licence (>=2). Read the file COPYING that comes with GRASS
Added: grass-addons/grass7/imagery/i.lmf/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.lmf/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/Makefile 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.lmf
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/imagery/i.lmf/fitting.c
===================================================================
--- grass-addons/grass7/imagery/i.lmf/fitting.c (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/fitting.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,88 @@
+
+/*************************************************************
+******Triangular Function Fitting********************
+*************************************************************/
+#include<stdio.h>
+#include<stdlib.h>
+
+#define NMAX 200
+#define MAXF 50
+int invert_matrix(double a[][MAXF], int order);
+
+int fitting(int npoint, int nfunc, double *dat, int *idx1, double f[][MAXF],
+ double *c, double *vfit)
+{
+ int i, j, k, k1, k2, nn;
+
+ double tf[MAXF][MAXF];
+
+ double tinv[MAXF][MAXF];
+
+ double vec[MAXF];
+
+ double tfij;
+
+ double sum;
+
+ nn = nfunc * 2 + 2;
+
+ /*Value Initialization */
+ for (i = 0; i < nn; i++) {
+ c[i] = 0.0;
+ }
+ for (i = 0; i < npoint; i++) {
+ vfit[i] = 0.0;
+ }
+
+ /*Matrix Initialization */
+ for (i = 0; i < nn; i++) {
+ for (j = 0; j < nn; j++) {
+ tf[i][j] = 0.0;
+ tf[j][i] = 0.0;
+ }
+ vec[i] = 0.0;
+ }
+
+ /*Making Matrix */
+ for (i = 0; i < npoint; i++) {
+ if (idx1[i] == 1) {
+ for (k1 = 0; k1 < nn; k1++) {
+ for (k2 = 0; k2 < nn; k2++) {
+ tf[k1][k2] = tf[k1][k2] + f[i][k2] * f[i][k1];
+ }
+ vec[k1] = vec[k1] + dat[i] * f[i][k1];
+ }
+ }
+ }
+
+ /*Matrix Copy and Inversion */
+ for (i = 0; i < nn; i++) {
+ for (j = i; j < nn; j++) {
+ tfij = tf[i][j];
+ tinv[i][j] = tfij;
+ tinv[j][i] = tfij;
+ }
+ }
+ invert_matrix(tinv, nn);
+
+ /*Calculation of Coefficients */
+ for (i = 0; i < nn; i++) {
+ sum = 0.0;
+ for (j = 0; j < nn; j++) {
+ sum = sum + tinv[i][j] * vec[j];
+ }
+ c[i] = sum;
+ }
+
+ /*Calculating theoretical value */
+ for (i = 0; i < npoint; i++) {
+ sum = 0.0;
+ for (k = 0; k < nn; k++) {
+ sum = sum + c[k] * f[i][k];
+ }
+ vfit[i] = sum;
+ }
+ return;
+}
+
+
Added: grass-addons/grass7/imagery/i.lmf/i.lmf.html
===================================================================
--- grass-addons/grass7/imagery/i.lmf/i.lmf.html (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/i.lmf.html 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,36 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.lmf</em> calculates the Local maximum fitting of a temporal time-series, intially for vegetation indices, it also works for surface reflectance.
+
+This is a first level port, only a fast fitting is done, see TODO.
+
+<h2>NOTES</h2>
+Original links are found here
+
+SAWADA, 2001:
+http://www.affrc.go.jp/ANDES/sawady/index.html
+
+NAGATANI et al., 2002:
+http://www.gisdevelopment.net/aars/acrs/2002/pos2/184.pdf
+
+Yann Chemin and Kiyoshi Honda repaired it and ported it from SGI/OpenMP to Linux.
+http://www.rsgis.ait.ac.th/~honda/lmf/lmf.html
+
+
+<h2>TODO</h2>
+Port the full detailed algorithm from Fortran, and vastly unemcomber/clean it.
+It will make the algorithm must slower though, gaining only a marginal fitting strength, for my actual experience with VIs curves.
+
+<h2>SEE ALSO</h2>
+
+<em>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Yann Chemin, International Rice Research Institute, The Philippines<br>
+
+
+<p>
+<i>Last changed: $Date: 2011-11-09 03:26:45 +0530 (Wed, 09 Nov 2011) $</i>
Added: grass-addons/grass7/imagery/i.lmf/invert_matrix.c
===================================================================
--- grass-addons/grass7/imagery/i.lmf/invert_matrix.c (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/invert_matrix.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,126 @@
+
+/****************************************************************
+ * Inverse and Determinant of A by the Gauss-Jordan Method.
+ * m is the order of the square matrix, A.
+ * A-Inverse replaces A.
+ * Determinant of A is placed in DET.
+ * Cooley and Lohnes (1971:63)
+ ****************************************************************/
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+
+#define IDIM 50
+#define MAXF 50
+int invert_matrix(double a[][MAXF], int order)
+{
+ double ipvt[IDIM + 1];
+
+ double pvt[IDIM + 1];
+
+ double ind[IDIM + 1][1];
+
+ int i, j, k, l, l1;
+
+ double amax, swap;
+
+ int irow, icol;
+
+ int m = order;
+
+ for (j = 0; j < m; j++) {
+ ipvt[j] = 0.0;
+ }
+ for (i = 0; i < m; i++) {
+
+ /*SEARCH FOR THE PIVOT ELEMENT */
+ double temporary = 0;
+
+ amax = 0.0;
+ for (j = 0; j < m; j++) {
+ if (ipvt[j] != 1) {
+ for (k = 0; k < m; k++) {
+ if (ipvt[k] == 0 && abs(amax) < abs(a[j][k])) {
+ irow = j;
+ icol = k;
+ amax = a[j][k];
+ }
+ else if (ipvt[j] == 1) {
+ temporary = 1;
+ break;
+ }
+ else if (ipvt[j] > 1 || ipvt[j] < 0) {
+ temporary = 5;
+ break;
+ }
+ }
+ if (temporary == 5 || temporary == 1) {
+ break;
+ }
+ }
+ else {
+ break;
+ }
+ if (temporary == 5) {
+ break;
+ }
+ }
+ if (temporary != 5) {
+
+ /*INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL */
+ if (irow != icol) {
+ for (l = 0; l < m; l++) {
+ swap = a[irow][l];
+ a[irow][l] = a[icol][l];
+ a[icol][l] = swap;
+ }
+ }
+ ind[i][0] = irow;
+ ind[i][1] = icol;
+ pvt[i] = a[icol][icol];
+
+ /*DIVIDE THE PIVOT ROW BY THE PIVOT ELEMENT */
+ a[icol][icol] = 1.0;
+ for (l = 0; l < m; l++) {
+ a[icol][l] = a[icol][l] / pvt[i];
+ }
+
+ /*REDUCE THE NON-PIVOT ROWS */
+ for (l1 = 0; l1 < m; l1++) {
+ if (l1 != icol) {
+ swap = a[l1][icol];
+ a[l1][icol] = 0.0;
+ if (swap < pow(0.1, -30) && swap > pow(-0.1, -30)) {
+ swap = 0.0;
+ }
+ for (l = 0; l < m; l++) {
+ a[l1][l] = a[l1][l] - a[icol][l] * swap;
+ }
+ }
+ }
+ }
+ else {
+ break;
+ }
+ }
+
+ /*INTERCHANGE THE COLUMNS */
+ for (i = 0; i < m; i++) {
+ l = m - i - 1;
+ if (ind[l][0] != ind[l][1]) {
+ irow = ind[l][0];
+ icol = ind[l][1];
+ for (k = 0; k < m; k++) {
+ swap = a[k][irow];
+ a[k][irow] = a[k][icol];
+ a[k][icol] = swap;
+ }
+ }
+ else if (ind[l][0] == ind[l][1]) {
+ break;
+ }
+ }
+ return;
+}
+
+
Added: grass-addons/grass7/imagery/i.lmf/lmf.c
===================================================================
--- grass-addons/grass7/imagery/i.lmf/lmf.c (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/lmf.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,148 @@
+
+/********************************************************************
+* LMF SUBROUTINE: Trying to process one pixel at a time ****
+* in a single function ****
+********************************************************************
+********************************************************************
+* NBANDS = how many dates in this time-series ****
+* NPOINT = how many dates of this time-series cover one year ****
+* (NPOINT=46 for MODIS 8-days product) ****
+********************************************************************/
+#include<stdio.h>
+#include<stdlib.h>
+
+#define MAXB 200
+ /* NMAX=MAXB */
+#define NMAX 200
+#define MAXF 50
+
+#define TRUE 1
+#define FALSE 0
+
+#define PI 3.1415926535897932385
+int make_matrix(int n, int npoint, int nfunc, int *numk, double f[][MAXF]);
+
+int fitting(int npoint, int nfunc, double *dat, int *idx1, double f[][MAXF],
+ double *c, double *vfit);
+int minmax(int n, int nwin, double *dat);
+
+int maxmin(int n, int nwin, double *dat);
+
+int lmf(int nbands, int npoint, double *inpix, double *outpix)
+{
+ int i, j, k, idk, norder, nfunc;
+
+ int nwin, mmsw, nds, isum;
+
+ int fitsw, fitmd;
+
+ double tcld, thh, thl;
+
+ double vaic, delta;
+
+ int idx1[NMAX];
+
+
+ /*for fitting */
+ double dat[NMAX];
+
+ double f[NMAX][MAXF], c[MAXF];
+
+ double vfit[NMAX];
+
+ double ipoint[MAXB];
+
+
+ /*for iteration */
+ double dat0[NMAX], c0[MAXF];
+
+ double dat1[NMAX];
+
+
+ /*for Matrix Inversion */
+ int numk[MAXF] =
+ { 1, 2, 3, 4, 6, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0 };
+
+ /*Set Defaults */
+ nfunc = 6;
+ for (i = 0; i < NMAX; i++) {
+ for (j = 0; j < MAXF; j++) {
+ f[i][j] = 0.0;
+ }
+ }
+ make_matrix(nbands, npoint, nfunc, numk, f);
+ nwin = 3;
+ mmsw = 1;
+
+ /*nds=npoint; */
+ nds = 2 * nfunc + 2;
+
+ /*LMF Process */
+ /*Checking Pixel Value */
+ fitsw = 1;
+
+ /*fitting method */
+ tcld = 0.15;
+ fitmd = 2;
+
+ /*FITMD0---- */
+ if (fitmd == 0) {
+ thh = tcld * 2.0;
+ thl = -tcld;
+ }
+
+ /*FITMD1---- */
+ if (fitmd == 1) {
+ thh = tcld;
+ thl = -tcld * 2.0;
+ }
+
+ /*FITMD2---- */
+ if (fitmd == 2) {
+ thh = tcld;
+ thl = -tcld;
+ }
+
+ /*initialize arrays */
+ for (k = 0; k < nbands; k++) {
+ dat0[k] = 0.0;
+ outpix[k] = 0.0;
+ ipoint[k] = 1;
+ idx1[k] = TRUE;
+ }
+ for (k = 0; k < nds; k++) {
+ c[k] = 0.0;
+ }
+ for (k = 0; k < nbands; k++) {
+ if (inpix[k] > 0.375) {
+ ipoint[k] = 0;
+ }
+ if (k >= 2) {
+ idk = (inpix[k] - inpix[k - 1]);
+ if (idk > 0.125) {
+ ipoint[k] = 0;
+ }
+ }
+ dat0[k] = inpix[k];
+ dat[k] = dat0[k];
+ }
+
+ /*Minimax Filter */
+ if (mmsw >= 1 && fitmd == 0) {
+ minmax(nbands, nwin, dat);
+ }
+ else if (mmsw >= 1 && fitmd == 1) {
+ maxmin(nbands, nwin, dat);
+ }
+
+ /*Fitting */
+ fitting(nbands, 3, dat, idx1, f, c, vfit);
+ for (k = 0; k < nbands; k++) {
+ outpix[k] = vfit[k];
+ }
+ return;
+}
+
+
Added: grass-addons/grass7/imagery/i.lmf/main.c
===================================================================
--- grass-addons/grass7/imagery/i.lmf/main.c (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/main.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,189 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.lmf
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Clean temporal signature using what is called
+ * Local Maximum Fitting (LMF)
+ * Initially created for Vegetation Indices from AVHRR.
+ *
+ * Translated from Fortran, unconfused, removed unused and
+ * other redundant variables. Removed BSQ Binary loading.
+ *
+ * SHORT HISTORY OF THE CODE
+ * -------------------------
+ * Original Fortran Beta-level code was written
+ * by Yoshito Sawada (2000), in OpenMP Fortran
+ * for SGI Workstation. Recovered broken code from
+ * their website in 2003, regenerated missing code and
+ * made it work in Linux.
+ * ORIGINAL WEBSITE
+ * ----------------
+ * http://www.affrc.go.jp/ANDES/sawady/index.html
+ * ENGLISH WEBSITE
+ * ---------------
+ * http://www.rsgis.ait.ac.th/~honda/lmf/lmf.html
+ *
+ * COPYRIGHT: (C) 2008 by the GRASS Development Team
+ *
+ * This program is free software under the GNU Lesser General Public
+ * License. Read the file COPYING that comes with GRASS for details.
+ *
+ *****************************************************************************/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+#define MAXFILES 366
+
+int lmf(int nbands, int npoint, double *inpix, double *outpix);
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /*region+header info */
+ char *mapset; /*mapset name */
+ int nrows, ncols;
+ int row, col;
+ struct GModule *module;
+ struct Option *input, *ndate, *output;
+ struct History history; /*metadata */
+ char *name; /*input raster name */
+ char *result[MAXFILES]; /*output raster name */
+ /*File Descriptors */
+ int nfiles;
+ int infd[MAXFILES];
+ int outfd[MAXFILES];
+ char **names;
+ char **ptr;
+ int ok;
+ int i = 0, j = 0;
+ void *inrast[MAXFILES];
+ DCELL *outrast[MAXFILES];
+ int data_format; /* 0=double 1=float 2=32bit signed int 5=8bit unsigned int (ie text) */
+ RASTER_MAP_TYPE in_data_type[MAXFILES]; /* 0=numbers 1=text */
+ RASTER_MAP_TYPE out_data_type = DCELL_TYPE;
+ int ndates;
+ double inpix[MAXFILES] = { 0.0 };
+ double outpix[MAXFILES] = { 0.0 };
+ /************************************/
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("LMF"));
+ G_add_keyword(_("Vegetation Indices"));
+ G_add_keyword(_("Atmospheric correction"));
+ G_add_keyword(_("Temporal"));
+ module->description =
+ _("Temporal Local Maximum Fitting of vegetation indices, works also for surface reflectance data. Number of bands is potentially several years, nfiles and ndates are respectively the number of pixels and the number of pixels in a year\n");
+
+ /* Define the different options */
+
+ input = G_define_standard_option(G_OPT_R_INPUTS);
+ input->description = _("Names of input layers");
+
+ ndate = G_define_option();
+ ndate->key = _("ndate");
+ ndate->type = TYPE_INTEGER;
+ ndate->required = YES;
+ ndate->gisprompt = _("parameter, integer");
+ ndate->description = _("Number of map layers per year");
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output->description = _("Name of the output layer");
+
+ nfiles = 1;
+ /********************/
+ if (G_parser(argc, argv))
+ exit(-1);
+
+ ok = 1;
+ names = input->answers;
+ ptr = input->answers;
+
+ ndates = atoi(ndate->answer);
+
+ for (; *ptr != NULL; ptr++) {
+ if (nfiles >= MAXFILES)
+ G_fatal_error(_("%s - too many ETa files. Only %d allowed"),
+ G_program_name(), MAXFILES);
+ name = *ptr;
+ if (!ok) {
+ continue;
+ }
+ infd[nfiles] = Rast_open_old(name, "");
+ if (infd[nfiles] < 0) {
+ ok = 0;
+ continue;
+ }
+ /* Allocate input buffer */
+ in_data_type[nfiles] = Rast_map_type(name, "");
+ infd[nfiles] = Rast_open_old(name, "");
+ Rast_get_cellhd(name, "", &cellhd);
+ inrast[nfiles] = Rast_allocate_buf(in_data_type[nfiles]);
+ nfiles++;
+ }
+ nfiles--;
+ if (nfiles <= 10) {
+ G_fatal_error(_("The min specified input map is ten"));
+ }
+
+ /***************************************************/
+ /* Allocate output buffer, use input map data_type */
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ for (i = 0; i < nfiles; i++) {
+ sprintf(result[i], output->answer, ".", i + 1);
+ outrast[i] = Rast_allocate_buf(out_data_type);
+ /* Create New raster files */
+ outfd[i] = Rast_open_new(result[i], 1);
+ }
+ /*******************/
+ /* Process pixels */
+ for (row = 0; row < nrows; row++) {
+ DCELL de;
+ DCELL d[MAXFILES];
+
+ G_percent(row, nrows, 2);
+ /* read input map */
+ for (i = 1; i <= nfiles; i++) {
+ Rast_get_row(infd[i], inrast[i], row, in_data_type[i]);
+ }
+ /*process the data */
+ for (col = 0; col < ncols; col++) {
+ for (i = 1; i <= nfiles; i++) {
+ switch (in_data_type[i]) {
+ case CELL_TYPE:
+ inpix[i - 1] = (double)((CELL *) inrast[i])[col];
+ break;
+ case FCELL_TYPE:
+ inpix[i - 1] = (double)((FCELL *) inrast[i])[col];
+ break;
+ case DCELL_TYPE:
+ inpix[i - 1] = (double)((DCELL *) inrast[i])[col];
+ break;
+ }
+ }
+ lmf(nfiles, ndates, inpix, outpix);
+ /* Put the resulting temporal curve
+ * in the output layers */
+ for (i = 0; i < nfiles; i++) {
+ outrast[i][col] = outpix[i];
+ }
+ }
+ for (i = 0; i < nfiles; i++) {
+ Rast_put_row(outfd[i], outrast[i], out_data_type);
+ }
+ }
+ for (i = 1; i <= nfiles; i++) {
+ G_free(inrast[i]);
+ Rast_close(infd[i]);
+ G_free(outrast[i]);
+ Rast_close(outfd[i]);
+ }
+ return 0;
+}
Added: grass-addons/grass7/imagery/i.lmf/make_matrix.c
===================================================================
--- grass-addons/grass7/imagery/i.lmf/make_matrix.c (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/make_matrix.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,50 @@
+/* Making Triangular Function Matrix/Vector */
+
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+
+#define NMAX 200
+#define MAXF 50
+
+#define PI 3.1415926535897932385
+int make_matrix(int n, int npoint, int nfunc, int *numk, double f[][MAXF])
+{
+ int nn;
+
+ double eps, vn, angle, theta;
+
+ int i, j, j1, j2;
+
+ double di, dk;
+
+ nn = 2 * nfunc + 2;
+ eps = 1.0;
+
+ /*Set Constants */
+ vn = (double)npoint;
+ angle = 2.0 * PI / vn;
+
+ /*Matrix Clear */
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < nn; j++) {
+ f[i][j] = 0.0;
+ }
+ }
+
+ /*Making Matrix */
+ for (i = 0; i < n; i++) {
+ di = (double)(i + 1);
+ for (j = 0; j < nfunc; j++) {
+ j1 = (j + 1) * 2;
+ j2 = j1 + 1;
+ dk = (double)numk[j];
+ theta = angle * di * dk;
+ f[i][j1] = sin(theta);
+ f[i][j2] = cos(theta);
+ } f[i][0] = 1.0;
+ f[i][1] = 1.0E-1 * (double)(i + 1);
+ } return;
+}
+
+
Added: grass-addons/grass7/imagery/i.lmf/maxmin.c
===================================================================
--- grass-addons/grass7/imagery/i.lmf/maxmin.c (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/maxmin.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,51 @@
+/* Maxmin Filter */
+#include<stdio.h>
+#include<stdlib.h>
+
+#define NMAX 200
+int maxmin(int n, int nwin, double *dat)
+{
+ int i, j, jf, jr;
+
+ double da1[NMAX];
+
+ double vmin0 = 9.9E99;
+
+ double vminf, vminr;
+
+ for (i = 0; i < nwin; i++) {
+ da1[i] = dat[0];
+ }
+ for (i = 0; i < n; i++) {
+ da1[i + nwin] = dat[i];
+ }
+ for (i = 0; i < nwin; i++) {
+ da1[i + n + nwin] = dat[n];
+ }
+ for (i = 0; i < n; i++) {
+ vminf = vmin0;
+ vminr = vmin0;
+ for (j = 0; j < nwin; j++) {
+ jf = i + nwin + j;
+ jr = i + nwin - j;
+ if (da1[jf] < vminf) {
+ vminf = da1[jf];
+ }
+ if (da1[jr] < vminr) {
+ vminr = da1[jr];
+ }
+ }
+ if (vminf >= vminr) {
+ dat[i] = vminf;
+ }
+ else {
+ dat[i] = vminr;
+ }
+ if (dat[i] > 1.0E10) {
+ dat[i] = 0.0;
+ }
+ }
+ return;
+}
+
+
Added: grass-addons/grass7/imagery/i.lmf/minmax.c
===================================================================
--- grass-addons/grass7/imagery/i.lmf/minmax.c (rev 0)
+++ grass-addons/grass7/imagery/i.lmf/minmax.c 2012-03-15 10:16:18 UTC (rev 51067)
@@ -0,0 +1,51 @@
+/* Minmax Filter */
+#include<stdio.h>
+#include<stdlib.h>
+
+#define NMAX 200
+int minmax(int n, int nwin, double *dat)
+{
+ int i, j, jf, jr;
+
+ double da1[NMAX];
+
+ double vmax0 = -9.9E99;
+
+ double vmaxf, vmaxr;
+
+ for (i = 0; i < nwin; i++) {
+ da1[i] = dat[0];
+ }
+ for (i = 0; i < n; i++) {
+ da1[i + nwin] = dat[i];
+ }
+ for (i = 0; i < nwin; i++) {
+ da1[i + n + nwin] = dat[n];
+ }
+ for (i = 0; i < n; i++) {
+ vmaxf = vmax0;
+ vmaxr = vmax0;
+ for (j = 0; j < nwin; j++) {
+ jf = i + nwin + j;
+ jr = i + nwin - j;
+ if (da1[jf] > vmaxf) {
+ vmaxf = da1[jf];
+ }
+ if (da1[jr] > vmaxr) {
+ vmaxr = da1[jr];
+ }
+ }
+ if (vmaxf <= vmaxr) {
+ dat[i] = vmaxf;
+ }
+ else {
+ dat[i] = vmaxr;
+ }
+ if (dat[i] <= -1.0E10) {
+ dat[i] = 0.0;
+ }
+ }
+ return;
+}
+
+
More information about the grass-commit
mailing list