[GRASS-dev] Implementation of the High Pass Filter Additive Fusion technique (i.fusion.hpf)

GRASS GIS GRASS GIS developer
1

Dear GRASS GIS users,

together with Nikos Ves, we share the “i.fusion.hpf” idea/proof of concept. At the moment, we have a custom shell script named i.fusion.hpf (an attempt for a proper GRASS add-on), which implements the High Pass Filter Additive (HPFA) Fusion Technique for Pan-Sharpening [*]. Nikos V started already porting to Python. How can we proceed in sharing it?

We require some help in testing. Specifically, comparing the results between the steps implemented in GRASS GIS and the respective “HPF Resolution Merge” in ERDAS. It seems that GRASS’ outputs are slightly more smooth. I will upload some screenshots in the Wiki. Perhaps I have done something wrong in the order of applying the algorithm’s steps (in the shell script)!?

Note, however, that the script misses the (optional, as explained in the related publication) histogram matching step. In ERDAS’ implementation, it seems that histo-matching is performed by default. Though, it might be possible to extract the “model” and subtract this final optional step so as to get 1:1 comparable outputs.

Background

HPFA seems to outperform the well known Pan-Sharpening techniques (incl. Brovey, IHS, PCA). The algorithm comprises the following steps:

  1. Computing ratio of low (Multi-Spectral) to high (Panchromatic) resolutions

  2. High Pass Filtering the Panchromatic Image

  3. Resampling MSX image to the higher resolution

  4. Adding weighted High-Pass-Filetred image to the upsampled MSX image

  5. Optionally, matching histogram of Pansharpened image to the one of the original MSX image

A few more words on the script

  • Stunning!, Crisp and colorful images (- currently after applying color rebalancing manually)

  • Extremely easy to use, i.e.: “i.fusion.hpf pan=Pan_DNs msx=Band1[,Band2,Band3,…]”

  • Grasping and testing the various parameters that define the High-Pass filter’s kernel size and center value is also a matter of short time.

  • It will work with any kind of imagery (after minor modifications)

  • However, it can be easily adapted for GRASS 7 / converted to Python

  • The attached script, badly coded by a non-programer, is in Bash for G64.

  • Currently works only for integers (with minor tweaking it can work with r.mfilter.fp to crunch Floating Points as well)

  • Lacks of the histogram matching operation

  • Something I don’t understand about g.tempfile – how to use it?

Two questions

? Can someone confirm that the part of the existing “i.pansharpen” code that performs histogram matching (code lines 348 - 431), do so as “linearly stretching an image to match another image’s Mean and StdDev”?

? Would it be desired to get the HPFA algorithm integrated in i.pansharpen?

For the records, the replication of the HPFA fusion technique in GRASS-GIS, as well as the “filter creation” bash one-liner, were Nikos Ve’s ideas. I followed-up with a bash proof of concept.

Nikos

[*] Optimizing the High-Pass Filter Addition Technique for Image Fusion (2008) by Ute G. Gangkofner, Pushkar S. Pradhan, and Derrold W. Holcomb

(attachments)

i.fusion.hpf.tar.gz (8.27 KB)

2

I don't understand how to use g.tempfile. The manual simply refers to (for
bash): temp1=`g.tempfile pid=$$`

My attempts to use it like that, end up with: "Illegal filename. Character </>
not allowed." Checking if the map's "filename" exists, returns for example:

--%<---
ERROR: Raster map
       </mnemosyne/geo/grassdb/nc_spm_08/landsat/.tmp/Resilience/26834.0>
       not found in current mapset
--->%--

I have been through the mailing list and have found past references to the
same problem. For example <http://lists.osgeo.org/pipermail/grass-user/2012-March/063977.html&gt;\. However, I can't get past the problem.

I tried to use naming conventions such as "i.fusion.hpf.tmp.$$" but I still
have some problems. I would like to use g.tempfile as it seems to be a
"cleaner" way(?).

Thank you, Nikos

3

The script does not seem to work as posted. Working on it... (trying to
understand g.tempfile above all right now).

Nikos

4

This creates confusion all the time.

For normal files, you use system tools (posix tools) for temporary files, e.g. mktemp.

