i.atcorr: cache fix, manual update

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@46015 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.atcorr/6s.cpp

@@ -90,10 +90,13 @@ int init_6S(char* icnd_name)
 	c directional reflectances                                             c
 	c**********************************************************************/
 
+    /* NOTE: wlmoy is not affected by a height and/or vis change */
     float wlmoy;
     if(iwave.iwave != -1) wlmoy = iwave.equivwl();
     else wlmoy = iwave.wl;
 
+    iwave.wlmoy = wlmoy;
+
     discom(geom, atms, aero, aerocon, alt, iwave);
     float tamoy, tamoyp, pizmoy, pizmoyp;
     if(aero.iaer != 0) specinterp(wlmoy, tamoy, tamoyp, pizmoy, pizmoyp, aerocon, alt);
@@ -112,9 +115,7 @@ void pre_compute_hv(const float height, const float vis)
     alt.set_height(height);
     alt.init(atms, aerocon);
    
-    float wlmoy;
-    if(iwave.iwave != -1) wlmoy = iwave.equivwl();
-    else wlmoy = iwave.wl;
+    float wlmoy = iwave.wlmoy;
 
     discom(geom, atms, aero, aerocon, alt, iwave);
     float tamoy, tamoyp, pizmoy, pizmoyp;
@@ -128,9 +129,7 @@ void pre_compute_v(const float vis)
     aerocon.set_visibility(vis, atms);
     alt.init(atms, aerocon);
 
-    float wlmoy;
-    if(iwave.iwave != -1) wlmoy = iwave.equivwl();
-    else wlmoy = iwave.wl;
+    float wlmoy = iwave.wlmoy;
 
     discom(geom, atms, aero, aerocon, alt, iwave);
     float tamoy, tamoyp, pizmoy, pizmoyp;
@@ -144,9 +143,7 @@ void pre_compute_h(const float height)
     alt.set_height(height);
     alt.init(atms, aerocon);
 
-    float wlmoy;
-    if(iwave.iwave != -1) wlmoy = iwave.equivwl();
-    else wlmoy = iwave.wl;
+    float wlmoy = iwave.wlmoy;
 
     discom(geom, atms, aero, aerocon, alt, iwave);
     float tamoy, tamoyp, pizmoy, pizmoyp;

imagery/i.atcorr/Altitude.cpp

@@ -211,7 +211,7 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
 	palt = 0;
 	pps = atms.p[0];
 	idatmp = 0;
-	taer55p = 0;
+	original_taer55p = taer55p = 0;
 	puw = 0;
     }
     else if(xpp >= 100)
@@ -219,7 +219,7 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
 	/* satellite case of equivalent */
 	palt = 1000;
 	pps = 0;
-	taer55p = aerocon.taer55;
+	original_taer55p = taer55p = aerocon.taer55;
 	puw = 0;
 	ftray = 1;
 	idatmp = 4;
@@ -227,9 +227,12 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
     else
     {
 	/* "real" plane case */
-	cin >> puw;
-	cin >> puo3;
-	cin.ignore(numeric_limits<int>::max(),'\n');	/* ignore comments */
+	cin >> original_puw;
+	cin >> original_puo3;
+	cin.ignore(numeric_limits < int >::max(), '\n');	/* ignore comments */
+
+	puw = original_puw;
+	puo3 = original_puo3;
 	if ( puw < 0 )
 	{
 	    presplane(atms);
@@ -252,7 +255,8 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
 
 	palt = plane_sim.zpl[33] - atms.z[0];
 	pps = plane_sim.ppl[33];
-	cin >> taer55p;
+	cin >> original_taer55p;
+	taer55p = original_taer55p;
 
 	if ((taer55p > 0) || ((aerocon.taer55 - taer55p) < 1e-03))
 	{
@@ -274,6 +278,87 @@ void Altitude::init(AtmosModel &atms, const AerosolConcentration &aerocon)
     }
 }
 
+void Altitude::update_hv(AtmosModel & atms, const AerosolConcentration & aerocon)
+{
+    xps = original_xps;
+    xpp = original_xpp;
+
+    float uwus;
+    float uo3us;
+
+    if (xps <= 0) {
+	xps = 0;
+	uwus = 1.424f;
+	uo3us = 0.344f;
+    }
+    else if (atms.idatm != 8)
+	pressure(atms, atms.uw, atms.uo3);
+    else
+	pressure(atms, uwus, uo3us);
+
+    if (xpp <= 0) {
+	/* ground measurement option */
+	palt = 0;
+	pps = atms.p[0];
+	idatmp = 0;
+	taer55p = 0;
+	puw = 0;
+    }
+    else if (xpp >= 100) {
+	/* satellite case of equivalent */
+	palt = 1000;
+	pps = 0;
+	taer55p = aerocon.taer55;
+	puw = 0;
+	ftray = 1;
+	idatmp = 4;
+    }
+    else {
+	/* "real" plane case */
+
+	puw = original_puw;
+	puo3 = original_puo3;
+
+	if (puw < 0) {
+	    presplane(atms);
+	    idatmp = 2;
+
+	    if (atms.idatm == 8) {
+		puwus = puw;
+		puo3us = puo3;
+		puw *= atms.uw / uwus;
+		puo3 *= atms.uo3 / uo3us;
+		idatmp = 8;
+	    }
+	}
+	else {
+	    presplane(atms);
+	    idatmp = 8;
+	}
+
+	palt = plane_sim.zpl[33] - atms.z[0];
+	pps = plane_sim.ppl[33];
+	taer55p = original_taer55p;
+
+	if ((taer55p > 0) || ((aerocon.taer55 - taer55p) < 1e-03)) {
+	    /* a scale heigh of 2km is assumed in case no value is given for taer55p */
+	    taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 2)));
+	}
+	else {
+	    /* compute effective scale heigh */
+	    double sham = exp(-palt / 4);
+	    double sha = 1 - (taer55p / aerocon.taer55);
+
+	    if (sha >= sham)
+		taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / 4)));
+	    else {
+		sha = -palt / log(sha);
+		taer55p = (float)(aerocon.taer55 * (1 - exp(-palt / sha)));
+	    }
+	}
+    }
+}
+
 void Altitude::parse()
 {
     cin >> original_xps;

imagery/i.atcorr/Altitude.h

@@ -71,6 +71,9 @@ private:
      and used in subsequent calls to init to set xps and xpp */
     float original_xps;
     float original_xpp;
+    float original_taer55p;
+    float original_puw;
+    float original_puo3;
 
 	void pressure(AtmosModel& atms, float& uw, float& uo3);
 
@@ -84,8 +87,10 @@ public:
 
     /* Set the height to be used the next time init is called */
     void set_height(const float height) { original_xps = height; }
-    /* call init when ever xps changes */
+    /* call init only once: init parses input file */
     void init(AtmosModel& atms, const AerosolConcentration &aerocon);
+    /* call update_hv whenever xps changes */
+    void update_hv(AtmosModel& atms, const AerosolConcentration &aerocon);
     
 	static Altitude Parse();
 };

