[GRASS-SVN] r46016 - grass/branches/develbranch_6/imagery/i.atcorr
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 18 02:16:18 EDT 2011
Author: mmetz
Date: 2011-04-17 23:16:18 -0700 (Sun, 17 Apr 2011)
New Revision: 46016
Added:
grass/branches/develbranch_6/imagery/i.atcorr/rbtree.cpp
grass/branches/develbranch_6/imagery/i.atcorr/rbtree.h
Modified:
grass/branches/develbranch_6/imagery/i.atcorr/6s.cpp
grass/branches/develbranch_6/imagery/i.atcorr/Altitude.cpp
grass/branches/develbranch_6/imagery/i.atcorr/Altitude.h
grass/branches/develbranch_6/imagery/i.atcorr/common.h
grass/branches/develbranch_6/imagery/i.atcorr/description.html
grass/branches/develbranch_6/imagery/i.atcorr/main.cpp
Log:
i.atcorr: cache fix, manual update (backport from trunk r46015)
Modified: grass/branches/develbranch_6/imagery/i.atcorr/6s.cpp
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/6s.cpp 2011-04-18 06:13:33 UTC (rev 46015)
+++ grass/branches/develbranch_6/imagery/i.atcorr/6s.cpp 2011-04-18 06:16:18 UTC (rev 46016)
@@ -90,10 +90,13 @@
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 @@
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 @@
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 @@
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;
Modified: grass/branches/develbranch_6/imagery/i.atcorr/Altitude.cpp
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/Altitude.cpp 2011-04-18 06:13:33 UTC (rev 46015)
+++ grass/branches/develbranch_6/imagery/i.atcorr/Altitude.cpp 2011-04-18 06:16:18 UTC (rev 46016)
@@ -211,7 +211,7 @@
palt = 0;
pps = atms.p[0];
idatmp = 0;
- taer55p = 0;
+ original_taer55p = taer55p = 0;
puw = 0;
}
else if(xpp >= 100)
@@ -219,7 +219,7 @@
/* 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 @@
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 @@
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::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;
Modified: grass/branches/develbranch_6/imagery/i.atcorr/Altitude.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/Altitude.h 2011-04-18 06:13:33 UTC (rev 46015)
+++ grass/branches/develbranch_6/imagery/i.atcorr/Altitude.h 2011-04-18 06:16:18 UTC (rev 46016)
@@ -71,6 +71,9 @@
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 @@
/* 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();
};
Modified: grass/branches/develbranch_6/imagery/i.atcorr/common.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/common.h 2011-04-18 06:13:33 UTC (rev 46015)
+++ grass/branches/develbranch_6/imagery/i.atcorr/common.h 2011-04-18 06:16:18 UTC (rev 46016)
@@ -4,7 +4,7 @@
/* Includes */
#include <stdio.h>
#include <stdlib.h>
-#include <iostream>
+#include <iostream> /* ??? */
#include <fstream>
#include <string>
#include <cmath>
Modified: grass/branches/develbranch_6/imagery/i.atcorr/description.html
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/description.html 2011-04-18 06:13:33 UTC (rev 46015)
+++ grass/branches/develbranch_6/imagery/i.atcorr/description.html 2011-04-18 06:16:18 UTC (rev 46016)
@@ -16,12 +16,11 @@
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 @@
<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 >0 for northern hemisphere and <0 for
-southern. Longitude must be >0 for eastern hemisphere and <0 for
+center. Latitude must be > 0 for northern hemisphere and < 0 for
+southern. Longitude must be > 0 for eastern hemisphere and < 0 for
western.
</blockquote>
@@ -357,7 +356,7 @@
<h3>E. Target altitude (xps), sensor platform (xpp)</h3>
Target altitude (xps, in negative [km]):
-<blockquote>xps <= 0 means the target is at the sea level.
+<blockquote>xps >= 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 @@
Convert DN (digital number = pixel values) to Radiance at top-of-atmosphere (TOA), using the
formula
<div class="code"><pre>
- Lλ = ((LMAXλ - LMINλ)/(QCALMAX-QCALMIN)) * (QCAL-QCALMIN) + LMINλ
+ Lλ = ((LMAXλ - LMINλ)/(QCALMAX-QCALMIN)) * (QCAL-QCALMIN) + LMINλ
</pre></div>
Where:
<ul>
-<li> Lλ = Spectral Radiance at the sensor's aperture in Watt/(meter squared * ster * μm), the
+<li> Lλ = Spectral Radiance at the sensor's aperture in Watt/(meter squared * ster * µm), the
apparent radiance as seen by the satellite sensor;</li>
<li> QCAL = the quantized calibrated pixel value in DN;</li>
-<li> LMINλ = the spectral radiance that is scaled to QCALMIN in watts/(meter squared * ster * μm);</li>
-<li> LMAXλ = the spectral radiance that is scaled to QCALMAX in watts/(meter squared * ster * μm);</li>
-<li> QCALMIN = the minimum quantized calibrated pixel value (corresponding to LMINλ) in DN;</li>
-<li> QCALMAX = the maximum quantized calibrated pixel value (corresponding to LMAXλ) in DN=255.</li>
+<li> LMINλ = the spectral radiance that is scaled to QCALMIN in watts/(meter squared * ster * µm);</li>
+<li> LMAXλ = the spectral radiance that is scaled to QCALMAX in watts/(meter squared * ster * µm);</li>
+<li> QCALMIN = the minimum quantized calibrated pixel value (corresponding to LMINλ) in DN;</li>
+<li> QCALMAX = the maximum quantized calibrated pixel value (corresponding to LMAXλ) in DN=255.</li>
</ul>
-LMINλ and LMAXλ are the radiances related to the minimal and maximal DN value, and are reported
+LMINλ and LMAXλ 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 @@
<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 @@
Finally, run the atmospheric correction (-r for reflectance input map; -a for date >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 ialt=elevation icnd=icnd_lsat4.txt oimg=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 @@
<br>Daniel Victoria, Anne Ghisla
<p><em>RapidEye sensors addition 11/2010:</em>
-<br>Peter Löwe, Anne Ghisla
+<br>Peter Löwe, Anne Ghisla
<p>
<i>Last changed: $Date$</i>
Modified: grass/branches/develbranch_6/imagery/i.atcorr/main.cpp
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/main.cpp 2011-04-18 06:13:33 UTC (rev 46015)
+++ grass/branches/develbranch_6/imagery/i.atcorr/main.cpp 2011-04-18 06:16:18 UTC (rev 46016)
@@ -29,7 +29,9 @@
* 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>
@@ -42,14 +44,22 @@
#include "Transform.h"
#include "6s.h"
+#include "rbtree.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 */
+ 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) */
@@ -71,10 +81,13 @@
int max;
};
+struct RBitem
+{
+ int alt; /* elevation */
+ int vis; /* visibility */
+ TransformInput ti; /* transformation parameters */
+};
-int hit = 0;
-int mis = 0;
-
/* function prototypes */
static void adjust_region(char *, const char *);
static CELL round_c(FCELL);
@@ -126,138 +139,88 @@
G_put_raster_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;
+ struct RBitem *a, *b;
-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;
- }
+ a = (struct RBitem *)ti_a;
+ b = (struct RBitem *)ti_b;
- 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;
+ /* 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 for transformation input parameters */
+class TICache
+{
+ struct RB_TREE *RBTree;
+ unsigned int tree_size;
-/* 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
-{
- enum TIMapSize
+ private:
+ struct RBitem set_alt_vis(float alt, float vis)
{
- 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 */
- };
+ struct RBitem rbitem;
- TICache tic[MAX_TICs]; /* array of TICaches */
- float visi[MAX_TICs];
- int p;
+ /* 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);
+
+ return rbitem;
+ }
-public:
- struct Position
+ public:
+ TICache()
{
- int i, j;
- Position() : i(-1), j(-1) {}
- Position(int x, int y) : i(x), j(y) {}
- bool valid() { return i != -1 && j != -1; }
- };
-
- 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();
+ RBTree = rbtree_create(compare_hv, sizeof(struct RBitem));
+ tree_size = 0;
}
+ int search(float alt, float vis, TransformInput *ti)
+ {
+ struct RBitem search_ti = set_alt_vis(alt, vis);
+ struct RBitem *found_ti;
- TransformInput get(Position pos) { return tic[pos.i].get(pos.j); }
+ found_ti = (struct RBitem *)rbtree_find(RBTree, &search_ti);
+ if (found_ti) {
+ *ti = found_ti->ti;
+ return 1;
+ }
+ else
+ return 0;
- void add(TransformInput ti, float vis, float alt) {
- tic[p].add(ti, alt);
- visi[p] = vis;
- p++;
- if(p >= MAX_TICs) p = 0;
}
-};
+ void add(TransformInput ti, float alt, float vis)
+ {
+ struct RBitem insert_ti = set_alt_vis(alt, vis);
-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;
- }
-};
+ /* add safety check here */
+ tree_size++;
+ insert_ti.ti = ti;
-typedef std::map<IntPair, TransformInput> CacheMap;
-
-
-const TransformInput& optimize_va(const FCELL& vis, const FCELL& alt)
-{
- static CacheMap timap;
- static TransformInput ti;
-
- IntPair key(vis, alt);
- CacheMap::iterator it = timap.find(key);
-
- 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.
@@ -289,8 +252,10 @@
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*)G_allocate_raster_buf(FCELL_TYPE);
if(ialt_fd >= 0) alt = (FCELL*)G_allocate_raster_buf(FCELL_TYPE);
@@ -330,15 +295,40 @@
G_set_f_null_value(&buf[col], 1);
continue;
}
- if (ialt_fd >= 0)
- alt[col] /= 1000.0f; /* converting to km from input which should be in meter */
+ 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));
+ }
+
/* 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? */
+ 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 { /* no optimizations */
pre_compute_hv(alt[col], vis[col]);
ti = compute();
@@ -350,18 +340,16 @@
{
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(); /* ... */
+ if (optimize) {
+ int in_cache = ticache.search(0, vis[col], &ti);
- ticache.add(ti, vis[col]);
- }
- }
+ 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 */
@@ -373,18 +361,16 @@
{
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(); /* ... */
+ if (optimize) {
+ int in_cache = ticache.search(alt[col], 0, &ti);
- ticache.add(ti, alt[col]);
- }
- }
+ 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 */
@@ -525,7 +511,8 @@
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;
}
@@ -573,6 +560,7 @@
int ivis_fd = -1; /* input visibility map's file descriptor */
const char *iimg_mapset, *ialt_mapset, *iviz_mapset;
struct History hist;
+ struct Cell_head orig_window;
/* Define module */
define_module();
@@ -592,6 +580,7 @@
G_fatal_error(_("Unable to open raster map <%s>"),
G_fully_qualified_name(opts.iimg->answer, iimg_mapset));
+ G_get_set_window(&orig_window);
adjust_region(opts.iimg->answer, iimg_mapset);
if(opts.ialt->answer) {
@@ -638,7 +627,15 @@
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);
@@ -657,5 +654,8 @@
Scaling is ignored and color ranges might not be correct. */
copy_colors(opts.iimg->answer, iimg_mapset, opts.oimg->answer);
+ G_set_window(&orig_window);
+ G_message(_("Atmospheric correction complete."));
+
exit(EXIT_SUCCESS);
}
Added: grass/branches/develbranch_6/imagery/i.atcorr/rbtree.cpp
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/rbtree.cpp (rev 0)
+++ grass/branches/develbranch_6/imagery/i.atcorr/rbtree.cpp 2011-04-18 06:16:18 UTC (rev 46016)
@@ -0,0 +1,552 @@
+/*!
+ * \file rbtree.c
+ *
+ * \brief binary search tree
+ *
+ * Generic balanced binary search tree (Red Black Tree) implementation
+ *
+ * (C) 2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public License
+ * (>=v2). Read the file COPYING that comes with GRASS for details.
+ *
+ * \author Original author Julienne Walker 2003, 2008
+ * GRASS implementation Markus Metz, 2009
+ */
+
+/* balanced binary search tree implementation
+ *
+ * this one is a Red Black Tree, no parent pointers, no threads
+ * The core code comes from Julienne Walker's tutorials on binary search trees
+ * original license: public domain
+ * http://eternallyconfuzzled.com/tuts/datastructures/jsw_tut_rbtree.aspx
+ * some ideas come from libavl (GPL >= 2)
+ *
+ * Red Black Trees are used to maintain a data structure with
+ * search, insertion and deletion in O(log N) time
+ */
+
+#include <cassert>
+#include <cstdlib>
+#include <cstring>
+
+extern "C"
+{
+#include <grass/gis.h>
+#include <grass/glocale.h>
+}
+#include "rbtree.h"
+
+/* internal functions */
+struct RB_NODE *rbtree_single(struct RB_NODE *, int);
+struct RB_NODE *rbtree_double(struct RB_NODE *, int);
+void *rbtree_first(struct RB_TRAV *);
+void *rbtree_next(struct RB_TRAV *);
+struct RB_NODE *rbtree_make_node(size_t, void *);
+int is_red(struct RB_NODE *);
+
+
+/* create new tree and initialize
+ * returns pointer to new tree, NULL for memory allocation error
+ */
+struct RB_TREE *rbtree_create(rb_compare_fn *compare, size_t rb_datasize)
+{
+ struct RB_TREE *tree = (struct RB_TREE *)malloc(sizeof(struct RB_TREE));
+
+ if (tree == NULL) {
+ G_warning("RB tree: Out of memory!");
+ return NULL;
+ }
+
+ assert(compare);
+
+ tree->datasize = rb_datasize;
+ tree->rb_compare = compare;
+ tree->count = 0;
+ tree->root = NULL;
+
+ return tree;
+}
+
+/* add an item to a tree
+ * non-recursive top-down insertion
+ * the algorithm does not allow duplicates and also does not warn about a duplicate
+ * returns 1 on success, 0 on failure
+ */
+int rbtree_insert(struct RB_TREE *tree, void *data)
+{
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ /* create a new root node for tree */
+ tree->root = rbtree_make_node(tree->datasize, data);
+ if (tree->root == NULL)
+ return 0;
+ }
+ else {
+ struct RB_NODE head = { 0 }; /* False tree root */
+ struct RB_NODE *g, *t; /* Grandparent & parent */
+ struct RB_NODE *p, *q; /* Iterator & parent */
+ int dir = 0, last = 0;
+
+ /* Set up helpers */
+ t = &head;
+ g = p = NULL;
+ q = t->link[1] = tree->root;
+
+ /* Search down the tree */
+ for (;;) {
+ if (q == NULL) {
+ /* Insert new node at the bottom */
+ p->link[dir] = q = rbtree_make_node(tree->datasize, data);
+ if (q == NULL)
+ return 0;
+ }
+ else if (is_red(q->link[0]) && is_red(q->link[1])) {
+ /* Color flip */
+ q->red = 1;
+ q->link[0]->red = 0;
+ q->link[1]->red = 0;
+ }
+
+ /* Fix red violation */
+ if (is_red(q) && is_red(p)) {
+ int dir2 = t->link[1] == g;
+
+ if (q == p->link[last])
+ t->link[dir2] = rbtree_single(g, !last);
+ else
+ t->link[dir2] = rbtree_double(g, !last);
+ }
+
+ last = dir;
+ dir = tree->rb_compare(q->data, data);
+
+ /* Stop if found. This check also disallows duplicates in the tree */
+ if (dir == 0)
+ break;
+
+ dir = dir < 0;
+
+ /* Move the helpers down */
+ if (g != NULL)
+ t = g;
+
+ g = p, p = q;
+ q = q->link[dir];
+ }
+
+ /* Update root */
+ tree->root = head.link[1];
+ }
+
+ /* Make root black */
+ tree->root->red = 0;
+
+ tree->count++;
+
+ return 1;
+}
+
+/* remove an item from a tree that matches given data
+ * non-recursive top-down removal
+ * returns 1 on successful removal
+ * returns 0 if data item was not found
+ */
+int rbtree_remove(struct RB_TREE *tree, const void *data)
+{
+ struct RB_NODE head = { 0 }; /* False tree root */
+ struct RB_NODE *q, *p, *g; /* Helpers */
+ struct RB_NODE *f = NULL; /* Found item */
+ int dir = 1, removed = 0;
+
+ assert(tree && data);
+
+ if (tree->root == NULL) {
+ return 0; /* empty tree, nothing to remove */
+ }
+
+ /* Set up helpers */
+ q = &head;
+ g = p = NULL;
+ q->link[1] = tree->root;
+
+ /* Search and push a red down */
+ while (q->link[dir] != NULL) {
+ int last = dir;
+
+ /* Update helpers */
+ g = p, p = q;
+ q = q->link[dir];
+ dir = tree->rb_compare(q->data, data);
+
+ /* Save found node */
+ if (dir == 0)
+ f = q;
+
+ dir = dir < 0;
+
+ /* Push the red node down */
+ if (!is_red(q) && !is_red(q->link[dir])) {
+ if (is_red(q->link[!dir]))
+ p = p->link[last] = rbtree_single(q, dir);
+ else if (!is_red(q->link[!dir])) {
+ struct RB_NODE *s = p->link[!last];
+
+ if (s != NULL) {
+ if (!is_red(s->link[!last]) && !is_red(s->link[last])) {
+ /* Color flip */
+ p->red = 0;
+ s->red = 1;
+ q->red = 1;
+ }
+ else {
+ int dir2 = g->link[1] == p;
+
+ if (is_red(s->link[last]))
+ g->link[dir2] = rbtree_double(p, last);
+ else if (is_red(s->link[!last]))
+ g->link[dir2] = rbtree_single(p, last);
+
+ /* Ensure correct coloring */
+ q->red = g->link[dir2]->red = 1;
+ g->link[dir2]->link[0]->red = 0;
+ g->link[dir2]->link[1]->red = 0;
+ }
+ }
+ }
+ }
+ }
+
+ /* Replace and remove if found */
+ if (f != NULL) {
+ free(f->data);
+ f->data = q->data;
+ p->link[p->link[1] == q] = q->link[q->link[0] == NULL];
+ free(q);
+ q = NULL;
+ tree->count--;
+ removed = 1;
+ }
+ else
+ G_debug(2, "RB tree: data not found in search tree");
+
+ /* Update root and make it black */
+ tree->root = head.link[1];
+ if (tree->root != NULL)
+ tree->root->red = 0;
+
+ return removed;
+}
+
+/* find data item in tree
+ * returns pointer to data item if found else NULL
+ */
+void *rbtree_find(struct RB_TREE *tree, const void *data)
+{
+ struct RB_NODE *curr_node = tree->root;
+ int cmp;
+
+ assert(tree && data);
+
+ while (curr_node != NULL) {
+ cmp = tree->rb_compare(curr_node->data, data);
+ if (cmp == 0)
+ return curr_node->data; /* found */
+
+ curr_node = curr_node->link[cmp < 0];
+ }
+ return NULL;
+}
+
+/* initialize tree traversal
+ * (re-)sets trav structure
+ * returns 0
+ */
+int rbtree_init_trav(struct RB_TRAV *trav, struct RB_TREE *tree)
+{
+ assert(trav && tree);
+
+ trav->tree = tree;
+ trav->curr_node = tree->root;
+ trav->first = 1;
+ trav->top = 0;
+
+ return 0;
+}
+
+/* traverse the tree in ascending order
+ * useful to get all items in the tree non-recursively
+ * struct RB_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+void *rbtree_traverse(struct RB_TRAV *trav)
+{
+ assert(trav);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_debug(1, "RB tree: empty tree");
+ else
+ G_debug(1, "RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return rbtree_next(trav);
+ else {
+ trav->first = 0;
+ return rbtree_first(trav);
+ }
+}
+
+/* find start point to traverse the tree in ascending order
+ * useful to get a selection of items in the tree
+ * magnitudes faster than traversing the whole tree
+ * may return first item that's smaller or first item that's larger
+ * struct RB_TRAV *trav needs to be initialized first
+ * returns pointer to data, NULL when finished
+ */
+void *rbtree_traverse_start(struct RB_TRAV *trav, const void *data)
+{
+ int dir = 0;
+
+ assert(trav && data);
+
+ if (trav->curr_node == NULL) {
+ if (trav->first)
+ G_warning("RB tree: empty tree");
+ else
+ G_warning("RB tree: finished traversing");
+
+ return NULL;
+ }
+
+ if (!trav->first)
+ return rbtree_next(trav);
+
+ /* else first time, get start node */
+
+ trav->first = 0;
+ trav->top = 0;
+
+ while (trav->curr_node != NULL) {
+ dir = trav->tree->rb_compare(trav->curr_node->data, data);
+ /* exact match, great! */
+ if (dir == 0)
+ return trav->curr_node->data;
+ else {
+ dir = dir < 0;
+ /* end of branch, also reached if
+ * smallest item is larger than search template or
+ * largest item is smaller than search template */
+ if (trav->curr_node->link[dir] == NULL)
+ return trav->curr_node->data;
+
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[dir];
+ }
+ }
+
+ return NULL; /* should not happen */
+}
+
+/* two functions needed to fully traverse the tree: initialize and continue
+ * useful to get all items in the tree non-recursively
+ * this one here uses a stack
+ * parent pointers or threads would also be possible
+ * but these would need to be added to RB_NODE
+ * -> more memory needed for standard operations
+ */
+
+/* start traversing the tree
+ * returns pointer to smallest data item
+ */
+void *rbtree_first(struct RB_TRAV *trav)
+{
+ /* get smallest item */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+
+ return trav->curr_node->data; /* return smallest item */
+}
+
+/* continue traversing the tree in ascending order
+ * returns pointer to data item, NULL when finished
+ */
+void *rbtree_next(struct RB_TRAV *trav)
+{
+ if (trav->curr_node->link[1] != NULL) {
+ /* something on the right side: larger item */
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[1];
+
+ /* go down, find smallest item in this branch */
+ while (trav->curr_node->link[0] != NULL) {
+ trav->up[trav->top++] = trav->curr_node;
+ trav->curr_node = trav->curr_node->link[0];
+ }
+ }
+ else {
+ /* at smallest item in this branch, go back up */
+ struct RB_NODE *last;
+
+ do {
+ if (trav->top == 0) {
+ trav->curr_node = NULL;
+ break;
+ }
+ last = trav->curr_node;
+ trav->curr_node = trav->up[--trav->top];
+ } while (last == trav->curr_node->link[1]);
+ }
+
+ if (trav->curr_node != NULL) {
+ return trav->curr_node->data;
+ }
+ else
+ return NULL; /* finished traversing */
+}
+
+/* destroy the tree */
+void rbtree_destroy(struct RB_TREE *tree)
+{
+ struct RB_NODE *it;
+ struct RB_NODE *save = tree->root;
+
+ /*
+ Rotate away the left links so that
+ we can treat this like the destruction
+ of a linked list
+ */
+ while((it = save) != NULL) {
+ if (it->link[0] == NULL) {
+ /* No left links, just kill the node and move on */
+ save = it->link[1];
+ free(it->data);
+ it->data = NULL;
+ free(it);
+ it = NULL;
+ }
+ else {
+ /* Rotate away the left link and check again */
+ save = it->link[0];
+ it->link[0] = save->link[1];
+ save->link[1] = it;
+ }
+ }
+ free(tree);
+ tree = NULL;
+}
+
+/* used for debugging: check for errors in tree structure */
+int rbtree_debug(struct RB_TREE *tree, struct RB_NODE *root)
+{
+ int lh, rh;
+
+ if (root == NULL)
+ return 1;
+ else {
+ struct RB_NODE *ln = root->link[0];
+ struct RB_NODE *rn = root->link[1];
+ int lcmp = 0, rcmp = 0;
+
+ /* Consecutive red links */
+ if (is_red(root)) {
+ if (is_red(ln) || is_red(rn)) {
+ G_warning("Red Black Tree debugging: Red violation");
+ return 0;
+ }
+ }
+
+ lh = rbtree_debug(tree, ln);
+ rh = rbtree_debug(tree, rn);
+
+ if (ln) {
+ lcmp = tree->rb_compare(ln->data, root->data);
+ }
+
+ if (rn) {
+ rcmp = tree->rb_compare(rn->data, root->data);
+ }
+
+ /* Invalid binary search tree:
+ * left node >= parent or right node <= parent */
+ if ((ln != NULL && lcmp > -1)
+ || (rn != NULL && rcmp < 1)) {
+ G_warning("Red Black Tree debugging: Binary tree violation");
+ return 0;
+ }
+
+ /* Black height mismatch */
+ if (lh != 0 && rh != 0 && lh != rh) {
+ G_warning("Red Black Tree debugging: Black violation");
+ return 0;
+ }
+
+ /* Only count black links */
+ if (lh != 0 && rh != 0)
+ return is_red(root) ? lh : lh + 1;
+ else
+ return 0;
+ }
+}
+
+/*******************************************************
+ * *
+ * internal functions for Red Black Tree maintenance *
+ * *
+ *******************************************************/
+
+/* add a new node to the tree */
+struct RB_NODE *rbtree_make_node(size_t datasize, void *data)
+{
+ struct RB_NODE *new_node = (struct RB_NODE *)malloc(sizeof(*new_node));
+
+ if (new_node == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+
+ new_node->data = malloc(datasize);
+ if (new_node->data == NULL)
+ G_fatal_error("RB Search Tree: Out of memory!");
+
+ memcpy(new_node->data, data, datasize);
+ new_node->red = 1; /* 1 is red, 0 is black */
+ new_node->link[0] = NULL;
+ new_node->link[1] = NULL;
+
+ return new_node;
+}
+
+/* check for red violation */
+int is_red(struct RB_NODE *root)
+{
+ if (root)
+ return root->red == 1;
+
+ return 0;
+}
+
+/* single rotation */
+struct RB_NODE *rbtree_single(struct RB_NODE *root, int dir)
+{
+ struct RB_NODE *newroot = root->link[!dir];
+
+ root->link[!dir] = newroot->link[dir];
+ newroot->link[dir] = root;
+
+ root->red = 1;
+ newroot->red = 0;
+
+ return newroot;
+}
+
+/* double rotation */
+struct RB_NODE *rbtree_double(struct RB_NODE *root, int dir)
+{
+ root->link[!dir] = rbtree_single(root->link[!dir], !dir);
+ return rbtree_single(root, dir);
+}
Property changes on: grass/branches/develbranch_6/imagery/i.atcorr/rbtree.cpp
___________________________________________________________________
Added: svn:mime-type
+ text/x-c++src
Added: svn:eol-style
+ native
Added: grass/branches/develbranch_6/imagery/i.atcorr/rbtree.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.atcorr/rbtree.h (rev 0)
+++ grass/branches/develbranch_6/imagery/i.atcorr/rbtree.h 2011-04-18 06:16:18 UTC (rev 46016)
@@ -0,0 +1,112 @@
+/*************************************************************
+ * USAGE *
+ *************************************************************
+ *
+ * NOTE: duplicates are not supported
+ *
+ * custom compare function
+ * extern int my_compare_fn(const void *, const void *);
+ * int my_compare_fn(const void *a, const void *b) {
+ * if ((mydatastruct *) a < (mydatastruct *) b)
+ * return -1;
+ * else if ((mydatastruct *) a > (mydatastruct *) b)
+ * return 1;
+ * else if ((mydatastruct *) a == (mydatastruct *) b)
+ * return 0;
+ * }
+ *
+ * create and initialize tree:
+ * struct RB_TREE *mytree = rbtree_create(my_compare_fn, item_size);
+ *
+ * insert items to tree:
+ * struct mydatastruct data = <some data>;
+ * if (rbtree_insert(mytree, &data) == 0)
+ * G_warning("could not insert data");
+ *
+ * find item in tree:
+ * struct mydatastruct data = <some data>;
+ * if (rbtree_find(mytree, &data) == 0)
+ * G_message("data not found");
+ *
+ * delete item from tree:
+ * struct mydatastruct data = <some data>;
+ * if (rbtree_remove(mytree, &data) == 0)
+ * G_warning("could not find data in tree");
+ *
+ * traverse tree (get all items in tree in ascending order):
+ * struct RB_TRAV trav;
+ * rbtree_init_trav(&trav, tree);
+ * while ((data = rbtree_traverse(&trav)) != NULL) {
+ * if (my_compare_fn(data, threshold_data) == 0) break;
+ * <do something with data>;
+ * }
+ *
+ * get a selection of items: all data > data1 and < data2
+ * start in tree where data is last smaller or first larger compared to data1
+ * struct RB_TRAV trav;
+ * rbtree_init_trav(&trav, tree);
+ * data = rbtree_traverse_start(&trav, &data1);
+ * <do something with data>;
+ * while ((data = rbtree_traverse(&trav)) != NULL) {
+ * if (data > data2) break;
+ * <do something with data>;
+ * }
+ *
+ * destroy tree:
+ * rbtree_destroy(mytree);
+ *
+ * debug the whole tree with
+ * rbtree_debug(mytree, mytree->root);
+ *
+ *************************************************************/
+
+#include <stddef.h>
+
+/* maximum RB Tree height */
+#define RBTREE_MAX_HEIGHT 64 /* should be more than enough */
+
+/* routine to compare data items
+ * return -1 if rb_a < rb_b
+ * return 0 if rb_a == rb_b
+ * return 1 if rb_a > rb_b
+ */
+typedef int rb_compare_fn(const void *rb_a, const void *rb_b);
+
+struct RB_NODE
+{
+ unsigned char red; /* 0 = black, 1 = red */
+ void *data; /* any kind of data */
+ struct RB_NODE *link[2]; /* link to children: link[0] for smaller, link[1] for larger */
+};
+
+struct RB_TREE
+{
+ struct RB_NODE *root; /* root node */
+ size_t datasize; /* item size */
+ size_t count; /* number of items in tree. */
+ rb_compare_fn *rb_compare; /* function to compare data */
+};
+
+struct RB_TRAV
+{
+ struct RB_TREE *tree; /* tree being traversed */
+ struct RB_NODE *curr_node; /* current node */
+ struct RB_NODE *up[RBTREE_MAX_HEIGHT]; /* stack of parent nodes */
+ int top; /* index for stack */
+ int first; /* little helper flag */
+};
+
+/* tree functions */
+struct RB_TREE *rbtree_create(rb_compare_fn *, size_t);
+void rbtree_destroy(struct RB_TREE *);
+int rbtree_insert(struct RB_TREE *, void *);
+int rbtree_remove(struct RB_TREE *, const void *);
+void *rbtree_find(struct RB_TREE *, const void *);
+
+/* tree traversal functions */
+int rbtree_init_trav(struct RB_TRAV *, struct RB_TREE *);
+void* rbtree_traverse(struct RB_TRAV *);
+void *rbtree_traverse_start(struct RB_TRAV *, const void *);
+
+/* debug tree from given node downwards */
+int rbtree_debug(struct RB_TREE *, struct RB_NODE *);
Property changes on: grass/branches/develbranch_6/imagery/i.atcorr/rbtree.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list