For large files (large means large as a map), use g.tempfile which does not work so well as posix tools in terms of safety but it uses grass database because there is an assumption that there should be enough space.

For temporary maps there are two way. You can create temporary vector map in C (and in Python using C types). You cannot do this for rasters. The universal method is to use unique (enough) map name. Function for this should be placed to Python library, so we have the standardized method. It should be noted that g.tempfile -d is not usable for this.

I hope that this helps and that it is correct.

Vaclav

···

On Thu, Nov 14, 2013 at 12:05 PM, Nikos Alexandris <nik@nikosalexandris.net> wrote:

I don’t understand how to use g.tempfile. The manual simply refers to (for
bash): temp1=g.tempfile pid=$$

My attempts to use it like that, end up with: “Illegal filename. Character </>
not allowed.” Checking if the map’s “filename” exists, returns for example:

–%<—
ERROR: Raster map
</mnemosyne/geo/grassdb/nc_spm_08/landsat/.tmp/Resilience/26834.0>
not found in current mapset
—>%–

I have been through the mailing list and have found past references to the
same problem. For example <http://lists.osgeo.org/pipermail/grass-user/2012-March/063977.html>. However, I can’t get past the problem.

I tried to use naming conventions such as “i.fusion.hpf.tmp.$$” but I still
have some problems. I would like to use g.tempfile as it seems to be a
“cleaner” way(?).

Thank you, Nikos

grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

5

Nikos wrote:
> I don't understand how to use g.tempfile. The manual simply refers to (for

bash): temp1=`g.tempfile pid=$$`

My attempts to use it like that, end up with: "Illegal filename. Character
</> not allowed." Checking if the map's "filename" exists,
returns for example:

--%<---
ERROR: Raster map
</mnemosyne/geo/grassdb/nc_spm_08/landsat/.tmp/Resilience/26834.0>
not found in current mapset
--->%--

I have been through the mailing list and have found past references to the
same problem. For example
<http://lists.osgeo.org/pipermail/grass-user/2012-March/063977.html&gt;\.
However, I can't get past the problem.

I tried to use naming conventions such as "i.fusion.hpf.tmp.$$" but I
still have some problems. I would like to use g.tempfile as it seems
to be a "cleaner" way(?).

g.tempfile is for actual files like ascii text files, not for map names.

the convention used in the grass scripts for temporary map names is arbitrary,
but e.g. the g.module command might make one called: tmp_gmodule.$$.elev,
for vector maps replace '.' with '_' and bring out your curly brackets. :slight_smile:

then cleanup is as simple as g.mremove 'tmp_gmodule.$$.*', and if maps
are left behind it's obvious later that they are temporary and where they
came from.

there's no automatic handling, but I don't see where that would save any time.

what problems did you have with "i.fusion.hpf.tmp.$$" ?

Hamish

6

Thank you both Hamish and Vaclav.

Hamish:

g.tempfile is for actual files like ascii text files, not for map names.

got it!

the convention used in the grass scripts for temporary map names is

arbitrary, but e.g. the g.module command might make one called:

tmp_gmodule.$$.elev, for vector maps replace ‘.’ with ‘_’ and bring out

your curly brackets. :slight_smile:

ok. I’ll simply use “i.fusion.hpf.tmp.SomeString” and take care to have unique "SomeString"s. I’ll try to clean up the “curly” hair :slight_smile:

then cleanup is as simple as g.mremove ‘tmp_gmodule.$$.*’, and if maps

are left behind it’s obvious later that they are temporary and where they

came from.

there’s no automatic handling, but I don’t see where that would save any

time.

what problems did you have with “i.fusion.hpf.tmp.$$” ?

I think the problem was (among my misunderstanding) that $$ was the same for a bunch of names fed as variables here and there. All fixed now, it works :slight_smile:

Nikos

7

Here it goes, as Nikos Ves (depending on free time) works on the Python version (and I’ll try to help), the bash proof of concept works. Attached here.

After some initial testing, we think that GRASS’ results are just slightly more smooth for when applying the exact same parameters to get the HPFA Pan-Sharpened material.

To run a test, get the attached tarred file, untar, compile as a grass-addon against G64. Then, the simplest execution form will be:

note, the algorithm processes each low resolution image separately

i.fusion.hpf pan=lsat7_2000_80 msx=lsat7_2000_10

or, let’s get the usual landsat bands 1, 2 and 3

i.fusion.hpf pan=lsat7_2000_80 msx=lsat7_2000_10,lsat7_2000_20,lsat7_2000_30

some manual work to rebalance colors (required)

i.landsat.rgb r=hpf_lsat7_2000_30 g=hpf_lsat7_2000_20 b=hpf_lsat7_2000_10

d.mon x0

d.rgb r=hpf_lsat7_2000_30 g=hpf_lsat7_2000_20 b=hpf_lsat7_2000_10

increase crispness (not very obvious)

i.fusion.hpf pan=lsat7_2000_80 msx=lsat7_2000_10,lsat7_2000_20,lsat7_2000_30 --o modulator=max

i.landsat.rgb r=hpf_lsat7_2000_30 g=hpf_lsat7_2000_20 b=hpf_lsat7_2000_10

d.mon x1

d.rgb r=hpf_lsat7_2000_30 g=hpf_lsat7_2000_20 b=hpf_lsat7_2000_10

hmmm… let’s override the algorithms ratio computation

i.fusion.hpf pan=lsat7_2000_80 msx=lsat7_2000_10,lsat7_2000_20,lsat7_2000_30 --o modulator=max ratio=6

i.landsat.rgb r=hpf_lsat7_2000_30 g=hpf_lsat7_2000_20 b=hpf_lsat7_2000_10

d.mon x2

d.rgb r=hpf_lsat7_2000_30 g=hpf_lsat7_2000_20 b=hpf_lsat7_2000_10

BTW, it works really great with QuickBird imagery. Will upload in the next days screenshots with freely available HR imagery (i.e. QB2, IKONOS) and Landsat5, 7 and 8. Just a few Landsat 8 screenshots uploaded in the Wiki:

http://grasswiki.osgeo.org/wiki/File:Bengkalis_Island_Landsat8_RGB_Red_NIR_Blue.jpg

http://grasswiki.osgeo.org/wiki/File:Bengkalis_Island_Landsat8_Panchromatic.jpg

http://grasswiki.osgeo.org/wiki/File:Bengkalis_Island_Landsat8_HPFA_Sharpened_RGB_Red_NIR_Blue_Default_Parameters.jpg

http://grasswiki.osgeo.org/wiki/File:Bengkalis_Island_Landsat8_HPFA_Sharpened_RGB_Red_NIR_Blue_Max_Crispness.jpg

Nikos

(attachments)

i.fusion.tar.gz (8.26 KB)

8

Hamish:

what problems did you have with "i.fusion.hpf.tmp.$$" ?

Nikos:

I think the problem was (among my misunderstanding) that $$
was the same for a bunch of names fed as variables here and there.

'$$' gets replaced with the PID of the calling process, so within a single script it will be consistent, and it should be unique among processes currently running on the system.

Hamish

9

On 14/11/13 02:25, Nikos Alexandris wrote:

Dear GRASS GIS users,

together with Nikos Ves, we share the "i.fusion.hpf" idea/proof of
concept. At the moment, we have a custom shell script named
`i.fusion.hpf` (an attempt for a proper GRASS add-on), which implements
the High Pass Filter Additive (HPFA) Fusion Technique for Pan-Sharpening
[*]. Nikos V started already porting to Python. How can we proceed in
sharing it?

[...]

Two questions

? Can someone confirm that the part of the existing "i.pansharpen" code
that performs histogram matching (code lines 348 - 431), do so as
"linearly stretching an image to match another image's Mean and StdDev"?

AFAICT, it applies the method described in [1]. I don't know (and don't have the time to think about) what this method does in terms of mean and stddev.

? Would it be desired to get the HPFA algorithm integrated in i.pansharpen?

Yes. I think that if we have a generic module such as i.pansharpen, it would be preferable to have all pansharpening methods in that one module.

Mortiz

10

Nikos Alexandris wrote:

together with Nikos Ves, we share the “i.fusion.hpf” idea/proof of

concept. At the moment, we have a custom shell script named

