[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