added new channel stretch option and updated manual

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

scripts/i.pansharpen/i.pansharpen.html

@@ -5,78 +5,86 @@ multispectral image to sharpen 3 lower resolution bands. The 3
 lower resolution bands can then be combined into an RGB color image at a
 higher (more detailed) resolution than is possible using the original 3
 bands. For example, Landsat ETM has low resolution spectral bands 1 (blue),
-2 (green), 3 (red), 4 (near IR), 5 (mid-IR), and 7 (mid-IR) at 30m resolution, 
+2 (green), 3 (red), 4 (near IR), 5 (mid-IR), and 7 (mid-IR) at 30m resolution,
 and a high resolution panchromatic band 8 at 15m resolution. Pan sharpening
-allows bands 3-2-1 (or other combinations of 30m resolution bands like 4-3-2 
+allows bands 3-2-1 (or other combinations of 30m resolution bands like 4-3-2
 or 5-4-2) to be combined into a 15m resolution color image.
 <br><br>
-i.pansharpen offers a choice of three different 'pan sharpening' 
+i.pansharpen offers a choice of three different 'pan sharpening'
 algorithms: IHS, Brovey, and PCA.
 <br><br>
-For <em>IHS pan sharpening</em>, the original 3 lower resolution bands, selected 
-as red, green and blue channels for creating an RGB composite image, are 
-transformed into IHS (intensity, hue, and saturation) color space. The 
-panchromatic band is then substituted for the intensity channel (I), combined 
-with the original hue (H) and saturation (S) channels, and transformed back to 
-RGB color space at the higher resolution of the panchromatic band. The 
+For <em>IHS pan sharpening</em>, the original 3 lower resolution bands, selected
+as red, green and blue channels for creating an RGB composite image, are
+transformed into IHS (intensity, hue, and saturation) color space. The
+panchromatic band is then substituted for the intensity channel (I), combined
+with the original hue (H) and saturation (S) channels, and transformed back to
+RGB color space at the higher resolution of the panchromatic band. The
 algorithm for this can be represented as: RGB -> IHS -> [pan]HS -> RGB.
 <br><br>
-With a <em>Brovey pan sharpening</em>, each of the 3 lower resolution bands and 
-panchromatic band are combined using the following algorithm to calculate 
+With a <em>Brovey pan sharpening</em>, each of the 3 lower resolution bands and
+panchromatic band are combined using the following algorithm to calculate
 3 new bands at the higher resolution (example for band 1):
 <pre>
-                         band1 
+                         band1
     new band1 = ----------------------- * panband
                  band1 + band2 + band3
 </pre>
-In <em>PCA pan sharpening</em>, a principal component analysis is performed on the 
+In <em>PCA pan sharpening</em>, a principal component analysis is performed on the
 original 3 lower resolution bands to create 3 principal component images
 (PC1, PC2, and PC3) and their associated eigenvectors (EV), such that:
 <pre>
-    
+
      band1  band2  band3
 PC1: EV1-1  EV1-2  EV1-3
 PC2: EV2-1  EV2-2  EV2-3
-PC3: EV3-1  EV3-2  EV3-3  
+PC3: EV3-1  EV3-2  EV3-3
 
 and
 
 PC1 = EV1-1 * band1 + EV1-2 * band2 + EV1-3 * band3 - mean(bands 1,2,3)
 
 </pre>
-An inverse PCA is then performed, substituting the panchromatic band for PC1. 
-To do this, the eigenvectors matrix is inverted (in this case transposed), the 
-PC images are multiplied by the eigenvectors with the panchromatic band 
-substituted for PC1, and mean of each band is added to each transformed image 
+An inverse PCA is then performed, substituting the panchromatic band for PC1.
+To do this, the eigenvectors matrix is inverted (in this case transposed), the
+PC images are multiplied by the eigenvectors with the panchromatic band
+substituted for PC1, and mean of each band is added to each transformed image
 band using the following algorithm (example for band 1):
 <pre>
 
-band1' = pan * EV1-1 + PC2 * EV2-1 + PC3 * EV3-1 + mean(band1)
-   
+band1 = pan * EV1-1 + PC2 * EV1-2 + PC3 * EV1-3 + mean(band1)
+
 </pre>
-The assignment of the channels depends on the satellite. Examples of satellite 
-imagery with high resolution panchromatic bands, and lower resolution spectral 
+The assignment of the channels depends on the satellite. Examples of satellite
+imagery with high resolution panchromatic bands, and lower resolution spectral
 bands include Landsat 7 ETM, QuickBird, and SPOT.
 <br>
 <h2>NOTES</h2>
-The module currently only works for 8-bit images.
+The module works for 2-bit to 30-bit images. All images are rescaled to 8-bit
+for processing. By default, the entire possible range for the selected bit depth is
+rescaled to 8-bit. For example, the range of 0-65535 for a 16-bit image is
+rescaled to 0-255). The 'r' flag allows the range of pixel values actually 
+present in an image rescaled to a full 8-bit range. For example, a 16 bit image
+might only have pixels that range from 70 to 35000; this range of 70-35000 would
+be rescaled to 0-255. This can give better visual distinction to features,
+especially when the range of actual values in an image only occupies a
+relatively limited portion of the possible range.
 <br><br>
