[GRASS-SVN] r52632 - in grass/trunk/scripts: . i.pansharpen
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Aug 11 23:14:49 PDT 2012
Author: cmbarton
Date: 2012-08-11 23:14:49 -0700 (Sat, 11 Aug 2012)
New Revision: 52632
Added:
grass/trunk/scripts/i.pansharpen/
grass/trunk/scripts/i.pansharpen/Makefile
grass/trunk/scripts/i.pansharpen/i.pansharpen.html
grass/trunk/scripts/i.pansharpen/i.pansharpen.py
grass/trunk/scripts/i.pansharpen/rgb_brovey321.jpg
grass/trunk/scripts/i.pansharpen/rgb_ihs321.jpg
grass/trunk/scripts/i.pansharpen/rgb_landsat321.jpg
grass/trunk/scripts/i.pansharpen/rgb_pca321.jpg
Log:
Added new pan sharpening script
Added: grass/trunk/scripts/i.pansharpen/Makefile
===================================================================
--- grass/trunk/scripts/i.pansharpen/Makefile (rev 0)
+++ grass/trunk/scripts/i.pansharpen/Makefile 2012-08-12 06:14:49 UTC (rev 52632)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=i.pansharpen
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass/trunk/scripts/i.pansharpen/i.pansharpen.html
===================================================================
--- grass/trunk/scripts/i.pansharpen/i.pansharpen.html (rev 0)
+++ grass/trunk/scripts/i.pansharpen/i.pansharpen.html 2012-08-12 06:14:49 UTC (rev 52632)
@@ -0,0 +1,189 @@
+<h2>DESCRIPTION</h2>
+
+<em><b>i.pansharpen</b></em> uses a high resolution panchromatic band from a
+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,
+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
+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'
+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
+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
+3 new bands at the higher resolution (example for band 1):
+<pre>
+ band1
+ new band1 = ----------------------- * panband
+ band1 + band2 + band3
+</pre>
+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
+
+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
+band using the following algorithm (example for band 1):
+<pre>
+
+band1' = pan * EV1-1 + PC2 * EV2-1 + PC3 * EV3-1 + 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
+bands include Landsat 7 ETM, QuickBird, and SPOT.
+<br>
+<h2>NOTES</h2>
+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
+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,
+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
+<em>r.composite</em>. Colors may be optionally optimized with <em>i.landsat.rgb</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.
+
+<h2>EXAMPLE</h2>
+
+Pan sharpening of a Landsat image from Boulder, Colorado, USA:
+
+<div class="code"><pre>
+# R, G, B composite at 30m
+d.rgb b=p034r032_7dt20010924_z13_10 g=lp034r032_7dt20010924_z13_20
+ r=p034r032_7dt20010924_z13_30
+
+# i.pansharpen with IHS algorithm
+i.pansharpen ms3=p034r032_7dt20010924_z13_30 ms2=p034r032_7dt20010924_z13_20
+ ms1=p034r032_7dt20010924_z13_10 pan=p034r032_7dp20010924_z13_80
+ output_prefix=ihs321 sharpen=ihs
+
+# display at 15m
+d.rgb b=ihs321_blue g=ihs321_green r=ihs321_red
+</pre></div>
+
+
+<b><em>Results:</em></b>
+
+<p><center>
+ <table border=1>
+ <tr>
+ <td align=center>
+ <img src="rgb_landsat321.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>
+ <img src="rgb_brovey321.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>
+ </font>
+ </td>
+ </tr>
+ <tr>
+ <td align=center>
+ <img src="rgb_ihs321.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>
+ <img src="rgb_pca321.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>
+ </font>
+ </td>
+ </tr>
+ </table>
+</center>
+<br>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.his.rgb.html">i.his.rgb</a>,
+<a href="i.rgb.his.html">i.rgb.his</a>,
+<a href="i.pca">i.pca</a>,
+<a href="d.rgb.html">d.rgb</a>,
+<a href="r.composite.html">r.composite</a>
+</em>
+
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li>Original Brovey formula reference unknown, probably... <br>
+ Roller, N.E.G. and Cox, S., (1980). Comparison of Landsat MSS
+ and merged MSS/RBV data for analysis of natural vegetation.
+ 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
+ classification. International Journal of Remote Sensing, 25(17), 3529-3539.
+
+<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
+ 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
+ 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
+ classification. International Journal of Geoinformatics, 1(1):51-61
+ (<a href="http://www.grassbook.org/neteler/papers/neteler2005_IJG_051-061_draft.pdf">PDF</a>)
+
+<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>
+
+
+
+<h2>AUTHORS</h2>
+
+Michael Barton (Arizona State University, USA)--with contributions from Markus
+Neteler (ITC-irst, Italy); Glynn Clements; Luca Delucchi (Fondazione E. Mach,
+Italy); Markus Metz; and Hamish Bowman.
+<p><i>Last changed: $Date: $</i>
Added: grass/trunk/scripts/i.pansharpen/i.pansharpen.py
===================================================================
--- grass/trunk/scripts/i.pansharpen/i.pansharpen.py (rev 0)
+++ grass/trunk/scripts/i.pansharpen/i.pansharpen.py 2012-08-12 06:14:49 UTC (rev 52632)
@@ -0,0 +1,435 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: i.pansharpen
+#
+# AUTHOR(S): Overall script by Michael Barton (ASU)
+# Brovey transformation in i.fusion.brovey by Markus Neteler <<neteler at osgeo org>>
+# i.fusion brovey converted to Python by Glynn Clements
+# IHS and PCA transformation added by Michael Barton (ASU)
+# histogram matching algorithm by Michael Barton and Luca Delucchi, Fondazione E. Mach (Italy)
+# Thanks to Markus Metz for help with PCA inversion
+# Thanks to Hamish Bowman for parallel processing algorithm
+#
+# PURPOSE: Sharpening of 3 RGB channels using a high-resolution panchromatic channel
+#
+# COPYRIGHT: (C) 2002-2012 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.
+#
+# REFERENCES:
+# Roller, N.E.G. and Cox, S., 1980. Comparison of Landsat MSS and merged MSS/RBV
+# data for analysis of natural vegetation. Proc. of the 14th International
+# Symposium on Remote Sensing of Environment, San Jose, Costa Rica, 23-30 April, pp. 1001-1007.
+#
+# Amarsaikhan, D., & Douglas, T. (2004). Data fusion and multisource image classification.
+# International Journal of Remote Sensing, 25(17), 3529-3539.
+#
+# Behnia, P. (2005). Comparison between four methods for data fusion of ETM+
+# multispectral and pan images. Geo-spatial Information Science, 8(2), 98-103
+#
+# for LANDSAT 5: see Pohl, C 1996 and others
+#
+#############################################################################
+
+#%Module
+#% description: Image fusion algorithms to sharpen multispectral with high-res panchromatic channels
+#% keywords: imagery
+#% keywords: fusion
+#% keywords: sharpen
+#% keywords: Brovey
+#% keywords: IHS
+#% keywords: PCA
+#% overwrite: yes
+#%End
+#%option
+#% key: sharpen
+#% description: Choose pan sharpening method
+#% options: brovey,ihs,pca
+#% answer: ihs
+#% required: yes
+#%end
+#%option G_OPT_R_INPUT
+#% key: ms3
+#% description: Input raster map for red channel
+#%end
+#%option G_OPT_R_INPUT
+#% key: ms2
+#% description: Input raster map for green channel
+#%end
+#%option G_OPT_R_INPUT
+#% key: ms1
+#% description: Input raster map for blue channel
+#%end
+#% option G_OPT_R_INPUT
+#% key: pan
+#% description: Input raster map for high resolution panchromatic channel
+#%end
+#%option
+#% key: output_prefix
+#% type: string
+#% description: Prefix for output raster maps
+#% required : yes
+#%end
+#%flag
+#% key: s
+#% description: Serial processing rather than parallel processing
+#%end
+#%flag
+#% key: l
+#% description: Rebalance blue channel for landsat maps
+#%end
+
+import sys
+import os
+import numpy as np
+import grass.script as grass
+
+def main():
+ sharpen = options['sharpen'] # sharpening algorithm
+ ms1 = options['ms1'] # blue channel
+ ms2 = options['ms2'] # green channel
+ ms3 = options['ms3'] # red channel
+ pan = options['pan'] # high res pan channel
+ out = options['output_prefix'] # prefix for output RGB maps
+ bladjust = flags['l'] # adjust blue channel
+ sproc = flags['s'] # serial processing
+
+ outb = grass.core.find_file('%s_blue' % out)
+ outg = grass.core.find_file('%s_green' % out)
+ outr = grass.core.find_file('%s_red' % out)
+
+ if (outb['name'] != '' or outg['name'] != '' or outr['name'] != '') and \
+ (not grass.overwrite() and not flags['o']):
+ grass.warning(_('Maps with selected output prefix names already exist. \
+ Delete them or use overwrite flag'))
+ return
+
+ pid = str(os.getpid())
+
+ #get PAN resolution:
+ kv = grass.raster_info(map = pan)
+ nsres = kv['nsres']
+ ewres = kv['ewres']
+ panres = (nsres + ewres) / 2
+
+ # clone current region
+ grass.use_temp_region()
+
+ grass.run_command('g.region', res = panres, align = pan)
+
+ grass.message('\n ')
+ grass.message(_("Performing pan sharpening with hi res pan image: %f" % panres))
+
+ if sharpen == "brovey":
+ grass.message(_("Using Brovey algorithm"))
+
+ #pan/intensity histogram matching using linear regression
+ outname = 'tmp%s_pan1' % pid
+ panmatch1 = matchhist(pan, ms1, outname)
+
+ outname = 'tmp%s_pan2' % pid
+ panmatch2 = matchhist(pan, ms2, outname)
+
+ outname = 'tmp%s_pan3' % pid
+ panmatch3 = matchhist(pan, ms3, outname)
+
+ outr = '%s_red' % out
+ outg = '%s_green' % out
+ outb = '%s_blue' % out
+
+ #calculate brovey transformation
+ grass.message('\n ')
+ grass.message(_("Calculating Brovey transformation..."))
+
+ if sproc:
+ # serial processing
+ e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
+ "$outr" = "$ms3" * "$panmatch3" / k
+ "$outg" = "$ms2" * "$panmatch2" / k
+ "$outb" = "$ms1" * "$panmatch1" / k'''
+ grass.mapcalc(e, outr=outr, outg=outg, outb=outb,
+ panmatch1=panmatch1, panmatch2=panmatch2, panmatch3=panmatch3,
+ ms1=ms1, ms2=ms2, ms3=ms3, overwrite=True)
+ else:
+ # parallel processing
+ pb = grass.mapcalc_start('%s_blue = (%s * %s) / (%s + %s + %s)' %
+ (out, ms1, panmatch1, ms1, ms2, ms3),
+ overwrite=True)
+ pg = grass.mapcalc_start('%s_green = (%s * %s) / (%s + %s + %s)' %
+ (out, ms2, panmatch2, ms1, ms2, ms3),
+ overwrite=True)
+ pr = grass.mapcalc_start('%s_red = (%s * %s) / (%s + %s + %s)' %
+ (out, ms3, panmatch3, ms1, ms2, ms3),
+ overwrite=True)
+
+ pb.wait()
+ pg.wait()
+ pr.wait()
+
+
+ # Cleanup
+ grass.run_command('g.remove', quiet=True, rast='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
+
+ elif sharpen == "ihs":
+ grass.message(_("Using IHS<->RGB algorithm"))
+ #transform RGB channels into IHS color space
+ grass.message('\n ')
+ grass.message(_("Transforming to IHS color space..."))
+ grass.run_command('i.rgb.his', overwrite=True,
+ red_input=ms3,
+ green_input=ms2,
+ blue_input=ms1,
+ hue_output="tmp%s_hue" % pid,
+ intensity_output="tmp%s_int" % pid,
+ saturation_output="tmp%s_sat" % pid)
+
+ #pan/intensity histogram matching using linear regression
+ target = "tmp%s_int" % pid
+ outname = "tmp%s_pan_int" % pid
+ panmatch = matchhist(pan, target, outname)
+
+ #substitute pan for intensity channel and transform back to RGB color space
+ grass.message('\n ')
+ grass.message(_("Transforming back to RGB color space and sharpening..."))
+ grass.run_command('i.his.rgb', overwrite=True,
+ hue_input="tmp%s_hue" % pid,
+ intensity_input="%s" % panmatch,
+ saturation_input="tmp%s_sat" % pid,
+ red_output="%s_red" % out,
+ green_output="%s_green" % out,
+ blue_output="%s_blue" % out)
+
+ # Cleanup
+ grass.run_command('g.remove', quiet=True, rast=panmatch)
+
+ elif sharpen == "pca":
+ grass.message(_("Using PCA/inverse PCA algorithm"))
+ grass.message('\n ')
+ grass.message(_("Creating PCA images and calculating eigenvectors..."))
+
+ #initial PCA with RGB channels
+ pca_out = grass.read_command('i.pca', quiet=True, rescale='0,0', input='%s,%s,%s' % (ms1, ms2, ms3), output_prefix='tmp%s.pca' % pid)
+ b1evect = []
+ b2evect = []
+ b3evect = []
+ for l in pca_out.replace('(',',').replace(')',',').splitlines():
+ b1evect.append(float(l.split(',')[1]))
+ b2evect.append(float(l.split(',')[2]))
+ b3evect.append(float(l.split(',')[3]))
+
+ #inverse PCA with hi res pan channel substituted for principal component 1
+ pca1 = 'tmp%s.pca.1' % pid
+ pca2 = 'tmp%s.pca.2' % pid
+ pca3 = 'tmp%s.pca.3' % pid
+ b1evect1 = b1evect[0]
+ b1evect2 = b1evect[1]
+ b1evect3 = b1evect[2]
+ b2evect1 = b2evect[0]
+ b2evect2 = b2evect[1]
+ b2evect3 = b2evect[2]
+ b3evect1 = b3evect[0]
+ b3evect2 = b3evect[1]
+ b3evect3 = b3evect[2]
+
+ outname = 'tmp%s_pan' % pid
+ panmatch = matchhist(pan, ms1, outname)
+
+ grass.message('\n ')
+ grass.message(_("Performing inverse PCA ..."))
+
+ stats1 = grass.parse_command("r.univar", map=ms1, flags='g',
+ parse=(grass.parse_key_val, { 'sep' : '=' }))
+ stats2 = grass.parse_command("r.univar", map=ms2, flags='g',
+ parse=(grass.parse_key_val, { 'sep' : '=' }))
+ stats3 = grass.parse_command("r.univar", map=ms3, flags='g',
+ parse=(grass.parse_key_val, { 'sep' : '=' }))
+
+ b1mean = float(stats1['mean'])
+ b2mean = float(stats2['mean'])
+ b3mean = float(stats3['mean'])
+
+
+ if sproc:
+ # serial processing
+ e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
+ "$outr" = "$ms3" * "$panmatch3" / k
+ "$outg" = "$ms2" * "$panmatch2" / k
+ "$outb" = "$ms1" * "$panmatch1" / k'''
+
+ outr = '%s_red' % out
+ outg = '%s_green' % out
+ outb = '%s_blue' % out
+
+ cmd1 = "$outb = ($panmatch * $b1evect1) + ($pca2 * $b2evect1) + ($pca3 * $b3evect1) + $b1mean"
+ cmd2 = "$outg = ($panmatch * $b1evect2) + ($pca2 * $b2evect1) + ($pca3 * $b3evect2) + $b2mean"
+ cmd3 = "$outr = ($panmatch * $b1evect3) + ($pca2 * $b2evect3) + ($pca3 * $b3evect3) + $b3mean"
+
+ cmd = '\n'.join([cmd1, cmd2, cmd3])
+
+ grass.mapcalc(cmd, outb=outb, outg=outg, outr=outr,
+ panmatch=panmatch, pca2=pca2, pca3=pca3,
+ b1evect1=b1evect1, b2evect1=b2evect1, b3evect1=b3evect1,
+ b1evect2=b1evect2, b2evect2=b2evect2, b3evect2=b3evect2,
+ b1evect3=b1evect3, b2evect3=b2evect3, b3evect3=b3evect3,
+ b1mean=b1mean, b2mean=b2mean, b3mean=b3mean,
+ overwrite=True)
+ else:
+ # parallel processing
+ pb = grass.mapcalc_start('%s_blue = (%s * %f) + (%s * %f) + (%s * %f) + %f'
+ % (out, panmatch, b1evect1, pca2, b2evect1, pca3, b3evect1, b1mean),
+ overwrite=True)
+
+ pg = grass.mapcalc_start('%s_green = (%s * %f) + (%s * %f) + (%s * %f) + %f'
+ % (out, panmatch, b1evect2, pca2, b2evect2, pca3, b3evect2, b2mean),
+ overwrite=True)
+
+ pr = grass.mapcalc_start('%s_red = (%s * %f) + (%s * %f) + (%s * %f) + %f'
+ % (out, panmatch, b1evect3, pca2, b2evect3, pca3, b3evect3, b3mean),
+ overwrite=True)
+
+ pr.wait()
+ pg.wait()
+ pb.wait()
+
+
+ # Cleanup
+ grass.run_command('g.mremove', flags='f', quiet=True, rast='tmp%s*,%s' % (pid,panmatch))
+
+ #Could add other sharpening algorithms here, e.g. wavelet transformation
+
+ grass.message('\n ')
+ grass.message(_("Assigning grey equalized color tables to output images..."))
+ #equalized grey scales give best contrast
+ for ch in ['red', 'green', 'blue']:
+ grass.run_command('r.colors', quiet=True, map = "%s_%s" % (out, ch), flags="e", col = 'grey')
+
+ #Landsat too blue-ish because panchromatic band less sensitive to blue light,
+ # so output blue channed can be modified
+ if bladjust:
+ grass.message('\n ')
+ 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)
+ os.remove(rules)
+
+ #output notice
+ grass.message('\n ')
+ grass.message(_("The following pan-sharpened output maps have been generated:"))
+ for ch in ['red', 'green', 'blue']:
+ grass.message(_("%s_%s") % (out, ch))
+
+ grass.message('\n ')
+ grass.message(_("To visualize output, run: g.region -p rast=%s.red" % out))
+ grass.message(_("d.rgb r=%s_red g=%s_green b=%s_blue" % (out, out, out)))
+ grass.message('\n ')
+ grass.message(_("If desired, combine channels into a single RGB map with 'r.composite'."))
+ grass.message(_("Channel colors can be rebalanced using i.landsat.rgb."))
+
+ # write cmd history:
+ for ch in ['red', 'green', 'blue']:
+ grass.raster_history("%s_%s" % (out, ch))
+
+ # Cleanup
+ grass.run_command('g.mremove', flags="f", rast="tmp%s*" % pid, quiet=True)
+
+
+def matchhist(original, target, matched):
+ #pan/intensity histogram matching using numpy arrays
+ grass.message('\n ')
+ grass.message(_("Histogram matching..."))
+
+ # input images
+ original = original.split('@')[0]
+ target = target.split('@')[0]
+ images = [original, target]
+
+ # create a dictionary to hold arrays for each image
+ arrays = {}
+
+ for i 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, fs=':')
+ stats = stats_out.communicate()[0].split('\n')[:-1]
+ stats_dict = dict( s.split(':', 1) for s in stats)
+ total_cells = 0 # total non-null cells
+ for j in stats_dict:
+ stats_dict[j] = int(stats_dict[j])
+ if j != '*':
+ total_cells += stats_dict[j]
+
+ # Make a 2x256 structured array for each image with a
+ # 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'))
+ cum_cells = 0 # cumulative total of cells for sum of current and all lower grey values
+
+ for n in range(0,256):
+ if str(n) in stats_dict:
+ num_cells = stats_dict[str(n)]
+ else:
+ num_cells = 0
+
+ cum_cells += num_cells
+
+ # cdf is the the number of cells at or below a given grey value
+ # divided by the total number of cells
+ cdf = float(cum_cells) / float(total_cells)
+
+ # insert values into array
+ arrays[i][n] = (n, cdf)
+
+ # open file for reclass rules
+ outfile = open(grass.tempfile(), 'w')
+ new_grey = 0
+
+ for i in arrays[original]:
+ # for each grey value and corresponding cdf value in original, find the
+ # cdf value in target that is closest to the target cdf value
+ difference_list = []
+ for j in arrays[target]:
+ # make a list of the difference between each original cdf value and
+ # the target cdf value
+ difference_list.append(abs(i[1] - j[1]))
+
+ # get the smallest difference in the list
+ min_difference = min(difference_list)
+
+ for j in arrays[target]:
+ # find the grey value in target that correspondes to the cdf
+ # closest to the original cdf
+ if j[1] == i[1] + min_difference or j[1] == i[1] - min_difference:
+ # build a reclass rules file from the original grey value and
+ # corresponding grey value from target
+ out_line = "%d = %d\n" % (i[0], j[0])
+ outfile.write(out_line)
+ break
+
+ outfile.close()
+
+ # create reclass of target from reclass rules file
+ result = grass.core.find_file(matched, element = 'cell')
+ if result['fullname']:
+ grass.run_command('g.remove', quiet=True, rast=matched)
+ grass.run_command('r.reclass', input=original, out=matched, rules=outfile.name)
+ else:
+ grass.run_command('r.reclass', input=original, out=matched, rules=outfile.name)
+
+ # Cleanup
+ # remove the rules file
+ grass.try_remove(outfile.name)
+
+ # return reclass of target with histogram that matches original
+ return matched
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
Property changes on: grass/trunk/scripts/i.pansharpen/i.pansharpen.py
___________________________________________________________________
Added: svn:executable
+ *
Added: grass/trunk/scripts/i.pansharpen/rgb_brovey321.jpg
===================================================================
(Binary files differ)
Property changes on: grass/trunk/scripts/i.pansharpen/rgb_brovey321.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass/trunk/scripts/i.pansharpen/rgb_ihs321.jpg
===================================================================
(Binary files differ)
Property changes on: grass/trunk/scripts/i.pansharpen/rgb_ihs321.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass/trunk/scripts/i.pansharpen/rgb_landsat321.jpg
===================================================================
(Binary files differ)
Property changes on: grass/trunk/scripts/i.pansharpen/rgb_landsat321.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass/trunk/scripts/i.pansharpen/rgb_pca321.jpg
===================================================================
(Binary files differ)
Property changes on: grass/trunk/scripts/i.pansharpen/rgb_pca321.jpg
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
More information about the grass-commit
mailing list