| 1 | #!/usr/bin/env python
|
|---|
| 2 |
|
|---|
| 3 | ############################################################################
|
|---|
| 4 | #
|
|---|
| 5 | # MODULE: i.pansharpen
|
|---|
| 6 | #
|
|---|
| 7 | # AUTHOR(S): Overall script by Michael Barton (ASU)
|
|---|
| 8 | # Brovey transformation in i.fusion.brovey by Markus Neteler <<neteler at osgeo org>>
|
|---|
| 9 | # i.fusion brovey converted to Python by Glynn Clements
|
|---|
| 10 | # IHS and PCA transformation added by Michael Barton (ASU)
|
|---|
| 11 | # histogram matching algorithm by Michael Barton and Luca Delucchi, Fondazione E. Mach (Italy)
|
|---|
| 12 | # Thanks to Markus Metz for help with PCA inversion
|
|---|
| 13 | # Thanks to Hamish Bowman for parallel processing algorithm
|
|---|
| 14 | #
|
|---|
| 15 | # PURPOSE: Sharpening of 3 RGB channels using a high-resolution panchromatic channel
|
|---|
| 16 | #
|
|---|
| 17 | # COPYRIGHT: (C) 2002-2019 by the GRASS Development Team
|
|---|
| 18 | #
|
|---|
| 19 | # This program is free software under the GNU General Public
|
|---|
| 20 | # License (>=v2). Read the file COPYING that comes with GRASS
|
|---|
| 21 | # for details.
|
|---|
| 22 | #
|
|---|
| 23 | # REFERENCES:
|
|---|
| 24 | # Roller, N.E.G. and Cox, S., 1980. Comparison of Landsat MSS and merged MSS/RBV
|
|---|
| 25 | # data for analysis of natural vegetation. Proc. of the 14th International
|
|---|
| 26 | # Symposium on Remote Sensing of Environment, San Jose, Costa Rica, 23-30 April, pp. 1001-1007.
|
|---|
| 27 | #
|
|---|
| 28 | # Amarsaikhan, D., & Douglas, T. (2004). Data fusion and multisource image classification.
|
|---|
| 29 | # International Journal of Remote Sensing, 25(17), 3529-3539.
|
|---|
| 30 | #
|
|---|
| 31 | # Behnia, P. (2005). Comparison between four methods for data fusion of ETM+
|
|---|
| 32 | # multispectral and pan images. Geo-spatial Information Science, 8(2), 98-103
|
|---|
| 33 | #
|
|---|
| 34 | # for LANDSAT 5: see Pohl, C 1996 and others
|
|---|
| 35 | #
|
|---|
| 36 | #############################################################################
|
|---|
| 37 |
|
|---|
| 38 | #%Module
|
|---|
| 39 | #% description: Image fusion algorithms to sharpen multispectral with high-res panchromatic channels
|
|---|
| 40 | #% keyword: imagery
|
|---|
| 41 | #% keyword: fusion
|
|---|
| 42 | #% keyword: sharpen
|
|---|
| 43 | #% keyword: Brovey
|
|---|
| 44 | #% keyword: IHS
|
|---|
| 45 | #% keyword: HIS
|
|---|
| 46 | #% keyword: PCA
|
|---|
| 47 | #% overwrite: yes
|
|---|
| 48 | #%End
|
|---|
| 49 | #%option G_OPT_R_INPUT
|
|---|
| 50 | #% key: red
|
|---|
| 51 | #% description: Name of raster map to be used for <red>
|
|---|
| 52 | #%end
|
|---|
| 53 | #%option G_OPT_R_INPUT
|
|---|
| 54 | #% key: green
|
|---|
| 55 | #% description: Name of raster map to be used for <green>
|
|---|
| 56 | #%end
|
|---|
| 57 | #%option G_OPT_R_INPUT
|
|---|
| 58 | #% key: blue
|
|---|
| 59 | #% description: Name of raster map to be used for <blue>
|
|---|
| 60 | #%end
|
|---|
| 61 | #% option G_OPT_R_INPUT
|
|---|
| 62 | #% key: pan
|
|---|
| 63 | #% description: Name of raster map to be used for high resolution panchromatic channel
|
|---|
| 64 | #%end
|
|---|
| 65 | #%option G_OPT_R_BASENAME_OUTPUT
|
|---|
| 66 | #%end
|
|---|
| 67 | #%option
|
|---|
| 68 | #% key: method
|
|---|
| 69 | #% description: Method for pan sharpening
|
|---|
| 70 | #% options: brovey,ihs,pca
|
|---|
| 71 | #% answer: ihs
|
|---|
| 72 | #% required: yes
|
|---|
| 73 | #%end
|
|---|
| 74 | #%option
|
|---|
| 75 | #% key: bitdepth
|
|---|
| 76 | #% type: integer
|
|---|
| 77 | #% description: Bit depth of image (must be in range of 2-30)
|
|---|
| 78 | #% options: 2-32
|
|---|
| 79 | #% answer: 8
|
|---|
| 80 | #% required: yes
|
|---|
| 81 | #%end
|
|---|
| 82 | #%flag
|
|---|
| 83 | #% key: s
|
|---|
| 84 | #% description: Serial processing rather than parallel processing
|
|---|
| 85 | #%end
|
|---|
| 86 | #%flag
|
|---|
| 87 | #% key: l
|
|---|
| 88 | #% description: Rebalance blue channel for LANDSAT
|
|---|
| 89 | #%end
|
|---|
| 90 | #%flag
|
|---|
| 91 | #% key: r
|
|---|
| 92 | #% description: Rescale (stretch) the range of pixel values in each channel to the entire 0-255 8-bit range for processing (see notes)
|
|---|
| 93 | #%end
|
|---|
| 94 |
|
|---|
| 95 | import os
|
|---|
| 96 |
|
|---|
| 97 | try:
|
|---|
| 98 | import numpy as np
|
|---|
| 99 | hasNumPy = True
|
|---|
| 100 | except ImportError:
|
|---|
| 101 | hasNumPy = False
|
|---|
| 102 |
|
|---|
| 103 | import grass.script as grass
|
|---|
| 104 |
|
|---|
| 105 |
|
|---|
| 106 | def main():
|
|---|
| 107 | if not hasNumPy:
|
|---|
| 108 | grass.fatal(_("Required dependency NumPy not found. Exiting."))
|
|---|
| 109 |
|
|---|
| 110 | sharpen = options['method'] # sharpening algorithm
|
|---|
| 111 | ms1_orig = options['blue'] # blue channel
|
|---|
| 112 | ms2_orig = options['green'] # green channel
|
|---|
| 113 | ms3_orig = options['red'] # red channel
|
|---|
| 114 | pan_orig = options['pan'] # high res pan channel
|
|---|
| 115 | out = options['output'] # prefix for output RGB maps
|
|---|
| 116 | bits = options['bitdepth'] # bit depth of image channels
|
|---|
| 117 | bladjust = flags['l'] # adjust blue channel
|
|---|
| 118 | sproc = flags['s'] # serial processing
|
|---|
| 119 | rescale = flags['r'] # rescale to spread pixel values to entire 0-255 range
|
|---|
| 120 |
|
|---|
| 121 | # Checking bit depth
|
|---|
| 122 | bits = float(bits)
|
|---|
| 123 | if bits < 2 or bits > 30:
|
|---|
| 124 | grass.warning(_("Bit depth is outside acceptable range"))
|
|---|
| 125 | return
|
|---|
| 126 |
|
|---|
| 127 | outb = grass.core.find_file('%s_blue' % out)
|
|---|
| 128 | outg = grass.core.find_file('%s_green' % out)
|
|---|
| 129 | outr = grass.core.find_file('%s_red' % out)
|
|---|
| 130 |
|
|---|
| 131 | if (outb['name'] != '' or outg['name'] != '' or outr['name'] != '') and not grass.overwrite():
|
|---|
| 132 | grass.warning(_('Maps with selected output prefix names already exist.'
|
|---|
| 133 | ' Delete them or use overwrite flag'))
|
|---|
| 134 | return
|
|---|
| 135 |
|
|---|
| 136 | pid = str(os.getpid())
|
|---|
| 137 |
|
|---|
| 138 | # convert input image channels to 8 bit for processing
|
|---|
| 139 | ms1 = 'tmp%s_ms1' % pid
|
|---|
| 140 | ms2 = 'tmp%s_ms2' % pid
|
|---|
| 141 | ms3 = 'tmp%s_ms3' % pid
|
|---|
| 142 | pan = 'tmp%s_pan' % pid
|
|---|
| 143 |
|
|---|
| 144 | if rescale == False:
|
|---|
| 145 | if bits == 8:
|
|---|
| 146 | grass.message(_("Using 8bit image channels"))
|
|---|
| 147 | if sproc:
|
|---|
| 148 | # serial processing
|
|---|
| 149 | grass.run_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
|
|---|
| 150 | quiet=True, overwrite=True)
|
|---|
| 151 | grass.run_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
|
|---|
| 152 | quiet=True, overwrite=True)
|
|---|
| 153 | grass.run_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
|
|---|
| 154 | quiet=True, overwrite=True)
|
|---|
| 155 | grass.run_command('g.copy', raster='%s,%s' % (pan_orig, pan),
|
|---|
| 156 | quiet=True, overwrite=True)
|
|---|
| 157 | else:
|
|---|
| 158 | # parallel processing
|
|---|
| 159 | pb = grass.start_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
|
|---|
| 160 | quiet=True, overwrite=True)
|
|---|
| 161 | pg = grass.start_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
|
|---|
| 162 | quiet=True, overwrite=True)
|
|---|
| 163 | pr = grass.start_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
|
|---|
| 164 | quiet=True, overwrite=True)
|
|---|
| 165 | pp = grass.start_command('g.copy', raster='%s,%s' % (pan_orig, pan),
|
|---|
| 166 | quiet=True, overwrite=True)
|
|---|
| 167 |
|
|---|
| 168 | pb.wait()
|
|---|
| 169 | pg.wait()
|
|---|
| 170 | pr.wait()
|
|---|
| 171 | pp.wait()
|
|---|
| 172 |
|
|---|
| 173 | else:
|
|---|
| 174 | grass.message(_("Converting image chanels to 8bit for processing"))
|
|---|
| 175 | maxval = pow(2, bits) - 1
|
|---|
| 176 | if sproc:
|
|---|
| 177 | # serial processing
|
|---|
| 178 | grass.run_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
|
|---|
| 179 | output=ms1, to='0,255', quiet=True, overwrite=True)
|
|---|
| 180 | grass.run_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
|
|---|
| 181 | output=ms2, to='0,255', quiet=True, overwrite=True)
|
|---|
| 182 | grass.run_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
|
|---|
| 183 | output=ms3, to='0,255', quiet=True, overwrite=True)
|
|---|
| 184 | grass.run_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
|
|---|
| 185 | output=pan, to='0,255', quiet=True, overwrite=True)
|
|---|
| 186 |
|
|---|
| 187 | else:
|
|---|
| 188 | # parallel processing
|
|---|
| 189 | pb = grass.start_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
|
|---|
| 190 | output=ms1, to='0,255', quiet=True, overwrite=True)
|
|---|
| 191 | pg = grass.start_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
|
|---|
| 192 | output=ms2, to='0,255', quiet=True, overwrite=True)
|
|---|
| 193 | pr = grass.start_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
|
|---|
| 194 | output=ms3, to='0,255', quiet=True, overwrite=True)
|
|---|
| 195 | pp = grass.start_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
|
|---|
| 196 | output=pan, to='0,255', quiet=True, overwrite=True)
|
|---|
| 197 |
|
|---|
| 198 | pb.wait()
|
|---|
| 199 | pg.wait()
|
|---|
| 200 | pr.wait()
|
|---|
| 201 | pp.wait()
|
|---|
| 202 |
|
|---|
| 203 | else:
|
|---|
| 204 | grass.message(_("Rescaling image chanels to 8bit for processing"))
|
|---|
| 205 |
|
|---|
| 206 | min_ms1 = int(grass.raster_info(ms1_orig)['min'])
|
|---|
| 207 | max_ms1 = int(grass.raster_info(ms1_orig)['max'])
|
|---|
| 208 | min_ms2 = int(grass.raster_info(ms2_orig)['min'])
|
|---|
| 209 | max_ms2 = int(grass.raster_info(ms2_orig)['max'])
|
|---|
| 210 | min_ms3 = int(grass.raster_info(ms3_orig)['min'])
|
|---|
| 211 | max_ms3 = int(grass.raster_info(ms3_orig)['max'])
|
|---|
| 212 | min_pan = int(grass.raster_info(pan_orig)['min'])
|
|---|
| 213 | max_pan = int(grass.raster_info(pan_orig)['max'])
|
|---|
| 214 |
|
|---|
| 215 | maxval = pow(2, bits) - 1
|
|---|
| 216 | if sproc:
|
|---|
| 217 | # serial processing
|
|---|
| 218 | grass.run_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1),
|
|---|
| 219 | output=ms1, to='0,255', quiet=True, overwrite=True)
|
|---|
| 220 | grass.run_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2),
|
|---|
| 221 | output=ms2, to='0,255', quiet=True, overwrite=True)
|
|---|
| 222 | grass.run_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3),
|
|---|
| 223 | output=ms3, to='0,255', quiet=True, overwrite=True)
|
|---|
| 224 | grass.run_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan),
|
|---|
| 225 | output=pan, to='0,255', quiet=True, overwrite=True)
|
|---|
| 226 |
|
|---|
| 227 | else:
|
|---|
| 228 | # parallel processing
|
|---|
| 229 | pb = grass.start_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1),
|
|---|
| 230 | output=ms1, to='0,255', quiet=True, overwrite=True)
|
|---|
| 231 | pg = grass.start_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2),
|
|---|
| 232 | output=ms2, to='0,255', quiet=True, overwrite=True)
|
|---|
| 233 | pr = grass.start_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3),
|
|---|
| 234 | output=ms3, to='0,255', quiet=True, overwrite=True)
|
|---|
| 235 | pp = grass.start_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan),
|
|---|
| 236 | output=pan, to='0,255', quiet=True, overwrite=True)
|
|---|
| 237 |
|
|---|
| 238 | pb.wait()
|
|---|
| 239 | pg.wait()
|
|---|
| 240 | pr.wait()
|
|---|
| 241 | pp.wait()
|
|---|
| 242 |
|
|---|
| 243 |
|
|---|
| 244 | # get PAN resolution:
|
|---|
| 245 | kv = grass.raster_info(map=pan)
|
|---|
| 246 | nsres = kv['nsres']
|
|---|
| 247 | ewres = kv['ewres']
|
|---|
| 248 | panres = (nsres + ewres) / 2
|
|---|
| 249 |
|
|---|
| 250 | # clone current region
|
|---|
| 251 | grass.use_temp_region()
|
|---|
| 252 | grass.run_command('g.region', res=panres, align=pan)
|
|---|
| 253 |
|
|---|
| 254 | # Select sharpening method
|
|---|
| 255 | grass.message(_("Performing pan sharpening with hi res pan image: %f" % panres))
|
|---|
| 256 | if sharpen == "brovey":
|
|---|
| 257 | brovey(pan, ms1, ms2, ms3, out, pid, sproc)
|
|---|
| 258 | elif sharpen == "ihs":
|
|---|
| 259 | ihs(pan, ms1, ms2, ms3, out, pid, sproc)
|
|---|
| 260 | elif sharpen == "pca":
|
|---|
| 261 | pca(pan, ms1, ms2, ms3, out, pid, sproc)
|
|---|
| 262 | # Could add other sharpening algorithms here, e.g. wavelet transformation
|
|---|
| 263 |
|
|---|
| 264 | grass.message(_("Assigning grey equalized color tables to output images..."))
|
|---|
| 265 |
|
|---|
| 266 | # equalized grey scales give best contrast
|
|---|
| 267 | grass.message(_("setting pan-sharpened channels to equalized grey scale"))
|
|---|
| 268 | for ch in ['red', 'green', 'blue']:
|
|---|
| 269 | grass.run_command('r.colors', quiet=True, map="%s_%s" % (out, ch),
|
|---|
| 270 | flags="e", color='grey')
|
|---|
| 271 |
|
|---|
| 272 | # Landsat too blue-ish because panchromatic band less sensitive to blue
|
|---|
| 273 | # light, so output blue channed can be modified
|
|---|
| 274 | if bladjust:
|
|---|
| 275 | grass.message(_("Adjusting blue channel color table..."))
|
|---|
| 276 | blue_colors = ['0 0 0 0\n5% 0 0 0\n67% 255 255 255\n100% 255 255 255']
|
|---|
| 277 | # these previous colors are way too blue for landsat
|
|---|
| 278 | # 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']
|
|---|
| 279 | bc = grass.feed_command('r.colors', quiet = True, map = "%s_blue" % out, rules = "-")
|
|---|
| 280 | bc.stdin.write('\n'.join(blue_colors))
|
|---|
| 281 | bc.stdin.close()
|
|---|
| 282 |
|
|---|
| 283 | # output notice
|
|---|
| 284 | grass.verbose(_("The following pan-sharpened output maps have been generated:"))
|
|---|
| 285 | for ch in ['red', 'green', 'blue']:
|
|---|
| 286 | grass.verbose(_("%s_%s") % (out, ch))
|
|---|
| 287 |
|
|---|
| 288 | grass.verbose(_("To visualize output, run: g.region -p raster=%s_red" % out))
|
|---|
| 289 | grass.verbose(_("d.rgb r=%s_red g=%s_green b=%s_blue" % (out, out, out)))
|
|---|
| 290 | grass.verbose(_("If desired, combine channels into a single RGB map with 'r.composite'."))
|
|---|
| 291 | grass.verbose(_("Channel colors can be rebalanced using i.colors.enhance."))
|
|---|
| 292 |
|
|---|
| 293 | # write cmd history:
|
|---|
| 294 | for ch in ['red', 'green', 'blue']:
|
|---|
| 295 | grass.raster_history("%s_%s" % (out, ch))
|
|---|
| 296 |
|
|---|
| 297 | # create a group with the three outputs
|
|---|
| 298 | #grass.run_command('i.group', group=out,
|
|---|
| 299 | # input="{n}_red,{n}_blue,{n}_green".format(n=out))
|
|---|
| 300 |
|
|---|
| 301 | # Cleanup
|
|---|
| 302 | grass.message(_("cleaning up temp files"))
|
|---|
| 303 | try:
|
|---|
| 304 | grass.run_command('g.remove', flags="f", type="raster",
|
|---|
| 305 | pattern="tmp%s*" % pid, quiet=True)
|
|---|
| 306 | except:
|
|---|
| 307 | ""
|
|---|
| 308 |
|
|---|
| 309 | def brovey(pan, ms1, ms2, ms3, out, pid, sproc):
|
|---|
| 310 | grass.verbose(_("Using Brovey algorithm"))
|
|---|
| 311 |
|
|---|
| 312 | # pan/intensity histogram matching using linear regression
|
|---|
| 313 | grass.message(_("Pan channel/intensity histogram matching using linear regression"))
|
|---|
| 314 | outname = 'tmp%s_pan1' % pid
|
|---|
| 315 | panmatch1 = matchhist(pan, ms1, outname)
|
|---|
| 316 |
|
|---|
| 317 | outname = 'tmp%s_pan2' % pid
|
|---|
| 318 | panmatch2 = matchhist(pan, ms2, outname)
|
|---|
| 319 |
|
|---|
| 320 | outname = 'tmp%s_pan3' % pid
|
|---|
| 321 | panmatch3 = matchhist(pan, ms3, outname)
|
|---|
| 322 |
|
|---|
| 323 | outr = '%s_red' % out
|
|---|
| 324 | outg = '%s_green' % out
|
|---|
| 325 | outb = '%s_blue' % out
|
|---|
| 326 |
|
|---|
| 327 | # calculate brovey transformation
|
|---|
| 328 | grass.message(_("Calculating Brovey transformation..."))
|
|---|
| 329 |
|
|---|
| 330 | if sproc:
|
|---|
| 331 | # serial processing
|
|---|
| 332 | e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
|
|---|
| 333 | "$outr" = 1 * round("$ms3" * "$panmatch3" / k)
|
|---|
| 334 | "$outg" = 1 * round("$ms2" * "$panmatch2" / k)
|
|---|
| 335 | "$outb" = 1 * round("$ms1" * "$panmatch1" / k)'''
|
|---|
| 336 | grass.mapcalc(e, outr=outr, outg=outg, outb=outb,
|
|---|
| 337 | panmatch1=panmatch1, panmatch2=panmatch2,
|
|---|
| 338 | panmatch3=panmatch3, ms1=ms1, ms2=ms2, ms3=ms3,
|
|---|
| 339 | overwrite=True)
|
|---|
| 340 | else:
|
|---|
| 341 | # parallel processing
|
|---|
| 342 | pb = grass.mapcalc_start('%s_blue = 1 * round((%s * %s) / (%s + %s + %s))' %
|
|---|
| 343 | (out, ms1, panmatch1, ms1, ms2, ms3),
|
|---|
| 344 | overwrite=True)
|
|---|
| 345 | pg = grass.mapcalc_start('%s_green = 1 * round((%s * %s) / (%s + %s + %s))' %
|
|---|
| 346 | (out, ms2, panmatch2, ms1, ms2, ms3),
|
|---|
| 347 | overwrite=True)
|
|---|
| 348 | pr = grass.mapcalc_start('%s_red = 1 * round((%s * %s) / (%s + %s + %s))' %
|
|---|
| 349 | (out, ms3, panmatch3, ms1, ms2, ms3),
|
|---|
| 350 | overwrite=True)
|
|---|
| 351 |
|
|---|
| 352 | pb.wait(), pg.wait(), pr.wait()
|
|---|
| 353 | try:
|
|---|
| 354 | pb.terminate(), pg.terminate(), pr.terminate()
|
|---|
| 355 | except:
|
|---|
| 356 | ""
|
|---|
| 357 |
|
|---|
| 358 | # Cleanup
|
|---|
| 359 | try:
|
|---|
| 360 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 361 | name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
|
|---|
| 362 | except:
|
|---|
| 363 | ""
|
|---|
| 364 |
|
|---|
| 365 | def ihs(pan, ms1, ms2, ms3, out, pid, sproc):
|
|---|
| 366 | grass.verbose(_("Using IHS<->RGB algorithm"))
|
|---|
| 367 | # transform RGB channels into IHS color space
|
|---|
| 368 | grass.message(_("Transforming to IHS color space..."))
|
|---|
| 369 | grass.run_command('i.rgb.his', overwrite=True,
|
|---|
| 370 | red=ms3,
|
|---|
| 371 | green=ms2,
|
|---|
| 372 | blue=ms1,
|
|---|
| 373 | hue="tmp%s_hue" % pid,
|
|---|
| 374 | intensity="tmp%s_int" % pid,
|
|---|
| 375 | saturation="tmp%s_sat" % pid)
|
|---|
| 376 |
|
|---|
| 377 | # pan/intensity histogram matching using linear regression
|
|---|
| 378 | target = "tmp%s_int" % pid
|
|---|
| 379 | outname = "tmp%s_pan_int" % pid
|
|---|
| 380 | panmatch = matchhist(pan, target, outname)
|
|---|
| 381 |
|
|---|
| 382 | # substitute pan for intensity channel and transform back to RGB color space
|
|---|
| 383 | grass.message(_("Transforming back to RGB color space and sharpening..."))
|
|---|
| 384 | grass.run_command('i.his.rgb', overwrite=True,
|
|---|
| 385 | hue="tmp%s_hue" % pid,
|
|---|
| 386 | intensity="%s" % panmatch,
|
|---|
| 387 | saturation="tmp%s_sat" % pid,
|
|---|
| 388 | red="%s_red" % out,
|
|---|
| 389 | green="%s_green" % out,
|
|---|
| 390 | blue="%s_blue" % out)
|
|---|
| 391 |
|
|---|
| 392 | # Cleanup
|
|---|
| 393 | try:
|
|---|
| 394 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 395 | name=panmatch)
|
|---|
| 396 | except:
|
|---|
| 397 | ""
|
|---|
| 398 |
|
|---|
| 399 | def pca(pan, ms1, ms2, ms3, out, pid, sproc):
|
|---|
| 400 |
|
|---|
| 401 | grass.verbose(_("Using PCA/inverse PCA algorithm"))
|
|---|
| 402 | grass.message(_("Creating PCA images and calculating eigenvectors..."))
|
|---|
| 403 |
|
|---|
| 404 | # initial PCA with RGB channels
|
|---|
| 405 | pca_out = grass.read_command('i.pca', quiet=True, rescale='0,0',
|
|---|
| 406 | input='%s,%s,%s' % (ms1, ms2, ms3),
|
|---|
| 407 | output='tmp%s.pca' % pid)
|
|---|
| 408 | if len(pca_out) < 1:
|
|---|
| 409 | grass.fatal(_("Input has no data. Check region settings."))
|
|---|
| 410 |
|
|---|
| 411 | b1evect = []
|
|---|
| 412 | b2evect = []
|
|---|
| 413 | b3evect = []
|
|---|
| 414 | for l in pca_out.replace('(', ',').replace(')', ',').splitlines():
|
|---|
| 415 | b1evect.append(float(l.split(',')[1]))
|
|---|
| 416 | b2evect.append(float(l.split(',')[2]))
|
|---|
| 417 | b3evect.append(float(l.split(',')[3]))
|
|---|
| 418 |
|
|---|
| 419 | # inverse PCA with hi res pan channel substituted for principal component 1
|
|---|
| 420 | pca1 = 'tmp%s.pca.1' % pid
|
|---|
| 421 | pca2 = 'tmp%s.pca.2' % pid
|
|---|
| 422 | pca3 = 'tmp%s.pca.3' % pid
|
|---|
| 423 | b1evect1 = b1evect[0]
|
|---|
| 424 | b1evect2 = b1evect[1]
|
|---|
| 425 | b1evect3 = b1evect[2]
|
|---|
| 426 | b2evect1 = b2evect[0]
|
|---|
| 427 | b2evect2 = b2evect[1]
|
|---|
| 428 | b2evect3 = b2evect[2]
|
|---|
| 429 | b3evect1 = b3evect[0]
|
|---|
| 430 | b3evect2 = b3evect[1]
|
|---|
| 431 | b3evect3 = b3evect[2]
|
|---|
| 432 |
|
|---|
| 433 | # Histogram matching
|
|---|
| 434 | outname = 'tmp%s_pan1' % pid
|
|---|
| 435 | panmatch1 = matchhist(pan, ms1, outname)
|
|---|
| 436 |
|
|---|
| 437 | outname = 'tmp%s_pan2' % pid
|
|---|
| 438 | panmatch2 = matchhist(pan, ms2, outname)
|
|---|
| 439 |
|
|---|
| 440 | outname = 'tmp%s_pan3' % pid
|
|---|
| 441 | panmatch3 = matchhist(pan, ms3, outname)
|
|---|
| 442 |
|
|---|
| 443 | grass.message(_("Performing inverse PCA ..."))
|
|---|
| 444 |
|
|---|
| 445 | # Get mean value of each channel
|
|---|
| 446 | stats1 = grass.parse_command("r.univar", map=ms1, flags='g',
|
|---|
| 447 | parse=(grass.parse_key_val,
|
|---|
| 448 | {'sep': '='}))
|
|---|
| 449 | stats2 = grass.parse_command("r.univar", map=ms2, flags='g',
|
|---|
| 450 | parse=(grass.parse_key_val,
|
|---|
| 451 | {'sep': '='}))
|
|---|
| 452 | stats3 = grass.parse_command("r.univar", map=ms3, flags='g',
|
|---|
| 453 | parse=(grass.parse_key_val,
|
|---|
| 454 | {'sep': '='}))
|
|---|
| 455 |
|
|---|
| 456 | b1mean = float(stats1['mean'])
|
|---|
| 457 | b2mean = float(stats2['mean'])
|
|---|
| 458 | b3mean = float(stats3['mean'])
|
|---|
| 459 |
|
|---|
| 460 | if sproc:
|
|---|
| 461 | # serial processing
|
|---|
| 462 | outr = '%s_red' % out
|
|---|
| 463 | outg = '%s_green' % out
|
|---|
| 464 | outb = '%s_blue' % out
|
|---|
| 465 |
|
|---|
| 466 | cmd1 = "$outb = 1 * round(($panmatch1 * $b1evect1) + ($pca2 * $b1evect2) + ($pca3 * $b1evect3) + $b1mean)"
|
|---|
| 467 | cmd2 = "$outg = 1 * round(($panmatch2 * $b2evect1) + ($pca2 * $b2evect2) + ($pca3 * $b2evect3) + $b2mean)"
|
|---|
| 468 | cmd3 = "$outr = 1 * round(($panmatch3 * $b3evect1) + ($pca2 * $b3evect2) + ($pca3 * $b3evect3) + $b3mean)"
|
|---|
| 469 |
|
|---|
| 470 | cmd = '\n'.join([cmd1, cmd2, cmd3])
|
|---|
| 471 |
|
|---|
| 472 | grass.mapcalc(cmd, outb=outb, outg=outg, outr=outr,
|
|---|
| 473 | panmatch1=panmatch1,
|
|---|
| 474 | panmatch2=panmatch2,
|
|---|
| 475 | panmatch3=panmatch3,
|
|---|
| 476 | pca2=pca2,
|
|---|
| 477 | pca3=pca3,
|
|---|
| 478 | b1evect1=b1evect1,
|
|---|
| 479 | b2evect1=b2evect1,
|
|---|
| 480 | b3evect1=b3evect1,
|
|---|
| 481 | b1evect2=b1evect2,
|
|---|
| 482 | b2evect2=b2evect2,
|
|---|
| 483 | b3evect2=b3evect2,
|
|---|
| 484 | b1evect3=b1evect3,
|
|---|
| 485 | b2evect3=b2evect3,
|
|---|
| 486 | b3evect3=b3evect3,
|
|---|
| 487 | b1mean=b1mean,
|
|---|
| 488 | b2mean=b2mean,
|
|---|
| 489 | b3mean=b3mean,
|
|---|
| 490 | overwrite=True)
|
|---|
| 491 | else:
|
|---|
| 492 | # parallel processing
|
|---|
| 493 | pb = grass.mapcalc_start('%s_blue = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)'
|
|---|
| 494 | % (out,
|
|---|
| 495 | panmatch1,
|
|---|
| 496 | b1evect1,
|
|---|
| 497 | pca2,
|
|---|
| 498 | b1evect2,
|
|---|
| 499 | pca3,
|
|---|
| 500 | b1evect3,
|
|---|
| 501 | b1mean),
|
|---|
| 502 | overwrite=True)
|
|---|
| 503 |
|
|---|
| 504 | pg = grass.mapcalc_start('%s_green = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)'
|
|---|
| 505 | % (out,
|
|---|
| 506 | panmatch2,
|
|---|
| 507 | b2evect1,
|
|---|
| 508 | pca2,
|
|---|
| 509 | b2evect2,
|
|---|
| 510 | pca3,
|
|---|
| 511 | b2evect3,
|
|---|
| 512 | b2mean),
|
|---|
| 513 | overwrite=True)
|
|---|
| 514 |
|
|---|
| 515 | pr = grass.mapcalc_start('%s_red = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)'
|
|---|
| 516 | % (out,
|
|---|
| 517 | panmatch3,
|
|---|
| 518 | b3evect1,
|
|---|
| 519 | pca2,
|
|---|
| 520 | b3evect2,
|
|---|
| 521 | pca3,
|
|---|
| 522 | b3evect3,
|
|---|
| 523 | b3mean),
|
|---|
| 524 | overwrite=True)
|
|---|
| 525 |
|
|---|
| 526 | pb.wait(), pg.wait(), pr.wait()
|
|---|
| 527 | try:
|
|---|
| 528 | pb.terminate(), pg.terminate(), pr.terminate()
|
|---|
| 529 | except:
|
|---|
| 530 | ""
|
|---|
| 531 |
|
|---|
| 532 | # Cleanup
|
|---|
| 533 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 534 | name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
|
|---|
| 535 |
|
|---|
| 536 | def matchhist(original, target, matched):
|
|---|
| 537 | # pan/intensity histogram matching using numpy arrays
|
|---|
| 538 | grass.message(_("Histogram matching..."))
|
|---|
| 539 |
|
|---|
| 540 | # input images
|
|---|
| 541 | original = original.split('@')[0]
|
|---|
| 542 | target = target.split('@')[0]
|
|---|
| 543 | images = [original, target]
|
|---|
| 544 |
|
|---|
| 545 | # create a dictionary to hold arrays for each image
|
|---|
| 546 | arrays = {}
|
|---|
| 547 |
|
|---|
| 548 | for img in images:
|
|---|
| 549 | # calculate number of cells for each grey value for for each image
|
|---|
| 550 | stats_out = grass.pipe_command('r.stats', flags='cin', input=img,
|
|---|
| 551 | sep=':')
|
|---|
| 552 | stats = stats_out.communicate()[0].split('\n')[:-1]
|
|---|
| 553 | stats_dict = dict(s.split(':', 1) for s in stats)
|
|---|
| 554 | total_cells = 0 # total non-null cells
|
|---|
| 555 | for j in stats_dict:
|
|---|
| 556 | stats_dict[j] = int(stats_dict[j])
|
|---|
| 557 | if j != '*':
|
|---|
| 558 | total_cells += stats_dict[j]
|
|---|
| 559 |
|
|---|
| 560 | if total_cells < 1:
|
|---|
| 561 | grass.fatal(_("Input has no data. Check region settings."))
|
|---|
| 562 |
|
|---|
| 563 | # Make a 2x256 structured array for each image with a
|
|---|
| 564 | # cumulative distribution function (CDF) for each grey value.
|
|---|
| 565 | # Grey value is the integer (i4) and cdf is float (f4).
|
|---|
| 566 |
|
|---|
| 567 | arrays[img] = np.zeros((256, ), dtype=('i4,f4'))
|
|---|
| 568 | cum_cells = 0 # cumulative total of cells for sum of current and all lower grey values
|
|---|
| 569 |
|
|---|
| 570 | for n in range(0, 256):
|
|---|
| 571 | if str(n) in stats_dict:
|
|---|
| 572 | num_cells = stats_dict[str(n)]
|
|---|
| 573 | else:
|
|---|
| 574 | num_cells = 0
|
|---|
| 575 |
|
|---|
| 576 | cum_cells += num_cells
|
|---|
| 577 |
|
|---|
| 578 | # cdf is the the number of cells at or below a given grey value
|
|---|
| 579 | # divided by the total number of cells
|
|---|
| 580 | cdf = float(cum_cells) / float(total_cells)
|
|---|
| 581 |
|
|---|
| 582 | # insert values into array
|
|---|
| 583 | arrays[img][n] = (n, cdf)
|
|---|
| 584 |
|
|---|
| 585 | # open file for reclass rules
|
|---|
| 586 | outfile = open(grass.tempfile(), 'w')
|
|---|
| 587 |
|
|---|
| 588 | for i in arrays[original]:
|
|---|
| 589 | # for each grey value and corresponding cdf value in original, find the
|
|---|
| 590 | # cdf value in target that is closest to the target cdf value
|
|---|
| 591 | difference_list = []
|
|---|
| 592 | for j in arrays[target]:
|
|---|
| 593 | # make a list of the difference between each original cdf value and
|
|---|
| 594 | # the target cdf value
|
|---|
| 595 | difference_list.append(abs(i[1] - j[1]))
|
|---|
| 596 |
|
|---|
| 597 | # get the smallest difference in the list
|
|---|
| 598 | min_difference = min(difference_list)
|
|---|
| 599 |
|
|---|
| 600 | for j in arrays[target]:
|
|---|
| 601 | # find the grey value in target that corresponds to the cdf
|
|---|
| 602 | # closest to the original cdf
|
|---|
| 603 | if j[1] <= i[1] + min_difference and j[1] >= i[1] - min_difference:
|
|---|
| 604 | # build a reclass rules file from the original grey value and
|
|---|
| 605 | # corresponding grey value from target
|
|---|
| 606 | out_line = "%d = %d\n" % (i[0], j[0])
|
|---|
| 607 | outfile.write(out_line)
|
|---|
| 608 | break
|
|---|
| 609 |
|
|---|
| 610 | outfile.close()
|
|---|
| 611 |
|
|---|
| 612 | # create reclass of target from reclass rules file
|
|---|
| 613 | result = grass.core.find_file(matched, element='cell')
|
|---|
| 614 | if result['fullname']:
|
|---|
| 615 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 616 | name=matched)
|
|---|
| 617 | grass.run_command('r.reclass', input=original, out=matched,
|
|---|
| 618 | rules=outfile.name)
|
|---|
| 619 | else:
|
|---|
| 620 | grass.run_command('r.reclass', input=original, out=matched,
|
|---|
| 621 | rules=outfile.name)
|
|---|
| 622 |
|
|---|
| 623 | # Cleanup
|
|---|
| 624 | # remove the rules file
|
|---|
| 625 | grass.try_remove(outfile.name)
|
|---|
| 626 |
|
|---|
| 627 | # return reclass of target with histogram that matches original
|
|---|
| 628 | return matched
|
|---|
| 629 |
|
|---|
| 630 | if __name__ == "__main__":
|
|---|
| 631 | options, flags = grass.parser()
|
|---|
| 632 | main()
|
|---|