i.fusion.hpf (an attempt for a proper GRASS add-on), which implements

the High Pass Filter Additive (HPFA) Fusion Technique for Pan-Sharpening

[*]. Nikos V started already porting to Python. How can we proceed in

sharing it?

[…]

Two questions

? Can someone confirm that the part of the existing “i.pansharpen” code

that performs histogram matching (code lines 348 - 431), do so as

“linearly stretching an image to match another image’s Mean and StdDev”?

Moritz Lennert:

AFAICT, it applies the method described in [1].

Is that a reference also indicated in i.pansharpen’s manual?

I don’t know (and don’t have the time to think about) what this method does

in terms of mean and stddev.

My guess was/is that it is not the same, i.e. it does not match Mean and StdDev. As a quick test, I tried the identical (me thinks) tool in WhiteboxGIS “Histogram Matching (Two Images)”, does not give identical Means and StdDevs after the operation – which is the case with i.pansharpen too if I am not wrong.

? Would it be desired to get the HPFA algorithm integrated in

i.pansharpen?

Yes. I think that if we have a generic module such as i.pansharpen, it

would be preferable to have all pansharpening methods in that one module.

There is one “difference” in that HPFA treats all bands to be sharpened separately. And, in this manner, it can be (mis-)used to sharpen any low-res band. For example, WorldView-2 products have 8 multi-spectral bands. Hence the “not red= green= blue=” design so far from my side.

@Michael, is this of interest for you?

Thank you Moritz, Nikos

11

On 15/11/13 10:50, Nikos Alexandris wrote:

Nikos Alexandris wrote:

> > together with Nikos Ves, we share the "i.fusion.hpf" idea/proof of

> > concept. At the moment, we have a custom shell script named

> > `i.fusion.hpf` (an attempt for a proper GRASS add-on), which implements

> > the High Pass Filter Additive (HPFA) Fusion Technique for
Pan-Sharpening

> > [*]. Nikos V started already porting to Python. How can we proceed in

> > sharing it?

>

> [...]

>

> > Two questions

> >

> > ? Can someone confirm that the part of the existing "i.pansharpen" code

> > that performs histogram matching (code lines 348 - 431), do so as

> > "linearly stretching an image to match another image's Mean and
StdDev"?

Moritz Lennert:

> AFAICT, it applies the method described in [1].

Is that a reference also indicated in i.pansharpen's manual?

It's not in the manual, but there's a long list of other references. But Michael is the one who knows where the inspiration came from. AFAIK, this is the classical, generic method of histogram matching.

> I don't know (and don't have the time to think about) what this
method does

> in terms of mean and stddev.

My guess was/is that it is not the same, i.e. it does not match Mean and
StdDev. As a quick test, I tried the identical (me thinks) tool in
WhiteboxGIS "Histogram Matching (Two Images)", does not give identical
Means and StdDevs after the operation -- which is the case with
i.pansharpen too if I am not wrong.

I just did a quick test:

pan in:

mean: 31.813
standard deviation: 3.75447

ms in:

mean: 15.2307
standard deviation: 3.55858

pan out:

mean: 15.6117
standard deviation: 3.23408

So for this example, mean seems to have been adjusted, but stddev not.

> > ? Would it be desired to get the HPFA algorithm integrated in

> > i.pansharpen?

>

> Yes. I think that if we have a generic module such as i.pansharpen, it

> would be preferable to have all pansharpening methods in that one module.

There is one "difference" in that HPFA treats all bands to be sharpened
separately. And, in this manner, it can be (mis-)used to sharpen any
low-res band. For example, WorldView-2 products have 8 multi-spectral
bands. Hence the "not red= green= blue=" design so far from my side.

