i.pansharpen limited to 8-bit imagery

[Nikos Alexandris]

Getting 11-bit IKONOS bands pan-sharpened by using i.pansharpen works only if the bands are rescaled to 8-bit, i.e. to range inside [0, 255]. In this case, color re-balancing (i.e. of Blue, Green and Red fusion-ed with Pan -IKONOS spectral bands) after the sharpening appears acceptable.

Any other attempts using out of [0, 255] range values, for example using Top-of-Atmosphere Reflectances (floating point values < 1.0) or the raw IKONOS bands that range between [0, 2047], are not successful.

Rescaling, however, might lead to some loss of spectral information (?). Hence, it may not be desired as pan-sharpening should not be limited in creating nice looking (pseudo-)color composites for viewing tasks (i.e. image interpretation) or printing. Sharpened/Fusioned images might as well be used as an input for segmentation/classification etc.

Looking at the i.pansharpen source code (​as suggested by Moritz Lennert in the grass-user mailing list) I see that in all algorithms a function for histogram matching is called and IIUC, in that function, values seem to be limited to [0,255].

Related code(lines): ​def matchhist and ​arrays[i = np.zeros((256,),dtype=('i4,f4'))]

Not sure if this is a defect or an enhancement request. Many of the High Resolution images are frequently delivered in >8-bit formats.

[cmbarton]

The question is whether the histogram matching code can handle resolutions of >8 bit.

And I guess I don't know. Have to make some modifications and see.

That said, the HIS pan sharpening only works with integers, because that is a limit of i.rgb.his and i.his.rgb. I don't know if they can handle integer ranges outside 0-255, however. You might give it a try to transform one of your images with one of these and see what it converts it to.

Michael

[mlennert]

Replying to cmbarton:

The question is whether the histogram matching code can handle resolutions of >8 bit.

And I guess I don't know. Have to make some modifications and see.

In the long run, it would be interesting to integrate i.histo.match [1] from addons into trunk. Then you could call that from i.pansharpen. However, the code of i.histo.match could definitely use some optimisation for performance.

[1] ​https://trac.osgeo.org/grass/browser/grass-addons/grass7/imagery/i.histo.match

[Nikos Alexandris]

Replying to cmbarton:

The question is whether the histogram matching code can handle resolutions of >8 bit. And I guess I don't know. Have to make some modifications and see.

Can we alter the manual somehow? To show how it can be done (now) with >8-bit images? I have already posted some example at ​http://grasswiki.osgeo.org/wiki/IKONOS#Pan_Sharpening.

That said, the HIS pan sharpening only works with integers, because that is a limit of i.rgb.his and i.his.rgb. I don't know if they can handle integer ranges outside 0-255, however. You might give it a try to transform one of your images with one of these and see what it converts it to.

I can/will try to check. First look in the manual of i.rgb.his does not mention anything about bit-depth(s).

[cmbarton]

If you have an image set that is more than 8bit, I can use it to test some things.

i.histo.match is a nice module. But its objective is different from histogram matching in i.pan.sharpen. So it would need modification to be used in this context. When I was writing i.pan.sharpen, I looked at the i.histo.match code but it was easier to use a much simpler algorithm. But since you know i.histo.match maybe you can see where the code could be modified to be used in i.pan.sharpen.

Thanks Michael

[Nikos Alexandris]

Replying to cmbarton:

If you have an image set that is more than 8bit, I can use it to test some things.

Unfortunately I cannot share the data I work with. However, there are publicly available 11-bit IKONOS images: ​ftp://ftp.glcf.umd.edu/glcf/China_earthquake_May_2008/IKONOS/

i.histo.match is a nice module. But its objective is different from histogram matching in i.pan.sharpen.

Unfortunately i.histo.match has its own limitation(s): it understands only integers - no floating points. See ​http://lists.osgeo.org/pipermail/grass-dev/2013-February/062161.html. Though, not limited to 2⁸.

So it would need modification to be used in this context. When I was writing i.pan.sharpen, I looked at the i.histo.match code but it was easier to use a much simpler algorithm. But since you know i.histo.match maybe you can see where the code could be modified to be used in i.pan.sharpen.

Maybe... It certainly would be a nice quick-"fix" if i.pansharpen can be released from being restricted to 2⁸? I mean, to expand in the "integer" planet as a first step?

Will try to check further in the coming days/week. However, I'd like to kindly ask again if it is possible to alter the documentation? It is not a nice surprise trying to work with i.pansharpen and getting strange colors or empty images as a result. See also ​http://gis.stackexchange.com/q/39174/5256. At least mention the current limitation or... other ideas? Then again, it is in G7, the development version :-)

[cmbarton]

Actually, looking at the code again, I'm pretty sure that it works as well with floating point values and ranges beyond 256 as it does with 8bit images--for the Brovey and PCA sharpening. The histogram matching routine accepts floating points numbers. It runs them through r.stats to create a 0-255 histogram. So it converts them to 8bit before outputting them.

But I don't think this is what is causing your color problems. Some images I ran through this in testing came out great and other did not, using the same resolution of image (all 8 bit in my tests). I think that the issue is preprocessing. IIRC, the best results come if you can estimate radiance (i.e., dealing with atmospheric corrections and shadows due to topography) at the surface and the worst come from uncorrected images.

You can do some color correction post facto using Markus' landsat color correction script. But mostly you need to do a significant amount of preprocessing to get the best results. It would be good if someone with more image processing experience could weigh in and offer some insight.

The ideal script--or perhaps better the ideal workflow--would combine a set of preprocessing steps with the pan sharpening. The latter is a set of processing algorithms that is completely ignorant of color.

Michael

[Nikos Alexandris]

Replying to cmbarton:

Actually, looking at the code again, I'm pretty sure that it works as well with floating point values and ranges beyond 256 as it does with 8bit images--for the Brovey and PCA sharpening. The histogram matching routine accepts floating points numbers. It runs them through r.stats to create a 0-255 histogram. So it converts them to 8bit before outputting them.

OK, maybe this is the "crux" then?

But I don't think this is what is causing your color problems. Some images I ran through this in testing came out great and other did not, using the same resolution of image (all 8 bit in my tests).

No problems when using 8-bit images here as well. The point is how to pan-sharpen 11-bit images, or 16-bit, or else? Any such examples available?

FWIW, I might not have described well the "request". Apologies -- honestly. In effect, it's not about color-balancing per se. It is about getting meaningful output from i.pansharpen.

Michael,

I tested various combinations here: "?CELL" types & "sharpen" methods: raw DNs 11-bit [0, 2047], DNs converted to 8-bit [0, 255], Radiances (random example of a Blue band ranging in [0, 394.363700815314]), Reflectances in [0, 1.0]. The only way to get a meaningful output is to rescale the material to [0, 255]. In all other cases without exceptions, colors are at best fancy, at worse "empty"; and I can't get them, no matter what I try, re-balanced (tried grey, grey255, grey1.0, grey.eq, with or without "-e", etc.) pre- and post- sharpening.

Among the various attempts, many of the outputs of i.pansharpen when _not_ using 8-bit as an input, range between 0 and 4 or 1!

I think that the issue is preprocessing. IIRC, the best results come if you can estimate radiance (i.e., dealing with atmospheric corrections and shadows due to topography) at the surface and the worst come from uncorrected images.

I disagree that this has to do with the degree of "fancyness" I am facing :-). It might, and will, help to get more visually appealing images/composites. "Precision" is, however, not the problem here.

You can do some color correction post facto using Markus' landsat color correction script.

Sure, it works for 8-bit inputs -- and very nice as I can see. Even more, no need to balance before sharpening actually. Why should one do that anyway?

But mostly you need to do a significant amount of preprocessing to get the best results. It would be good if someone with more image processing experience could weigh in and offer some insight.

With or without pre-processing, I am not able to get meaningful images here. I am probably doing something wrong. I am available to test and share anything (but the images I work with). But we can work with the publicly available material.

The ideal script--or perhaps better the ideal workflow--would combine a set of preprocessing steps with the pan sharpening. The latter is a set of processing algorithms that is completely ignorant of color.

That would be really great! But this would be something for High(est) Precision hunting. My problem, as well as the problems described by Haziq in the respective Q&A is to get at least something that approaches a real-RGB image after pan-sharpening and using the R, G and B channels to draw an RGB composite.

Thanks

[mlennert]

Replying to cmbarton:

If you have an image set that is more than 8bit, I can use it to test some things.

i.histo.match is a nice module. But its objective is different from histogram matching in i.pan.sharpen. So it would need modification to be used in this context. When I was writing i.pan.sharpen, I looked at the i.histo.match code but it was easier to use a much simpler algorithm. But since you know i.histo.match maybe you can see where the code could be modified to be used in i.pan.sharpen.

Your histogram matching code matches the histogram of a source image to that of a target image, whereas i.histo.match matches the histogram of each given image to the cumulative histogram of all images. Both approaches are valid, and both should be available in a histogram matching module.

I think in both modules the routines can probably be improved in terms of performance.

Even though it shouldn't be too complicated to integrate source-target matching in i.histo.match, I don't have the time to look into that, now. If anyone has a student that needs a project, this would be a nice exercise ;-)

Moritz

[cmbarton]

For pan sharpening, the reason for histogram matching is that the higher resolution pan image needs to replace one of the lower resolution images as part of the sharpening algorithm. So the pan band needs to be histogram matched to the one it is replacing. This is why I don't think that matching it to an average of all bands is a good idea. I could be mistaken, but it seems incorrect conceptually.

The routine I used is the standard one discussed for this. There may well be a more accurate one and a better way to deal with floating point and large integer data. Using numpy arrays makes this one quite fast. It might benefit from parallelization in some way. I parallelized other parts of the code but not this one. So, yes, if there is a student out there with not enough to do... ;-)

Michael

[mlennert]

Replying to cmbarton:

For pan sharpening, the reason for histogram matching is that the higher resolution pan image needs to replace one of the lower resolution images as part of the sharpening algorithm. So the pan band needs to be histogram matched to the one it is replacing. This is why I don't think that matching it to an average of all bands is a good idea. I could be mistaken, but it seems incorrect conceptually.

I agree with you that for i.pansharpen that's the way to go. What I said is that I think i.histo.match should ideally allow the use to chose between the different techniques.

The routine I used is the standard one discussed for this. There may well be a more accurate one and a better way to deal with floating point and large integer data. Using numpy arrays makes this one quite fast. It might benefit from parallelization in some way. I parallelized other parts of the code but not this one. So, yes, if there is a student out there with not enough to do... ;-)

In terms of optimisation I was mostly thinking about the double loop through arrays[target] within the 'for i in arrays[original]' loop. I think that there are more efficient search algorithms, but yours definitely works.

Moritz

[Nikos Alexandris]

Some (minor) things to improve i.pansharpen:

• make "sharpen=brovey" progress notifications more informative, for example instead of solely "Histogram matching...", inform the user what is being histo-matched to what. I.e., "Matching histogram of Input_Map to histogram of Reference_Map"

• it seems that the module (at least when "sharp=brovey" ?) does not like dashes in the input images names:

i.pansharpen sharpen=brovey ms1=M1BS-01_P001.1 ms2=M1BS-01_P001.2 ms3=M1BS-01_P001.3 pan=P1BS-01_P001 output_prefix=sharpen_brovey_01_P001

Performing pan sharpening with hi res pan image: 1.000000
Using Brovey algorithm

Histogram matching...
 100%
 100%

Histogram matching...
 100%
 100%

Histogram matching...
 100%
 100%

Calculating Brovey transformation...
Invalid map <M1BS>
Invalid map <01_P001.1>
Parse error
ERROR: parse error
Invalid map <M1BS>
Invalid map <01_P001.2>
Parse error
ERROR: parse error
Invalid map <M1BS>
Invalid map <01_P001.3>
Parse error
ERROR: parse error

Assigning grey equalized color tables to output images...
ERROR: Raster map <sharpen_brovey_01_P001_red> not found
ERROR: Raster map <sharpen_brovey_01_P001_green> not found
ERROR: Raster map <sharpen_brovey_01_P001_blue> not found

The following pan-sharpened output maps have been generated:
sharpen_brovey_01_P001_red
sharpen_brovey_01_P001_green
sharpen_brovey_01_P001_blue

To visualize output, run: g.region -p rast=sharpen_brovey_01_P001.red
d.rgb r=sharpen_brovey_01_P001_red g=sharpen_brovey_01_P001_green
b=sharpen_brovey_01_P001_blue

If desired, combine channels into a single RGB map with 'r.composite'.
Channel colors can be rebalanced using i.landsat.rgb.
WARNING: Unable to write history for <sharpen_brovey_01_P001_red>. Raster
         map <sharpen_brovey_01_P001_red> not found in current mapset.
WARNING: Unable to write history for <sharpen_brovey_01_P001_green>. Raster
         map <sharpen_brovey_01_P001_green> not found in current mapset.
WARNING: Unable to write history for <sharpen_brovey_01_P001_blue>. Raster
         map <sharpen_brovey_01_P001_blue> not found in current mapset.

[Nikos Alexandris]

Attaching two overview screenshots (uploaded in GRASS-Wiki, links below) with all "sharpen= ihs, brovey, pca" methods in 11-bit IKONOS DNs, and double precision both Radiance and Reflectance values. (Unfortunately, it's not possible to draw d.rgb directly in a d.frame.split-ed d.mon, thus I created composites)

None of them appears ok to me. Even more, NULL cells have been introduced.

• ​http://grasswiki.osgeo.org/wiki/File:IKONOS_PanSharpening_Three_Inputs_Greyscaled_Three_Methods.png

• ​http://grasswiki.osgeo.org/wiki/File:IKONOS_PanSharpening_Three_Inputs_Rebalanced_Three_Methods.png

Working with rescaled 8-bit versions, produced, so far for me, nice outputs.

[mlennert]

Replying to nikosa:

However, I'd like to kindly ask again if it is possible to alter the documentation? It is not a nice surprise trying to work with i.pansharpen and getting strange colors or empty images as a result. See also ​http://gis.stackexchange.com/q/39174/5256. At least mention the current limitation or... other ideas? Then again, it is in G7, the development version :-)

I just (r57332) added a sentence to the NOTES section about the module currently working only for 8-bit images.

Several remarks to justify that:

• The histogram matching algorithm in the module is limited to 8-bit as it reclasses only sources values between 0 and 255 to target value between 0 and 255 (it also calls r.stats with the -i flag, and so rounds everything to integer, IIUC). In other words, any values above 255 silently ignored ! To check, just comment out line where the temporary histogram matched (reclassed) images are g.removed. As the histogram matching is used in all three algorithms, they are all limited to 8-bit.

• the i.rgb.his/i.his.rgb modules are also limited to 8-bit, AFAIK. There is an old ticket with a patch proposed by Hamish on that: #774. See also [1]

So, in conclusion, to be able to use i.pansharpen with >8-bit images, the histogram matching algorithm needs to be rewritten. You could try to just skip the histogram matching part (just comment out the relevant call to the function and replace the panmatch images by the pan image and see what happens (I'll have to read up on the theory again to understand why histogram matching is necessary here).

You can also apply the brovey or pca method yourself using r.mapcalc and i.pca.

[1] ​http://lists.osgeo.org/pipermail/grass-user/2009-September/052573.html

[cmbarton]

Does someone have a small sample set of images to test in order to work out alternatives to 8bit image matching?

A key issue is that it will be necessary to know the potential value range for the images. I don't see a way around this being a user-entered value. Maybe there also needs to be a flag for floating point.

[Nikos Alexandris]

Replying to cmbarton:

Does someone have a small sample set of images to test in order to work out alternatives to 8bit image matching?

See, for example, publicly available IKONOS images: ​ftp://ftp.glcf.umd.edu/glcf/China_earthquake_May_2008/IKONOS/. I did all the series of tests (as postes in the mailing list) with the first one ​ftp://ftp.glcf.umd.edu/glcf/China_earthquake_May_2008/IKONOS/po_58204_0000000.20001116.China-Sichuan/

I also have a small script to pan-sharpen DNs, Radiance and Reflectancs (given all of the pre-processing has been done and some naming conventions are followed):

for Method_Input in "ihs DNs" "ihs Radiance" "ihs ToAR" \
"brovey DNs" "brovey Radiance" "brovey ToAR" \
"pca DNs" "pca Radiance" "pca ToAR"

do

  # parse "${Method_STR}" and set positional parameters
  set -- $Method_Input ; echo $1 $2

  # some echo
  echo "Pan-Sharpening the <${2}> images based on the <${1}> method"

  i.pansharpen \
  sharpen=${1} \
  pan=Pan_${2} \
  ms1=Blue_${2} \
  ms2=Green_${2} \
  ms3=Red_${2} \
  output=sharp_${1}_${2}

done

A key issue is that it will be necessary to know the potential value range for the images.

If it's an 11-bit image, the the input can be fixed to [0,2047]? Actually, most of the data come in either 8-bit or 16-bit form. The latter for not loosing the fine detail of >8-bit sensors (me thinks).

I don't see a way around this being a user-entered value. Maybe there also needs to be a flag for floating point.

Can't this be automated in Python? Right, RTFM for myself :D

[cmbarton]

Replying to nikosa:

Replying to cmbarton:

Does someone have a small sample set of images to test in order to work out alternatives to 8bit image matching?

See, for example, publicly available IKONOS images: ​ftp://ftp.glcf.umd.edu/glcf/China_earthquake_May_2008/IKONOS/. I did all the series of tests (as postes in the mailing list) with the first one ​ftp://ftp.glcf.umd.edu/glcf/China_earthquake_May_2008/IKONOS/po_58204_0000000.20001116.China-Sichuan/

I also have a small script to pan-sharpen DNs, Radiance and Reflectancs (given all of the pre-processing has been done and some naming conventions are followed):

for Method_Input in "ihs DNs" "ihs Radiance" "ihs ToAR" \
"brovey DNs" "brovey Radiance" "brovey ToAR" \
"pca DNs" "pca Radiance" "pca ToAR"

do

  # parse "${Method_STR}" and set positional parameters
  set -- $Method_Input ; echo $1 $2

  # some echo
  echo "Pan-Sharpening the <${2}> images based on the <${1}> method"

  i.pansharpen \
  sharpen=${1} \
  pan=Pan_${2} \
  ms1=Blue_${2} \
  ms2=Green_${2} \
  ms3=Red_${2} \
  output=sharp_${1}_${2}

done

A key issue is that it will be necessary to know the potential value range for the images.

If it's an 11-bit image, the the input can be fixed to [0,2047]? Actually, most of the data come in either 8-bit or 16-bit form. The latter for not loosing the fine detail of >8-bit sensors (me thinks).

I don't see a way around this being a user-entered value. Maybe there also needs to be a flag for floating point.

Can't this be automated in Python? Right, RTFM for myself :D

Actually, Python doesn't care if it is floating point or not (although you can specify). It's the GRASS routines that matter most, especially r.stats and r.mapcalc (can't can't change this for r.his.rgb/r.rgb.his or i.pca).

Michael

[Nikos Alexandris]

Replying to cmbarton:

The question is whether the histogram matching code can handle resolutions of >8 bit.

Michael, for viewing purposes we all agree (don't we?) that 8-bit per "channel" is more than enough. For analysis, my guess is that it would be interesting to have the "fine digits" preserved. But, with the histogram matching on the way, I don't know really if and at what extent the data are irreparably re-touched.

Anyway,

did you folks notice that this ticket #2048 just hit the upper... 2^{11 limit? :D}

[Nikos Alexandris]

Without the intention to add "noise", besides the many publications out there who perform segmentation and image classification on pan-sharpened 11(16)-bit data (need for fine radiometry), there are also consumer- & professional-grade displays that are capable of 11-bit. They are costly. Nevertheless, if the main concern is work and eye-health when spending countless hours in front of a display, one might consider getting one if possible.

GRASS will become more attractive with such enhancements.

Thanks, N

[neteler]

See also #774

[martinl]

What is the state of this issue?

[cmbarton]

I haven't made any modifications, although I think some other people have. I could look into this again when I'm on sabbatical next year.

[cmbarton]

Maybe while on sabbatical, I'll FINALLY have time to dig into this and see what can be done.

[cmbarton]

Before I start on this, I have a couple questions:

1. #774 is a discussion about support for >8bit values in i.rgb.his and i.his.rgb. I cannot tell from this discussion or from the module manual if this happened or not. I'm guessing it did not. Can anyone verify?

2. The main place where i.pansharpen is explicitly 8bit is in the histogram matching method. The main goal of pan sharpening, AFAICT, is to produce an RGB image that is visually pleasing at spatial resolutions equivalent to a high res pan band. The original 'data' are lost regardless of the bit depth. Given this, the most straightforward solution to #1 and #2 is to transform the original 4 bands (RGB and pan) to 8bit before doing any processing. Any reason not to pursue that strategy?

[Nikos Alexandris]

Replying to cmbarton:

Before I start on this, I have a couple questions:

1. #774 is a discussion about support for >8bit values in i.rgb.his and i.his.rgb. I cannot tell from this discussion or from the module manual if this happened or not. I'm guessing it did not. Can anyone verify?

Please read (also) comment ​https://trac.osgeo.org/grass/ticket/774#comment:21. Would be nice if anyone advanced/core dev would look into it.

[cmbarton]

To follow up with the proposal in comment 28 above, the easiest approach would be to add the option for the user to select bit depth (following one of the suggestions in the comments). I would add options for 8, 11, and 16 bit. Any reason to include 32 bit or other?

[Nikos Alexandris]

Replying to cmbarton:

To follow up with the proposal in comment 28 above, the easiest approach would be to add the option for the user to select bit depth (following one of the suggestions in the comments). I would add options for 8, 11, and 16 bit. Any reason to include 32 bit or other?

Dear Michael, it's great you work on this! I would leave all the range [2, 16] or [2, 32] open to the user. One never knows what needs there might be.

However,

• for bitness < 2, there I have a comment here ​https://github.com/NikosAlexandris/test_color_space_conversions/blob/master/test_color_space_roundtrips.py#L31 about an r.rescale flaw which means, my test failed for bitness = 2.

• for bitness >= 15 the test will fail for precision precision to 0.001

[cmbarton]

Niko, I'm not sure that much of the code--or assumptions underlying the pan sharpening--would work or make sense for bit depth <8.

I can check of course. The main reason to have user input for bit depth is to be able to calculate the potential raster value range to rescale to 8 bit for processing. This still assumes an integer raster rather than a floating point one.

[cmbarton]

I've attached an updated version of i.pansharpen (i.pansharpen2 here) that works with 8-32bit imagery. I've tested with artificially inflated bit depth images (8bit landsat rescaled to 11bit) and it seems to work fine. I'd appreciate it if others could test with imagery that is >8bit.

[Nikos Alexandris]

Thanks Michael, I've skimmed through.

1. What justifies the assumption at ​https://trac.osgeo.org/grass/attachment/ticket/2048/i.pansharpen2#L138 and L139?

Instead, what about returning an error, if the input bitness of an image is incompatible with what the program can handle, or at least, issue a big Warning to the user.

2. ​https://trac.osgeo.org/grass/attachment/ticket/2048/i.pansharpen2#L152 forces 8-bit range.

Hopefully we can address this too by using something like the modified versions of i.rgb.his and i.his.rgb (as demonstrated here ​https://trac.osgeo.org/grass/ticket/774#comment:13). There is no mathematical constrain to force an 8-bit range. At least not for IHS and PCA. Is there? Regarding "Brovey", the only assumption I can find by quickly searching on the internet, is "It assumes that the spectral range spanned by the panchromatic image is the same as that covered by the multispectral channels." (source: ​https://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/fundamentals-of-panchromatic-sharpening.htm).

[Nikos Alexandris]

As for the histogram matching routine, here two related open source examples (written in Rust):

​https://github.com/jblindsay/whitebox-tools/blob/266c939ad236c16503b088d2faa4500362912fef/src/tools/image_analysis/histogram_matching.rs

and

​https://github.com/jblindsay/whitebox-tools/blob/266c939ad236c16503b088d2faa4500362912fef/src/tools/image_analysis/histogram_matching_two_images.rs

My understanding is that the implementation bases upon (64-bit) floating point.

[cmbarton]

I am hoping that someone can test the revised code. I'll work on updating the docs a bit to match the new options.

The algorithms used in the script were devised with assumed 8bit image channels. The math might work with other bit depths, but other assumptions might be violated. Maybe someone else can clarify. Beyond that, the histogram matching method would need to be redesigned to accommodate imagery other than 8bit. Allowing for floating point image value would require yet more reworking, including new versions of i.rgb.his and i.his.rgb. Hence, rescaling the original images to 8 bit seems like a good solution at the present.

While I'm comfortable with downscaling from higher to lower bit depth (e.g., 16 down to 8), I'm less comfortable with upscaling imagery at less than 8 bits. I tried this on 3bit imagery (rescaled from 8bit landsat) and it produced visually poor results that do not show the level of feature differentiation seen in the 3 bit pan channel alone. So keeping it at ≥8bit seems advisable unless there is a reason pan "sharpening" of lower bit depth is actually useful for something. The cutoff at 32 bits is arbitrary and could easily be increased to higher values. Any reason to do this or not do this if the images are rescaled to 8bits for processing anyway?

If you run i.pansharpen in the terminal, anything outside the 8-32 bit range will already generate an error with a message that bit depth is outside the allowed range. If you do it in the GUI, a message will only be generated in the console. I could add a wxPython pop up, but not sure it is worth it. Maybe the best is to have the module fail (i.e., with a return) if bit depth is outside the acceptable range rather than set it back to 8.

[cmbarton]

Replying to Nikos Alexandris:

As for the histogram matching routine, here two related open source examples (written in Rust):

​https://github.com/jblindsay/whitebox-tools/blob/266c939ad236c16503b088d2faa4500362912fef/src/tools/image_analysis/histogram_matching.rs

and

​https://github.com/jblindsay/whitebox-tools/blob/266c939ad236c16503b088d2faa4500362912fef/src/tools/image_analysis/histogram_matching_two_images.rs

My understanding is that the implementation bases upon (64-bit) floating point.

I would prefer to get this upgrade to i.pansharpen tested and into GRASS for use before any complete rewrite of the module (maybe a task for someone else).

[Nikos Alexandris]

Replying to cmbarton:

I am hoping that someone can test the revised code. I'll work on updating the docs a bit to match the new options.

I will, as soon as I find some time.

The algorithms used in the script were devised with assumed 8bit image channels. The math might work with other bit depths, but other assumptions might be violated. Maybe someone else can clarify. Beyond that, the histogram matching method would need to be redesigned to accommodate imagery other than 8bit. Allowing for floating point image value would require yet more reworking, including new versions of i.rgb.his and i.his.rgb. Hence, rescaling the original images to 8 bit seems like a good solution at the present.

I understand. Let's go step by step then.

While I'm comfortable with downscaling from higher to lower bit depth (e.g., 16 down to 8), I'm less comfortable with upscaling imagery at less than 8 bits.

I agree. Yet, my understanding is, no rescaling should be performed. What goes in, should go out. Some loss of precision due to rounding applies. But no other loss should be justified. And this should be our target. Histogram matching doesn't or shouldn't care either about the input range. It should respect the given range, be it 4-bit or 64-bit. Some histogram matching routines, though, will introduce values to fill-in gaps while trying to match a "real world" histogram with the one of a function. But this is not rescaling.

I tried this on 3bit imagery (rescaled from 8bit landsat) and it produced visually poor results that do not show the level of feature differentiation seen in the 3 bit pan channel alone. So keeping it at ≥8bit seems advisable unless there is a reason pan "sharpening" of lower bit depth is actually useful for something.

Pan-sharpening fuses low and high spatial resolution images. These images are expected to capture the same object, in different dimensions. If we have 4-bit images, one high-res at 1m and one, or more, low-res at 4m, then these will be a valid subject for pan-sharpening. Please correct me if I am wrong.

The cutoff at 32 bits is arbitrary and could easily be increased to higher values. Any reason to do this or not do this if the images are rescaled to 8bits for processing anyway?

What about we support the maximum possible bitness that the GRASS raster engine can handle?

If you run i.pansharpen in the terminal, anything outside the 8-32 bit range will already generate an error with a message that bit depth is outside the allowed range.

I did and I did not read any error. Unless I did not try your latest update. In any case, thanks for taking the time to work on it.

[cmbarton]

Here is another limit I didn't realize.

r.rescale bombs with 32bit imagery, but can handle 30bit. So that needs to be the upper cut off for now.

[cmbarton]

Also, if you try to go beyond the range limit in the GUI, the value will be reset to the minimum (if you try to go below) or the maximum (if you try to go above. No error message is generated because the values are automatically corrected back to the limited range. In the terminal, you do get an error if you try to go outside the range.

[cmbarton]

updated pan sharpening module

[cmbarton]

I just uploaded a slightly revised version of i.pansharpen with bit depth in the 2-30 range and faster parallel processing of image channel rescaling.

[Nikos Alexandris]

Michael,

it works for me, i.e. 16-bit Landsat 8 bands used as input, are processed using the default method, and output 8-bit pan-sharpened bands.

Here,

# input
for RASTER in B8 B4 B5 B1 ;do echo $RASTER ;r.info -gr $RASTER ;echo ;done

B8
north=4423507.5
south=4188892.5
east=704707.5
west=474292.5
nsres=15
ewres=15
rows=15641
cols=15361
cells=240261401
datatype=CELL
ncats=0
min=0
max=65535

B4
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

B5
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

B1
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

and

#process
i.pansharpen pan=B8 red=B4 green=B5 blue=B1 bitdepth=16 output=pan_ihs method=ihs --o
i.pansharpen pan=B8 red=B4 green=B5 blue=B1 bitdepth=16 output=pan_brovey method=brovey --o
i.pansharpen pan=B8 red=B4 green=B5 blue=B1 bitdepth=16 output=pan_pca method=pca --o

and

# output
for RASTER in $(glr pattern=pan*) ;do echo $RASTER ;r.info -gr $RASTER ;echo ;done

pan_brovey_blue
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=0
max=21.8425196850394

pan_brovey_green
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=0
max=65.8663101604278

pan_brovey_red
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=0
max=28.5573122529644

pan_ihs_blue
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=CELL
ncats=0
min=0
max=77

pan_ihs_green
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=CELL
ncats=0
min=0
max=109

pan_ihs_red
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=CELL
ncats=0
min=0
max=87

pan_pca_blue
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=31.40914130045
max=96.916420771218

pan_pca_green
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=44.8984732158751
max=71.6383695727307

pan_pca_red
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=-41.5285779366908
max=33.7977309887422

[Nikos Alexandris]

Visually, after an i.colors.enhance, only PCA is a bit weird (​https://i.imgur.com/VsASM4V.png). Brovey (​https://i.imgur.com/d6QTSPK.png) and IHS (​https://i.imgur.com/CXTN1bc.png) are more as epxected.

[cmbarton]

I found a problem with the inverse PCA calculation and have fixed it. I'm uploading the revised version.

[Nikos Alexandris]

Update(d) with latest ​https://trac.osgeo.org/grass/attachment/ticket/2048/i.pansharpen.py and repeat for PCA. Here ​https://i.imgur.com/lkEoBRw.png, I think it is still not what we should expect color-wise. If someone else sees the same fuzzy colors, maybe we should go again through the steps for the PCA method.

# process
i.pansharpen pan=B8 red=B4 green=B5 blue=B1 bitdepth=16 output=pan_pca method=pca --o

# in & out
for RASTER in $(echo B8 B4 B5 B1 $(glr pattern=pan_pca*)) ;do echo $RASTER ;r.info -gr $RASTER ;echo ;done

B8
north=4423507.5
south=4188892.5
east=704707.5
west=474292.5
nsres=15
ewres=15
rows=15641
cols=15361
cells=240261401
datatype=CELL
ncats=0
min=0
max=65535

B4
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

B5
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

B1
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

pan_pca_blue
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=37.4635679507703
max=73.8477809701631

pan_pca_green
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=-60.815057189378
max=71.127075165179

pan_pca_red
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=20.3850355262278
max=72.9376049039095

[cmbarton]

new i.pansharpen examples

[cmbarton]

I just attached images generated from the newest version of i.pansharpen using the Landsat7 imagery.

Brovey: i_pansharpen_rgb_brovey542.jpg
IHS: i_pansharpen_rgb_ihs542.jpg​
PCA: i_pansharpen_rgb_pca542.jpg​
Unsharpened RGB: i_pansharpen_rgb_landsat542.jpg​

All look equally sharp (and equivalent to the pan band) and the colors seem closely similar.

[veroandreo]

Replying to Nikos Alexandris:

Update(d) with latest ​https://trac.osgeo.org/grass/attachment/ticket/2048/i.pansharpen.py and repeat for PCA. Here ​https://i.imgur.com/lkEoBRw.png, I think it is still not what we should expect color-wise. If someone else sees the same fuzzy colors, maybe we should go again through the steps for the PCA method.

# process
i.pansharpen pan=B8 red=B4 green=B5 blue=B1 bitdepth=16 output=pan_pca method=pca --o

# in & out
for RASTER in $(echo B8 B4 B5 B1 $(glr pattern=pan_pca*)) ;do echo $RASTER ;r.info -gr $RASTER ;echo ;done

B8
north=4423507.5
south=4188892.5
east=704707.5
west=474292.5
nsres=15
ewres=15
rows=15641
cols=15361
cells=240261401
datatype=CELL
ncats=0
min=0
max=65535

B4
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

B5
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

B1
north=4423515
south=4188885
east=704715
west=474285
nsres=30
ewres=30
rows=7821
cols=7681
cells=60073101
datatype=CELL
ncats=0
min=0
max=65535

pan_pca_blue
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=37.4635679507703
max=73.8477809701631

pan_pca_green
north=4315237.5
south=4302412.5
east=524647.5
west=512302.5
nsres=15
ewres=15
rows=855
cols=823
cells=703665
datatype=DCELL
ncats=0
min=-60.815057189378
max=71.127075165179

Maybe a silly question, but why does pan_pca_green has a minimum of -60? Shouldn't that be >= 0?

[cmbarton]

Can you give me access to these files so I can test?

[Nikos Alexandris]

Replying to cmbarton:

Can you give me access to these files so I can test?

Sent off-public (just to not waste space here) a sample from the Landsat 8 scene I use(d).

[cmbarton]

I tested i.pansharpen with these files. I was not able to reproduce the reported problem (no negative values in sharpened channels and all methods produce similar (though not identical) results. I did find an unrelated bug that can happen in rare cases and fixed it, and did a few other updates. Updated in trunk and merged to RB 7.6.

[martinl]

Ticket retargeted after milestone closed