[GRASS-SVN] r46015 - grass/trunk/imagery/i.atcorr
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 18 02:13:33 EDT 2011
Author: mmetz
Date: 2011-04-17 23:13:33 -0700 (Sun, 17 Apr 2011)
New Revision: 46015
Modified:
grass/trunk/imagery/i.atcorr/6s.cpp
grass/trunk/imagery/i.atcorr/Altitude.cpp
grass/trunk/imagery/i.atcorr/Altitude.h
grass/trunk/imagery/i.atcorr/Makefile
grass/trunk/imagery/i.atcorr/common.h
grass/trunk/imagery/i.atcorr/i.atcorr.html
grass/trunk/imagery/i.atcorr/main.cpp
Log:
i.atcorr: cache fix, manual update
Modified: grass/trunk/imagery/i.atcorr/6s.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/6s.cpp 2011-04-17 18:19:14 UTC (rev 46014)
+++ grass/trunk/imagery/i.atcorr/6s.cpp 2011-04-18 06:13:33 UTC (rev 46015)
@@ -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/trunk/imagery/i.atcorr/Altitude.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/Altitude.cpp 2011-04-17 18:19:14 UTC (rev 46014)
+++ grass/trunk/imagery/i.atcorr/Altitude.cpp 2011-04-18 06:13:33 UTC (rev 46015)
@@ -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/trunk/imagery/i.atcorr/Altitude.h
===================================================================
--- grass/trunk/imagery/i.atcorr/Altitude.h 2011-04-17 18:19:14 UTC (rev 46014)
+++ grass/trunk/imagery/i.atcorr/Altitude.h 2011-04-18 06:13:33 UTC (rev 46015)
@@ -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/trunk/imagery/i.atcorr/Makefile
===================================================================
--- grass/trunk/imagery/i.atcorr/Makefile 2011-04-17 18:19:14 UTC (rev 46014)
+++ grass/trunk/imagery/i.atcorr/Makefile 2011-04-18 06:13:33 UTC (rev 46015)
@@ -4,7 +4,7 @@
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 @@
ifneq ($(CXX),)
default: cmd
endif
+
Modified: grass/trunk/imagery/i.atcorr/common.h
===================================================================
--- grass/trunk/imagery/i.atcorr/common.h 2011-04-17 18:19:14 UTC (rev 46014)
+++ grass/trunk/imagery/i.atcorr/common.h 2011-04-18 06:13:33 UTC (rev 46015)
@@ -4,7 +4,7 @@
/* Includes */
#include <stdio.h>
#include <stdlib.h>
-#include <iostream>
+#include <iostream> /* ??? */
#include <fstream>
#include <string>
#include <cmath>
Modified: grass/trunk/imagery/i.atcorr/i.atcorr.html
===================================================================
--- grass/trunk/imagery/i.atcorr/i.atcorr.html 2011-04-17 18:19:14 UTC (rev 46014)
+++ grass/trunk/imagery/i.atcorr/i.atcorr.html 2011-04-18 06:13:33 UTC (rev 46015)
@@ -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 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 @@
<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/trunk/imagery/i.atcorr/main.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/main.cpp 2011-04-17 18:19:14 UTC (rev 46014)
+++ grass/trunk/imagery/i.atcorr/main.cpp 2011-04-18 06:13:33 UTC (rev 46015)
@@ -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 @@
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(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 @@
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 @@
/* 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;
+ 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.
@@ -271,140 +238,167 @@
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)
- 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 */
- /* 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(); /* ... */
+ /* 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 ? */
- ticache.add(ti, vis[col]);
- }
- }
- else
- {
- pre_compute_v(vis[col]); /* re-compute transformation inputs */
- ti = compute(); /* ... */
- }
- }
+ /* round to nearest visibility bin */
+ /* rounding result: watch out for fp representation error */
+ vis[col] = ((int) (vis[col] + 0.5));
+ }
- if((ialt_fd >= 0) && (prev_alt != alt[col]))
- {
- prev_alt = alt[col]; /* keep track of previous altitude */
+ /* 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(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 (!in_cache) {
+ pre_compute_hv(alt[col], vis[col]); /* re-compute transformation inputs */
+ ti = compute(); /* ... */
- ticache.add(ti, alt[col]);
- }
- }
- else
- {
- pre_compute_h(alt[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;
+ /* 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))
+ 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 @@
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 @@
" 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 @@
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 @@
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 @@
/* 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 @@
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);
}
More information about the grass-commit
mailing list