[GRASS-SVN] r65193 - in sandbox/alexandris: . i.fusion.hpf
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon May 4 12:29:42 PDT 2015
Author: nikosa
Date: 2015-05-04 12:29:42 -0700 (Mon, 04 May 2015)
New Revision: 65193
Added:
sandbox/alexandris/i.fusion.hpf/
sandbox/alexandris/i.fusion.hpf/Makefile
sandbox/alexandris/i.fusion.hpf/constants.py
sandbox/alexandris/i.fusion.hpf/high_pass_filter.py
sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.html
sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.py
Log:
Test commit of i.fusion.hpf
Added: sandbox/alexandris/i.fusion.hpf/Makefile
===================================================================
--- sandbox/alexandris/i.fusion.hpf/Makefile (rev 0)
+++ sandbox/alexandris/i.fusion.hpf/Makefile 2015-05-04 19:29:42 UTC (rev 65193)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.fusion.hpf
+
+ETCFILES = constants high_pass_filter
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+include $(MODULE_TOPDIR)/include/Make/Python.make
+
+default: script
Property changes on: sandbox/alexandris/i.fusion.hpf/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: sandbox/alexandris/i.fusion.hpf/constants.py
===================================================================
--- sandbox/alexandris/i.fusion.hpf/constants.py (rev 0)
+++ sandbox/alexandris/i.fusion.hpf/constants.py 2015-05-04 19:29:42 UTC (rev 65193)
@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+#-*- coding:utf-8 -*-
+
+"""
+ at author: Nikos Alexandris | October 2014
+"""
+
+"""
+Constants for the HPFA Image Fusion Technique:
+Kernel Size, Center Value, Modulation Factor (all depend on Resolution Ratio).
+
+Sources:
+
+- "Optimizing the High-Pass Filter Addition Technique for Image Fusion",
+Ute G. Gangkofner, Pushkar S. Pradhan, and Derrold W. Holcomb (2008).
+
+- “ERDAS IMAGINE.” Accessed March 19, 2015.
+http://doc.hexagongeospatial.com/ERDAS%20IMAGINE/ERDAS_IMAGINE_Help/#ii_hpfmerge_mergedialog.htm.
+
+"""
+
+RATIO_RANGES = (
+ (1, 2.5),
+ (2.5, 3.5),
+ (3.5, 5.5),
+ (5.5, 7.5),
+ (7.5, 9.5),
+ (9.5, float('inf')))
+
+KERNEL_SIZES = (5, 7, 9, 11, 13, 15)
+
+MATRIX_PROPERTIES = zip(RATIO_RANGES, KERNEL_SIZES)
+
+
+# Replicating ERDAS' Imagine parameters -------------------------------------
+CENTER_CELL = {
+ 'Low': [24, 48, 80, 120, 168, 336],
+ 'Mid': [28, 56, 93, 150, 210, 392],
+ 'High': [32, 64, 106, 180, 252, 448]}
+
+# Can't find a unique sequence pattern, so... python fun below: =============
+#ks = (5, 7, 9, 11, 13, 15)
+#lo = list([(s**2-1) for s in ks[0:-1]])
+#lo.append(int((ks[-1]**2-1)*1.5))
+#inc = list([v/6 for v in lo[0:3]])
+#inc.extend([v/4 for v in lo[3:5]])
+#inc.append(lo[-1]/6)
+#mi = list([a + b for a, b in zip(lo, inc)])
+#hi = list([a + b for a, b in zip(mi, inc)])
+#CENTER_CELL = {'High': hi, 'Low': lo, 'Mid': mi}
+# ===========================================================================
+
+
+MODULATOR = {
+ 'Min': [0.20, 0.35, 0.35, 0.50, 0.65, 1.00],
+ 'Mid': [0.25, 0.50, 0.50, 0.65, 1.00, 1.35],
+ 'Max': [0.30, 0.65, 0.65, 1.00, 1.40, 2.00]}
+
+MODULATOR_2 = {'Min': 0.25, 'Mid': 0.35, 'Max': 0.50}
Property changes on: sandbox/alexandris/i.fusion.hpf/constants.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: sandbox/alexandris/i.fusion.hpf/high_pass_filter.py
===================================================================
--- sandbox/alexandris/i.fusion.hpf/high_pass_filter.py (rev 0)
+++ sandbox/alexandris/i.fusion.hpf/high_pass_filter.py 2015-05-04 19:29:42 UTC (rev 65193)
@@ -0,0 +1,135 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+ at author: Nikos Alexandris | Created on Mon Nov 3 13:02:59 2014
+"""
+
+import os
+from constants import MATRIX_PROPERTIES, CENTER_CELL, MODULATOR, MODULATOR_2
+
+
+def kernel_size(ratio):
+ """High Pass Filter Additive image fusion compatible kernel size.
+ Based on a float ratio, ranging in (1.0, 10.0).
+ Returns a single integer"""
+ ks = [k for ((lo, hi), k) in MATRIX_PROPERTIES if lo <= ratio < hi]
+ return ks[0] # ks: kernel size, as integer?
+
+
+def center_cell(level, ks):
+ """High Pass Filter Additive image fusion compatible kernel center
+ cell value."""
+ level = level.capitalize()
+ ks_idx = [k for ((lo, hi), k) in MATRIX_PROPERTIES].index(ks)
+ center = [cc for cc in CENTER_CELL[level]][ks_idx]
+ return center
+
+
+def modulator(modulation, modulation2, ks, second_pass):
+ """Returning a modulation factor determining image Cripsness"""
+ ks_idx = [k for ((lo, hi), k) in MATRIX_PROPERTIES].index(ks)
+ if second_pass:
+ modulation2 = modulation2.capitalize()
+ modfac = MODULATOR_2[modulation2]
+ else:
+ modulation = modulation.capitalize()
+ modfac = [mf for mf in MODULATOR[modulation]][ks_idx]
+ return modfac
+
+
+class Kernel:
+ """HPF compatible Kernel (size),
+ where size is odd | Returns multi-line string"""
+ def __init__(self, size, level):
+ self.size = int(size)
+ self.center = center_cell(level, self.size)
+ self.kernel = ''
+
+ # loop over columns, return value for column, row
+ for row in range(self.size):
+
+ # middle row
+ midrow = (self.size/2)
+
+ # fill rows
+ if row != self.size/2:
+ self.kernel += "-1 " * self.size + "\n"
+
+ # fill mid row
+ else:
+
+ # single-digit center?
+ if len(str(self.center)) == 1:
+ self.center = " " + str(self.center)
+
+ # prettier output for double-digit or larger center
+ self.kernel += "-1 " * midrow + str(self.center) + \
+ " " + "-1 " * midrow + "\n"
+
+ # remove trailing spaces
+ self.kernel = os.linesep.join([s.rstrip()
+ for s in self.kernel.splitlines()
+ if s])
+
+ def size(self):
+ return self.size
+
+ def center(self):
+ return self.center
+
+ def __str__(self):
+ return "Kernel:\n" + self.kernel
+
+
+class High_Pass_Filter:
+ """Based on a suitable Kernel string, this class creates a
+ filter suitable for GRASS-GIS' r.mfilter module.
+ Returns a *NIX ASCII multi-line string whose contents is
+ a matrix defining the way in which raster data will be filtered
+ by r.mfilter. The format of this file is described in r.mfilter's
+ manual."""
+ def __init__(self,
+ ratio,
+ level='Low',
+ modulation='Mid',
+ second_pass=False,
+ modulation2='Min',
+ divisor=1,
+ type='P'):
+
+ # parameters
+ self.ratio = ratio
+ self.size = kernel_size(self.ratio)
+
+ if second_pass:
+ self.modulator_2 = modulator(None, modulation2, self.size, True)
+ else:
+ self.modulator = modulator(modulation, None, self.size, False)
+
+ # build kernel
+ self.kernel = Kernel(self.size, level).kernel
+ self.header = 'MATRIX ' + str(self.size)
+ self.divisor = 'DIVISOR ' + str(divisor)
+ self.type = 'TYPE ' + str(type)
+ self.footer = str(self.divisor) + '\n' + self.type
+
+ # build filter
+ self.filter = ''
+ self.filter += self.header + '\n'
+ self.filter += self.kernel + '\n'
+ self.filter += self.footer
+
+ def set_divisor(self, divisor):
+ self.divisor = divisor
+
+ def set_type_(self, type):
+ self.type = type
+
+ def __str__(self):
+ return "Filter:\n" + self.filter
+
+# reusable & stand-alone
+if __name__ == "__main__":
+ print "Constructing a Filter for the HPFA Image Fusion Technique"
+ print " (Running as stand-alone tool)\n"
Property changes on: sandbox/alexandris/i.fusion.hpf/high_pass_filter.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.html
===================================================================
--- sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.html (rev 0)
+++ sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.html 2015-05-04 19:29:42 UTC (rev 65193)
@@ -0,0 +1,111 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.fusion.hpf</em> is an implementation of the High Pass Filter Additive
+(HPFA) Fusion Technique. It combines high-resolution panchromatic data with
+lower resolution multispectral data, resulting in an output with both excellent
+detail and a realistic representation of original multispectral scene colors.
+
+The process involves a convolution using a High Pass Filter (HPF) on the high
+resolution data, then combining this with the lower resolution multispectral
+data. Optionally, a linear histogram matching technique is performed in a way that
+matches the resulting Pan-Sharpened image to the statistical mean and standard
+deviation of the original multi-spectral image.<br>
+
+<h3>Background</h3>
+<ol type="1">
+ <li>Computing ratio of low (Multi-Spectral) to high (Panchromatic)
+ resolutions</li>
+ <li>High Pass Filtering the Panchromatic Image</li>
+ <li>Resampling MSX image to the higher resolution</li>
+ <li>Adding weighted High-Pass-Filetred image to the upsampled MSX
+ image</li>
+ <li>Optionally, matching histogram of Pansharpened image to the one of
+ the original MSX image</li>
+</ol>
+
+<pre>
+Figure:
+ ____________________________________________________________________________
++ +
+| Pan Img -> High Pass Filter -> HP Img |
+| | |
+| v |
+| MSx Img -> Weighting Factors -> Weighted HP Img |
+| | | |
+| | v |
+| +------------------------> Addition to MSx Img => Fused MSx Image |
++____________________________________________________________________________+
+
+</pre>
+Source: Gangkofner, 2008
+
+<h2>NOTES</h2>
+<ul>
+ <li> Grasping and testing the various parameters that define the High-Pass
+ filter's kernel size and center value is a matter of short time.</li>
+ <li> Works with any number and type of raster imagery (8-bit, 16-bit)</li>
+ <li> The "black border" effect, possibly caused due to a non-perfect match of the high vs. the low
+ resolution of the input images, can be trimmed out by using the <code>trim</code>
+ option --a floating point "trimming factor" with which to multiply the
+ pixel size of the low resolution image-- and shrink the extent of the
+ output image.</li>
+</ul>
+
+<h2>EXAMPLE</h2>
+The module is fairly easy to use. Arbitrary examples:
+<ul>
+ <li>for one band
+<div class="code"><pre>
+ i.fusion.hpf pan=Panchromatic msx=Red
+</pre></div></li>
+
+ <li>for multiple bands
+<div class="code"><pre>
+i.fusion.hpf pan=Panchromatic msx=Red,Green,Blue,NIR
+</pre></div></li>
+</ul>
+
+Various illustrated examples detailed in the document <a
+ href="https://github.com/NikosAlexandris/i.fusion.hpf/blob/master/Documentation.pdf">i.fusion.hpf, implementation of the High Pass Filter Additive (HPFA) Image Fusion Technique</a> (PDF)
+
+<h2>TODO</h2>
+<ul>
+ <li> Go through <a
+ href="http://trac.osgeo.org/grass/wiki/Submitting/Python">Submitting
+ Python</a></li>
+ <li> Access input raster by row I/O</li>
+ <li> Support for parallel processing</li>
+ <li> Proper command history tracking.</li>
+ <li> Add timestamps (r.timestamp, temporal framework)</li>
+ <li> Deduplicate code where applicable</li>
+ <li> Make verbose level messages shorter, yet more informative (ie report center cell)</li>
+ <li> Test if it compiles in other systems</li>
+ <li> Check options to integrate in <a href="i.pansharpen.html">i.pansharpen</a>. Think of FFM methods vs. Others?</li>
+ <li> Improve Documentation.lyx</li>
+</ul>
+
+<h2>REFERENCES</h2>
+<ul>
+ <li>Gangkofner, U. G., Pradhan, P. S., and Holcomb, D. W. (2008). Optimizing
+the high-pass filter addition technique for image fusion.
+PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING, 74(9):1107--1118.</li>
+
+ <li>"ERDAS IMAGINE." Accessed March 19, 2015. <a
+ href="http://doc.hexagongeospatial.com/ERDAS%20IMAGINE/ERDAS_IMAGINE_Help/#ii_hpfmerge_mergedialog.htm">ERDAS
+ IMAGINE Help</a>.</li>
+
+ <li>Aniruddha Ghosh & P.K. Joshi (2013) Assessment of pan-sharpened very
+ high-resolution WorldView-2 images, International Journal of Remote
+ Sensing, 34:23, 8336-8359</li>
+</ul>
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="i.pansharpen.html">i.pansharpen</a>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Nikos Alexandris<br>
+
Property changes on: sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.py
===================================================================
--- sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.py (rev 0)
+++ sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.py 2015-05-04 19:29:42 UTC (rev 65193)
@@ -0,0 +1,595 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+ MODULE: i.fusion.hpf
+
+ AUTHOR(S): Nikos Alexandris <nik at nikosalexandris.net>
+ Converted from a bash shell script | Trikala, Nov. 2014
+
+
+ PURPOSE: HPF Resolution Merge -- Algorithm Replication in GRASS GIS
+
+ Module to combine high-resolution panchromatic data with
+ lower resolution multispectral data, resulting in an output
+ with both excellent detail and a realistic representation of
+ original multispectral scene colors.
+
+ The process involves a convolution using a High Pass Filter
+ (HPF) on the high resolution data, then combining this with
+ the lower resolution multispectral data.
+
+ Optionally, a linear histogram matching technique is performed
+ in a way that matches the resulting Pan-Sharpened imaged to
+ them statistical mean and standard deviation of the original
+ multi-spectral image. Credits for how to implement this
+ technique go to GRASS-GIS developer Moritz Lennert.
+
+
+ Source: "Optimizing the High-Pass Filter Addition Technique for
+ Image Fusion", Ute G. Gangkofner, Pushkar S. Pradhan,
+ and Derrold W. Holcomb (2008)
+
+ Figure 1:
+
++-----------------------------------------------------------------------------+
+| Pan Img -> High Pass Filter -> HP Img |
+| | |
+| v |
+| MSx Img -> Weighting Factors -> Weighted HP Img |
+| | | |
+| | v |
+| +------------------------> Addition to MSx Img => Fused MSx Image |
++-----------------------------------------------------------------------------+
+
+ COPYRIGHT: (C) 2014 - 2015 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.
+"""
+
+#%Module
+#% description: Fusing high resolution Panchromatic and low resolution Multi-Spectral data based on the High-Pass Filter Addition technique (Gangkofner, 2008)
+#% keywords: imagery
+#% keywords: fusion
+#% keywords: sharpening
+#% keywords: high pass filter
+#% keywords: HPFA
+#%End
+
+#%flag
+#% key: l
+#% description: Linearly match histogram of Pan-sharpened output to Multi-Spectral input
+#%end
+
+#%flag
+#% key: 2
+#% description: 2-Pass Processing (recommended) for large resolution ratio (>=5.5)
+#%end
+
+#%flag
+#% key: c
+#% description: Match color table of Pan-Sharpened output to Multi-Spectral input
+#%end
+
+#%option G_OPT_R_INPUT
+#% key: pan
+#% key_desc: filename
+#% description: High resolution Panchromatic image
+#% required : yes
+#%end
+
+#%option G_OPT_R_INPUTS
+#% key: msx
+#% key_desc: filename(s)
+#% description: Low resolution Multi-Spectral image(s)
+#% required: yes
+#% multiple: yes
+#%end
+
+#%option G_OPT_R_BASENAME_OUTPUT
+#% key: suffix
+#% key_desc: suffix string
+#% type: string
+#% label: Suffix for output image(s)
+#% description: Names of Pan-Sharpened image(s) will end with this suffix
+#% required: yes
+#% answer: hpf
+#%end
+
+#%option
+#% key: ratio
+#% key_desc: rational number
+#% type: double
+#% label: Custom ratio
+#% description: Custom ratio overriding standard calculation
+#% options: 1.0-10.0
+#% guisection: High Pass Filter
+#% required: no
+#%end
+
+#%option
+#% key: center
+#% key_desc: string
+#% type: string
+#% label: Center cell value
+#% description: Center cell value of the High-Pass-Filter
+#% descriptions: Level of center value (low, mid, high)
+#% options: low,mid,high
+#% required: no
+#% answer: low
+#% guisection: High Pass Filter
+#% multiple : no
+#%end
+
+#%option
+#% key: center2
+#% key_desc: string
+#% type: string
+#% label: 2nd Pass center cell value
+#% description: Center cell value for the second High-Pass-Filter (use -2 flag)
+#% descriptions: Level of center value for second pass
+#% options: low,mid,high
+#% required: no
+#% answer: low
+#% guisection: High Pass Filter
+#% multiple : no
+#%end
+
+#%option
+#% key: modulation
+#% key_desc: string
+#% type: string
+#% label: Modulation level
+#% description: Modulation level weighting the HPF image determining crispness
+#% descriptions: Levels of modulating factors
+#% options: min,mid,max
+#% required: no
+#% answer: mid
+#% guisection: Crispness
+#% multiple : no
+#%end
+
+#%option
+#% key: modulation2
+#% key_desc: string
+#% type: string
+#% label: 2nd Pass modulation level (use -2 flag)
+#% description: Modulation level weighting the second HPF image determining crispness (use -2 flag)
+#% descriptions: mid;Mid: 0.35;min;Minimum: 0.25;max;Maximum: 0.5;
+#% options: min,mid,max
+#% required: no
+#% answer: mid
+#% guisection: Crispness
+#% multiple : no
+#%end
+
+#%option
+#% key: trim
+#% key_desc: rational number
+#% type: double
+#% label: Trimming factor
+#% description: Trim output border pixels by a factor of the pixel size of the low resolution image. A factor of 1.0 may suffice.
+#% guisection: High Pass Filter
+#% required: no
+#%end
+
+# required librairies
+
+import os
+import sys
+sys.path.insert(1, os.path.join(os.path.dirname(sys.path[0]),
+ 'etc', 'i.fusion.hpf'))
+import atexit
+
+import grass.script as grass
+from grass.pygrass.modules.shortcuts import general as g
+from grass.pygrass.raster.abstract import Info
+
+from high_pass_filter import High_Pass_Filter
+
+if "GISBASE" not in os.environ:
+ print "You must be in GRASS GIS to run this program."
+ sys.exit(1)
+
+
+# globals
+
+ratio = float()
+tmp = ''
+tmp_hpf_matrix = ''
+modulator = float()
+modulator_2 = float()
+
+
+# helper functions
+
+def cleanup():
+ """Clean up temporary maps"""
+ grass.run_command('g.remove', flags='f', type="raster",
+ pattern='tmp.{pid}*'.format(pid=os.getpid(), quiet=True))
+
+
+def run(cmd, **kwargs):
+ """Pass arbitrary number of key-word arguments to grass commands and the
+ "quiet" flag by default."""
+ grass.run_command(cmd, quiet=True, **kwargs)
+
+
+def avg(img):
+ """Retrieving Average of input image"""
+ uni = grass.parse_command("r.univar", map=img, flags='g')
+ avg = float(uni['mean'])
+ return avg
+
+
+def stddev(img):
+ """Retrieving Standard Deviation of input image"""
+ uni = grass.parse_command("r.univar", map=img, flags='g')
+ sd = float(uni['stddev'])
+ return sd
+
+
+def hpf_weight(low_sd, hpf_sd, mod, pss):
+ """Returning an appropriate weighting value for the
+ High Pass Filtered image. The required inputs are:
+ - low_sd: StdDev of Low resolution image
+ - hpf_sd: StdDev of High Pass Filtered image
+ - mod: Appropriate Modulating Factor determining image crispness
+ - pss: Number of Pass (1st or 2nd)"""
+ wgt = low_sd / hpf_sd * mod # mod: modulator
+ msg = ' >> '
+ if pss == 2:
+ msg += '2nd Pass '
+ msg += 'Weighting = {l:.{dec}f} / {h:.{dec}f} * {m:.{dec}f} = {w:.{dec}f}'
+ msg = msg.format(l=low_sd, h=hpf_sd, m=mod, w=wgt, dec=3)
+ g.message(msg, flags='v')
+ return wgt
+
+
+def hpf_ascii(center, filter, tmpfile, pss):
+ """Exporting a High Pass Filter in a temporary ASCII file"""
+ if pss == 1:
+ global modulator
+ modulator = filter.modulator
+ msg_pass2 = ''
+
+ elif pss == 2:
+ global modulator_2
+ modulator_2 = filter.modulator_2
+ msg_pass2 = '2nd Pass '
+
+ # structure informative message
+ msg = " > {m}Filter Properties: size: {f}, center: {c}"
+ msg = msg.format(m=msg_pass2, f=filter.size, c=center)
+ g.message(msg, flags='v')
+
+ # open, write and close file
+ asciif = open(tmpfile, 'w')
+ asciif.write(filter.filter)
+ asciif.close()
+
+
+# main program
+
+def main():
+
+ global tmp_hpf_matrix
+
+ pan = options['pan']
+ msxlst = options['msx'].split(',')
+ outputsuffix = options['suffix']
+ custom_ratio = options['ratio']
+ center = options['center']
+ center2 = options['center2']
+ modulation = options['modulation']
+ modulation2 = options['modulation2']
+
+ if options['trim']:
+ trimming_factor = float(options['trim'])
+ else:
+ trimming_factor = False
+
+ histogram_match = flags['l']
+ second_pass = flags['2']
+ color_match = flags['c']
+
+# # Check & warn user about "ns == ew" resolution of current region ======
+# region = grass.region()
+# nsr = region['nsres']
+# ewr = region['ewres']
+#
+# if nsr != ewr:
+# msg = ('>>> Region's North:South ({ns}) and East:West ({ew}) '
+# 'resolutions do not match!')
+# msg = msg.format(ns=nsr, ew=ewr)
+# g.message(msg, flags='w')
+
+ mapset = grass.gisenv()['MAPSET'] # Current Mapset?
+ region = grass.region() # and region settings
+
+ # List images and their properties
+
+ imglst = [pan]
+ imglst.extend(msxlst) # List of input imagery
+
+ images = {}
+ for img in imglst: # Retrieving Image Info
+ images[img] = Info(img, mapset)
+ images[img].read()
+
+ panres = images[pan].nsres # Panchromatic resolution
+
+ grass.use_temp_region() # to safely modify the region
+ run('g.region', res=panres) # Respect extent, change resolution
+ g.message("|! Region's resolution matched to Pan's ({p})".format(p=panres))
+
+ # Loop Algorithm over Multi-Spectral images
+
+ for msx in msxlst:
+
+ global tmp
+
+ # Inform
+ g.message("\nProcessing image: {m}".format(m=msx))
+
+ # Tracking command history -- Why don't do this all r.* modules?
+ cmd_history = ''
+
+ #
+ # 1. Compute Ratio
+ #
+
+ g.message("\n|1 Determining ratio of low to high resolution")
+
+ # Custom Ratio? Skip standard computation method.
+ if custom_ratio:
+ global ratio
+ ratio = float(custom_ratio)
+ g.message('Using custom ratio, overriding standard method!',
+ flags='w')
+
+ # Multi-Spectral resolution(s), multiple
+ else:
+ # Image resolutions
+ g.message(" > Retrieving image resolutions")
+
+ msxres = images[msx].nsres
+
+ # check
+ if panres == msxres:
+ msg = ("The Panchromatic's image resolution ({pr}) "
+ "equals to the Multi-Spectral's one ({mr}). "
+ "Something is probably not right! "
+ "Please check your input images.")
+ msg = msg.format(pr=panres, mr=msxres)
+ grass.fatal(_(msg))
+
+ # compute ratio
+ ratio = msxres / panres
+ msg_ratio = (' >> Resolution ratio '
+ 'low ({m:.{dec}f}) to high ({p:.{dec}f}): {r:.1f}')
+ msg_ratio = msg_ratio.format(m=msxres, p=panres, r=ratio, dec=3)
+ g.message(msg_ratio)
+
+ # 2nd Pass requested, yet Ratio < 5.5
+ if second_pass and ratio < 5.5:
+ g.message(" >>> Resolution ratio < 5.5, skipping 2nd pass.\n"
+ " >>> If you insist, force it via the <ratio> option!",
+ flags='i')
+ second_pass = bool(0)
+
+ #
+ # 2. High Pass Filtering
+ #
+
+ g.message('\n|2 High Pass Filtering the Panchromatic Image')
+
+ # Temporary files
+ tmpfile = grass.tempfile() # Temporary file - replace with os.getpid?
+ tmp = 'tmp.' + grass.basename(tmpfile) # use its basenam
+ tmp_pan_hpf = '{tmp}_pan_hpf'.format(tmp=tmp) # HPF image
+ tmp_msx_blnr = '{tmp}_msx_blnr'.format(tmp=tmp) # Upsampled MSx
+ tmp_msx_hpf = '{tmp}_msx_hpf'.format(tmp=tmp) # Fused image
+
+ tmp_hpf_matrix = grass.tempfile() # ASCII filter
+
+ if second_pass and ratio > 5.5: # 2nd Pass?
+ tmp_pan_hpf_2 = '{tmp}_pan_hpf_2'.format(tmp=tmp) # 2nd Pass HPF image
+ tmp_hpf_matrix_2 = grass.tempfile() # 2nd Pass ASCII filter
+
+ # Construct Filter
+ hpf = High_Pass_Filter(ratio, center, modulation, False, None)
+ hpf_ascii(center, hpf, tmp_hpf_matrix, 1)
+
+ # Construct 2nd Filter
+ if second_pass and ratio > 5.5:
+ hpf_2 = High_Pass_Filter(ratio, center2, None, True, modulation2)
+ hpf_ascii(center2, hpf_2, tmp_hpf_matrix_2, 2)
+
+ # Filtering
+ run('r.mfilter', input=pan, filter=tmp_hpf_matrix,
+ output=tmp_pan_hpf,
+ title='High Pass Filtered Panchromatic image',
+ overwrite=True)
+
+ # 2nd Filtering
+ if second_pass and ratio > 5.5:
+ run('r.mfilter', input=pan, filter=tmp_hpf_matrix_2,
+ output=tmp_pan_hpf_2,
+ title='2-High-Pass Filtered Panchromatic Image',
+ overwrite=True)
+
+ #
+ # 3. Upsampling low resolution image
+ #
+
+ g.message("\n|3 Upsampling (bilinearly) low resolution image")
+
+ run('r.resamp.interp',
+ method='bilinear', input=msx, output=tmp_msx_blnr, overwrite=True)
+
+ #
+ # 4. Weighting the High Pass Filtered image(s)
+ #
+
+ g.message("\n|4 Weighting the High-Pass-Filtered image (HPFi)")
+
+ # Compute (1st Pass) Weighting
+ msg_w = " > Weighting = StdDev(MSx) / StdDev(HPFi) * " \
+ "Modulating Factor"
+ g.message(msg_w)
+
+ # StdDev of Multi-Spectral Image(s)
+ msx_avg = avg(msx)
+ msx_sd = stddev(msx)
+ g.message(" >> StdDev of <{m}>: {sd:.3f}".format(m=msx, sd=msx_sd))
+
+ # StdDev of HPF Image
+ hpf_sd = stddev(tmp_pan_hpf)
+ g.message(" >> StdDev of HPFi: {sd:.3f}".format(sd=hpf_sd))
+
+ # Modulating factor
+ g.message(" >> Modulating Factor: {m:.2f}".format(m=modulator))
+
+ # weighting HPFi
+ weighting = hpf_weight(msx_sd, hpf_sd, modulator, 1)
+
+ #
+ # 5. Adding weighted HPF image to upsampled Multi-Spectral band
+ #
+
+ g.message("\n|5 Adding weighted HPFi to upsampled image")
+ fusion = '{hpf} = {msx} + {pan} * {wgt}'
+ fusion = fusion.format(hpf=tmp_msx_hpf, msx=tmp_msx_blnr,
+ pan=tmp_pan_hpf, wgt=weighting)
+ grass.mapcalc(fusion)
+
+ # command history
+ hst = 'Weigthing applied: {msd:.3f} / {hsd:.3f} * {mod:.3f} | '
+ cmd_history += hst.format(msd=msx_sd, hsd=hpf_sd, mod=modulator)
+
+ if second_pass and ratio > 5.5:
+
+ #
+ # 4+ 2nd Pass Weighting the High Pass Filtered image
+ #
+
+ g.message("\n|4+ 2nd Pass Weighting the HPFi")
+
+ # StdDev of HPF Image #2
+ hpf_2_sd = stddev(tmp_pan_hpf_2)
+ g.message(" >> StdDev of 2nd HPFi: {h:.3f}".format(h=hpf_2_sd))
+
+ # Modulating factor #2
+ msg = ' >> 2nd Pass Modulating Factor: {m:.2f}'
+ g.message(msg.format(m=modulator_2))
+
+ # 2nd Pass weighting
+ weighting_2 = hpf_weight(msx_sd, hpf_2_sd, modulator_2, 2)
+
+ #
+ # 5+ Adding weighted HPF image to upsampled Multi-Spectral band
+ #
+
+ g.message("\n|5+ Adding small-kernel-based weighted 2nd HPFi "
+ "back to fused image")
+
+ add_back = '{final} = {msx_hpf} + {pan_hpf} * {wgt}'
+ add_back = add_back.format(final=tmp_msx_hpf, msx_hpf=tmp_msx_hpf,
+ pan_hpf=tmp_pan_hpf_2, wgt=weighting_2)
+ grass.mapcalc(add_back)
+
+ # 2nd Pass history entry
+ hst = "2nd Pass Weighting: {m:.3f} / {h:.3f} * {mod:.3f} | "
+ cmd_history += hst.format(m=msx_sd, h=hpf_2_sd, mod=modulator_2)
+
+ if color_match:
+ g.message("\n|* Matching output to input color table")
+ run('r.colors',
+ map=tmp_msx_hpf, raster=msx)
+
+ #
+ # 6. Stretching linearly the HPF-Sharpened image(s) to match the Mean
+ # and Standard Deviation of the input Multi-Sectral image(s)
+ #
+
+ if histogram_match:
+
+ # adapt output StdDev and Mean to the input(ted) ones
+ g.message("\n|+ Matching histogram of Pansharpened image "
+ "to %s" % (msx), flags='v')
+
+ # Collect stats for linear histogram matching
+ msx_hpf_avg = avg(tmp_msx_hpf)
+ msx_hpf_sd = stddev(tmp_msx_hpf)
+
+ # expression for mapcalc
+ lhm = '{out} = ({hpf} - {hpfavg}) / {hpfsd} * {msxsd} + {msxavg}'
+ lhm = lhm.format(out=tmp_msx_hpf, hpf=tmp_msx_hpf,
+ hpfavg=msx_hpf_avg, hpfsd=msx_hpf_sd,
+ msxsd=msx_sd, msxavg=msx_avg)
+
+ # compute
+ grass.mapcalc(lhm, quiet=True, overwrite=True)
+
+ # update history string
+ cmd_history += "Linear Histogram Matching: %s |" % lhm
+
+ #
+ # Optional. Trim to remove black border effect (rectangular only)
+ #
+
+ if trimming_factor:
+
+ tf = trimming_factor
+
+ # communicate
+ msg = '\n|* Trimming output image border pixels by '
+ msg += '{factor} times the low resolution\n'.format(factor=tf)
+ nsew = ' > Input extent: n: {n}, s: {s}, e: {e}, w: {w}'
+ nsew = nsew.format(n=region.n, s=region.s, e=region.e, w=region.w)
+ msg += nsew
+
+ g.message(msg)
+
+ # re-set borders
+ region.n -= tf * images[msx].nsres
+ region.s += tf * images[msx].nsres
+ region.e -= tf * images[msx].ewres
+ region.w += tf * images[msx].ewres
+
+ # communicate and act
+ msg = ' > Output extent: n: {n}, s: {s}, e: {e}, w: {w}'
+ msg = msg.format(n=region.n, s=region.s, e=region.e, w=region.w)
+ g.message(msg)
+
+ # modify only the extent
+ run('g.region',
+ n=region.n, s=region.s, e=region.e, w=region.w)
+ trim = "{out} = {input}".format(out=tmp_msx_hpf, input=tmp_msx_hpf)
+ grass.mapcalc(trim)
+
+ #
+ # End of Algorithm
+
+ # history entry
+ run("r.support", map=tmp_msx_hpf, history=cmd_history)
+
+ # add suffix to basename & rename end product
+ msx_name = "{base}.{suffix}"
+ msx_name = msx_name.format(base=msx.split('@')[0], suffix=outputsuffix)
+ run("g.rename", raster=(tmp_msx_hpf, msx_name))
+
+ # visualising-related information
+ grass.del_temp_region() # restoring previous region settings
+ g.message("\n|! Original Region restored")
+ g.message("\n>>> Hint, rebalancing colors (via i.colors.enhance) "
+ "may improve appearance of RGB composites!",
+ flags='i')
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ sys.exit(main())
Property changes on: sandbox/alexandris/i.fusion.hpf/i.fusion.hpf.py
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list