imagery/i.atcorr/Makefile

@@ -4,7 +4,7 @@ PGM = i.atcorr
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 
-LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(BTREE2LIB) 
 DEPENDENCIES = $(RASTERDEP) $(GISDEP)
 
 LINK = $(CXX)
@@ -12,3 +12,4 @@ LINK = $(CXX)
 ifneq ($(CXX),)
 default: cmd
 endif
+

imagery/i.atcorr/common.h

@@ -4,7 +4,7 @@
 /* Includes */
 #include <stdio.h>
 #include <stdlib.h>
-#include <iostream>
+#include <iostream> /* ??? */
 #include <fstream>
 #include <string>
 #include <cmath>

imagery/i.atcorr/i.atcorr.html

@@ -16,12 +16,11 @@ to cover the input raster map before the atmospheric correction is
 performed. The previous settings are restored afterwards.
 
 <p>
-Because using a continuous <b>elevation</b> or <b>visibility</b>
+Because using a <b>elevation</b> and/or <b>visibility</b>
 raster map makes execution time much longer, it is advised to use
-categorized raster maps instead, in conjuction with flag <b>-o</b>.
-This flag tells <em>i.atcorr</em> to try and speedup
-calculations. However, this option under certain conditions can make
-execution time longer.
+the optimization flag <b>-o</b>.
+This flag tells <em>i.atcorr</em> to try and speedup calculations. 
+However, this option will increase memory requirements.
 
 <p>
 If flag <b>-r</b> is used, the input raster data are treated as
@@ -164,8 +163,8 @@ the ascendant node at equator</td>
 <blockquote>
 * <em>NOTE</em>: for HRV, TM, ETM+, LISS and ASTER experiments,
 longitude and latitude are the coordinates of the scene
-center. Latitude must be &gt;0 for northern hemisphere and &lt;0 for
-southern. Longitude must be &gt;0 for eastern hemisphere and &lt;0 for
+center. Latitude must be &gt; 0 for northern hemisphere and &lt; 0 for
+southern. Longitude must be &gt; 0 for eastern hemisphere and &lt; 0 for
 western.
 </blockquote>
 
@@ -357,7 +356,7 @@ and 'aer' for aerosol).
 <h3>E. Target altitude (xps), sensor platform (xpp)</h3>
 
 Target altitude (xps, in negative [km]):
-<blockquote>xps &lt;= 0 means the target is at the sea level.
+<blockquote>xps &gt;= 0 means the target is at the sea level.
 <br>otherwise xps expresses the altitude of the target (e.g., mean elevation)
 in [km], given as negative value
 </blockquote>
@@ -583,20 +582,20 @@ If the overpass time is unknown, use the <a href="http://www-air.larc.nasa.gov/t
 Convert DN (digital number = pixel values) to Radiance at top-of-atmosphere (TOA), using the
 formula
 <div class="code"><pre>
-   L&#955; = ((LMAX&#955; - LMIN&#955;)/(QCALMAX-QCALMIN)) * (QCAL-QCALMIN) + LMIN&#955;
+   L&lambda; = ((LMAX&lambda; - LMIN&lambda;)/(QCALMAX-QCALMIN)) * (QCAL-QCALMIN) + LMIN&lambda;
 </pre></div>
 Where:
 <ul>
-<li> L&#955; = Spectral Radiance at the sensor's aperture in Watt/(meter squared * ster * &#956;m), the
+<li> L&lambda; = Spectral Radiance at the sensor's aperture in Watt/(meter squared * ster * &micro;m), the
       apparent radiance as seen by the satellite sensor;</li>
 <li> QCAL = the quantized calibrated pixel value in DN;</li>