i.pansharpen does not imply rgb either (although the description of the ms* parameters does suggest that. You can obviously use any ms bands you want.

Moritz

12

Moritz Lennert wrote:

…

I just did a quick test:

pan in:

mean: 31.813

standard deviation: 3.75447

ms in:

mean: 15.2307

standard deviation: 3.55858

pan out:

mean: 15.6117

standard deviation: 3.23408

So for this example, mean seems to have been adjusted, but stddev not.

Thanks. Looks not to be the exact same here too.

? Would it be desired to get the HPFA algorithm integrated in

i.pansharpen?

Yes. I think that if we have a generic module such as i.pansharpen, it

would be preferable to have all pansharpening methods in that one

module.

There is one “difference” in that HPFA treats all bands to be sharpened

separately. And, in this manner, it can be (mis-)used to sharpen any

low-res band. For example, WorldView-2 products have 8 multi-spectral

bands. Hence the “not red= green= blue=” design so far from my side.

i.pansharpen does not imply rgb either (although the description of the

ms* parameters does suggest that. You can obviously use any ms bands you

want.

Yes, I know. I just wanted to avoid this r= g= b= logic. Maybe I am wrong and perhaps it is better, indeed, to follow the common path.

13

Nikos Alexandris wrote:

BTW, it works really great with QuickBird imagery.

For the matter,

And, finally, slightly tweaked (other parameters offered by the algorithm/script): http://grasswiki.osgeo.org/wiki/File:RGB_HPF_Sharpened_Center_Low_Modulator_Min.jpg

Nikos

14

On Fri, Nov 15, 2013 at 10:05 AM, Moritz Lennert
<mlennert@club.worldonline.be> wrote:

On 14/11/13 02:25, Nikos Alexandris wrote:

Dear GRASS GIS users,

together with Nikos Ves, we share the "i.fusion.hpf" idea/proof of
concept. At the moment, we have a custom shell script named
`i.fusion.hpf` (an attempt for a proper GRASS add-on), which implements
the High Pass Filter Additive (HPFA) Fusion Technique for Pan-Sharpening
[*]. Nikos V started already porting to Python.

Great!

? Would it be desired to get the HPFA algorithm integrated in
i.pansharpen?

Yes. I think that if we have a generic module such as i.pansharpen, it would
be preferable to have all pansharpening methods in that one module.

+1
please integrate the method there...

Markus

15

I haven't had time today to look into this. I used a fairly standard histogram matching algorithm, for which I cited the reference. But I don't remember exactly how it went. Matching the SD was not part of it however.

Michael
______________________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University
Tempe, AZ 85287-2402
USA

voice: 480-965-6262 (SHESC), 480-965-8130/727-9746 (CSDC)
fax: 480-965-7671(SHESC), 480-727-0709 (CSDC)
www: http://csdc.asu.edu, http://shesc.asu.edu
    http://www.public.asu.edu/~cmbarton

On Nov 15, 2013, at 3:20 AM, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 15/11/13 10:50, Nikos Alexandris wrote:

Nikos Alexandris wrote:

> > together with Nikos Ves, we share the "i.fusion.hpf" idea/proof of

> > concept. At the moment, we have a custom shell script named

> > `i.fusion.hpf` (an attempt for a proper GRASS add-on), which implements

> > the High Pass Filter Additive (HPFA) Fusion Technique for
Pan-Sharpening

> > [*]. Nikos V started already porting to Python. How can we proceed in

> > sharing it?

>

> [...]

>

> > Two questions

> >

> > ? Can someone confirm that the part of the existing "i.pansharpen" code

> > that performs histogram matching (code lines 348 - 431), do so as

> > "linearly stretching an image to match another image's Mean and
StdDev"?

Moritz Lennert:

> AFAICT, it applies the method described in [1].

Is that a reference also indicated in i.pansharpen's manual?

It's not in the manual, but there's a long list of other references. But Michael is the one who knows where the inspiration came from. AFAIK, this is the classical, generic method of histogram matching.

> I don't know (and don't have the time to think about) what this
method does

> in terms of mean and stddev.

My guess was/is that it is not the same, i.e. it does not match Mean and
StdDev. As a quick test, I tried the identical (me thinks) tool in
WhiteboxGIS "Histogram Matching (Two Images)", does not give identical
Means and StdDevs after the operation -- which is the case with
i.pansharpen too if I am not wrong.

I just did a quick test:

pan in:

mean: 31.813
standard deviation: 3.75447

ms in:

mean: 15.2307
standard deviation: 3.55858

pan out:

mean: 15.6117
standard deviation: 3.23408

So for this example, mean seems to have been adjusted, but stddev not.

> > ? Would it be desired to get the HPFA algorithm integrated in

> > i.pansharpen?

>

> Yes. I think that if we have a generic module such as i.pansharpen, it

> would be preferable to have all pansharpening methods in that one module.

There is one "difference" in that HPFA treats all bands to be sharpened
separately. And, in this manner, it can be (mis-)used to sharpen any
low-res band. For example, WorldView-2 products have 8 multi-spectral
bands. Hence the "not red= green= blue=" design so far from my side.

i.pansharpen does not imply rgb either (although the description of the ms* parameters does suggest that. You can obviously use any ms bands you want.

Moritz

16

I agree that it would be best to have all pan sharpening algorithms together if possible. New ones should be addable as new classes or methods in the existing module. Note that I did employ parallel processing to the extent possible. It might be possible to apply this kind of processing to other sharpening algorithms.

Michael

C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University
Tempe, AZ 85287-2402
USA

voice: 480-965-6262 (SHESC), 480-965-8130/727-9746 (CSDC)
fax: 480-965-7671(SHESC), 480-727-0709 (CSDC)
www: http://csdc.asu.edu, http://shesc.asu.edu
http://www.public.asu.edu/~cmbarton

On Nov 15, 2013, at 8:20 AM, Nikos Alexandris <nik@nikosalexandris.net> wrote:

Moritz Lennert wrote:
…

I just did a quick test:

pan in:
mean: 31.813
standard deviation: 3.75447

ms in:
mean: 15.2307
standard deviation: 3.55858

pan out:
mean: 15.6117
standard deviation: 3.23408

So for this example, mean seems to have been adjusted, but stddev not.

Thanks. Looks not to be the exact same here too.

? Would it be desired to get the HPFA algorithm integrated in
i.pansharpen?

Yes. I think that if we have a generic module such as i.pansharpen, it
would be preferable to have all pansharpening methods in that one
module.

There is one “difference” in that HPFA treats all bands to be sharpened
separately. And, in this manner, it can be (mis-)used to sharpen any
low-res band. For example, WorldView-2 products have 8 multi-spectral
bands. Hence the “not red= green= blue=” design so far from my side.

i.pansharpen does not imply rgb either (although the description of the
ms* parameters does suggest that. You can obviously use any ms bands you
want.

Yes, I know. I just wanted to avoid this r= g= b= logic. Maybe I am wrong and perhaps it is better, indeed, to follow the common path.

17

I took a look at the i.pansharpen code. The method matchhist does the histogram matching. It creates cumulative distribution functions (CDF) of the source and target histograms and then finds the closest values to match at each point on the CDF. It is pretty thoroughly documented in the code. There are other methods of histogram matching, but IIRC, this was the most basic and widespread. As some others have commented, it assumes that images have 256 integer grey values. A more sophisticated histogram matching algorithm could utilize floating point values and a wider range of values. Hope this helps

Michael
____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

voice: 480-965-6262 (SHESC), 480-965-8130/727-9746 (CSDC)
fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC)
www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

On Nov 15, 2013, at 3:20 AM, Moritz Lennert <mlennert@club.worldonline.be> wrote:

On 15/11/13 10:50, Nikos Alexandris wrote:

Nikos Alexandris wrote:

together with Nikos Ves, we share the "i.fusion.hpf" idea/proof of

concept. At the moment, we have a custom shell script named

`i.fusion.hpf` (an attempt for a proper GRASS add-on), which implements

the High Pass Filter Additive (HPFA) Fusion Technique for

Pan-Sharpening

[*]. Nikos V started already porting to Python. How can we proceed in

sharing it?

[...]

Two questions

? Can someone confirm that the part of the existing "i.pansharpen" code

that performs histogram matching (code lines 348 - 431), do so as

"linearly stretching an image to match another image's Mean and

StdDev"?

Moritz Lennert:

AFAICT, it applies the method described in [1].

Is that a reference also indicated in i.pansharpen's manual?

It's not in the manual, but there's a long list of other references. But
Michael is the one who knows where the inspiration came from. AFAIK,
this is the classical, generic method of histogram matching.

I don't know (and don't have the time to think about) what this

method does

in terms of mean and stddev.

My guess was/is that it is not the same, i.e. it does not match Mean and
StdDev. As a quick test, I tried the identical (me thinks) tool in
WhiteboxGIS "Histogram Matching (Two Images)", does not give identical
Means and StdDevs after the operation -- which is the case with
i.pansharpen too if I am not wrong.

I just did a quick test:

pan in:

mean: 31.813
standard deviation: 3.75447

ms in:

mean: 15.2307
standard deviation: 3.55858

pan out:

mean: 15.6117
standard deviation: 3.23408

So for this example, mean seems to have been adjusted, but stddev not.

? Would it be desired to get the HPFA algorithm integrated in

i.pansharpen?

Yes. I think that if we have a generic module such as i.pansharpen, it

would be preferable to have all pansharpening methods in that one module.

There is one "difference" in that HPFA treats all bands to be sharpened
separately. And, in this manner, it can be (mis-)used to sharpen any
low-res band. For example, WorldView-2 products have 8 multi-spectral
bands. Hence the "not red= green= blue=" design so far from my side.

i.pansharpen does not imply rgb either (although the description of the
ms* parameters does suggest that. You can obviously use any ms bands you
want.

Moritz

18

On Sat, Nov 16, 2013 at 9:27 PM, Michael Barton <Michael.Barton@asu.edu> wrote:

I took a look at the i.pansharpen code. The method matchhist does the histogram matching. It creates cumulative distribution functions (CDF) of the source and target histograms and then finds the closest values to match at each point on the CDF. It is pretty thoroughly documented in the code. There are other methods of histogram matching, but IIRC, this was the most basic and widespread. As some others have commented, it assumes that images have 256 integer grey values. A more sophisticated histogram matching algorithm could utilize floating point values and a wider range of values. Hope this helps

How does it compare to the one used in the Addon "i.hist.match"?
(grass-addons/grass7/imagery/i.histo.match/i.histo.match.py)

There I don't see a 8bit limitation (I may be wrong). This might solve
ticket #2048.

Markus

19

Michael Barton wrote:

> I took a look at the i.pansharpen code. The method matchhist does the
> histogram matching. It creates cumulative distribution functions (CDF) of
> the source and target histograms and then finds the closest values to
> match at each point on the CDF. It is pretty thoroughly documented in the
> code. There are other methods of histogram matching, but IIRC, this was
> the most basic and widespread. As some others have commented, it assumes
> that images have 256 integer grey values. A more sophisticated histogram
> matching algorithm could utilize floating point values and a wider range
> of values. Hope this helps

Markus Neteler wrote:

How does it compare to the one used in the Addon "i.hist.match"?
(grass-addons/grass7/imagery/i.histo.match/i.histo.match.py)

My rough guess is that it is about the same logic. In the case of pansharpen
there is one reference cdf. In the case of i.histo.match there are some
assumptions I guess and some averaged values are used as a reference. See also
Moritz' comment:

" 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. "

There I don't see a 8bit limitation (I may be wrong). This might solve
ticket #2048.

Note, the IHS method for example depends on the respective modules which are
8-bit based too if I am not wrong.

Nikos

20

On Nov 16, 2013, at 5:46 PM, Nikos Alexandris <nik@nikosalexandris.net>
wrote:

" 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 also looked at i.histo.match when I was writing i.pansharpen. IIRC, its original goal was to create an average histogram that could be applied to multiple images tiled together, so that they would have a nice consistent appearance without patching them into a single image. It is not the kind of histogram matching needed for pan sharpening, however.

I do not have an image set that is >8bit and floating point for testing. But someone on the list did not too long ago.

Any histogram matching algorithm (including a variant of i.histo.match) could be added as a drop-in method to i.pansharpen. All it needs is to take an "original" and "target" input files (the target is the file whose histogram is modified to match the original), and a "matched" output file.

Michael

There I don't see a 8bit limitation (I may be wrong). This might solve
ticket #2048.

Note, the IHS method for example depends on the respective modules which are
8-bit based too if I am not wrong.

Nikos

____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

voice: 480-965-6262 (SHESC), 480-965-8130/727-9746 (CSDC)
fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC)
www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu