[GRASS-SVN] r52624 - grass-addons/grass6/imagery/i.despeckle
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Aug 10 15:29:01 PDT 2012
Author: neteler
Date: 2012-08-10 15:29:01 -0700 (Fri, 10 Aug 2012)
New Revision: 52624
Modified:
grass-addons/grass6/imagery/i.despeckle/i.despeckle
Log:
Maning Sambale: more methods implemented
Modified: grass-addons/grass6/imagery/i.despeckle/i.despeckle
===================================================================
--- grass-addons/grass6/imagery/i.despeckle/i.despeckle 2012-08-10 18:33:57 UTC (rev 52623)
+++ grass-addons/grass6/imagery/i.despeckle/i.despeckle 2012-08-10 22:29:01 UTC (rev 52624)
@@ -3,12 +3,18 @@
############################################################################
#
# MODULE: i.despeckle
-# AUTHOR(S): Markus Neteler
+# AUTHOR(S): Markus Neteler, Maning Sambale
# PURPOSE: FGAMMA -- Gamma SAR Speckle Filter, based on
# http://web.archive.org/web/20050826071553/http://www.pcigeomatics.com/cgi-bin/pcihlp/M_FGAMMA
+# http://web.archive.org/web/20070122173456/http://www.pcigeomatics.com/cgi-bin/pcihlp/FLE|ALGORITHM
+# http://web.archive.org/web/20071216181051/http://www.pcigeomatics.com/cgi-bin/pcihlp/FELEE
+# http://web.archive.org/web/20071216081728/http://www.pcigeomatics.com/cgi-bin/pcihlp/FFROST
+# http://web.archive.org/web/20071216181046/http://www.pcigeomatics.com/cgi-bin/pcihlp/FEFROST
#
-# COPYRIGHT: (C) 2012 GRASS Development Team/Markus Neteler
+# List: http://web.archive.org/web/20060721084141/http://www.pcigeomatics.com/cgi-bin/searchhlp.pl
#
+# COPYRIGHT: (C) 2012 GRASS Development Team/Markus Neteler/Maning Sambale
+#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
@@ -21,7 +27,7 @@
#
############################################################################
#%Module
-#% description: Applies GAMMA SAR Speckle Filter to map.
+#% description: Applies SAR Speckle Filter to a raster power map.
#% keywords: imagery, filter, despeckle
#%End
#%Option
@@ -30,7 +36,7 @@
#% required: yes
#% multiple: no
#% key_desc: name
-#% description: Name of input raster map
+#% description: Name of input raster power map
#% gisprompt: old,cell,raster
#%End
#%Option
@@ -45,9 +51,9 @@
#%Option
#% key: method
#% type: string
-#% required: no
+#% required: yes
#% multiple: no
-#% options: gamma
+#% options: gamma,lee,kuan
#% description: Despeckle operation
#% answer: gamma
#%End
@@ -84,53 +90,105 @@
fi
-
MAP=$GIS_OPT_INPUT
WIN=$GIS_OPT_SIZE
-NumberLooks=$GIS_OPT_LOOKS
+NumberLooks=$GIS_OPT_LOOKS # controls the amount of smoothing applied to image
OUTPUT=$GIS_OPT_OUTPUT
+METHOD=$GIS_OPT_METHOD
-######################
-# S = standard deviation of intensity within window
+#########################
+
+# S = standard deviation of intensity within window
r.neighbors in=$MAP out=$MAP.stddev method=stddev size=$WIN
-# Im = mean value of intensity within window
+# Im = mean value of intensity within window
r.neighbors in=$MAP out=$MAP.mean method=average size=$WIN
-# Ci = S / Im
-r.mapcalc "$MAP.ci = $MAP.stddev / $MAP.mean"
+# Ci = S / Im
+r.mapcalc "$MAP.ci = float($MAP.stddev) / float($MAP.mean)"
-# Cu = SQRT(1/NumberLooks)
-Cu="`echo $NumberLooks | awk '{printf "%f", sqrt(1/$1)}'`"
+# Cu = SQRT(1/NumberLooks)
+Cu="`echo $NumberLooks | awk '{printf "%f", sqrt(1.0/$1)}'`"
-# Cmax = SQRT(2)*Cu
-Cmax="`echo $Cu | awk '{printf "%f", sqrt(2) * $1}'`"
+######################
-# A = (1+Cu^2)/(Ci^2-Cu^2)
-r.mapcalc "$MAP.A = (1 + $Cu^2)/($MAP.ci^2 - $Cu^2)"
+if [ $METHOD == "lee" ]
+then
+ # A 7x7 filter usually gives the best results
+ #
+ # Note: Lee filter model requires that the signal represents power.
+ # If the input image is in amplitude format, each grey level must be
+ # squared to derive power and finally square root must be applied to
+ # the filtered result.
+ echo "Processing with LEE filter"
-# B = A-NumberLooks-1
-r.mapcalc "$MAP.B = $MAP.A - $NumberLooks - 1"
+ # W = 1 - Cu^2/Ci^2
+ r.mapcalc "$MAP.W = 1.0 - ($Cu^2/$MAP.ci^2)"
-# D = Im*Im*B*B + 4*A*NumberLooks*Im*Ic
-# Ic = center pixel in filter window
-r.mapcalc "$MAP.D = $MAP.mean * $MAP.mean * $MAP.B * $MAP.B + 4 * $MAP.A * $NumberLooks * $MAP.mean * $MAP"
+ # R = Ic * W + Im * (1 - W)
+ # Ic = center pixel in filter window
+ r.mapcalc "$OUTPUT = $MAP * $MAP.W + $MAP.mean * (1.0 - $MAP.W)"
-# Rf = (B*Im + SQRT(D))/(2*A)
-r.mapcalc "$MAP.rf = ($MAP.B * $MAP.mean + sqrt($MAP.D)) / (2 * $MAP.A)"
+ g.remove --q rast=$MAP.mean,$MAP.stddev,$MAP.rf,$MAP.ci,$MAP.W
+fi
-# R1 = Im for Ci <= Cu
-# R2 = Rf for Cu < Ci < Cmax
-# R3 = Ic for Ci >= Cmax
+if [ $METHOD == "kuan" ]
+then
+ echo "Processing with KUAN filter"
+ # A 7x7 filter usually gives the best results
+ #
+ # Note: KUAN filter model requires that the signal represents power.
-r.mapcalc "$OUTPUT = eval(r1 = $MAP.mean, r2 = $MAP.rf, r3 = $MAP, \\
- if($MAP.ci < $Cu, r1, \\
- if($MAP.ci > $Cu && $MAP.ci < $Cmax, r2, \\
- if($MAP.ci > $Cmax, r3, null() ) \\
- ) \\
- ))"
+ # W = (1 - Cu^2/Ci^2)/(1+Cu^2)
+ r.mapcalc "$MAP.W = (1.0 - $Cu^2/$MAP.ci^2) / (1.0 + $Cu^2)"
+ # R = Ic * W + Im * (1 - W)
+ # Ic = center pixel in filter window
+ r.mapcalc "$OUTPUT = $MAP * $MAP.W + $MAP.mean * (1.0 - $MAP.W)"
+
+ g.remove --q rast=$MAP.mean,$MAP.stddev,$MAP.ci,$MAP.W
+
+fi
+
+if [ $METHOD == "gamma" ]
+then
+ echo "Processing with GAMMA filter"
+ # A 7x7 filter usually gives the best results
+ #
+ # Note: GAMMA filter model requires that the signal represents power.
+
+ # Cmax = SQRT(2)*Cu
+ Cmax="`echo $Cu | awk '{printf "%f", sqrt(2.0) * $1}'`"
+
+ # A = (1+Cu^2)/(Ci^2-Cu^2)
+ r.mapcalc "$MAP.A = (1.0 + $Cu^2)/($MAP.ci^2 - $Cu^2)"
+
+ # B = A-NumberLooks-1
+ r.mapcalc "$MAP.B = $MAP.A - $NumberLooks - 1"
+
+ # D = Im*Im*B*B + 4*A*NumberLooks*Im*Ic
+ # Ic = center pixel in filter window
+ r.mapcalc "$MAP.D = $MAP.mean * $MAP.mean * $MAP.B * $MAP.B + 4 * $MAP.A * \\
+ $NumberLooks * $MAP.mean * $MAP"
+
+ # Rf = (B*Im + SQRT(D))/(2*A)
+ r.mapcalc "$MAP.rf = ($MAP.B * $MAP.mean + sqrt($MAP.D)) / (2 * $MAP.A)"
+
+ # R1 = Im for Ci <= Cu
+ # R2 = Rf for Cu < Ci < Cmax
+ # R3 = Ic for Ci >= Cmax
+ r.mapcalc "$OUTPUT = eval(r1 = $MAP.mean, r2 = $MAP.rf, r3 = $MAP, \\
+ if($MAP.ci < $Cu, r1, \\
+ if($MAP.ci > $Cu && $MAP.ci < $Cmax, r2, \\
+ if($MAP.ci > $Cmax, r3, null() ) \\
+ ) \\
+ ))"
+
+ g.remove --q rast=$MAP.mean,$MAP.stddev,$MAP.rf,$MAP.ci,$MAP.A,$MAP.B,$MAP.D,$MAP.W
+
+fi
+####################
+
r.colors $OUTPUT col=grey -e
-
-g.remove --q rast=$MAP.mean,$MAP.stddev,$MAP.rf,$MAP.ci,$MAP.A,$MAP.B,$MAP.D
+r.support map=$OUTPUT description="generated by i.descpeckle using $METHOD filter."
echo "Written <$OUTPUT>"
More information about the grass-commit
mailing list