-The command temporarily changes the computational region to the high 
-resolution of the panchromatic band during sharpening calculations, then 
-restores the previous region settings. The current region coordinates (and 
-null values) are respected. The high resolution panchromatic image is 
-histogram matched to the band it is replaces prior to substitution (i.e., the 
-intensity channel for IHS sharpening, the low res band selected for each color 
+i.pansharpen temporarily changes the computational region to the high
+resolution of the panchromatic band during sharpening calculations, then
+restores the previous region settings. The current region coordinates (and
+null values) are respected. The high resolution panchromatic image is
+histogram matched to the band it is replaces prior to substitution (i.e., the
+intensity channel for IHS sharpening, the low res band selected for each color
 channel with Brovey sharpening, and the PC1 image for PCA sharpening).
 <br><br>
-By default, the command will attempt to employ parallel processing, using 
-up to 3 cores simultaneously. The -s flag will disable parallel processing, 
+By default, the command will attempt to employ parallel processing, using
+up to 3 cores simultaneously. The -s flag will disable parallel processing,
 but does use an optimized r.mapcalc expression to reduce disk I/O.
 <br><br>
-The three pan-sharpened output channels may be combined with <em>d.rgb</em> or 
+The three pan-sharpened output channels may be combined with <em>d.rgb</em> or
 <em>r.composite</em>. Colors may be optionally optimized with <em>i.colors.enhance</em>.
-While the resulting color image will be at the higher resolution in all cases, 
-the 3 pan sharpening algorithms differ in terms of spectral response.  
+While the resulting color image will be at the higher resolution in all cases,
+the 3 pan sharpening algorithms differ in terms of spectral response.
 
 <h2>EXAMPLES</h2>
 
@@ -85,14 +93,14 @@ the 3 pan sharpening algorithms differ in terms of spectral response.
 Pan sharpening of a Landsat image from Boulder, Colorado, USA:
 
 <div class="code"><pre>
-# R, G, B composite at 30m 
+# R, G, B composite at 30m
 g.region raster=p034r032_7dt20010924_z13_10 -p
-d.rgb b=p034r032_7dt20010924_z13_10 g=lp034r032_7dt20010924_z13_20 
+d.rgb b=p034r032_7dt20010924_z13_10 g=lp034r032_7dt20010924_z13_20
     r=p034r032_7dt20010924_z13_30
 
 # i.pansharpen with IHS algorithm
-i.pansharpen red=p034r032_7dt20010924_z13_30 green=p034r032_7dt20010924_z13_20 
-    blue=p034r032_7dt20010924_z13_10 pan=p034r032_7dp20010924_z13_80 
+i.pansharpen red=p034r032_7dt20010924_z13_30 green=p034r032_7dt20010924_z13_20
+    blue=p034r032_7dt20010924_z13_10 pan=p034r032_7dp20010924_z13_80
     output=ihs321 method=ihs
 
 # ... likewise with method=brovey and method=pca
@@ -109,14 +117,14 @@ d.rgb b=ihs321_blue g=ihs321_green r=ihs321_red
   <table border=1>
   <tr>
     <td align=center>
-      &nbsp;<img src="i_pansharpen_rgb_landsat321.jpg" alt="R, G, B composite of Landsat at 30m">
+      &nbsp;<img src="i_pansharpen_rgb_landsat542.jpg" alt="R, G, B composite of Landsat at 30m">
       <br>
       <font size="-1">
       <i>R, G, B composite of Landsat at 30m</i>
       </font>
     </td>
     <td align=center>
-      &nbsp;<img src="i_pansharpen_rgb_brovey321.jpg" alt="R, G, B composite of Brovey sharpened image at 15m">
+      &nbsp;<img src="i_pansharpen_rgb_brovey542.jpg" alt="R, G, B composite of Brovey sharpened image at 15m">
       <br>
       <font size="-1">
       <i>R, G, B composite of Brovey sharpened image at 15m</i>
@@ -125,14 +133,14 @@ d.rgb b=ihs321_blue g=ihs321_green r=ihs321_red
   </tr>
   <tr>
     <td align=center>
-      &nbsp;<img src="i_pansharpen_rgb_ihs321.jpg" alt="R, G, B composite of IHS sharpened image at 15m">
+      &nbsp;<img src="i_pansharpen_rgb_ihs542.jpg" alt="R, G, B composite of IHS sharpened image at 15m">
       <br>
       <font size="-1">
       <i>R, G, B composite of IHS sharpened image at 15m</i>
       </font>
     </td>
     <td align=center>
-      &nbsp;<img src="i_pansharpen_rgb_pca321.jpg" alt="R, G, B composite of PCA sharpened image at 15m">
+      &nbsp;<img src="i_pansharpen_rgb_pca542.jpg" alt="R, G, B composite of PCA sharpened image at 15m">
       <br>
       <font size="-1">
       <i>R, G, B composite of PCA sharpened image at 15m"</i>
@@ -175,7 +183,6 @@ i.colors.enhance r=lsat7_2002_ihs_red g=lsat7_2002_ihs_green b=lsat7_2002_ihs_bl
 <h2>SEE ALSO</h2>
 
 <em>
-<a href="i.colors.enhance.html">i.colors.enhance</a>,
 <a href="i.his.rgb.html">i.his.rgb</a>,
 <a href="i.rgb.his.html">i.rgb.his</a>,
 <a href="i.pca.html">i.pca</a>,
@@ -193,26 +200,26 @@ i.colors.enhance r=lsat7_2002_ihs_red g=lsat7_2002_ihs_green b=lsat7_2002_ihs_bl
    Proc. of the 14th International Symposium on Remote Sensing
    of Environment, San Jose, Costa Rica, 23-30 April, pp. 1001-1007
 
-<li>Amarsaikhan, D., Douglas, T. (2004). Data fusion and multisource image 
+<li>Amarsaikhan, D., Douglas, T. (2004). Data fusion and multisource image
    classification. International Journal of Remote Sensing, 25(17), 3529-3539.
 
-<li>Behnia, P. (2005). Comparison between four methods for data fusion of ETM+ 
+<li>Behnia, P. (2005). Comparison between four methods for data fusion of ETM+
    multispectral and pan images. Geo-spatial Information Science, 8(2), 98-103.
-   
-<li>Du, Q., Younan, N. H., King, R., Shah, V. P. (2007). On the Performance 
-   Evaluation of Pan-Sharpening Techniques. Geoscience and Remote Sensing 
+
+<li>Du, Q., Younan, N. H., King, R., Shah, V. P. (2007). On the Performance
+   Evaluation of Pan-Sharpening Techniques. Geoscience and Remote Sensing
    Letters, IEEE, 4(4), 518-522.
 
-<li>Karathanassi, V., Kolokousis, P., Ioannidou, S. (2007). A comparison 
-   study on fusion methods using evaluation indicators. International Journal 
+<li>Karathanassi, V., Kolokousis, P., Ioannidou, S. (2007). A comparison
+   study on fusion methods using evaluation indicators. International Journal
    of Remote Sensing, 28(10), 2309-2341.
 
 <li>Neteler, M, D. Grasso, I. Michelazzi, L. Miori, S. Merler, and C.
-   Furlanello (2005). An integrated toolbox for image registration, fusion and 
+   Furlanello (2005). An integrated toolbox for image registration, fusion and
    classification. International Journal of Geoinformatics, 1(1):51-61
    (<a href="http://www.grassbook.org/wp-content/uploads/neteler/papers/neteler2005_IJG_051-061_draft.pdf">PDF</a>)
 
-<li>Pohl, C, and J.L van Genderen (1998). Multisensor image fusion in remote 
+<li>Pohl, C, and J.L van Genderen (1998). Multisensor image fusion in remote
     sensing: concepts, methods and application. Int. J. of Rem. Sens., 19, 823-854.
 </ul>
 
@@ -222,6 +229,6 @@ i.colors.enhance r=lsat7_2002_ihs_red g=lsat7_2002_ihs_green b=lsat7_2002_ihs_bl
 
 Michael Barton (Arizona State University, USA)<br>
 with contributions from Markus Neteler (ITC-irst, Italy); Glynn Clements;
-Luca Delucchi (Fondazione E. Mach, Italy); Markus Metz; and Hamish Bowman. 
+Luca Delucchi (Fondazione E. Mach, Italy); Markus Metz; and Hamish Bowman.
 
 <p><i>Last changed: $Date$</i>

scripts/i.pansharpen/i.pansharpen.py

@@ -87,6 +87,10 @@
 #% key: l
 #% description: Rebalance blue channel for LANDSAT
 #%end
+#%flag
+#% key: r
+#% description: Rescale (stretch) the range of pixel values in each channel to the entire 0-255 8-bit range for processing (see notes)
+#%end
 
 import os
 
@@ -112,6 +116,7 @@ def main():
     bits      = options['bitdepth'] # bit depth of image channels
     bladjust  = flags['l']  # adjust blue channel
     sproc     = flags['s']  # serial processing
+    rescale   = flags['r'] # rescale to spread pixel values to entire 0-255 range
 
     # Checking bit depth
     bits = float(bits)
@@ -136,57 +141,98 @@ def main():
     ms3 = 'tmp%s_ms3' % pid
     pan = 'tmp%s_pan' % pid
 
-    if bits == 8:
-        grass.message(_("Using 8bit image channels"))
-        if sproc:
-            # serial processing
-            grass.run_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
-                               quiet=True, overwrite=True)
-            grass.run_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
-                               quiet=True, overwrite=True)
-            grass.run_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
-                               quiet=True, overwrite=True)
-            grass.run_command('g.copy', raster='%s,%s' % (pan_orig, pan),
-                               quiet=True, overwrite=True)
-        else:
-            # parallel processing
-            pb = grass.start_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
+    if rescale == False:
+        if bits == 8:
+            grass.message(_("Using 8bit image channels"))
+            if sproc:
+                # serial processing
+                grass.run_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
                                    quiet=True, overwrite=True)
-            pg = grass.start_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
+                grass.run_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
                                    quiet=True, overwrite=True)
-            pr = grass.start_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
+                grass.run_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
                                    quiet=True, overwrite=True)
-            pp = grass.start_command('g.copy', raster='%s,%s' % (pan_orig, pan),
+                grass.run_command('g.copy', raster='%s,%s' % (pan_orig, pan),
                                    quiet=True, overwrite=True)
+            else:
+                # parallel processing
+                pb = grass.start_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
+                                       quiet=True, overwrite=True)
+                pg = grass.start_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
+                                       quiet=True, overwrite=True)
+                pr = grass.start_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
+                                       quiet=True, overwrite=True)
+                pp = grass.start_command('g.copy', raster='%s,%s' % (pan_orig, pan),
+                                       quiet=True, overwrite=True)
+
+                pb.wait()
+                pg.wait()
+                pr.wait()
+                pp.wait()
 
-            pb.wait()
-            pg.wait()
-            pr.wait()
-            pp.wait()
+        else:
+            grass.message(_("Converting image chanels to 8bit for processing"))
+            maxval = pow(2, bits) - 1
+            if sproc:
+                # serial processing
+                grass.run_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+                                   output=ms1, to='0,255', quiet=True, overwrite=True)
+                grass.run_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+                                   output=ms2, to='0,255', quiet=True, overwrite=True)
+                grass.run_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+                                   output=ms3, to='0,255', quiet=True, overwrite=True)
+                grass.run_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+                                   output=pan, to='0,255', quiet=True, overwrite=True)
+
+            else:
+                # parallel processing
+                pb = grass.start_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+                                       output=ms1, to='0,255', quiet=True, overwrite=True)
+                pg = grass.start_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+                                       output=ms2, to='0,255', quiet=True, overwrite=True)
+                pr = grass.start_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+                                       output=ms3, to='0,255', quiet=True, overwrite=True)
+                pp = grass.start_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+                                       output=pan, to='0,255', quiet=True, overwrite=True)
+
+                pb.wait()
+                pg.wait()
+                pr.wait()
+                pp.wait()
 
     else:
-        grass.message(_("Converting image chanels to 8bit for processing"))
+        grass.message(_("Rescaling image chanels to 8bit for processing"))
+
+        min_ms1 = int(grass.raster_info(ms1_orig)['min'])
+        max_ms1 = int(grass.raster_info(ms1_orig)['max'])
+        min_ms2 = int(grass.raster_info(ms2_orig)['min'])
+        max_ms2 = int(grass.raster_info(ms2_orig)['max'])
+        min_ms3 = int(grass.raster_info(ms3_orig)['min'])
+        max_ms3 = int(grass.raster_info(ms3_orig)['max'])
+        min_pan = int(grass.raster_info(pan_orig)['min'])
+        max_pan = int(grass.raster_info(pan_orig)['max'])
+
         maxval = pow(2, bits) - 1
         if sproc:
             # serial processing
-            grass.run_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+            grass.run_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1),
                                output=ms1, to='0,255', quiet=True, overwrite=True)
-            grass.run_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+            grass.run_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2),
                                output=ms2, to='0,255', quiet=True, overwrite=True)
-            grass.run_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+            grass.run_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3),
                                output=ms3, to='0,255', quiet=True, overwrite=True)
-            grass.run_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+            grass.run_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan),
                                output=pan, to='0,255', quiet=True, overwrite=True)
 
         else:
             # parallel processing
