[GRASS-SVN] r31308 - in grass-addons/gipe: . i.lmf

svn_grass at osgeo.org svn_grass at osgeo.org
Sat May 10 05:03:18 EDT 2008


Author: ychemin
Date: 2008-05-10 05:03:18 -0400 (Sat, 10 May 2008)
New Revision: 31308

Added:
   grass-addons/gipe/i.lmf/
   grass-addons/gipe/i.lmf/Makefile
   grass-addons/gipe/i.lmf/description.html
   grass-addons/gipe/i.lmf/fitting.c
   grass-addons/gipe/i.lmf/invert_matrix.c
   grass-addons/gipe/i.lmf/lmf.c
   grass-addons/gipe/i.lmf/main.c
   grass-addons/gipe/i.lmf/make_matrix.c
   grass-addons/gipe/i.lmf/maxmin.c
   grass-addons/gipe/i.lmf/minmax.c
Modified:
   grass-addons/gipe/imagery_Makefile
Log:
Added i.lmf temporal splining with maximum filter (initially made for vegetation indices)

Added: grass-addons/gipe/i.lmf/Makefile
===================================================================
--- grass-addons/gipe/i.lmf/Makefile	                        (rev 0)
+++ grass-addons/gipe/i.lmf/Makefile	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.lmf
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+	EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd

Added: grass-addons/gipe/i.lmf/description.html
===================================================================
--- grass-addons/gipe/i.lmf/description.html	                        (rev 0)
+++ grass-addons/gipe/i.lmf/description.html	2008-05-10 09:03:18 UTC (rev 31308)
@@ -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: 2008/05/10 15:51:43 $</i>

Added: grass-addons/gipe/i.lmf/fitting.c
===================================================================
--- grass-addons/gipe/i.lmf/fitting.c	                        (rev 0)
+++ grass-addons/gipe/i.lmf/fitting.c	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,73 @@
+/*************************************************************
+******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/gipe/i.lmf/invert_matrix.c
===================================================================
--- grass-addons/gipe/i.lmf/invert_matrix.c	                        (rev 0)
+++ grass-addons/gipe/i.lmf/invert_matrix.c	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,108 @@
+/****************************************************************
+ * 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/gipe/i.lmf/lmf.c
===================================================================
--- grass-addons/gipe/i.lmf/lmf.c	                        (rev 0)
+++ grass-addons/gipe/i.lmf/lmf.c	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,113 @@
+/********************************************************************
+* 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/gipe/i.lmf/main.c
===================================================================
--- grass-addons/gipe/i.lmf/main.c	                        (rev 0)
+++ grass-addons/gipe/i.lmf/main.c	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,228 @@
+/****************************************************************************
+ *
+ * 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/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
+	/************************************/
+	/* FMEO Declarations*****************/
+	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;
+
+	char	*fileName;
+#define fileNameLe 8
+#define fileNamePosition 3
+
+	int	ndates;
+	double	inpix[MAXFILES]={0.0};
+	double	outpix[MAXFILES]={0.0};
+	/************************************/
+
+	G_gisinit(argv[0]);
+
+	module = G_define_module();
+	module->keywords = _("LMF,Vegetation Indices,Atmospheric temporal correction");
+	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_INPUT) ;
+	input->multiple   = YES;
+	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");
+
+	/* FMEO init nfiles */
+	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;
+		/* find map in mapset */
+		mapset = G_find_cell2 (name, "");
+	        if (mapset == NULL)
+		{
+			G_fatal_error (_("cell file [%s] not found"), name);
+			ok = 0;
+		}
+		if (!ok){
+			continue;
+		}
+		infd[nfiles] = G_open_cell_old (name, mapset);
+		if (infd[nfiles] < 0)
+		{
+			ok = 0;
+			continue;
+		}
+		/* Allocate input buffer */
+		in_data_type[nfiles] = G_raster_map_type(name, mapset);
+		if( (infd[nfiles] = G_open_cell_old(name,mapset)) < 0){
+			G_fatal_error(_("Cannot open cell file [%s]"), name);
+		}
+		if( (G_get_cellhd(name,mapset,&cellhd)) < 0){
+			G_fatal_error(_("Cannot read file header of [%s]"), name);
+		}
+		inrast[nfiles] = G_allocate_raster_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 = G_window_rows();
+	ncols = G_window_cols();
+	for(i=0;i<nfiles;i++){
+		sprintf(result[i],output->answer,".",i+1);
+		if (G_legal_filename (result[i]) < 0)
+		{
+			G_fatal_error (_("[%s] is an illegal name"), result[i]);
+			ok = 0;
+		}
+		outrast[i] = G_allocate_raster_buf(out_data_type);
+		/* Create New raster files */
+		if ( (outfd[i] = G_open_raster_new (result[i],1)) < 0){
+			G_fatal_error (_("Could not open <%s>"),result[i]);
+		}
+	}
+	/*******************/
+	/* 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++)
+		{
+			if( (G_get_raster_row(infd[i],inrast[i],row,in_data_type[i])) < 0){
+				G_fatal_error (_("Could not read from <%s>"),name);
+			}
+		}
+		/*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++){
+			if (G_put_raster_row (outfd[i], outrast[i], out_data_type) < 0)
+				G_fatal_error (_("Cannot write to <%s>"),result[i]);
+		}
+	}
+	for (i=1;i<=nfiles;i++)
+	{
+		G_free (inrast[i]);
+		G_close_cell (infd[i]);
+		G_free (outrast[i]);
+		G_close_cell (outfd[i]);
+	}
+	return 0;
+}

Added: grass-addons/gipe/i.lmf/make_matrix.c
===================================================================
--- grass-addons/gipe/i.lmf/make_matrix.c	                        (rev 0)
+++ grass-addons/gipe/i.lmf/make_matrix.c	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,48 @@
+// Making Triangular Function Matrix/Vector
+
+// double  f[NMAX][MAXF];
+
+#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;
+	
+//	printf("make_matrix: n=%d\tnpoint=%d\tnfunc=%d\n",n,npoint,nfunc);
+        nn=2*nfunc+2;
+//	printf("make_matrix: nn=%d\n",nn);
+        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/gipe/i.lmf/maxmin.c
===================================================================
--- grass-addons/gipe/i.lmf/maxmin.c	                        (rev 0)
+++ grass-addons/gipe/i.lmf/maxmin.c	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,45 @@
+/* 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/gipe/i.lmf/minmax.c
===================================================================
--- grass-addons/gipe/i.lmf/minmax.c	                        (rev 0)
+++ grass-addons/gipe/i.lmf/minmax.c	2008-05-10 09:03:18 UTC (rev 31308)
@@ -0,0 +1,45 @@
+/* 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;
+}

Modified: grass-addons/gipe/imagery_Makefile
===================================================================
--- grass-addons/gipe/imagery_Makefile	2008-05-10 06:16:03 UTC (rev 31307)
+++ grass-addons/gipe/imagery_Makefile	2008-05-10 09:03:18 UTC (rev 31308)
@@ -50,6 +50,7 @@
 	i.group \
 	i.his.rgb \
 	i.latitude \
+	i.lmf \
 	i.longitude \
 	i.landsat.dehaze \
 	i.maxlik \



More information about the grass-commit mailing list