Ticket #2137 (assigned enhancement)

Opened 7 months ago

Last modified 3 weeks ago

Gaussian pyramid in gdal

Reported by: warmerdam Assigned to: warmerdam (accepted)
Priority: normal Milestone: 1.6.0
Component: GDAL_Raster Version: 1.4.2
Severity: normal Keywords: overview
Cc: krzystek@hm.edu, rouault

Description

From Peter Krzystek:

Frank,

it took me a while to find time to do something on that subject. Meanwhile I am working with gdal version 1.4.2 and have successfully added the 3x3 Gaussian kernel. I have tested it with a large number of sample Tiff files (tiled/untiled) covering the following formats. 

1x8 (uncompressed, lzw, packbits)
1x8 Pallete (uncompressed, lzw, packbits)
3x8 (uncompressed, lzw, packbits)
4x8 (uncompressed, lzw, packbits)
1x16 (uncompressed, lzw, packbits)
3x16 (uncompressed, lzw, packbits)
4x16 (uncompressed, lzw, packbits)

All the uncompressed formats - either tiled or untiled - show no problems. However, as you probably know, lzw and packbits tiff formats do not work properly. The only workaround in this case - as suggested in bug report 1205 (

http://trac.osgeo.org/gdal/query?status=%21closed&type=defect&order=id&desc=1

) is to convert the image to an uncompressed image, do the image pyramid, and to use geotiffcp to convert it back to either lzw or packbits.

Please find attached the piece of code that includes the 3x3 Gaussian. The change affects only the file "overview.cpp" and the functions "GDALRegenerateOverviews" and "GDALDownsampleChunk32RAs". Since I am not an experienced contributor to gdal code I send you the file as a compressed document. I hope this is fine. If not, could you please give me a short hint how I can get a user ID etc. and should proceed.

If you can accept the code change please let me know.

Best regards

Peter

P.S. The implementation of the Gaussian kernel follows the code for the option "AVER" that does the filtering in each scanline. Thus, no decomposition of the linear kernel is applied that would reduce the number of instructions (i.e. multiplications). However, since the kernel is only 3x3 this little drawback will not cause a severe time delay. If the size of the kernel (e.g. 3,5, etc.) should be an option, another parameter needs to be passed to the function "GDALDownsampleChunk32RAs".

Attachments

overview.patch (6.5 kB) - added by warmerdam on 01/07/08 12:41:14.
Updated patch version against trunk.

Change History

01/07/08 12:40:29 changed by warmerdam

  • status changed from new to assigned.

Patch replaced with one against trunk. This version of the patch also changes nCount to a double precision value (dfCountWeight) for summing the gaussian weights which I think is necessary for correct operation.

01/07/08 12:41:14 changed by warmerdam

  • attachment overview.patch added.

Updated patch version against trunk.

01/07/08 12:44:11 changed by warmerdam

  • cc changed from krzystek@hm.edu to krzystek@hm.edu, rouault.

Patch applied as r in trunk which includes some white space adjustments.

We still need user documentation (at the api and gdaladdo level), and some sort of automated testing.

I also think the gaussian case is missing nodata tests added to the average case by Even since 1.4.2. Adding Even as a cc in case he would like to add this.

01/07/08 12:44:50 changed by warmerdam

The last change was r13485.

01/07/08 18:06:09 changed by rouault

Copy&paste from a discussion on IRC 

[23:36] <even__> FrankW: about #2137, I don't understand why it's necessary to use double for the gauss weights. The weight matrix only contains integers
[23:36] <gdaltrac> Ticket #2137: Gaussian pyramid in gdal, http://trac.osgeo.org/gdal/ticket/2137
[23:36] *** cgs_bob_ est maintenant connu sous le nom de cgs_bob.
[23:37] <-- _adam_gfx a quitté ce serveur (" HydraIRC -> http://www.hydrairc.com <- IRC has never been so good").
[23:40] <FrankW> Hmm, why is gaussMatrix of type double?  I hadn't inspected that, but I now think the count should be int, and the gaussmatrix should be int.
[23:40] <FrankW> I made the change without careful inspection based on the fact that the array was double.
[23:40] <FrankW> feel free to "int-ify" things if you would like.
[23:40] <even__> no problem
[23:41] <even__> but I also think that there are issues with the indexes
[23:41] <even__> See http://rafb.net/p/7ialX277.html
[23:42] <FrankW> even__:  Hmm, that looks bad.
[23:43] <even__> something related to nSrcYOff, nSrcYOff2, nSrcXOff, nSrcXOff2 I think
[23:45] <even__> hum, apart from that problem, I'm not even sure that it can work
[23:45] <FrankW> Feel free to note the conditions in the ticket (or in a new one) and we can bounce it back to the author.
[23:45] <FrankW> Unless you are inclined to work on it.
[23:45] <even__> Gauss kernel are used for blurring effects, aren't they ?
[23:46] <even__> but I'm not aware we can use them, as simply as proposed, for reducing the size of a raster
[23:46] <FrankW> I'm not clear on why a gaussian kernel would be good for overviews.
[23:46] <FrankW> This is why "patches welcome" is not my mantra.  Patches often involve a lot of maintainer work to consider, fix and incorporate. :-)
[23:51] --> mloskot a rejoint ce canal (n=mloskot@osgeo/member/mloskot).
[23:53] <even__> FrankW: I don't know how the author wanted it to work. I see two possibilities : either you take the nearest pixel and it's neighbours and apply the 3x3 kernel on them
[23:53] <even__> either you consider the source area we use for average and compute a bigger kernel, and apply on that area
[23:53] <even__> (hope I'm clear)
[23:53] <even__> the current implementation is somewhere in the middle...
[23:55] <even__> solution 1 should be easy to fix, but I doubt it makes much sense, except for 1/3 reduction
[23:55] <even__> solution 2 is probably more interesting (I'm not a signal processing expert...), but harder to implement

The warnings mentionned : 

$ apps/gdaladdo -r gauss small_gauss.tif 2 4 8

0...10...20...ERROR 5: Access window out of range in RasterIO().  Requested
(0,128) of size 256x129 on raster of 256x256.
30ERROR 5: Access window out of range in RasterIO().  Requested
(0,0) of size 128x129 on raster of 128x128.
...40...50...60ERROR 5: Access window out of range in RasterIO().  Requested
(0,128) of size 256x129 on raster of 256x256.
..ERROR 5: Access window out of range in RasterIO().  Requested
(0,0) of size 128x129 on raster of 128x128.
.70...80...90..ERROR 5: Access window out of range in RasterIO().  Requested
(0,128) of size 256x129 on raster of 256x256.
.ERROR 5: Access window out of range in RasterIO().  Requested
(0,0) of size 128x129 on raster of 128x128.
100 - done.

Same command with average resampling works just fine. GeoTIFF is 512x512.

07/05/08 04:54:24 changed by rouault

Frank, I see that you have documented the GAUSS method in r14815, but I'm wondering if it is the right time as AFAIK it is still broken.

07/05/08 13:21:05 changed by warmerdam

Good point, I have removed the docs again in r14827.