add scale option to i.topo.corr

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

imagery/i.topo.corr/correction.c

@@ -11,6 +11,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <float.h>
 #include <unistd.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
@@ -19,20 +20,25 @@
 #include "local_proto.h"
 
 void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
-	       double zenith)
+	       double zenith, int do_scale)
 {
     int row, col, nrows, ncols;
     void *pref, *pcos;
 
-    double cos_z, cos_i, ref_i;
+    double cos_z, cos_i, ref_i, result;
     double n, sx, sxx, sy, sxy, tx, ty;
     double a, m, cka, ckb, kk;
 
+    double imin, imax, omin, omax, factor;    /* for scaling to input */
 
     nrows = Rast_window_rows();
     ncols = Rast_window_cols();
 
     cos_z = cos(D2R * zenith);
+    
+    imin = omin = DBL_MAX;
+    imax = omax = -DBL_MAX;
+    factor = 1;
 
     /* Calculating regression */
     if (method > NON_LAMBERTIAN) {
@@ -53,6 +59,12 @@ void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
 		if (!Rast_is_null_value(pref, band->type) &&
 		    !Rast_is_null_value(pcos, cosi->type)) {
 		    ref_i = Rast_get_d_value(pref, band->type);
+		    
+		    if (imin > ref_i)
+			imin = ref_i;
+		    if (imax < ref_i)
+			imax = ref_i;
+		    
 		    switch (method) {
 		    case MINNAERT:
 			if (cos_i > 0. && cos_z > 0. && ref_i > 0.) {
@@ -106,6 +118,45 @@ void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
 	cka = ckb = 0.;
 	kk = 1.;
     }
+
+    if (do_scale) {
+	/* Topographic correction, pass 1 */
+	for (row = 0; row < nrows; row++) {
+	    G_percent(row, nrows, 2);
+
+	    Rast_get_row(band->fd, band->rast, row, band->type);
+	    Rast_get_row(cosi->fd, cosi->rast, row, cosi->type);
+
+	    pref = band->rast;
+	    pcos = cosi->rast;
+	    
+	    for (col = 0; col < ncols; col++) {
+
+		cos_i = Rast_get_d_value(pcos, cosi->type);
+
+		if (!Rast_is_null_value(pref, band->type) &&
+		    !Rast_is_null_value(pcos, cosi->type)) {
+
+		    ref_i = Rast_get_d_value(pref, band->type);
+		    result = (DCELL) (ref_i * pow((cos_z + cka) /
+		                      (cos_i + ckb), kk));
+		    G_debug(3,
+			    "Old val: %f, cka: %f, cos_i: %f, ckb: %f, kk: %f, New val: %f",
+			    ref_i, cka, cos_i, ckb, kk, result);
+
+		    if (omin > result)
+			omin = result;
+		    if (omax < result)
+			omax = result;
+
+		}
+		pref = G_incr_void_ptr(pref, Rast_cell_size(band->type));
+		pcos = G_incr_void_ptr(pcos, Rast_cell_size(cosi->type));
+	    }
+	}
+	G_percent(1, 1, 1);
+	factor = (imax - imin) / (omax - omin);
+    }
     /* Topographic correction */
     for (row = 0; row < nrows; row++) {
 	G_percent(row, nrows, 2);
@@ -126,8 +177,13 @@ void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
 		!Rast_is_null_value(pcos, cosi->type)) {
 
 		ref_i = Rast_get_d_value(pref, band->type);
-		((DCELL *) out->rast)[col] =
-		    (DCELL) (ref_i * pow((cos_z + cka) / (cos_i + ckb), kk));
+		result = (DCELL) (ref_i * pow((cos_z + cka) /
+		                  (cos_i + ckb), kk));
+
+		if (do_scale)
+		    result = (result - omin) * factor + imin;
+
+		((DCELL *) out->rast)[col] = result;
 		G_debug(3,
 			"Old val: %f, cka: %f, cos_i: %f, ckb: %f, kk: %f, New val: %f",
 			ref_i, cka, cos_i, ckb, kk,
@@ -138,4 +194,5 @@ void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
 	}
 	Rast_put_row(out->fd, out->rast, out->type);
     }
+    G_percent(1, 1, 1);
 }

imagery/i.topo.corr/local_proto.h

@@ -23,6 +23,6 @@ typedef struct
 #define C_CORRECT			12
 
 void eval_cosi(Gfile *, Gfile *, double, double);
-void eval_tcor(int, Gfile *, Gfile *, Gfile *, double);
+void eval_tcor(int, Gfile *, Gfile *, Gfile *, double, int);
 
 #endif

imagery/i.topo.corr/main.c

@@ -32,11 +32,12 @@ int main(int argc, char *argv[])
 
     int i;
     struct Option *base, *output, *input, *zeni, *azim, *metho;
-    struct Flag *ilum;
+    struct Flag *ilum, *scl;
 
     Gfile dem, out, band;
     double zenith, azimuth;
     int method = COSINE;
+    int do_scale;
 
     /* initialize GIS environment */
     G_gisinit(argv[0]);
@@ -88,6 +89,10 @@ int main(int argc, char *argv[])
     ilum->key = 'i';
     ilum->description = _("Output sun illumination terrain model");
     
+    scl = G_define_flag();
+    scl->key = 's';
+    scl->description = _("Scale output to input and copy color rules");
+    
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -100,6 +105,7 @@ int main(int argc, char *argv[])
 
     zenith = atof(zeni->answer);
     out.type = DCELL_TYPE;
+    do_scale = scl->answer;
 
     /* Evaluate only cos_i raster file */
     /* i.topo.corr -i out=cosi.on07 base=SRTM_v2 zenith=33.3631 azimuth=59.8897 */
@@ -180,7 +186,7 @@ int main(int argc, char *argv[])
 	    band.rast = Rast_allocate_buf(band.type);
 	    dem.rast = Rast_allocate_buf(dem.type);
 	    /* ----- */
-	    eval_tcor(method, &out, &dem, &band, zenith);
+	    eval_tcor(method, &out, &dem, &band, zenith, do_scale);
 	    /* ----- */
 	    G_free(dem.rast);
 	    Rast_close(dem.fd);
@@ -196,10 +202,18 @@ int main(int argc, char *argv[])
 		struct FPRange range;
 		DCELL min, max;
 		struct Colors grey;
-
-		Rast_read_fp_range(out.name, G_mapset(), &range);
-		Rast_get_fp_range_min_max(&range, &min, &max);
-		Rast_make_grey_scale_colors(&grey, min, max);
+		int make_colors = 1;
+
+		if (do_scale) {
+		    if (Rast_read_colors(band.name, "", &grey) >= 0)
+			make_colors = 0;
+		}
+		
+		if (make_colors) {
+		    Rast_read_fp_range(out.name, G_mapset(), &range);
+		    Rast_get_fp_range_min_max(&range, &min, &max);
+		    Rast_make_grey_scale_colors(&grey, min, max);
+		}
 		Rast_write_colors(out.name, G_mapset(), &grey);
 	    }
 	}