[GRASS-SVN] r56148 - in grass-addons/grass7/imagery: . i.wavelet
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue May 7 04:09:23 PDT 2013
Author: ychemin
Date: 2013-05-07 04:09:23 -0700 (Tue, 07 May 2013)
New Revision: 56148
Added:
grass-addons/grass7/imagery/i.wavelet/
grass-addons/grass7/imagery/i.wavelet/Makefile
grass-addons/grass7/imagery/i.wavelet/arrays.c
grass-addons/grass7/imagery/i.wavelet/arrays.h
grass-addons/grass7/imagery/i.wavelet/i.wavelet.html
grass-addons/grass7/imagery/i.wavelet/main.c
grass-addons/grass7/imagery/i.wavelet/w_daubechies.h
grass-addons/grass7/imagery/i.wavelet/wt.c
grass-addons/grass7/imagery/i.wavelet/wt.h
grass-addons/grass7/imagery/i.wavelet/wt_haar.c
grass-addons/grass7/imagery/i.wavelet/wt_haar.h
Modified:
grass-addons/grass7/imagery/Makefile
Log:
wavelet decomposition, and multi-source recomposition in spectral/temporal dimension
Modified: grass-addons/grass7/imagery/Makefile
===================================================================
--- grass-addons/grass7/imagery/Makefile 2013-05-07 10:31:11 UTC (rev 56147)
+++ grass-addons/grass7/imagery/Makefile 2013-05-07 11:09:23 UTC (rev 56148)
@@ -14,7 +14,8 @@
i.points.auto \
i.rotate \
i.segment \
- i.vi.mpi
+ i.vi.mpi \
+ i.wavelet
include $(MODULE_TOPDIR)/include/Make/Dir.make
Added: grass-addons/grass7/imagery/i.wavelet/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/Makefile 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.wavelet
+
+LIBES = $(RASTERLIB) $(GISLIB) $(IMAGERYLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP) $(IMAGERYDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+
+default: cmd
Added: grass-addons/grass7/imagery/i.wavelet/arrays.c
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/arrays.c (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/arrays.c 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,42 @@
+#include<stdio.h>
+#include<stdlib.h>
+/*Create 1D arrays*/
+int* ai1d(int X) {
+ int* A = (int*) malloc(X*sizeof(int*));
+ return A;
+}
+
+float* af1d(int X) {
+ float* A = (float*) malloc(X*sizeof(float*));
+ return A;
+}
+
+double* ad1d(int X) {
+ double* A = (double*) malloc(X*sizeof(double*));
+ return A;
+}
+
+/*Create 2D arrays*/
+int** ai2d(int X, int Y) {
+ int i;
+ int** A = (int**) malloc(X*sizeof(int*));
+ for (i = 0; i < X; i++)
+ A[i] = (int*) malloc(Y*sizeof(int));
+ return A;
+}
+
+float** af2d(int X, int Y) {
+ int i;
+ float** A = (float**) malloc(X*sizeof(float*));
+ for (i = 0; i < X; i++)
+ A[i] = (float*) malloc(Y*sizeof(float));
+ return A;
+}
+
+double** ad2d(int X, int Y) {
+ int i;
+ double** A = (double**) malloc(X*sizeof(double*));
+ for (i = 0; i < X; i++)
+ A[i] = (double*) malloc(Y*sizeof(double));
+ return A;
+}
Added: grass-addons/grass7/imagery/i.wavelet/arrays.h
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/arrays.h (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/arrays.h 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,11 @@
+#include<stdio.h>
+/** Get size of an array **/
+#define N_ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))
+/*Create 1D arrays*/
+int* ai1d(int X);
+float* af1d(int X);
+double* ad1d(int X);
+/*Create 2D arrays*/
+int** ai2d(int X, int Y);
+float** af2d(int X, int Y);
+double** ad2d(int X, int Y);
Added: grass-addons/grass7/imagery/i.wavelet/i.wavelet.html
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/i.wavelet.html (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/i.wavelet.html 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,24 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.wavelet</em> Decomposes the time-series with the requested wavelet.
+The output will have 4 components, High-Pass Level 1 and 2, Low-Pass
+Level 1 and 2. Low Pass Level 1 can be recomposed using HP2 and LP2
+using the recomposition mode.<br>
+This module is designed to make a temporal fusion from multi-sources.
+<h2>NOTES</h2>
+
+<h2>TODO</h2>
+More wavelet families to be included.
+
+<h2>SEE ALSO</h2>
+
+<em>
+ <a href="r.hants.html">r.hants</a>,
+ <a href="i.lmf.html">i.lmf</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Yann Chemin, International Water Management Institute
+
+<p><i>Last changed: $Date: 2013-04-29 02:54:20 +0530 (Mon, 29 May 2013) $</i>
Added: grass-addons/grass7/imagery/i.wavelet/main.c
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/main.c (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/main.c 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,372 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.wavelet
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: A (multi-source) temporal/spectral (de)recomposition.
+ * Decompose 2 levels of a raster in temporal/spectral dimension
+ * producing Hig-Pass 1 and 2 + Low-Pass 1 and 2 (4 outputs)
+ * (-i): Recompose temporal/spectral dimensions
+ * using as input [HP2 + LP2] to recreate LP1
+ * and using as input HP1 along with the recreated LP1
+ * to produced the new data set.
+ *
+ * COPYRIGHT: (C) 2013 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>
+#include <grass/imagery.h>
+#include "wt.h"
+#include "wt_haar.h"
+#include "w_daubechies.h"
+
+#define MAXFILES 10000
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /*region+header info */
+ int nrows, ncols;
+ int row, col;
+ struct GModule *module;
+ struct Option *input, *output;
+ struct Option *ihp1, *ihp2, *ilp2;/*Recompose*/
+ struct Option *olp1, *ohp1, *ohp2, *olp2;/*Decompose*/
+ struct Option *resolution;/*wavelet resolution*/
+ struct Flag *flag1, *flag2, *flag3;
+ struct Flag *flag4, *flag5, *flag6;
+ struct History history; /*metadata */
+ struct Colors colors; /*Color rules */
+
+ /************************************/
+ char *temp; /*input raster name */
+ char *result; /*output raster name */
+ /*File Descriptors */
+ int nfiles;
+ int infd[MAXFILES];
+ int outfd[MAXFILES];
+ int outfd1[MAXFILES],outfd2[MAXFILES];
+ int outfd3[MAXFILES],outfd4[MAXFILES];
+ char **names;
+ char **ptr;
+ int i = 0, n = 0;
+ void *outrast[MAXFILES];
+ unsigned char *lp1[MAXFILES];
+ unsigned char *lp2[MAXFILES];
+ unsigned char *hp1[MAXFILES];
+ unsigned char *hp2[MAXFILES];
+
+ RASTER_MAP_TYPE in_data_type[MAXFILES]; /* 0=numbers 1=text */
+ DCELL **outlp1,**outhp1,**outlp2,**outhp2;
+ CELL val1, val2;
+ int res=0;/*wavelet sample rate option*/
+ char buffer[16];/*create file names with number extensions*/
+
+ /************************************/
+ struct Ref ref; /*group handling Decompose*/
+ struct Ref reflp2, refhp2, refhp1;/*group handling Recompose*/
+ int *fd0,*fd1,*fd2,*fd3,*fd4;/*file descriptors group rasters*/
+ /*DCELL **buf0,*buf1,*buf2,*buf3,*buf4;*//*buffers lines group rasters */
+ /************************************/
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("imagery"));
+ G_add_keyword(_("wavelet"));
+ G_add_keyword(_("fusion"));
+ module->description = _("Decompostion/Recomposition in temporal dimension using wavelets");
+
+ /* Define the different options for decomposition */
+ input = G_define_standard_option(G_OPT_I_GROUP);
+ input->key = "input_group_to_decompose";
+ input->required = NO;
+ input->guisection = "Decomposition";
+
+ olp1 = G_define_standard_option(G_OPT_R_OUTPUT);
+ olp1->key = "output_lp1_from_decomposition";
+ olp1->required = NO;
+ olp1->guisection = "Decomposition";
+ olp2 = G_define_standard_option(G_OPT_R_OUTPUT);
+ olp2->key = "output_lp2_from_decomposition";
+ olp2->required = NO;
+ olp2->guisection = "Decomposition";
+ ohp1 = G_define_standard_option(G_OPT_R_OUTPUT);
+ ohp1->key = "output_hp1_from_decomposition";
+ ohp1->required = NO;
+ ohp1->guisection = "Decomposition";
+ ohp2 = G_define_standard_option(G_OPT_R_OUTPUT);
+ ohp2->key = "output_hp2_from_decomposition";
+ ohp2->required = NO;
+ ohp2->guisection = "Decomposition";
+
+ /* Define the different options for recomposition */
+ ilp2 = G_define_standard_option(G_OPT_I_GROUP);
+ ilp2->key = _("input_lp2_group_for_recomposition");
+ ilp2->required = NO;
+ ilp2->guisection = _("Recomposition");
+ ihp1 = G_define_standard_option(G_OPT_I_GROUP);
+ ihp1->key = _("input_hp1_group_for_recomposition");
+ ihp1->required = NO;
+ ihp1->guisection = _("Recomposition");
+ ihp2 = G_define_standard_option(G_OPT_I_GROUP);
+ ihp2->key = _("input_hp2_group_for_recomposition");
+ ihp2->required = NO;
+ ihp2->guisection = _("Recomposition");
+
+ output = G_define_standard_option(G_OPT_R_INPUT);
+ output->key = _("output_from_recomposition");
+ output->required = NO;
+ output->guisection = _("Recomposition");
+
+ /* Define the different flags */
+ flag1 = G_define_flag();
+ flag1->key = 'i';
+ flag1->description = _("Recomposition (Default: Decomposition)");
+ flag1->guisection = _("Required");
+
+ flag2 = G_define_flag();
+ flag2->key = 'H';
+ flag2->description = _("Use Haar wavelets");
+ flag2->guisection = _("Wavelets");
+
+ flag3 = G_define_flag();
+ flag3->key = 'D';
+ flag3->description = _("Use Daubechies wavelets (specify resolution=4,6,8,10,12,14,16,18,20");
+ flag3->guisection = _("Wavelets");
+
+ /* Define the different values required */
+ resolution = G_define_option();
+ resolution->key = _("wavelet_sample_rate");
+ resolution->key_desc = _("integer");
+ resolution->type = TYPE_INTEGER;
+ resolution->multiple = NO;
+ resolution->description = _("4,6,8,10,12,14,16,18,20");
+ resolution->required = NO;
+ resolution->guisection = _("Optional");
+
+ nfiles = 1;
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ res=atoi(resolution->answer);
+ if((flag3->answer) && (res==4||res==6||res==8||res==10||res==12
+ ||res==14||res==16||res==18||res==20)){
+ /** Good to go with Flag3 => Daubechies **/
+ }else{ G_fatal_error(_("To use Daubechies, you need a valid resolution"));}
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ DCELL *dc;
+ if (!(flag1->answer)){
+ /* ****** DECOMPOSITION ******* */
+ if (!I_get_group_ref(input->answer, &ref))
+ G_fatal_error(_("Unable to read REF file for group <%s>"),
+ input->answer);
+ if (ref.nfiles <= 0)
+ G_fatal_error(_("Group <%s> contains no raster maps"),
+ input->answer);
+ /* Read Imagery Group */
+ fd0 = G_malloc(ref.nfiles * sizeof(int));
+ DCELL **buf0 = (DCELL **) G_malloc(ref.nfiles * sizeof(DCELL *));
+ for (n = 0; n < ref.nfiles; n++) {
+ buf0[n] = Rast_allocate_d_buf();
+ fd0[n] = Rast_open_old(ref.file[n].name, ref.file[n].mapset);
+ }
+ /* create temporal array */
+ DCELL *dc = G_malloc(ref.nfiles * sizeof(DCELL *));
+ DCELL *buf1 = G_malloc(ref.nfiles * sizeof(DCELL *));
+ DCELL *buf2 = G_malloc(ref.nfiles * sizeof(DCELL *));
+ DCELL *buf3 = G_malloc(ref.nfiles * sizeof(DCELL *));
+ DCELL *buf4 = G_malloc(ref.nfiles * sizeof(DCELL *));
+ /* Create New output raster files */
+ for (n = 0; n < ref.nfiles; n++) {
+ snprintf(buffer, sizeof(buffer), "%d", n);
+ temp=strcat(olp1->answer,".");
+ result=strcat(temp,buffer);
+ outfd1[n] = Rast_open_new(result, 1);
+ outlp1[n] = Rast_allocate_d_buf();
+ temp=strcat(ohp1->answer,".");
+ result=strcat(temp,buffer);
+ outfd2[n] = Rast_open_new(result, 1);
+ outhp1[n] = Rast_allocate_d_buf();
+ temp=strcat(olp2->answer,".");
+ result=strcat(temp,buffer);
+ outfd3[n] = Rast_open_new(result, 1);
+ outlp2[n] = Rast_allocate_d_buf();
+ temp=strcat(ohp2->answer,".");
+ result=strcat(temp,buffer);
+ outfd4[n] = Rast_open_new(result, 1);
+ outhp2[n] = Rast_allocate_d_buf();
+ }
+ /* read row */
+ for (row = 0; row < nrows; row++) {
+ for (n = 0; n < ref.nfiles; n++) {
+ /* Read one row of the signal input images */
+ Rast_get_d_row(fd0[n], buf0[n], row);
+ }
+ /* Process pixels */
+ /* #pragma parallel default(shared) private(col)*/
+ for (col = 0; col < ncols; col++) {
+ /* Extract temporal array */
+ for (n = 0; n < ref.nfiles; n++) {
+ dc[n] = buf0[n][col];
+ }
+ if (flag2->answer) {
+ dwt_haar_l2(dc, ref.nfiles, buf1, buf2, buf3, buf4);
+ }
+ else /*if (flag3->answer) which is daubechies only so far*/ {
+ dwt_l2(dc, ref.nfiles, buf1, buf2, buf3, buf4, d[res-4], d[res-3], res);
+ }
+ for (n = 0; n < ref.nfiles; n++) {
+ ((DCELL **) outlp1)[n][col] = buf1[n];
+ ((DCELL **) outhp1)[n][col] = buf2[n];
+ ((DCELL **) outlp2)[n][col] = buf3[n];
+ ((DCELL **) outhp2)[n][col] = buf4[n];
+ }
+ }
+ for (n = 0; n < ref.nfiles; n++) {
+ Rast_put_d_row(fd1[n], outlp1[n]);
+ Rast_put_d_row(fd2[n], outhp1[n]);
+ Rast_put_d_row(fd3[n], outlp2[n]);
+ Rast_put_d_row(fd4[n], outhp2[n]);
+ }
+ }
+ for (n = 0; n < ref.nfiles; n++) {
+ G_free(outlp1[n]);
+ G_free(outlp2[n]);
+ G_free(outhp1[n]);
+ G_free(outhp2[n]);
+ Rast_close(fd1[n]);
+ Rast_close(fd2[n]);
+ Rast_close(fd3[n]);
+ Rast_close(fd4[n]);
+ }
+ G_free(dc);
+ G_free(buf1);
+ G_free(buf2);
+ G_free(buf3);
+ G_free(buf4);
+ } else {
+ /* ****** RECOMPOSITION ******* */
+ if (!I_get_group_ref(ilp2->answer, &reflp2))
+ G_fatal_error(_("Unable to read REF file for group <%s>"),
+ ilp2->answer);
+ if (reflp2.nfiles <= 0)
+ G_fatal_error(_("Group <%s> contains no raster maps"),
+ ilp2->answer);
+ if (!I_get_group_ref(ihp2->answer, &refhp2))
+ G_fatal_error(_("Unable to read REF file for group <%s>"),
+ ihp2->answer);
+ if (refhp2.nfiles <= 0)
+ G_fatal_error(_("Group <%s> contains no raster maps"),
+ ihp2->answer);
+ if (!I_get_group_ref(ihp1->answer, &refhp1))
+ G_fatal_error(_("Unable to read REF file for group <%s>"),
+ ihp1->answer);
+ if (refhp1.nfiles <= 0)
+ G_fatal_error(_("Group <%s> contains no raster maps"),
+ ihp1->answer);
+
+ /* Read LP2 Imagery Group */
+ fd0 = G_malloc(reflp2.nfiles * sizeof(int));
+ DCELL **ibuf0 = (DCELL **) G_malloc(reflp2.nfiles * sizeof(DCELL *));
+ for (n = 0; n < reflp2.nfiles; n++) {
+ ibuf0[n] = Rast_allocate_d_buf();
+ fd0[n] = Rast_open_old(reflp2.file[n].name, reflp2.file[n].mapset);
+ }
+ /* Read HP2 Imagery Group */
+ fd1 = G_malloc(refhp2.nfiles * sizeof(int));
+ DCELL **ibuf1 = (DCELL **) G_malloc(refhp2.nfiles * sizeof(DCELL *));
+ for (n = 0; n < refhp2.nfiles; n++) {
+ ibuf1[n] = Rast_allocate_d_buf();
+ fd1[n] = Rast_open_old(refhp2.file[n].name, refhp2.file[n].mapset);
+ }
+ /* Read HP1 Imagery Group */
+ fd2 = G_malloc(refhp1.nfiles * sizeof(int));
+ DCELL **ibuf2 = (DCELL **) G_malloc(refhp1.nfiles * sizeof(DCELL *));
+ for (n = 0; n < refhp1.nfiles; n++) {
+ ibuf2[n] = Rast_allocate_d_buf();
+ fd2[n] = Rast_open_old(refhp1.file[n].name, refhp1.file[n].mapset);
+ }
+ /* create temporal array */
+ DCELL *rc = G_malloc(refhp1.nfiles * sizeof(DCELL *));
+ DCELL *buf0 = G_malloc(reflp2.nfiles * sizeof(DCELL *));
+ DCELL *buf1 = G_malloc(refhp2.nfiles * sizeof(DCELL *));
+ DCELL *buf2 = G_malloc(refhp1.nfiles * sizeof(DCELL *));
+ DCELL *buf3 = G_malloc(refhp1.nfiles * sizeof(DCELL *));/*LP1 in functions*/
+ /* Create New output raster files */
+ DCELL **outbuf = (DCELL **) G_malloc(refhp1.nfiles * sizeof(DCELL *));
+ for (n = 0; n < refhp1.nfiles; n++) {
+ snprintf(buffer, sizeof(buffer), "%d", n);
+ temp=strcat(output->answer,".");
+ result=strcat(temp,buffer);
+ outfd[n] = Rast_open_new(result, 1);
+ outrast[n] = Rast_allocate_d_buf();
+ }
+ /* read row */
+ for (row = 0; row < nrows; row++) {
+ /* Read one row of the input images */
+ for (n = 0; n < reflp2.nfiles; n++)
+ Rast_get_d_row(fd0[n], ibuf0[n], row);
+ for (n = 0; n < refhp2.nfiles; n++)
+ Rast_get_d_row(fd1[n], ibuf1[n], row);
+ for (n = 0; n < refhp1.nfiles; n++)
+ Rast_get_d_row(fd2[n], ibuf2[n], row);
+ /* Process pixels */
+ /* #pragma parallel default(shared) private(col)*/
+ for (col = 0; col < ncols; col++) {
+ /* Extract temporal array */
+ for (n = 0; n < reflp2.nfiles; n++)
+ buf0[n] = ibuf0[n][col];
+ for (n = 0; n < refhp2.nfiles; n++)
+ buf1[n] = ibuf1[n][col];
+ for (n = 0; n < refhp1.nfiles; n++)
+ buf2[n] = ibuf2[n][col];
+ if (flag2->answer) {
+ idwt_haar_l2(buf3, buf2, buf0, buf1, refhp1.nfiles, rc);
+ }
+ else /*if (flag3->answer) which is daubechies only so far*/ {
+ idwt_l2(buf3, buf2, buf0, buf1, refhp1.nfiles, rc, d[res-4], d[res-3], res);
+ }
+ for (n = 0; n < refhp1.nfiles; n++)
+ ((DCELL **) outbuf)[n][col] = rc[n];
+ }
+ for (n = 0; n < refhp1.nfiles; n++)
+ Rast_put_d_row(outfd[n], outbuf[n]);
+ }
+ for (n = 0; n < ref.nfiles; n++) {
+ G_free(outlp1[n]);
+ G_free(outlp2[n]);
+ G_free(outhp1[n]);
+ G_free(outhp2[n]);
+ Rast_close(fd1[n]);
+ Rast_close(fd2[n]);
+ Rast_close(fd3[n]);
+ Rast_close(fd4[n]);
+ }
+ G_free(buf0);
+ G_free(buf1);
+ G_free(buf2);
+ G_free(buf3);
+ }
+ /* Color table from 0.0 to 1.0 */
+ Rast_init_colors(&colors);
+ val1 = 0;
+ val2 = 1;
+ Rast_add_c_color_rule(&val1, 0, 0, 0, &val2, 255, 255, 255, &colors);
+ /* Metadata */
+ Rast_short_history(result, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result, &history);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/imagery/i.wavelet/w_daubechies.h
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/w_daubechies.h (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/w_daubechies.h 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,156 @@
+#include<stdio.h>
+
+/** Define the Daubechies discrete wavelets **/
+/** rows: h_4,g_4,h_6,g_6,h_8,g_8,h_10,g_10,h_12,g_12,h_14,g_14,h_16,g_16,h_20,g_20 **/
+static double d[18][20]={
+/*static double h_4[4] =*/
+{ 0.48296291314453414337487159986, 0.83651630373780790557529378092,
+ 0.22414386804201338102597276224, -0.12940952255126038117444941881,
+ 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 },
+/*static double g_4[4] =*/
+{ -0.12940952255126038117444941881, -0.22414386804201338102597276224,
+ 0.83651630373780790557529378092, -0.48296291314453414337487159986,
+ 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 },
+
+/*static double h_6[6] = */
+{ 0.33267055295008261599851158914, 0.80689150931109257649449360409,
+ 0.45987750211849157009515194215, -0.13501102001025458869638990670,
+ -0.08544127388202666169281916918, 0.03522629188570953660274066472,
+ 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 },
+/*static double g_6[6] = */
+{ 0.03522629188570953660274066472, 0.08544127388202666169281916918,
+ -0.13501102001025458869638990670, -0.45987750211849157009515194215,
+ 0.80689150931109257649449360409, -0.33267055295008261599851158914,
+ 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 },
+
+/*static double h_8[8] = */
+{ 0.23037781330889650086329118304, 0.71484657055291564708992195527,
+ 0.63088076792985890788171633830, -0.02798376941685985421141374718,
+ -0.18703481171909308407957067279, 0.03084138183556076362721936253,
+ 0.03288301166688519973540751355, -0.01059740178506903210488320852,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 },
+/*static double g_8[8] = */
+{ -0.01059740178506903210488320852, -0.03288301166688519973540751355,
+ 0.03084138183556076362721936253, 0.18703481171909308407957067279,
+ -0.02798376941685985421141374718, -0.63088076792985890788171633830,
+ 0.71484657055291564708992195527, -0.23037781330889650086329118304,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 },
+
+/*static double h_10[10] = */
+{ 0.16010239797419291448072374802, 0.60382926979718967054011930653,
+ 0.72430852843777292772807124410, 0.13842814590132073150539714634,
+ -0.24229488706638203186257137947, -0.03224486958463837464847975506,
+ 0.07757149384004571352313048939, -0.00624149021279827427419051911,
+ -0.01258075199908199946850973993, 0.00333572528547377127799818342,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 },
+/*static double g_10[10] = */
+{ 0.00333572528547377127799818342, 0.01258075199908199946850973993,
+ -0.00624149021279827427419051911, -0.07757149384004571352313048939,
+ -0.03224486958463837464847975506, 0.24229488706638203186257137947,
+ 0.13842814590132073150539714634, -0.72430852843777292772807124410,
+ 0.60382926979718967054011930653, -0.16010239797419291448072374802,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 },
+
+/*static double h_12[12] = */
+{ 0.11154074335010946362132391724, 0.49462389039845308567720417688,
+ 0.75113390802109535067893449844, 0.31525035170919762908598965481,
+ -0.22626469396543982007631450066, -0.12976686756726193556228960588,
+ 0.09750160558732304910234355254, 0.02752286553030572862554083950,
+ -0.03158203931748602956507908070, 0.00055384220116149613925191840,
+ 0.00477725751094551063963597525, -0.00107730108530847956485262161,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 },
+/*static double g_12[12] = */
+{ -0.00107730108530847956485262161, -0.00477725751094551063963597525,
+ 0.00055384220116149613925191840, 0.03158203931748602956507908070,
+ 0.02752286553030572862554083950, -0.09750160558732304910234355254,
+ -0.12976686756726193556228960588, 0.22626469396543982007631450066,
+ 0.31525035170919762908598965481, -0.75113390802109535067893449844,
+ 0.49462389039845308567720417688, -0.11154074335010946362132391724,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 },
+
+/*static double h_14[14] = */
+{ 0.07785205408500917901996352196, 0.39653931948191730653900039094,
+ 0.72913209084623511991694307034, 0.46978228740519312247159116097,
+ -0.14390600392856497540506836221, -0.22403618499387498263814042023,
+ 0.07130921926683026475087657050, 0.08061260915108307191292248036,
+ -0.03802993693501441357959206160, -0.01657454163066688065410767489,
+ 0.01255099855609984061298988603, 0.00042957797292136652113212912,
+ -0.00180164070404749091526826291, 0.00035371379997452024844629584,
+ 0.0,0.0,0.0,0.0,0.0,0.0 },
+/*static double g_14[14] = */
+{ 0.00035371379997452024844629584, 0.00180164070404749091526826291,
+ 0.00042957797292136652113212912, -0.01255099855609984061298988603,
+ -0.01657454163066688065410767489, 0.03802993693501441357959206160,
+ 0.08061260915108307191292248036, -0.07130921926683026475087657050,
+ -0.22403618499387498263814042023, 0.14390600392856497540506836221,
+ 0.46978228740519312247159116097, -0.72913209084623511991694307034,
+ 0.39653931948191730653900039094, -0.07785205408500917901996352196,
+ 0.0,0.0,0.0,0.0,0.0,0.0 },
+
+/*static double h_16[16] = */
+{ 0.05441584224310400995500940520, 0.31287159091429997065916237551,
+ 0.67563073629728980680780076705, 0.58535468365420671277126552005,
+ -0.01582910525634930566738054788, -0.28401554296154692651620313237,
+ 0.00047248457391328277036059001, 0.12874742662047845885702928751,
+ -0.01736930100180754616961614887, -0.04408825393079475150676372324,
+ 0.01398102791739828164872293057, 0.00874609404740577671638274325,
+ -0.00487035299345157431042218156, -0.00039174037337694704629808036,
+ 0.00067544940645056936636954757, -0.00011747678412476953373062823,
+ 0.0,0.0,0.0,0.0 },
+/*static double g_16[16] = */
+{ -0.00011747678412476953373062823, -0.00067544940645056936636954757,
+ -0.00039174037337694704629808036, 0.00487035299345157431042218156,
+ 0.00874609404740577671638274325, -0.01398102791739828164872293057,
+ -0.04408825393079475150676372324, 0.01736930100180754616961614887,
+ 0.12874742662047845885702928751, -0.00047248457391328277036059001,
+ -0.28401554296154692651620313237, 0.01582910525634930566738054788,
+ 0.58535468365420671277126552005, -0.67563073629728980680780076705,
+ 0.31287159091429997065916237551, -0.05441584224310400995500940520,
+ 0.0,0.0,0.0,0.0 },
+
+/*static double h_18[18] = */
+{ 0.03807794736387834658869765888, 0.24383467461259035373204158165,
+ 0.60482312369011111190307686743, 0.65728807805130053807821263905,
+ 0.13319738582500757619095494590, -0.29327378327917490880640319524,
+ -0.09684078322297646051350813354, 0.14854074933810638013507271751,
+ 0.03072568147933337921231740072, -0.06763282906132997367564227483,
+ 0.00025094711483145195758718975, 0.02236166212367909720537378270,
+ -0.00472320475775139727792570785, -0.00428150368246342983449679500,
+ 0.00184764688305622647661912949, 0.00023038576352319596720521639,
+ -0.00025196318894271013697498868, 0.00003934732031627159948068988,
+ 0.0,0.0 },
+/*static double g_18[18] = */
+{ 0.00003934732031627159948068988, 0.00025196318894271013697498868,
+ 0.00023038576352319596720521639, -0.00184764688305622647661912949,
+ -0.00428150368246342983449679500, 0.00472320475775139727792570785,
+ 0.02236166212367909720537378270, -0.00025094711483145195758718975,
+ -0.06763282906132997367564227483, -0.03072568147933337921231740072,
+ 0.14854074933810638013507271751, 0.09684078322297646051350813354,
+ -0.29327378327917490880640319524, -0.13319738582500757619095494590,
+ 0.65728807805130053807821263905, -0.60482312369011111190307686743,
+ 0.24383467461259035373204158165, -0.03807794736387834658869765888,
+ 0.0,0.0 },
+
+/*static double h_20[20] = */
+{ 0.02667005790055555358661744877, 0.18817680007769148902089297368,
+ 0.52720118893172558648174482796, 0.68845903945360356574187178255,
+ 0.28117234366057746074872699845, -0.24984642432731537941610189792,
+ -0.19594627437737704350429925432, 0.12736934033579326008267723320,
+ 0.09305736460357235116035228984, -0.07139414716639708714533609308,
+ -0.02945753682187581285828323760, 0.03321267405934100173976365318,
+ 0.00360655356695616965542329142, -0.01073317548333057504431811411,
+ 0.00139535174705290116578931845, 0.00199240529518505611715874224,
+ -0.00068585669495971162656137098, -0.00011646685512928545095148097,
+ 0.00009358867032006959133405013, -0.00001326420289452124481243668 },
+/*static double g_20[20] = */
+{ -0.00001326420289452124481243668, -0.00009358867032006959133405013,
+ -0.00011646685512928545095148097, 0.00068585669495971162656137098,
+ 0.00199240529518505611715874224, -0.00139535174705290116578931845,
+ -0.01073317548333057504431811411, -0.00360655356695616965542329142,
+ 0.03321267405934100173976365318, 0.02945753682187581285828323760,
+ -0.07139414716639708714533609308, -0.09305736460357235116035228984,
+ 0.12736934033579326008267723320, 0.19594627437737704350429925432,
+ -0.24984642432731537941610189792, -0.28117234366057746074872699845,
+ 0.68845903945360356574187178255, -0.52720118893172558648174482796,
+ 0.18817680007769148902089297368, -0.02667005790055555358661744877 }};
+
Added: grass-addons/grass7/imagery/i.wavelet/wt.c
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/wt.c (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/wt.c 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,76 @@
+#include<stdio.h>
+#include<math.h>
+int dwt_l2(double *signal, int length, double *LP1, double *HP1, double *HP2, double *LP2, double *h, double *g, int l){
+ int n,i,j,k;
+ double Hilbert_lp, Hilbert_hp;
+ for (n=0;n<2;n++){
+ printf("Decomposition level %d\tData length = %d\t", n+1, length);
+ printf("& Scaling length = %f\n", length/pow(2,n+1));
+ for (i=0;i<length/pow(2,n+1);i++){
+ if (n==0){
+ for (j=0;j<l;j++){
+ k=i*2+j;
+ while (k>=length/pow(2,n)){
+ k -= length/pow(2,n);
+ }
+ Hilbert_lp = signal[k] * h[j];
+ Hilbert_hp = signal[k] * g[j];
+ LP1[i] += Hilbert_lp;
+ HP1[i] += Hilbert_hp;
+ }
+ }
+ if (n==1){
+ for (j=0;j<l;j++){
+ k=i*2+j;
+ while (k>= length/pow(2,n)){
+ k -= length/pow(2,n);
+ }
+ Hilbert_lp = LP1[k] * h[j];
+ Hilbert_hp = LP1[k] * g[j];
+ LP2[i] += Hilbert_lp;
+ HP2[i] += Hilbert_hp;
+ }
+ }
+ }
+ }
+ return (1);
+}
+
+
+
+
+/** REVERSE MODE FUNCTION **/
+int idwt_l2(double *LP1, double *HP1, double *LP2, double *HP2, int length, double *out, double *h, double *g, int l){
+ double Hilbert_lp, Hilbert_hp;
+ int n,i,j,k;
+ for (n=0;n<2;n++){
+ printf("Recomposition level %d\tData length = %d\t", n+1, length);
+ printf("& Scaling length = %f\n", length*pow(2,n));
+ for (i=0;i<length*pow(2,n);i++){
+ if (n==0){
+ for (j=0;j<l;j++){
+ k=i*2+j;
+ while (k>=length*pow(2,n+1)){
+ k -= length*pow(2,n+1);
+ }
+ Hilbert_lp = LP2[i] * h[j];
+ Hilbert_hp = HP2[i] * g[j];
+ LP1[k] += Hilbert_lp + Hilbert_hp;
+ }
+ }
+ if (n==1){
+ for (j=0;j<l;j++){
+ k=i*2+j;
+ while (k>= length*pow(2,n+1)){
+ k -= length*pow(2,n+1);
+ }
+ Hilbert_lp = LP1[i] * h[j];
+ Hilbert_hp = HP1[i] * g[j];
+ out[k] += Hilbert_lp + Hilbert_hp;
+ }
+ }
+ }
+ }
+ return (1);
+}
+
Added: grass-addons/grass7/imagery/i.wavelet/wt.h
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/wt.h (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/wt.h 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,5 @@
+int dwt_l2(double *signal, int length, double *LP1, double *HP1, double *HP2, double *LP2, double *h, double *g, int l);
+
+/** REVERSE MODE FUNCTION **/
+int idwt_l2(double *LP1, double *HP1, double *LP2, double *HP2, int length, double *out, double *h, double *g, int l);
+
Added: grass-addons/grass7/imagery/i.wavelet/wt_haar.c
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/wt_haar.c (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/wt_haar.c 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,57 @@
+#include<stdio.h>
+#include<math.h>
+/*Public Domain*/
+
+/*Descomposition of a signal following the Haar wavelet method for 2 levels*/
+int dwt_haar_l2(double *signal, int length, double *LP1, double *HP1, double *LP2, double *HP2){
+ int n, i;
+ double summation, difference;
+ for (n=0; n<2; n++){
+ printf ("Decomposition level %d\t", n+1);
+ printf ("Data length = %d\t", length);
+ printf ("& Scaling length = %f\n", length/pow(2,n+1));
+ for (i=0;i<length/pow(2,n+1);i++){
+ if (n==0){
+ summation = signal[i * 2] + signal[i * 2 + 1];
+ difference = signal[i * 2] - signal[i * 2 + 1];
+ LP1[i] = summation;
+ HP1[i] = difference;
+ }
+ if (n==1){
+ summation = LP1[i * 2] + LP1[i * 2 + 1];
+ difference = LP1[i * 2] - LP1[i * 2 + 1];
+ LP2[i] = summation;
+ HP2[i] = difference;
+ }
+ }
+ }
+ return (1);
+}
+
+
+/*Recomposition of a signal from its three wavelet coefs from Level 1 & 2: HP1 available, LP1 made by LP2 & HP2*/
+int idwt_haar_l2(double *LP1, double *HP1, double *LP2, double *HP2, int length, double *out){
+ int n, i;
+ double summation, difference;
+ for (n=0;n<2;n++){
+ printf ("Recomposition level %d\t", (2-n));
+ printf ("Data length = %d\t", length);
+ printf ("& Scaling length = %f\n", length/pow(2,(2-n)));
+ for (i=0;i<length/pow(2,(2-n));i++){
+ if ((2-n-1)==1){
+ summation = (LP2[i] + HP2[i])/2;
+ difference = (LP2[i] - HP2[i])/2;
+ LP1[i*2] = summation;
+ LP1[i*2+1] = difference;
+ }
+ if ((2-n-1)==0){
+ summation = (LP1[i] + HP1[i])/2;
+ difference = (LP1[i] - HP1[i])/2;
+ out[i*2] = summation;
+ out[i*2+1] = difference;
+ }
+ }
+ }
+ return (1);
+}
+
Added: grass-addons/grass7/imagery/i.wavelet/wt_haar.h
===================================================================
--- grass-addons/grass7/imagery/i.wavelet/wt_haar.h (rev 0)
+++ grass-addons/grass7/imagery/i.wavelet/wt_haar.h 2013-05-07 11:09:23 UTC (rev 56148)
@@ -0,0 +1,6 @@
+/*Descomposition of a signal following the Haar wavelet method for 2 levels*/
+int dwt_haar_l2(double *signal, int length, double *LP1, double *HP1, double *LP2, double *HP2);
+
+/*Recomposition of a signal from its three wavelet coefs from Level 1 & 2: HP1 available, LP1 made by LP2 & HP2*/
+int idwt_haar_l2(double *LP1, double *HP1, double *LP2, double *HP2, int length, double *out);
+
More information about the grass-commit
mailing list