-            pb = grass.start_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+            pb = grass.start_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1),
                                    output=ms1, to='0,255', quiet=True, overwrite=True)
-            pg = grass.start_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+            pg = grass.start_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2),
                                    output=ms2, to='0,255', quiet=True, overwrite=True)
-            pr = grass.start_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+            pr = grass.start_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3),
                                    output=ms3, to='0,255', quiet=True, overwrite=True)
-            pp = grass.start_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+            pp = grass.start_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan),
                                    output=pan, to='0,255', quiet=True, overwrite=True)
 
             pb.wait()
@@ -194,6 +240,7 @@ def main():
             pr.wait()
             pp.wait()
 
+
     # get PAN resolution:
     kv = grass.raster_info(map=pan)
     nsres = kv['nsres']
@@ -226,13 +273,12 @@ def main():
     # light, so output blue channed can be modified
     if bladjust:
         grass.message(_("Adjusting blue channel color table..."))
-        rules = grass.tempfile()
-        colors = open(rules, 'w')
-        colors.write('5 0 0 0\n20 200 200 200\n40 230 230 230\n67 255 255 255 \n')
-        colors.close()
-
-        grass.run_command('r.colors', map="%s_blue" % out, rules=rules, quiet=True)
-        os.remove(rules)
+        blue_colors = ['0 0 0 0\n5% 0 0 0\n67% 255 255 255\n100% 255 255 255']
+        # these previous colors are way too blue for landsat
+        # blue_colors = ['0 0 0 0\n10% 0 0 0\n20% 200 200 200\n40% 230 230 230\n67% 255 255 255\n100% 255 255 255']
+        bc = grass.feed_command('r.colors', quiet = True, map = "%s_blue" % out, rules = "-")
+        bc.stdin.write('\n'.join(blue_colors))
+        bc.stdin.close()
 
     # output notice
     grass.verbose(_("The following pan-sharpened output maps have been generated:"))
@@ -499,9 +545,9 @@ def matchhist(original, target, matched):
     # create a dictionary to hold arrays for each image
     arrays = {}
 
-    for i in images:
+    for img in images:
         # calculate number of cells for each grey value for for each image
-        stats_out = grass.pipe_command('r.stats', flags='cin', input=i,
+        stats_out = grass.pipe_command('r.stats', flags='cin', input=img,
                                        sep=':')
         stats = stats_out.communicate()[0].split('\n')[:-1]
         stats_dict = dict(s.split(':', 1) for s in stats)
@@ -518,7 +564,7 @@ def matchhist(original, target, matched):
         #   cumulative distribution function (CDF) for each grey value.
         #   Grey value is the integer (i4) and cdf is float (f4).
 
-        arrays[i] = np.zeros((256, ), dtype=('i4,f4'))
+        arrays[img] = np.zeros((256, ), dtype=('i4,f4'))
         cum_cells = 0  # cumulative total of cells for sum of current and all lower grey values
 
         for n in range(0, 256):
@@ -534,7 +580,7 @@ def matchhist(original, target, matched):
             cdf = float(cum_cells) / float(total_cells)
 
             # insert values into array
-            arrays[i][n] = (n, cdf)
+            arrays[img][n] = (n, cdf)
 
     # open file for reclass rules
     outfile = open(grass.tempfile(), 'w')