[GRASS-SVN] r72153 - grass/trunk/imagery/i.atcorr
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jan 26 01:19:45 PST 2018
Author: mmetz
Date: 2018-01-26 01:19:45 -0800 (Fri, 26 Jan 2018)
New Revision: 72153
Modified:
grass/trunk/imagery/i.atcorr/main.cpp
Log:
i.atcorr: optimize cache with transformation parameters
Modified: grass/trunk/imagery/i.atcorr/main.cpp
===================================================================
--- grass/trunk/imagery/i.atcorr/main.cpp 2018-01-26 09:13:18 UTC (rev 72152)
+++ grass/trunk/imagery/i.atcorr/main.cpp 2018-01-26 09:19:45 UTC (rev 72153)
@@ -37,6 +37,7 @@
#include <cstdlib>
#include <map>
+#include <cmath>
extern "C"
{
@@ -49,13 +50,13 @@
#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.
+/* MODIS surface reflectance (MOD09) uses pre-computed 6S parameters
+ * for 10 aerosol optical depths. Here we use finer grain. */
+/* TICache: create bins for altitude and visibility */
+#define BIN_ALT 10 /* unit is meters */
+#define BIN_VIS 100 /* unit is meters */
+
/* uncomment to disable cache usage */
/* #define _NO_OPTIMIZE_ */
@@ -95,7 +96,6 @@
/* 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);
@@ -129,7 +129,26 @@
return (CELL) (-(-x + .5));
}
+/* round height, input and output unit is meters */
+static int round_h(float x)
+{
+ x /= BIN_ALT;
+ x = floor(x + 0.5);
+ x *= BIN_ALT;
+ return x;
+}
+
+/* round visibility, input unit is km, output unit is meters */
+static int round_v(float x)
+{
+ x *= 1000. / BIN_VIS;
+ x = floor(x + 0.5);
+ x *= BIN_VIS;
+
+ return x;
+}
+
/* Converts the buffer to cell and write it to disk */
static void write_fp_to_cell(int ofd, FCELL *buf)
{
@@ -182,11 +201,8 @@
{
struct RBitem rbitem;
- /* although alt and vis are already rounded,
- * the + 0.5 is needed for fp representation errors */
- /* alt has been converted to kilometers */
- rbitem.alt = (int) (alt * BIN_ALT + 0.5);
- /* vis should be kilometers */
+ /* alt and vis must be in meters */
+ rbitem.alt = (alt < 0 ? (int) (alt - 0.5) : (int) (alt + 0.5));
rbitem.vis = (int) (vis + 0.5);
return rbitem;
@@ -248,8 +264,8 @@
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 cur_alt = -100000;
+ int cur_vis = -100000;
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);
@@ -300,80 +316,77 @@
continue;
}
if (ialt_fd >= 0) {
- if (alt[col] < 0)
- alt[col] = 0; /* on or below sea level, all the same for 6S */
- else
- alt[col] /= 1000.0f; /* converting to km from input which should be in meter */
-
- /* round to nearest altitude bin */
- /* rounding result: watch out for fp representation error */
- alt[col] = ((int) (alt[col] * BIN_ALT + 0.5)) / BIN_ALT;
+ alt[col] = round_h(alt[col]); /* in meters */
}
if (ivis_fd >= 0) {
- if (vis[col] < 0)
- vis[col] = 0; /* negative visibility is invalid, print a WARNING ? */
-
- /* round to nearest visibility bin */
- /* rounding result: watch out for fp representation error */
- vis[col] = ((int) (vis[col] + 0.5));
+ if (vis[col] < 0) {
+ /* negative visibility is invalid */
+ G_warning(_("Negative visibility!"));
+ vis[col] = 0;
+ }
+ if (vis[col] < 5.0) {
+ /* too low visibility, text comes from 6S main.f L109-113 */
+ G_warning(_("The visibility must be better than 5.0km, "
+ "for smaller values calculations might be no more valid!"));
+ }
+ vis[col] = round_v(vis[col]); /* from km to meters */
}
- /* check if both maps are active and if whether any value has changed */
+ /* check if altitude or visibility 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];
+ ((cur_vis != vis[col]) || (cur_alt != alt[col]))) {
+ cur_alt = alt[col];
+ cur_vis = vis[col];
if (optimize) {
- int in_cache = ticache.search(alt[col], vis[col], &ti);
+ int in_cache = ticache.search(cur_alt, cur_vis, &ti);
if (!in_cache) {
- pre_compute_hv(alt[col], vis[col]); /* re-compute transformation inputs */
+ pre_compute_hv(alt[col] / 1000., vis[col] / 1000.); /* re-compute transformation inputs */
ti = compute(); /* ... */
- ticache.add(ti, alt[col], vis[col]);
+ ticache.add(ti, cur_alt, cur_vis);
}
}
else {
- pre_compute_hv(alt[col], vis[col]); /* re-compute transformation inputs */
+ pre_compute_hv(alt[col] / 1000., vis[col] / 1000.); /* re-compute transformation inputs */
ti = compute(); /* ... */
}
}
- else { /* only one of the maps is being used */
+ else {
+ if ((ivis_fd >= 0) && (cur_vis != vis[col])) {
+ cur_vis = vis[col];
- 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);
+ int in_cache = ticache.search(cur_alt, cur_vis, &ti);
if (!in_cache) {
- pre_compute_v(vis[col]); /* re-compute transformation inputs */
+ pre_compute_v(vis[col] / 1000.); /* re-compute transformation inputs */
ti = compute(); /* ... */
- ticache.add(ti, 0, vis[col]);
+ ticache.add(ti, cur_alt, cur_vis);
}
}
else {
- pre_compute_v(vis[col]); /* re-compute transformation inputs */
+ pre_compute_v(vis[col] / 1000.); /* re-compute transformation inputs */
ti = compute(); /* ... */
}
}
- if ((ialt_fd >= 0) && (prev_alt != alt[col])) {
- prev_alt = alt[col]; /* keep track of previous altitude */
+ if ((ialt_fd >= 0) && (cur_alt != alt[col])) {
+ cur_alt = alt[col];
if (optimize) {
- int in_cache = ticache.search(alt[col], 0, &ti);
+ int in_cache = ticache.search(alt[col], cur_vis, &ti);
if (!in_cache) {
- pre_compute_h(alt[col]); /* re-compute transformation inputs */
+ pre_compute_h(alt[col] / 1000.); /* re-compute transformation inputs */
ti = compute(); /* ... */
- ticache.add(ti, alt[col], 0);
+ ticache.add(ti, cur_alt, cur_vis);
}
}
else {
- pre_compute_h(alt[col]); /* re-compute transformation inputs */
+ pre_compute_h(alt[col] / 1000.); /* re-compute transformation inputs */
ti = compute(); /* ... */
}
}
More information about the grass-commit
mailing list