[GRASS-SVN] r46528 -
grass/branches/releasebranch_6_4/imagery/i.atcorr
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jun 3 02:51:30 EDT 2011
Author: mmetz
Date: 2011-06-02 23:51:30 -0700 (Thu, 02 Jun 2011)
New Revision: 46528
Added:
grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.cpp
grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.h
Modified:
grass/branches/releasebranch_6_4/imagery/i.atcorr/6s.cpp
grass/branches/releasebranch_6_4/imagery/i.atcorr/Altitude.cpp
grass/branches/releasebranch_6_4/imagery/i.atcorr/Altitude.h
grass/branches/releasebranch_6_4/imagery/i.atcorr/description.html
grass/branches/releasebranch_6_4/imagery/i.atcorr/main.cpp
Log:
backport i.atcorr
Modified: grass/branches/releasebranch_6_4/imagery/i.atcorr/6s.cpp
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.atcorr/6s.cpp 2011-06-03 06:46:11 UTC (rev 46527)
+++ grass/branches/releasebranch_6_4/imagery/i.atcorr/6s.cpp 2011-06-03 06:51:30 UTC (rev 46528)
@@ -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/releasebranch_6_4/imagery/i.atcorr/Altitude.cpp
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.atcorr/Altitude.cpp 2011-06-03 06:46:11 UTC (rev 46527)
+++ grass/branches/releasebranch_6_4/imagery/i.atcorr/Altitude.cpp 2011-06-03 06:51:30 UTC (rev 46528)
@@ -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/releasebranch_6_4/imagery/i.atcorr/Altitude.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.atcorr/Altitude.h 2011-06-03 06:46:11 UTC (rev 46527)
+++ grass/branches/releasebranch_6_4/imagery/i.atcorr/Altitude.h 2011-06-03 06:51:30 UTC (rev 46528)
@@ -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/releasebranch_6_4/imagery/i.atcorr/description.html
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.atcorr/description.html 2011-06-03 06:46:11 UTC (rev 46527)
+++ grass/branches/releasebranch_6_4/imagery/i.atcorr/description.html 2011-06-03 06:51:30 UTC (rev 46528)
@@ -13,12 +13,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
@@ -626,11 +625,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,
@@ -653,21 +649,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>
Modified: grass/branches/releasebranch_6_4/imagery/i.atcorr/main.cpp
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.atcorr/main.cpp 2011-06-03 06:46:11 UTC (rev 46527)
+++ grass/branches/releasebranch_6_4/imagery/i.atcorr/main.cpp 2011-06-03 06:51:30 UTC (rev 46528)
@@ -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,25 @@
#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.
+
+/* uncomment to disable cache usage */
+/* #define _NO_OPTIMIZE_ */
+
/* 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,19 +84,22 @@
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);
-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 define_module (void);
-static struct Options define_options (void);
-static void read_scale (Option *, ScaleRange &);
+static void adjust_region(char *, 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);
+static void copy_colors(char *, const char *, char *);
+static void define_module(void);
+static struct Options define_options(void);
+static void read_scale(Option *, ScaleRange &);
/*
@@ -91,7 +107,7 @@
Atmospheric corrections should be done on the whole
satelite image, not just portions.
*/
-static void adjust_region (char *name, const char *mapset)
+static void adjust_region(char *name, const char *mapset)
{
struct Cell_head iimg_head; /* the input image header file */
@@ -105,7 +121,7 @@
/* Rounds a floating point cell value */
-static CELL round_c (FCELL x)
+static CELL round_c(FCELL x)
{
if (x >= 0.0)
return (CELL)(x + .5);
@@ -115,7 +131,7 @@
/* 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;
int col;
@@ -126,138 +142,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.
@@ -274,9 +240,9 @@
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,
+static void process_raster(int ifd, InputMask imask, ScaleRange iscale,
int ialt_fd, int ivis_fd, int ofd, bool oflt,
- ScaleRange oscale, bool optimize)
+ ScaleRange oscale)
{
FCELL* buf; /* buffer for the input values */
FCELL* alt = NULL; /* buffer for the elevation values */
@@ -284,13 +250,21 @@
FCELL prev_alt = -1.f;
FCELL prev_vis = -1.f;
int row, col, nrows, ncols;
+ /* switch on optimization automatically if elevation and/or visibility map is given */
+ bool optimize = (ialt_fd >= 0 || ivis_fd >= 0);
+
+#ifdef _NO_OPTIMIZE_
+ optimize = false;
+#endif
/* 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*)G_allocate_raster_buf(FCELL_TYPE);
if(ialt_fd >= 0) alt = (FCELL*)G_allocate_raster_buf(FCELL_TYPE);
@@ -323,22 +297,47 @@
/* loop over all the values in the row */
for(col = 0; col < ncols; col++)
{
- if(vis && G_is_f_null_value(&vis[col]) ||
- alt && G_is_f_null_value(&alt[col]) ||
+ if ((vis && G_is_f_null_value(&vis[col])) ||
+ (alt && G_is_f_null_value(&alt[col])) ||
G_is_f_null_value(&buf[col]))
{
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 +349,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 +370,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 */
@@ -400,7 +395,7 @@
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%%"));
+ G_warning(_("The output data will overflow. Reflectance > 100%%"));
}
/* write output */
@@ -418,7 +413,7 @@
/* Copy the colors from map named iname to the map named oname */
-static void copy_colors (char *iname, const char *imapset, char *oname)
+static void copy_colors(char *iname, const char *imapset, char *oname)
{
struct Colors colors;
@@ -428,7 +423,7 @@
/* Define our module so that Grass can print it if the user wants to know more. */
-static void define_module (void)
+static void define_module(void)
{
struct GModule *module;
@@ -458,13 +453,13 @@
/* Define the options and flags */
-static struct Options define_options (void)
+static struct Options define_options(void)
{
struct Options opts;
- opts.iimg = G_define_standard_option (G_OPT_R_INPUT);
- opts.iimg->key = "iimg";
-
+ opts.iimg = G_define_standard_option(G_OPT_R_INPUT);
+ opts.iimg->key = "iimg";
+
opts.iscl = G_define_option();
opts.iscl->key = "iscl";
opts.iscl->type = TYPE_INTEGER;
@@ -526,13 +521,14 @@
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;
@@ -545,7 +541,7 @@
if(range.min==range.max)
{
- G_warning (_("Scale range length should be > 0; Using default values: [0,255]"));
+ G_warning(_("Scale range length should be > 0; Using default values: [0,255]"));
range.min = 0;
range.max = 255;
@@ -574,6 +570,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();
@@ -590,16 +587,17 @@
if ( (iimg_mapset = G_find_cell2 ( opts.iimg->answer, "") ) == NULL )
G_fatal_error ( _("Raster map <%s> not found"), opts.iimg->answer);
if((iimg_fd = G_open_cell_old(opts.iimg->answer, iimg_mapset)) < 0)
- G_fatal_error (_("Unable to open raster map <%s>"),
+ 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) {
if ( (ialt_mapset = G_find_cell2 ( opts.ialt->answer, "") ) == NULL )
G_fatal_error ( _("Raster map <%s> not found"), opts.ialt->answer);
if((ialt_fd = G_open_cell_old(opts.ialt->answer, ialt_mapset)) < 0)
- G_fatal_error (_("Unable to open raster map <%s>"),
+ G_fatal_error(_("Unable to open raster map <%s>"),
G_fully_qualified_name(opts.ialt->answer, ialt_mapset));
}
@@ -607,7 +605,7 @@
if ( (iviz_mapset = G_find_cell2 ( opts.ivis->answer, "") ) == NULL )
G_fatal_error ( _("Raster map <%s> not found"), opts.ivis->answer);
if((ivis_fd = G_open_cell_old(opts.ivis->answer, iviz_mapset)) < 0)
- G_fatal_error (_("Unable to open raster map <%s>"),
+ G_fatal_error(_("Unable to open raster map <%s>"),
G_fully_qualified_name(opts.ivis->answer, iviz_mapset));
}
@@ -615,13 +613,13 @@
if(opts.oflt->answer)
{
if((oimg_fd = G_open_fp_cell_new(opts.oimg->answer)) < 0)
- G_fatal_error (_("Unable to create raster map <%s>"),
+ G_fatal_error(_("Unable to create raster map <%s>"),
opts.oimg->answer);
}
else
{
if((oimg_fd = G_open_raster_new(opts.oimg->answer, CELL_TYPE)) < 0)
- G_fatal_error (_("Unable to create raster map <%s>"),
+ G_fatal_error(_("Unable to create raster map <%s>"),
opts.oimg->answer);
}
@@ -639,9 +637,14 @@
if(opts.etmbefore->answer) imask = (InputMask)(imask | ETM_BEFORE);
if(opts.etmafter->answer) imask = (InputMask)(imask | ETM_AFTER);
+ /* switch on optimization automatically if elevation and/or visibility map is given */
+ if (opts.optimize->answer)
+ G_important_message(_("Optimization is switched on automatically, the -o flag has no effect"));
+
/* 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);
/* Close the input and output file descriptors */
@@ -658,5 +661,8 @@
Scaling is ignored and color ranges might not be correct. */
copy_colors(opts.iimg->answer, iimg_mapset, opts.oimg->answer);
- exit (EXIT_SUCCESS);
+ G_set_window(&orig_window);
+ G_message(_("Atmospheric correction complete."));
+
+ exit(EXIT_SUCCESS);
}
Added: grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.cpp
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.cpp (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.cpp 2011-06-03 06:51:30 UTC (rev 46528)
@@ -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);
+}
Added: grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.h (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.atcorr/rbtree.h 2011-06-03 06:51:30 UTC (rev 46528)
@@ -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 *);
More information about the grass-commit
mailing list