-<li> LMIN&#955; = the spectral radiance that is scaled to QCALMIN in watts/(meter squared * ster * &#956;m);</li>
-<li> LMAX&#955; = the spectral radiance that is scaled to QCALMAX in watts/(meter squared * ster * &#956;m);</li>
-<li> QCALMIN = the minimum quantized calibrated pixel value (corresponding to LMIN&#955;) in DN;</li>
-<li> QCALMAX = the maximum quantized calibrated pixel value (corresponding to LMAX&#955;) in DN=255.</li>
+<li> LMIN&lambda; = the spectral radiance that is scaled to QCALMIN in watts/(meter squared * ster * &micro;m);</li>
+<li> LMAX&lambda; = the spectral radiance that is scaled to QCALMAX in watts/(meter squared * ster * &micro;m);</li>
+<li> QCALMIN = the minimum quantized calibrated pixel value (corresponding to LMIN&lambda;) in DN;</li>
+<li> QCALMAX = the maximum quantized calibrated pixel value (corresponding to LMAX&lambda;) in DN=255.</li>
 </ul>
 
-LMIN&#955; and LMAX&#955; are the radiances related to the minimal and maximal DN value, and are reported
+LMIN&lambda; and LMAX&lambda; are the radiances related to the minimal and maximal DN value, and are reported
 in the metadata file for each image, or in the table 1. High gain or low gain is also reported
 in the metadata file of each Landsat image. The minimal DN value (QCALMIN) is 1 for Landsat ETM+
 images (see
@@ -623,11 +622,8 @@ r.mapcalc "lsat7_2002_40_rad=((241.1 - (-5.1)) / (255.0 - 1.0)) * (lsat7_2002_40
 
 
 <div class="code"><pre>
-# using an integer DEM greatly accelerates the i.atcorr computations
-r.mapcalc "elev_int = round(elevation)"
-
 # find mean elevation (target above sea level, used as initialization value in control file)
-r.univar elev_int
+r.univar elevation
 </pre></div>
 
 Create a control file 'icnd.txt' for channel 4 (NIR), based on metadata. For the overpass time,
@@ -650,21 +646,21 @@ we need to define decimal hours:<br>
 Finally, run the atmospheric correction (-r for reflectance input map; -a for date &gt;July 2000;
 -o to use cache acceleration):
 <div class="code"><pre>
-i.atcorr -r -a -o lsat7_2002_40_rad ialt=elev_int icnd=icnd_lsat4.txt oimg=lsat7_2002_40_atcorr
+i.atcorr -r -a -o lsat7_2002_40_rad elev=elevation parameters=icnd_lsat4.txt output=lsat7_2002_40_atcorr
 </pre></div>
 
 Note that the altitude value from 'icnd_lsat4.txt' file is read at the beginning
 to compute the initial transform. It is necessary to give a value which could
-the the mean value of the elevation model. For the atmospheric correction then
+be the mean value of the elevation model. For the atmospheric correction then
 the raster elevation values are used from the map.
 <p>
 Note that the process is computationally intensive.<br>
 Note also, that <em>i.atcorr</em> reports solar elevation angle above horizon rather than solar zenith angle.
 
 <h2><font color="red">REMAINING DOCUMENTATION ISSUES</font></h2>
-1. It should be explained under what circumstances the use of categorized maps
-in conjuction with flag <em>-o</em> can slow down the calculations instead of
-speeding them up.
+1. The influence and importance of the visibility value or map should be 
+explained, also how to obtain an estimate for either visibility or aerosol 
+optical depth at 550nm.
 
 <h2>SEE ALSO</h2>
 
@@ -714,7 +710,7 @@ Not Found
 <br>Daniel Victoria, Anne Ghisla
 
 <p><em>RapidEye sensors addition 11/2010:</em>
-<br>Peter Löwe, Anne Ghisla
+<br>Peter L&ouml;we, Anne Ghisla
 
 <p>
 <i>Last changed: $Date$</i>

imagery/i.atcorr/main.cpp

@@ -1,3 +1,4 @@
+
 /***************************************************************************
                  atcorr - atmospheric correction for Grass GIS
 
@@ -29,40 +30,51 @@
 
 * Addition of IRS-1C LISS, Feb 2009: Markus Neteler
 
-TODO: use dynamic allocation for TiCache 
+* input elevation/visibility map: efficient cache with dynamic memory 
+* allocation: Markus Metz, Apr 2011
+
 ***************************************************************************/
 
 #include <cstdlib>
 #include <map>
 
-extern "C" {
+extern "C"
+{
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include <grass/glocale.h>
+#include <grass/rbtree.h>
 }
 
 #include "Transform.h"
 #include "6s.h"
 
+/* TICache: create 1 meter bins for altitude in km */
+/* 10m bins are also ok */
+/* BIN_ALT = 1000 / <bin size in meters> */
+#define BIN_ALT 1000.
+/* TICache: create 1 km bins for visibility */
+#define BIN_VIS 1.
+
 /* Input options and flags */
 struct Options
 {
     /* options */
-    struct Option *iimg;    /* input satellite image */
-    struct Option *iscl;    /* input data is scaled to this range */
-    struct Option *ialt;    /* an input elevation map in km used to increase */
-                            /* atmospheric correction accuracy, including this */
-                            /* will make computations take much, much longer */
-    struct Option *ivis;    /* an input visibility map in km (same purpose and effect as ialt) */
-    struct Option *icnd;    /* the input conditions file */
-    struct Option *oimg;    /* output image name */
-    struct Option *oscl;    /* scale the output data (reflectance values) to this range */
+    struct Option *iimg;	/* input satellite image */
+    struct Option *iscl;	/* input data is scaled to this range */
+    struct Option *ialt;	/* an input elevation map in meters used to increase */
+    /* atmospheric correction accuracy, including this */
+    /* will make computations take much, much longer */
+    struct Option *ivis;	/* an input visibility map in km (same purpose and effect as ialt) */
+    struct Option *icnd;	/* the input conditions file */
+    struct Option *oimg;	/* output image name */
+    struct Option *oscl;	/* scale the output data (reflectance values) to this range */
 
     /* flags */
-    struct Flag *oflt;      /* output data as floating point and do not round */
-    struct Flag *irad;      /* treat input values as reflectance instead of radiance values */
-    struct Flag *etmafter;  /* treat input data as a satelite image of type etm+ taken after July 1, 2000 */
-    struct Flag *etmbefore; /* treat input data as a satelite image of type etm+ taken before July 1, 2000 */
+    struct Flag *oflt;		/* output data as floating point and do not round */
+    struct Flag *irad;		/* treat input values as reflectance instead of radiance values */
+    struct Flag *etmafter;	/* treat input data as a satelite image of type etm+ taken after July 1, 2000 */
+    struct Flag *etmbefore;	/* treat input data as a satelite image of type etm+ taken before July 1, 2000 */
     struct Flag *optimize;
 };
 
@@ -72,16 +84,20 @@ struct ScaleRange
     int max;
 };
 
-
-int hit = 0;
-int mis = 0;
+struct RBitem
+{
+    int alt;		/* elevation */
+    int vis;		/* visibility */
+    TransformInput ti;	/* transformation parameters */
+};
 
 /* function prototypes */
 static void adjust_region(const char *);
 static CELL round_c(FCELL);
 static void write_fp_to_cell(int, FCELL *);
-static void process_raster(int, InputMask, ScaleRange, int, int, int, bool, ScaleRange, bool);
-static void copy_colors(char *, const char *, char *);
+static void process_raster(int, InputMask, ScaleRange, int, int, int, bool,
+			   ScaleRange, bool);
+static void copy_colors(const char *, char *);
 static void define_module(void);
 static struct Options define_options(void);
 static void read_scale(Option *, ScaleRange &);
@@ -91,7 +107,7 @@ static void read_scale(Option *, ScaleRange &);
    Adjust the region to that of the input raster.
    Atmospheric corrections should be done on the whole
    satelite image, not just portions.
-*/
+ */
 static void adjust_region(const char *name)
 {
     struct Cell_head iimg_head;	/* the input image header file */
@@ -105,157 +121,108 @@ static void adjust_region(const char *name)
 /* Rounds a floating point cell value */
 static CELL round_c(FCELL x)
 {
-    if(x >= 0.0)
-	return (CELL)(x + .5);
+    if (x >= 0.0)
+	return (CELL) (x + .5);
 
-    return (CELL)(-(-x + .5));
+    return (CELL) (-(-x + .5));
 }
 
 
 /* Converts the buffer to cell and write it to disk */
-static void write_fp_to_cell(int ofd, FCELL* buf)
+static void write_fp_to_cell(int ofd, FCELL *buf)
 {
-    CELL* cbuf;
+    CELL *cbuf;
     int col;
 
-    cbuf = (CELL*)Rast_allocate_buf(CELL_TYPE);
+    cbuf = (CELL *) Rast_allocate_buf(CELL_TYPE);
 
-    for(col = 0; col < Rast_window_cols(); col++) cbuf[col] = round_c(buf[col]);
+    for (col = 0; col < Rast_window_cols(); col++)
+	cbuf[col] = round_c(buf[col]);
     Rast_put_row(ofd, cbuf, CELL_TYPE);
 }
 
-
-/* See also Cache note below */
-class TICache
+/* compare function for RB tree */
+static int compare_hv(const void *ti_a, const void *ti_b)
 {
-    enum TICacheSize
-    {
-        MAX_TIs = 4096 /* this value is a guess, increase it if in general 
-                        * more categories are used. TODO: use dynamic allocation
-                        * since 4096 is the limit on 32bit */
-    };
-    TransformInput tis[MAX_TIs];
-    float alts[MAX_TIs];
-    int p;
-
-public:
-    TICache() { p = 0; for(int i = 0; i < MAX_TIs; i++) alts[i] = -1; }
-    int search(float alt) { 
-	for(int i = 0; i < MAX_TIs; i++) 
-	    if(alt == alts[i]) 
-	    {
-		hit++;
-		return i;
-	    } 
-	mis++;
-	return -1; 
-    }
-
-    TransformInput get(int i) { return tis[i]; }
-    void add(TransformInput ti, float alt) { 
-	tis[p] = ti; 
-	alts[p] = alt; 
-	p++; 
-	if(p >= MAX_TIs) p = 0; 
+    struct RBitem *a, *b;
+
+    a = (struct RBitem *)ti_a;
+    b = (struct RBitem *)ti_b;
+
+    /* check most common case first
+     * also faster if either an altitude or a visibility map is given,
+     * but not both */
+    if (a->alt == b->alt) {
+	if (a->vis == b->vis)
+	    return 0;
+	if (a->vis > b->vis)
+	    return 1;
+	else if (a->vis < b->vis)
+	    return -1;
     }
-};
+    else if (a->alt > b->alt)
+	return 1;
+    else if (a->alt < b->alt)
+	return -1;
 
+    /* should not be reached */
+    return 0;
+}
 
-/* the transform input map, is a array of ticaches.
-   The first key is the visibility which matches to a TICache for the altitudes.
-   This code is horrible, i just spent 20min writing and 5min debugging it. */
-
-/* Cache note:
-   The DEM cases are in range 0 < DEM < 8888 for the World in case of using an 
-   integer DEM values in meters. So the cache should ideally store 8888 different
-   cases for the World-type conditions if all happen in the same image. */
-
-class TIMap
+/* Cache for transformation input parameters */
+class TICache
 {
-    enum TIMapSize
-    {
-	MAX_TICs = 4096  /* this value is a guess. It means that <size> TI's will be 
-                          * the max combinations of vis/alt pairs. TODO: use dynamic allocation
-                          * since 4096 is the limit on 32bit */
-    };
-
-    TICache tic[MAX_TICs]; /* array of TICaches */
-    float visi[MAX_TICs];
-    int p;
+    struct RB_TREE *RBTree;
+    unsigned int tree_size;
 
-public:
-    struct Position
+  private:
+    struct RBitem set_alt_vis(float alt, float vis)
     {
-	int i, j;
-	Position() : i(-1), j(-1) {}
-	Position(int x, int y) : i(x), j(y) {}
-	bool valid() { return i != -1 && j != -1; }
-    };	
+	struct RBitem rbitem;
+
+	/* although alt and vis are already rounded,
+	 * the + 0.5 is needed for fp representation errors */
+	/* alt has been converted to kilometers */
+	rbitem.alt = (int) (alt * BIN_ALT + 0.5);
+	/* vis should be kilometers */
+	rbitem.vis = (int) (vis + 0.5);
 	
-    TIMap() { p = 0; for(int i = 0; i < MAX_TICs; i++) visi[i] = -1; }
-    Position search(float vis, float alt) { 
-	for(int i = 0; i < MAX_TICs; i++)
-	    if(vis == visi[i]) {
-		Position pos;
-		pos.i = i;
-		pos.j = tic[i].search(alt);
-		return pos;
-	    } 
-	return Position();
+	return rbitem;
     }
 
-    TransformInput get(Position pos) { return tic[pos.i].get(pos.j); }
-
-    void add(TransformInput ti, float vis, float alt) {
-	tic[p].add(ti, alt);
-	visi[p] = vis;
-	p++;
-	if(p >= MAX_TICs) p = 0;
+  public:
+      TICache()
+    {
+	RBTree = rbtree_create(compare_hv, sizeof(struct RBitem));
+	tree_size = 0;
     }
-};
-
-
-struct IntPair
-{
-    FCELL x;
-    FCELL y;
-	
-    IntPair(FCELL i, FCELL j) : x(i), y(j) {}
-	
-    bool operator<(const IntPair& b) const
-	{
-	    if(x < b.x) return true;
-	    else if(x > b.x) return false;
-	    else if(y < b.y) return true;
-	    return false;
-	}	
-};
+    int search(float alt, float vis, TransformInput *ti)
+    {
+	struct RBitem search_ti = set_alt_vis(alt, vis);
+	struct RBitem *found_ti;
 
+	found_ti = (struct RBitem *)rbtree_find(RBTree, &search_ti);
+	if (found_ti) {
+	    *ti = found_ti->ti;
+	    return 1;
+	}
+	else
+	    return 0;
 
-typedef std::map<IntPair, TransformInput> CacheMap;
+    }
 
+    void add(TransformInput ti, float alt, float vis)
+    {
+	struct RBitem insert_ti = set_alt_vis(alt, vis);
 
-const TransformInput& optimize_va(const FCELL& vis, const FCELL& alt)
-{
-    static CacheMap timap;
-    static TransformInput ti;
+	/* add safety check here */
+	tree_size++;
 
-    IntPair key(vis, alt);
-    CacheMap::iterator it = timap.find(key);
+	insert_ti.ti = ti;
 
-    if(it != timap.end()) /* search found key */
-    {
-	ti = (*it).second;
+	rbtree_insert(RBTree, &insert_ti);
     }
-    else
-    {
-	pre_compute_hv(alt, vis);
-	ti = compute();
-	timap.insert(std::make_pair(key, ti));
-    }
-	
-    return ti;
-}	
+};
 
 
 /* Process the raster and do atmospheric corrections.
@@ -271,140 +238,167 @@ const TransformInput& optimize_va(const FCELL& vis, const FCELL& alt)
    ofd: output file descriptor
    oflt: if true use FCELL_TYPE for output
    oscale: output file's range (default is min = 0, max = 255)
-*/
+ */
 static void process_raster(int ifd, InputMask imask, ScaleRange iscale,
-			    int ialt_fd, int ivis_fd, int ofd, bool oflt,
-			    ScaleRange oscale, bool optimize)
+			   int ialt_fd, int ivis_fd, int ofd, bool oflt,
+			   ScaleRange oscale, bool optimize)
 {
-    FCELL* buf;         /* buffer for the input values */
-    FCELL* alt = NULL;         /* buffer for the elevation values */
-    FCELL* vis = NULL;         /* buffer for the visibility values */
-    FCELL  prev_alt = -1.f;
-    FCELL  prev_vis = -1.f;
+    FCELL *buf;			/* buffer for the input values */
+    FCELL *alt = NULL;		/* buffer for the elevation values */
+    FCELL *vis = NULL;		/* buffer for the visibility values */
+    FCELL prev_alt = -1.f;
+    FCELL prev_vis = -1.f;
     int row, col, nrows, ncols;
 
     /* do initial computation with global elevation and visibility values */
     TransformInput ti;
+
     ti = compute();
 
-    TICache ticache;    /* use this to increase computation speed when an elevation map with categories are given */
-	
+    /* use a cache to increase computation speed when an elevation map 
+     * and/or a visibility map is given */
+    TICache ticache;
+
     /* allocate memory for buffers */
-    buf = (FCELL*)Rast_allocate_buf(FCELL_TYPE);
-    if(ialt_fd >= 0) alt = (FCELL*)Rast_allocate_buf(FCELL_TYPE);
-    if(ivis_fd >= 0) vis = (FCELL*)Rast_allocate_buf(FCELL_TYPE);
+    buf = (FCELL *) Rast_allocate_buf(FCELL_TYPE);
+    if (ialt_fd >= 0)
+	alt = (FCELL *) Rast_allocate_buf(FCELL_TYPE);
+    if (ivis_fd >= 0)
+	vis = (FCELL *) Rast_allocate_buf(FCELL_TYPE);
 
     nrows = Rast_window_rows();
     ncols = Rast_window_cols();
 
-    for(row = 0; row < nrows; row++)
-    {
-	G_percent(row, nrows, 1);     /* keep the user informed of our progress */
-		
-        /* read the next row */
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 1);	/* keep the user informed of our progress */
+
+	/* read the next row */
 	Rast_get_row(ifd, buf, row, FCELL_TYPE);
 
-        /* read the next row of elevation values */
-        if(ialt_fd >= 0)
+	/* read the next row of elevation values */
+	if (ialt_fd >= 0)
 	    Rast_get_row(ialt_fd, alt, row, FCELL_TYPE);
 
-        /* read the next row of elevation values */
-        if(ivis_fd >= 0)
+	/* read the next row of elevation values */
+	if (ivis_fd >= 0)
 	    Rast_get_row(ivis_fd, vis, row, FCELL_TYPE);
 
-        /* loop over all the values in the row */
-	for(col = 0; col < ncols; col++)
-	{
-	    if((vis && Rast_is_f_null_value(&vis[col])) || 
-	       (alt && Rast_is_f_null_value(&alt[col])) || 
-	              Rast_is_f_null_value(&buf[col]))
-	    {
-	        Rast_set_f_null_value(&buf[col], 1);
-	        continue;
+	/* loop over all the values in the row */
+	for (col = 0; col < ncols; col++) {
+	    if ((vis && Rast_is_f_null_value(&vis[col])) ||
+		(alt && Rast_is_f_null_value(&alt[col])) ||
+		Rast_is_f_null_value(&buf[col])) {
+		Rast_set_f_null_value(&buf[col], 1);
+		continue;
+	    }
+	    if (ialt_fd >= 0) {
+		if (alt[col] < 0)
+		    alt[col] = 0; /* on or below sea level, all the same for 6S */
+		else
+		    alt[col] /= 1000.0f;	/* converting to km from input which should be in meter */
+
+		/* round to nearest altitude bin */
+		/* rounding result: watch out for fp representation error */
+		alt[col] = ((int) (alt[col] * BIN_ALT + 0.5)) / BIN_ALT;
+	    }
+	    if (ivis_fd >= 0) {
+		if (vis[col] < 0)
+		    vis[col] = 0; /* negative visibility is invalid, print a WARNING ? */
+
+		/* round to nearest visibility bin */
+		/* rounding result: watch out for fp representation error */
+		vis[col] = ((int) (vis[col] + 0.5));
 	    }
-	    if (ialt_fd >= 0)
-		alt[col] /= 1000.0f; /* converting to km from input which should be in meter */
-
-            /* check if both maps are active and if whether any value has changed */
-            if((ialt_fd >= 0) && (ivis_fd >= 0) && ((prev_vis != vis[col]) || (prev_alt != alt[col])))
-            {
-               	prev_alt = alt[col]; /* update new values */
-               	prev_vis = vis[col];
- 		if(optimize) ti = optimize_va(vis[col], alt[col]); /* try to optimize? */
-		else { /* no optimizations */
-		    pre_compute_hv(alt[col], vis[col]);
-		    ti = compute();
-		}	
-            }
-            else    /* only one of the maps is being used */
-            {
-                if((ivis_fd >= 0) && (prev_vis != vis[col]))
-                {
-                    prev_vis = vis[col];        /* keep track of previous visibility */
-                    
-                    if(optimize)
-                    {
-                        int p = ticache.search(vis[col]);
-                        if(p >= 0) ti = ticache.get(p);
-                        else
-                        {
-                            pre_compute_v(vis[col]);    /* re-compute transformation inputs */
-                            ti = compute();             /* ... */
-
-                            ticache.add(ti, vis[col]);                        
-                        }
-                    }
-                    else
-                    {
-                        pre_compute_v(vis[col]);    /* re-compute transformation inputs */
-                        ti = compute();             /* ... */
-                    }
-                }
-
-                if((ialt_fd >= 0) && (prev_alt != alt[col]))
-                {
-                    prev_alt = alt[col];        /* keep track of previous altitude */
-
-                    if(optimize)
-                    {
-                        int p = ticache.search(alt[col]);
-                        if(p >= 0) ti = ticache.get(p);
-                        else
-                        {
-                            pre_compute_h(alt[col]);    /* re-compute transformation inputs */
-                            ti = compute();             /* ... */
-
-                            ticache.add(ti, alt[col]);
-                        }
-                    }
-                    else
-                    {
-                        pre_compute_h(alt[col]);    /* re-compute transformation inputs */
-                        ti = compute();             /* ... */
-                    }
-                }
-            }
-	    G_debug(3, "Computed r%d (%d), c%d (%d)", row, nrows, col, ncols);
-            /* transform from iscale.[min,max] to [0,1] */
-            buf[col] = (buf[col] - iscale.min) / ((float)iscale.max - (float)iscale.min);
-            buf[col] = transform(ti, imask, buf[col]);
-            /* transform from [0,1] to oscale.[min,max] */
-            buf[col] = buf[col] * ((float)oscale.max - (float)oscale.min) + oscale.min;
 
-            if(~oflt && (buf[col] > (float)oscale.max))
+	    /* check if both maps are active and if whether any value has changed */
+	    if ((ialt_fd >= 0) && (ivis_fd >= 0) &&
+		((prev_vis != vis[col]) || (prev_alt != alt[col]))) {
+		prev_alt = alt[col];	/* update new values */
+		prev_vis = vis[col];
+		if (optimize) {
+		    int in_cache = ticache.search(alt[col], vis[col], &ti);
+
+		    if (!in_cache) {
+			pre_compute_hv(alt[col], vis[col]);	/* re-compute transformation inputs */
+			ti = compute();	/* ... */
+
+			ticache.add(ti, alt[col], vis[col]);
+		    }
+		}
+		else {
+		    pre_compute_hv(alt[col], vis[col]);	/* re-compute transformation inputs */
+		    ti = compute();	/* ... */
+		}
+	    }
+	    else {		/* only one of the maps is being used */
+
+		if ((ivis_fd >= 0) && (prev_vis != vis[col])) {
+		    prev_vis = vis[col];	/* keep track of previous visibility */
+
+		    if (optimize) {
+			int in_cache = ticache.search(0, vis[col], &ti);
+
+			if (!in_cache) {
+			    pre_compute_v(vis[col]);	/* re-compute transformation inputs */
+			    ti = compute();	/* ... */
+
+			    ticache.add(ti, 0, vis[col]);
+			}
+		    }
+		    else {
+			pre_compute_v(vis[col]);	/* re-compute transformation inputs */
+			ti = compute();	/* ... */
+		    }
+		}
+
+		if ((ialt_fd >= 0) && (prev_alt != alt[col])) {
+		    prev_alt = alt[col];	/* keep track of previous altitude */
+
+		    if (optimize) {
+			int in_cache = ticache.search(alt[col], 0, &ti);
+
+			if (!in_cache) {
+			    pre_compute_h(alt[col]);	/* re-compute transformation inputs */
+			    ti = compute();	/* ... */
+
+			    ticache.add(ti, alt[col], 0);
+			}
+		    }
+		    else {
+			pre_compute_h(alt[col]);	/* re-compute transformation inputs */
+			ti = compute();	/* ... */
+		    }
+		}
+	    }
+	    G_debug(3, "Computed r%d (%d), c%d (%d)", row, nrows, col, ncols);
+	    /* transform from iscale.[min,max] to [0,1] */
+	    buf[col] =
+		(buf[col] - iscale.min) / ((float)iscale.max -
+					   (float)iscale.min);
+	    buf[col] = transform(ti, imask, buf[col]);
+	    /* transform from [0,1] to oscale.[min,max] */
+	    buf[col] =
+		buf[col] * ((float)oscale.max - (float)oscale.min) +
+		oscale.min;
+
+	    if (~oflt && (buf[col] > (float)oscale.max))
 		G_warning(_("The output data will overflow. Reflectance > 100%%"));
 	}
 
-        /* write output */
-	if(oflt) Rast_put_row(ofd, buf, FCELL_TYPE);
-	else write_fp_to_cell(ofd, buf);
+	/* write output */
+	if (oflt)
+	    Rast_put_row(ofd, buf, FCELL_TYPE);
+	else
+	    write_fp_to_cell(ofd, buf);
     }
     G_percent(1, 1, 1);
-    
+
     /* free allocated memory */
     G_free(buf);
-    if(ialt_fd >= 0) G_free(alt);
-    if(ivis_fd >= 0) G_free(vis);
+    if (ialt_fd >= 0)
+	G_free(alt);
+    if (ivis_fd >= 0)
+	G_free(vis);
 }
 
 
@@ -425,7 +419,8 @@ static void define_module(void)
     struct GModule *module;
 
     module = G_define_module();
-    module->label = _("Performs atmospheric correction using the 6S algorithm.");
+    module->label =
+	_("Performs atmospheric correction using the 6S algorithm.");
     module->description =
 	_("6S - Second Simulation of Satellite Signal in the Solar Spectrum.");
     G_add_keyword(_("imagery"));
@@ -446,7 +441,7 @@ static void define_module(void)
        " The code is provided as is and is not to be sold. See notes on\n"
        " http://loasys.univ-lille1.fr/informatique/sixs_gb.html\n"
        " http://www.ltid.inpe.br/dsr/mauro/6s/index.html\n"
-       " and on http://www.cs.sun.ac.za/~caz/index.html\n";*/
+       " and on http://www.cs.sun.ac.za/~caz/index.html\n"; */
 }
 
 
@@ -456,41 +451,41 @@ static struct Options define_options(void)
     struct Options opts;
 
     opts.iimg = G_define_standard_option(G_OPT_R_INPUT);
-    
+
     opts.iscl = G_define_option();
-    opts.iscl->key          = "range";
-    opts.iscl->type         = TYPE_INTEGER;
-    opts.iscl->key_desc     = "min,max";
-    opts.iscl->required     = NO;
-    opts.iscl->answer       = "0,255";
-    opts.iscl->description  = _("Input range");
+    opts.iscl->key = "range";
+    opts.iscl->type = TYPE_INTEGER;
+    opts.iscl->key_desc = "min,max";
+    opts.iscl->required = NO;
+    opts.iscl->answer = "0,255";
+    opts.iscl->description = _("Input range");
     opts.iscl->guisection = _("Input");
 
     opts.ialt = G_define_standard_option(G_OPT_R_ELEV);
-    opts.ialt->required	        = NO;
-    opts.ialt->description      = _("Name of input elevation raster map (in m)");
-    opts.ialt->guisection       = _("Input");
+    opts.ialt->required = NO;
+    opts.ialt->description = _("Name of input elevation raster map (in m)");
+    opts.ialt->guisection = _("Input");
 
     opts.ivis = G_define_standard_option(G_OPT_R_INPUT);
-    opts.ivis->key		= "visibility";
-    opts.ivis->required	        = NO;
-    opts.ivis->description	= _("Name of input visibility raster map (in km)");
-    opts.ivis->guisection       = _("Input");
+    opts.ivis->key = "visibility";
+    opts.ivis->required = NO;
+    opts.ivis->description = _("Name of input visibility raster map (in km)");
+    opts.ivis->guisection = _("Input");
 
     opts.icnd = G_define_standard_option(G_OPT_F_INPUT);
-    opts.icnd->key		= "parameters";
-    opts.icnd->required	        = YES;
-    opts.icnd->description	= _("Name of input text file with 6S parameters");
+    opts.icnd->key = "parameters";
+    opts.icnd->required = YES;
+    opts.icnd->description = _("Name of input text file with 6S parameters");
 
     opts.oimg = G_define_standard_option(G_OPT_R_OUTPUT);
 
     opts.oscl = G_define_option();
-    opts.oscl->key          = "rescale";
-    opts.oscl->type         = TYPE_INTEGER;
-    opts.oscl->key_desc     = "min,max";
-    opts.oscl->answer       = "0,255";
-    opts.oscl->required     = NO;
-    opts.oscl->description  = _("Rescale output raster map");
+    opts.oscl->key = "rescale";
+    opts.oscl->type = TYPE_INTEGER;
+    opts.oscl->key_desc = "min,max";
+    opts.oscl->answer = "0,255";
+    opts.oscl->required = NO;
+    opts.oscl->description = _("Rescale output raster map");
     opts.oscl->guisection = _("Output");
 
     opts.oflt = G_define_flag();
@@ -500,111 +495,112 @@ static struct Options define_options(void)
 
     opts.irad = G_define_flag();
     opts.irad->key = 'r';
-    opts.irad->description = _("Input raster map converted to reflectance (default is radiance)");
+    opts.irad->description =
+	_("Input raster map converted to reflectance (default is radiance)");
     opts.irad->guisection = _("Input");
 
     opts.etmafter = G_define_flag();
     opts.etmafter->key = 'a';
-    opts.etmafter->description = _("Input from ETM+ image taken after July 1, 2000");
+    opts.etmafter->description =
+	_("Input from ETM+ image taken after July 1, 2000");
     opts.etmafter->guisection = _("Input");
 
     opts.etmbefore = G_define_flag();
     opts.etmbefore->key = 'b';
-    opts.etmbefore->description = _("Input from ETM+ image taken before July 1, 2000");
+    opts.etmbefore->description =
+	_("Input from ETM+ image taken before July 1, 2000");
     opts.etmbefore->guisection = _("Input");
 
     opts.optimize = G_define_flag();
     opts.optimize->key = 'o';
-    opts.optimize->description = _("Try to increase computation speed when categorized altitude or/and visibility map is used");
+    opts.optimize->description =
+	_("Try to increase computation speed when altitude and/or visibility map is used");
 
     return opts;
 }
 
 /* Read the min and max values from the iscl and oscl options */
-void read_scale(Option *scl, ScaleRange &range)
+void read_scale(Option * scl, ScaleRange & range)
 {
     /* set default values */
     range.min = 0;
     range.max = 255;
 
-    if(scl->answer)
-    {
-        sscanf(scl->answers[0], "%d", &range.min);
-        sscanf(scl->answers[1], "%d", &range.max);
+    if (scl->answer) {
+	sscanf(scl->answers[0], "%d", &range.min);
+	sscanf(scl->answers[1], "%d", &range.max);
 
-        if(range.min==range.max)
-        {
-            G_warning(_("Scale range length should be > 0; Using default values: [0,255]"));
+	if (range.min == range.max) {
+	    G_warning(_("Scale range length should be > 0; Using default values: [0,255]"));
 
-            range.min = 0;
-            range.max = 255;
-        }
+	    range.min = 0;
+	    range.max = 255;
+	}
     }
 
     /* swap values if max is smaller than min */
-    if(range.max < range.min)
-    {
-        int temp;
-        temp = range.max;
-        range.max = range.min;
-        range.min = temp;
+    if (range.max < range.min) {
+	int temp;
+
+	temp = range.max;
+	range.max = range.min;
+	range.min = temp;
     }
 }
 
 
-int main(int argc, char* argv[])
+int main(int argc, char *argv[])
 {
-    struct Options opts;        
-    struct ScaleRange iscale;   /* input file's data is scaled to this interval */
-    struct ScaleRange oscale;   /* output file's scale */
-    int iimg_fd;	        /* input image's file descriptor */
-    int oimg_fd;	        /* output image's file descriptor */
-    int ialt_fd = -1;       /* input elevation map's file descriptor */
-    int ivis_fd = -1;       /* input visibility map's file descriptor */
+    struct Options opts;
+    struct ScaleRange iscale;	/* input file's data is scaled to this interval */
+    struct ScaleRange oscale;	/* output file's scale */
+    int iimg_fd;		/* input image's file descriptor */
+    int oimg_fd;		/* output image's file descriptor */
+    int ialt_fd = -1;		/* input elevation map's file descriptor */
+    int ivis_fd = -1;		/* input visibility map's file descriptor */
     struct History hist;
-    
+    struct Cell_head orig_window;
+
     /* Define module */
     define_module();
-  
+
     /* Define the different input options */
     opts = define_options();
 
     /**** Start ****/
     G_gisinit(argv[0]);
-    if(G_parser(argc, argv) < 0)
+    if (G_parser(argc, argv) < 0)
 	exit(EXIT_FAILURE);
 
+    G_get_set_window(&orig_window);
     adjust_region(opts.iimg->answer);
 
     /* open input raster */
-    if((iimg_fd = Rast_open_old(opts.iimg->answer, "")) < 0)
-	G_fatal_error(_("Unable to open raster map <%s>"),
-		       opts.iimg->answer);
-        
-    if(opts.ialt->answer) {
-	if((ialt_fd = Rast_open_old(opts.ialt->answer, "")) < 0)
-            G_fatal_error(_("Unable to open raster map <%s>"),
-			   opts.ialt->answer);
+    if ((iimg_fd = Rast_open_old(opts.iimg->answer, "")) < 0)
+	G_fatal_error(_("Unable to open raster map <%s>"), opts.iimg->answer);
+
+    if (opts.ialt->answer) {
+	if ((ialt_fd = Rast_open_old(opts.ialt->answer, "")) < 0)
+	    G_fatal_error(_("Unable to open raster map <%s>"),
+			  opts.ialt->answer);
     }
 
-    if(opts.ivis->answer) {
-	if((ivis_fd = Rast_open_old(opts.ivis->answer, "")) < 0)
-            G_fatal_error(_("Unable to open raster map <%s>"),
-			   opts.ivis->answer);
+    if (opts.ivis->answer) {
+	if ((ivis_fd = Rast_open_old(opts.ivis->answer, "")) < 0)
+	    G_fatal_error(_("Unable to open raster map <%s>"),
+			  opts.ivis->answer);
     }
-                
+
     /* open a floating point raster or not? */
-    if(opts.oflt->answer)
-    {
-	if((oimg_fd = Rast_open_fp_new(opts.oimg->answer)) < 0)
+    if (opts.oflt->answer) {
+	if ((oimg_fd = Rast_open_fp_new(opts.oimg->answer)) < 0)
 	    G_fatal_error(_("Unable to create raster map <%s>"),
-			   opts.oimg->answer);
+			  opts.oimg->answer);
     }
-    else
-    {
-	if((oimg_fd = Rast_open_new(opts.oimg->answer, CELL_TYPE)) < 0)
+    else {
+	if ((oimg_fd = Rast_open_new(opts.oimg->answer, CELL_TYPE)) < 0)
 	    G_fatal_error(_("Unable to create raster map <%s>"),
-			   opts.oimg->answer);
+			  opts.oimg->answer);
     }
 
     /* read the scale parameters */
@@ -613,24 +609,37 @@ int main(int argc, char* argv[])
 
     /* initialize this 6s computation and parse the input conditions file */
     init_6S(opts.icnd->answer);
-	
-    InputMask imask = RADIANCE;         /* the input mask tells us what transformations if any
-					   needs to be done to make our input values, reflectance
-					   values scaled between 0 and 1 */
-    if(opts.irad->answer) imask = REFLECTANCE;
-    if(opts.etmbefore->answer) imask = (InputMask)(imask | ETM_BEFORE);
-    if(opts.etmafter->answer) imask = (InputMask)(imask | ETM_AFTER);
+
+    InputMask imask = RADIANCE;	/* the input mask tells us what transformations if any
+				   needs to be done to make our input values, reflectance
+				   values scaled between 0 and 1 */
+    if (opts.irad->answer)
+	imask = REFLECTANCE;
+    if (opts.etmbefore->answer)
+	imask = (InputMask) (imask | ETM_BEFORE);
+    if (opts.etmafter->answer)
+	imask = (InputMask) (imask | ETM_AFTER);
+
+    if ((ialt_fd >= 0 || ivis_fd >= 0) && !opts.optimize->answer) {
+	G_message(_("An elevation and/or visibility map is given, but the optimization flag is not set."));
+	G_message(_("This can take some time."));
+    }
+
+    /* switch on optimization automatically if elevation and/or visibility map is given? */
 
     /* process the input raster and produce our atmospheric corrected output raster. */
+    G_message(_("Atmospheric correction..."));
     process_raster(iimg_fd, imask, iscale, ialt_fd, ivis_fd,
-                   oimg_fd, opts.oflt->answer, oscale, opts.optimize->answer);
+		   oimg_fd, opts.oflt->answer, oscale, opts.optimize->answer);
 
 
     /* Close the input and output file descriptors */
     Rast_short_history(opts.oimg->answer, "raster", &hist);
     Rast_close(iimg_fd);
-    if(opts.ialt->answer) Rast_close(ialt_fd);
-    if(opts.ivis->answer) Rast_close(ivis_fd);
+    if (opts.ialt->answer)
+	Rast_close(ialt_fd);
+    if (opts.ivis->answer)
+	Rast_close(ivis_fd);
     Rast_close(oimg_fd);
 
     Rast_command_history(&hist);
@@ -640,5 +649,8 @@ int main(int argc, char* argv[])
        Scaling is ignored and color ranges might not be correct. */
     copy_colors(opts.iimg->answer, opts.oimg->answer);
 
+    Rast_set_window(&orig_window);
+    G_message(_("Atmospheric correction complete."));
+
     exit(EXIT_SUCCESS);
 }