Opened 13 years ago

Closed 12 years ago

Last modified 12 years ago

#2137 closed enhancement (fixed)

Gaussian pyramid in gdal

Reported by: warmerdam Owned by: Even Rouault
Priority: normal Milestone: 1.6.0
Component: GDAL_Raster Version: 1.4.2
Severity: normal Keywords: overview
Cc: krzystek@…, Even Rouault, warmerdam@…

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 (7)

overview.patch (6.5 KB) - added by warmerdam 13 years ago.
Updated patch version against trunk.
overview_1.patch (6.9 KB) - added by twinwizard 12 years ago.
New patch of overview.cpp against trunk
tiff_ovr.patch (2.1 KB) - added by twinwizard 12 years ago.
patch of tiff_ovr.py against trunk
overview.cpp (60.9 KB) - added by twinwizard 12 years ago.
overview.cpp.gauss (60.3 KB) - added by Even Rouault 12 years ago.
geotiff.cpp (233.9 KB) - added by twinwizard 12 years ago.
gt_overview.cpp (23.7 KB) - added by twinwizard 12 years ago.

Download all attachments as: .zip

Change History (51)

comment:1 Changed 13 years ago by warmerdam

Status: newassigned

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.

Changed 13 years ago by warmerdam

Attachment: overview.patch added

Updated patch version against trunk.

comment:2 Changed 13 years ago by warmerdam

Cc: Even Rouault added

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.

comment:3 Changed 13 years ago by warmerdam

The last change was r13485.

comment:4 Changed 13 years ago by Even 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.

comment:5 Changed 12 years ago by Even 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.

comment:6 Changed 12 years ago by warmerdam

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

comment:7 Changed 12 years ago by Even Rouault

I've decided to revert in r15048 the code for the gaussian algorithm as I don't still see the point of it for overview building and as the issues I've raised have not been fixed.

Gaussian filtering could be interesting for image processing, but when you build overviews, I don't understand why you would want to give more weight to some pixels than to others : a small shift of 1 pixel in the source window can potentially lead to significant differences in the result overview. That seems really arbitrary. Well, nearest neighbour is arbitrary but it makes sense (it is implemented in a deterministic way, but we could as well implement it by selecting a random pixel in the source window. That would lead roughly to the same result). Average gives the same weight to all the pixels in the source window.

What could maybe make sense is to apply gaussian filtering on the full resolution dataset, and then compute overviews with averaging.

comment:8 in reply to:  7 Changed 12 years ago by twinwizard

Cc: warmerdam@… added

Dear rouault, thank you for asking. I have build a new version in which your remarks are fixed. I have also tested it with lots of sample images.

Replying to rouault:

I've decided to revert in r15048 the code for the gaussian algorithm as I don't still see the point of it for overview building and as the issues I've raised have not been fixed.

Would it be acceptable for you if I would add the Gaussian option again to the code? I would also contribute docs and sample images. I could do this and finalize it within the next 2 weeks. Please let me know your comments. See comments on your remarks below. Best regards Peter

Gaussian filtering could be interesting for image processing, but when you build overviews, I don't understand why you would want to give more weight to some pixels than to others : a small shift of 1 pixel in the source window can potentially lead to significant differences in the result overview. That seems really arbitrary. Well, nearest neighbour is arbitrary but it makes sense (it is implemented in a deterministic way, but we could as well implement it by selecting a random pixel in the source window. That would lead roughly to the same result). Average gives the same weight to all the pixels in the source window.

It can be clearly shown (any textbook from image processing would be fine) that a simple averaging leads to so called subsampling effects (=Moire effects). Instead, the Gaussian approach suppresses image frequencies which are beyond the so-called Nyquist frequency, so that a subsampling (which is the case when creating overviews) is not affected by these image frequencies. Furthermore, if GDAL is applied for photogrammetric and remote sensing purposes Gaussian overviews (as an option) are needed anyway. So, Gaussian overviews are meanwhile a standard since they are theoretically correct. Instead, a simple averaging is just an option when the overviews are solely used for visualization.

What could maybe make sense is to apply gaussian filtering on the full resolution dataset, and then compute overviews with averaging.

This would be an option to smooth the full resolution image with a Gaussian kernel. As far as the overviews are concerned see comments above.

comment:9 Changed 12 years ago by Even Rouault

Peter,

feel free to attach a new patch (against recent trunk please, as there were quite a few changes in overview.cpp recently) that adds again the option. I'll try to review it. You seem much better aware of signal processing than me, so I trust you on the usefullness of such a resampling.

But the size of the kernel should be somehow linked to the resampling factor, right ? This wasn't the case in the first version.

Even

comment:10 in reply to:  9 Changed 12 years ago by twinwizard

Even, thank you for cooperation. I will do my best in the next days. Usually, the size of the kernel is either 3x3 or 5x5 assuming a resampling factor of 2. If another resampling factor is desired the size of the kernel must change. Best regards Peter

Replying to rouault:

Peter,

feel free to attach a new patch (against recent trunk please, as there were quite a few changes in overview.cpp recently) that adds again the option. I'll try to review it. You seem much better aware of signal processing than me, so I trust you on the usefullness of such a resampling.

But the size of the kernel should be somehow linked to the resampling factor, right ? This wasn't the case in the first version.

Even

comment:11 Changed 12 years ago by twinwizard

Even, I have just finished my internal tests on the Gaussian filter. Which images should I use from gdal for testing? I can find a lot in the autotest directory. Furthermore, where should automated testing and documentation be placed? Best regards Peter

comment:12 Changed 12 years ago by Even Rouault

You can add tests in autotest/gcore/tiff_ovr.py

As far as the images, in addition to some "real-world images" that I guess you have tried, with GDALCreate(), you can generate simply images of different sizes, of different block sizes (tiled/stripped for a GeoTIFF), with different resampling factor, to test how your algorithm behaves on edges. Difficult to say exactly what to test...

comment:13 in reply to:  12 Changed 12 years ago by twinwizard

I am ready in setting up tests for Gaussian overviews based on sample images. How should I proceed? Should I attach a new patch for tiff_ovr.py? And how should I add the images. Just adding the images to data of gcore? Futhermore, I still have problems in building and installing python bindings from gdal trunk using cygwin. All recommendations I found did not work. Do you a know a good documentation? I used the python bindings for python25 I found on http://pypi.python.org/pypi/GDAL/. Best regards Peter

Replying to rouault:

You can add tests in autotest/gcore/tiff_ovr.py

As far as the images, in addition to some "real-world images" that I guess you have tried, with GDALCreate(), you can generate simply images of different sizes, of different block sizes (tiled/stripped for a GeoTIFF), with different resampling factor, to test how your algorithm behaves on edges. Difficult to say exactly what to test...

comment:14 Changed 12 years ago by Even Rouault

I'm not a specialist of python bindings setup. I've never tryied on cygwin. That's not the easiest environment to setup I guess! not sure this is supported in fact. But you could probably use the python bindings for GDAL 1.5.2 against GDAL trunk. That should work (at least, it works for me on Linux. Not all the autotest from GDAL trunk will work against bindings of 1.5.2 as there were some additions/fixes, but that should mostly work )

Yes a patch to tiff_ovr.py is fine.

It would be great if you could reuse existing images in gcore/data or gdrivers/data. We don't really want to add big images to keep the archive size reasonnably small. But as I said you can generate on-the-fly test images of various dimensions with the Create() method in the gcore/tmp directory.

Try to do your best but don't exhaust yourself on setting up the test part.

Changed 12 years ago by twinwizard

Attachment: overview_1.patch added

New patch of overview.cpp against trunk

Changed 12 years ago by twinwizard

Attachment: tiff_ovr.patch added

patch of tiff_ovr.py against trunk

comment:15 Changed 12 years ago by twinwizard

Even, I have added patches for overview.cpp and tiff_ovr.py. The latter ist located in autotest/gcore. I have tested the patch for overview.cpp again with geo images covering the 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) In all cases I got successful results meaning that the overviews where built without any errors. Since the images are to large they cannot be uploaded to gdal. If you want access to the images I can arrange a ftp access. The patch of tiff_ovr.py includes two tests on images located in gcore/data. Since this images are rather small they can only be used for formal testing. If one wants to compare AVERAGE with GAUSSIAN filtering one needs large images. I hope everything is fine. Let me know if something does not work. If documentation is needed I will do it. Let me know what to do. Regards Peter

comment:16 Changed 12 years ago by Even Rouault

Hum, I don't manage to apply your patch, probably because there's too bit of reformatting and 'patch' get lost.

Could you try to genereate it by doing "diff -u overview.cpp.ori overview.cpp > overview.cpp.patch", so that it can get applied by "patch -p1 < overview.cpp.patch". If that doesn't work in your tree, then attach the whole file.

Just by looking at it, I don't see any guarantee that the 'j' index in the loop doing the real stuff will not exceed 2. Correleted to that, I don't understand the nSrcYOff2 computation, and why you need a +1 on nFullResYChunk. For a resampling factor of 2, this probably works, but with other resampling factors, that will likely fail/crash.

The gaussMatrix[][] should be a static const array. No need to initialize it all the times.

comment:17 Changed 12 years ago by Even Rouault

The GDALRegenerateOverviewsMultiBand() should be supported too. But my previous comments must be addressed first.

comment:18 in reply to:  16 Changed 12 years ago by twinwizard

I have added as an "patch" the complete new overview.cpp, since diff -u etc. did not work for some reason. Replying to rouault:

Hum, I don't manage to apply your patch, probably because there's too bit of reformatting and 'patch' get lost.

Could you try to genereate it by doing "diff -u overview.cpp.ori overview.cpp > overview.cpp.patch", so that it can get applied by "patch -p1 < overview.cpp.patch". If that doesn't work in your tree, then attach the whole file.

Just by looking at it, I don't see any guarantee that the 'j' index in the loop doing the real stuff will not exceed 2. Correleted to that, I don't understand the nSrcYOff2 computation, and why you need a +1 on nFullResYChunk. For a resampling factor of 2, this probably works, but with other resampling factors, that will likely fail/crash.

According your remarks I have rearranged the code as follows:

  • nFullResYChunk is not enlarged by 1 anymore.
  • The index j is now controlled so that it does not exceed the dimension of the gaussian matrix dim. The dimension of the gaussian matrix is variable depending on the resampling factor.
  • The computation of nSrcYOff2 is now the same as with AVERAGING.
  • The algorithm works now with other resampling factors. However, as far as Gaussian pyramids are concerned, the factor 2 is usually applied. It would make sense to add this to the help text, for instance in gdaladdo.

    The gaussMatrix[][] should be a static const array. No need to initialize it all the times.

I have added 3 gaussian matrices as static const array.

comment:19 Changed 12 years ago by Even Rouault

Your updated file is not against svn trunk.

comment:20 Changed 12 years ago by Even Rouault

Code indentation should also be coherent with the indentation of the rest of the file, that is to say only spaces, not tabulations. It would be good to use UNIX end of lines, so that I can easily workaround that with dos2unix...

comment:21 in reply to:  19 Changed 12 years ago by twinwizard

Replying to rouault:

Your updated file is not against svn trunk.

I used the current overview.cpp which can be checked out as https://svn.osgeo.org/gdal/branches/1.5/gdal/gcore/overview.cpp. I assume that you mean that I did not update (commit) the new code using svn or - for instance- TortoiseSVN. If yes, I assumed that I have to attach at first a patch for overview.cpp.

comment:22 in reply to:  20 Changed 12 years ago by twinwizard

Replying to rouault:

Code indentation should also be coherent with the indentation of the rest of the file, that is to say only spaces, not tabulations. It would be good to use UNIX end of lines, so that I can easily workaround that with dos2unix...

I have removed all the tabulations which came from my VisualC++ editor.

comment:23 in reply to:  19 Changed 12 years ago by twinwizard

Replying to rouault:

Your updated file is not against svn trunk.

Forgot to say that I have attched a new overview.cpp

comment:24 Changed 12 years ago by Even Rouault

Ah sorry, I must clarify what I meant by "svn trunk"... It's the term we use to speak about the development version which lies at http://svn.osgeo.org/gdal/trunk, whereas http://svn.osgeo.org/gdal/branches/1.5 is not the trunk, but a maintenance branch. All new functionality must be added to developement version aka svn trunk. Maintenance branch is only for bug fixing. (We didn't really invent that term of "trunk" and "branches". If you read the subversion manual, it's the proposed structure of the repository)

So you must checkout http://svn.osgeo.org/gdal/trunk and adapt your patch to work with that version of overview.cpp.

We could certainly do better in explaining this process in this page of the wiki (http://trac.osgeo.org/gdal/wiki/HowToContribute) ...

comment:25 in reply to:  24 Changed 12 years ago by twinwizard

Replying to rouault:

Ah sorry, I must clarify what I meant by "svn trunk"... It's the term we use to speak about the development version which lies at http://svn.osgeo.org/gdal/trunk, whereas http://svn.osgeo.org/gdal/branches/1.5 is not the trunk, but a maintenance branch. All new functionality must be added to developement version aka svn trunk. Maintenance branch is only for bug fixing. (We didn't really invent that term of "trunk" and "branches". If you read the subversion manual, it's the proposed structure of the repository)

So you must checkout http://svn.osgeo.org/gdal/trunk and adapt your patch to work with that version of overview.cpp.

When compiling trunk under windows xp and using nmake.opt compiling fails and no gdal.dll is created. I used the same compiling options as for branches.

We could certainly do better in explaining this process in this page of the wiki (http://trac.osgeo.org/gdal/wiki/HowToContribute) ...

comment:26 Changed 12 years ago by Even Rouault

I don't use windows so I can hardly help you. But when you say "I used the same compiling options as for branches", do you mean that you used nmake.opt from branches/1.5 or that you somehow merged your changes in it. If you don't manage to build, you should probably explain your problem here and/or the gdal-dev mailing list with the actual error message so that someone can give some directions where to go...

comment:27 in reply to:  26 Changed 12 years ago by twinwizard

Replying to rouault:

I don't use windows so I can hardly help you. But when you say "I used the same compiling options as for branches", do you mean that you used nmake.opt from branches/1.5 or that you somehow merged your changes in it. If you don't manage to build, you should probably explain your problem here and/or the gdal-dev mailing list with the actual error message so that someone can give some directions where to go...

I just used the nmake.opt which comes along with trunk. Even when no changes are made to nmake.opt gdal is not compiling correctly. If I compile gdal from branches compiling works fine. The first error messages is:

D:\temp\testTrunc\test\frmts>cd gtiff && nmake /nologo /f makefile.vc && cd ..

cl /nologo /MD /EHsc /Ox /W3 /D_CRT_SECURE_NO_DEPRECATE /D_CRT_NONSTDC_NO_DEPRECATE /DNDEBUG -I..\..\port -I..\..\ogr -I..\..\gcore -I..\..\alg -I..\..\ogr\ogrsf_frmts -Ilibtiff -DBIGTIFF_SUPPORT -Ilibgeotiff /c geotiff.cpp

geotiff.cpp geotiff.cpp(2684) : error C2374: 'iBand' : Neudefinition; Mehrfachinitialisierung

geotiff.cpp(2634) : Siehe Deklaration von 'iBand'

which means translated that iBand is newly defined. This is strange since iBand is defined as an int in a for loop.

nmake is invoked as: nmake -f makefile.vc MSVC_VER=1400 The error occurs both with debug flaf on and off. I am using xp and VC2005. I have also VC6.0 on my machine. However, this seems to be no problem with gdal from branches.

If I am editing geotiff.cp and fixing the error the next compile ends with error messages coming from jpef stuff. Any idea? Sorry.

comment:28 Changed 12 years ago by Even Rouault

Well, ok, thanks to your precise report for the geotiff.cpp, I think I fixed the issue in r15302 for that point. But please be as precise with your jpeg issue to help me fixing stuff. (I repeat again : I have no VC compiler to experiment!). The root of the problem is that VC6 is not strictly compliant with the C++ standard and it doesn't support the 'int i' in the second for loop :

for(int i=0;i<N;i++)
{
}
for(int i=0;i<N;i++)
{
}

I'd suppose that the same pattern happens with the jpeg stuff

comment:29 in reply to:  28 Changed 12 years ago by twinwizard

Replying to rouault:

Well, ok, thanks to your precise report for the geotiff.cpp, I think I fixed the issue in r15302 for that point. But please be as precise with your jpeg issue to help me fixing stuff. (I repeat again : I have no VC compiler to experiment!). The root of the problem is that VC6 is not strictly compliant with the C++ standard and it doesn't support the 'int i' in the second for loop :

for(int i=0;i<N;i++)
{
}
for(int i=0;i<N;i++)
{
}

I'd suppose that the same pattern happens with the jpeg stuff

This is the error message which happens in tif_ojpeg.c:

tif_ojpeg.c(599) : error C2632: 'long' gefolgt von 'long' ist unzulaessig tif_ojpeg.c(601) : error C2632: 'long' gefolgt von 'long' ist unzulaessig tif_ojpeg.c(606) : error C2632: 'long' gefolgt von 'long' ist unzulaessig tif_ojpeg.c(613) : error C2632: 'long' gefolgt von 'long' ist unzulaessig tif_ojpeg.c(620) : error C2632: 'long' gefolgt von 'long' ist unzulaessig

which means that type (long long) is not accepted.

When editing the corresponding lines I end up in numerous linking errors (=unresolved external symbols). I have cut the entire list because of space reasons.

cd .. cd vb6 nmake /nologo /f makefile.vc cd .. if exist gdal.lib del gdal.lib lib /nologo /out:gdal.lib port\*.obj gcore\*.obj alg\*.obj frmts\o\*.obj ogr\ogrsf_frmts\ogrsf_frmts.lib ogr\ogr.lib vb6\vb6_support.obj link /nologo /dll /INCLUDE:_OGRFeatureStylePuller /INCLUDE:_OSRValidate /INCLUDE:_OPTGetProjectionMethods /INCLUDE:_OGR_G_GetPointCount /INCLUDE:_OGRRegisterAll /INCLUDE:_GDALSimpleImageWarp@36 /INCLUDE:_GDALReprojectImage@48 /INCLUDE:_GDALComputeMedianCutPCT@32 /INCLUDE:_GDALDitherRGB2PCT@28 /INCLUDE:_OCTNewCoordinateTransformation@8 /INCLUDE:_vbSafeArrayToPtr@16 port\*.obj gcore\*.obj alg\*.obj frmts\o\*.obj ogr\ogrsf_frmts\ogrsf_frmts.lib ogr\ogr.lib vb6\vb6_support.obj odbc32.lib odbccp32.lib user32.lib gcore\Version.res /out:gdal16dev.dll /implib:gdal_i.lib /debug

Bibliothek gdal_i.lib und Objekt gdal_i.exp wird erstellt

cpl_minizip_unzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflateInit2_ cpl_vsil_gzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflateInit2_ cpl_minizip_unzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _crc32 cpl_vsil_gzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _crc32 cpl_minizip_unzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflate cpl_vsil_gzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflate cpl_minizip_unzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflateEnd cpl_vsil_gzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflateEnd cpl_vsil_gzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflateCopy cpl_vsil_gzip.obj : error LNK2001: Nichtaufgeloestes externes Symbol _inflateReset gdaldefaultoverviews.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GTIFFBuildOverviews gdaljp2metadata.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GTIFWktFromMemBuf gdaljp2metadata.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GTIFMemBufFromWkt gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_BLX gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_ADRG gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_GXF gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_USGSDEM gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_RIK gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_LCP gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_DIPEx .... .... _GDALRegister_RPFTOC gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_NITF gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_GTiff gdalallregister.obj : error LNK2001: Nichtaufgeloestes externes Symbol _GDALRegister_VRT gdal16dev.dll : fatal error LNK1120: 78 unaufgeloeste externe Verweise

comment:30 Changed 12 years ago by Even Rouault

Owner: changed from warmerdam to Even Rouault
Status: assignednew

Hum, for your last VC6 issue, maybe report it on the mailing list so that someone more competent than me can help you. I'm skeptical about "long long" not being handled by VC6, as I can see "long long" was already being used at least in frmts/gtiff/libtiff/tif_read.c of 1.5 branch....

comment:31 in reply to:  30 Changed 12 years ago by twinwizard

Resolution: fixed
Status: newclosed

Replying to rouault:

Hum, for your last VC6 issue, maybe report it on the mailing list so that someone more competent than me can help you. I'm skeptical about "long long" not being handled by VC6, as I can see "long long" was already being used at least in frmts/gtiff/libtiff/tif_read.c of 1.5 branch....

I fortunately fixed my compiler problems, so that I could successfully compile trunk. The "long long" problem in tif_ojpeg.c disappeared when using VC2005. I have attched the new overview.cpp as a "patch". This new overview.cpp is updated against svn trunk.

Changed 12 years ago by twinwizard

Attachment: overview.cpp added

comment:32 Changed 12 years ago by Even Rouault

Resolution: fixed
Status: closedreopened

Reopening (the building problems do not really close the main subject of the ticket ;-))

comment:33 Changed 12 years ago by Even Rouault

ok, I tried your patch in my tree. Apart from stylistic changes&improvements, there were 2 or 3 bugs I've fixed. I attach my modified version as overview.cpp.

I just tested in on 2 or 3 examples and the only difference I could see (with my eyes) between gauss and average was that there's something like a 1 pixel shift in x direction and y direction between both algorithms. Is that intended ?

Can you really prove me with an example that, at least in some cases, your implementation of Gauss overview leads to better results than averaging ? I need to be really convinced that there's enough added value to justify the extra maintenance that any new code introduces.

Changed 12 years ago by Even Rouault

Attachment: overview.cpp.gauss added

comment:34 Changed 12 years ago by Even Rouault

Finally my modified version is called "overview.cpp.gauss" as I can't overwrite a file attached by someone else.

comment:35 in reply to:  33 Changed 12 years ago by twinwizard

Replying to rouault:

ok, I tried your patch in my tree. Apart from stylistic changes&improvements, there were 2 or 3 bugs I've fixed. I attach my modified version as overview.cpp.

I just tested in on 2 or 3 examples and the only difference I could see (with my eyes) between gauss and average was that there's something like a 1 pixel shift in x direction and y direction between both algorithms. Is that intended ?

Good observation, however, this is quite clear. When averaging an image with a 2x2 kernel, the pixel location of the new pixel in the overview refers to the center of the 2x2 window. Thus, if you think about a georeferenced image, the new overview is not only scaled by the factor 2 but als shifted by half of pixel (Probably not corrected in GDAL?). In case of the (current) Gaussian implementation the algorithm starts from the pixel location (1,1) (assuming (0,0) is the center of upper left pixel) and jumps with a jump interval of 2. Thus, the Gaussian overviews are shifted by half a pixel against the averaged overviews. Usually, the pixel shift of overviews should be corrected if precise measuring in the overviews is intended. This is independent from the way how overviews are created.

Can you really prove me with an example that, at least in some cases, your implementation of Gauss overview leads to better results than averaging? I need to be really convinced that there's enough added value to justify the extra maintenance that any new code introduces.

First theoretically. Each image contains frequencies. If a overviews is created from an image, this is equivalent to a subsampling. The Nyquist theorem says, that the subsampling frequency (that is the jump interval of 2) must be at least twice as large as the highest frequency in an image. However, usually an image contains high frequencies. Thus, the image must be smoothed (for instance with a Gaussian kernel) in order to reduce the high frequencies so that the subsampling is correct.

Practically, overviews are created in professional remote sensing software packages (the list is long) by smoothing at first the image with a 3x3 or 5x5 kernel. The second step is to take every second pixel for the overview. Nothing else is done in the current overview patch. If you like I can provide you sample images on a ftp server that show you that the Gaussian approach is superior to the averaging approach. The images are to large to upload to GDAL. Do you have another idea?

comment:36 in reply to:  34 Changed 12 years ago by twinwizard

Replying to rouault:

Finally my modified version is called "overview.cpp.gauss" as I can't overwrite a file attached by someone else.

How can I download the file? Just by copy and paste?

comment:37 Changed 12 years ago by Even Rouault

Resolution: fixed
Status: reopenedclosed

Peter,

Thanks for your explanations. I've finally commited the file in SVN in r15323. Someone on IRC looked at the code and expressed doubts about the correctness of the computation of the source pixel.

I would have personnaly think about the following :

nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
nSrcXOff = nSrcXOff - nGaussMatrixDim/2;
nSrcXOff2 = nSrcXOff + nGaussMatrixDim;

I'd be still interested to test it with sample images if you can upload some on a ftp/http server

As far as downloading attachments in Trac is concerned, you can notice at the bottom of the page of the attachment a link called "Download in other formats: original format".

comment:38 in reply to:  37 ; Changed 12 years ago by twinwizard

Replying to rouault:

Peter,

Thanks for your explanations. I've finally commited the file in SVN in r15323. Someone on IRC looked at the code and expressed doubts about the correctness of the computation of the source pixel.

I would have personnaly think about the following :

nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
nSrcXOff = nSrcXOff - nGaussMatrixDim/2;
nSrcXOff2 = nSrcXOff + nGaussMatrixDim;

I'd be still interested to test it with sample images if you can upload some on a ftp/http server

I will set the ftp sever next week. Hopefully, it will work on Monday. I will let you know.

As far as downloading attachments in Trac is concerned, you can notice at the bottom of the page of the attachment a link called "Download in other formats: original format".

comment:39 in reply to:  38 Changed 12 years ago by twinwizard

Replying to twinwizard:

Replying to rouault:

Peter,

Thanks for your explanations. I've finally commited the file in SVN in r15323. Someone on IRC looked at the code and expressed doubts about the correctness of the computation of the source pixel.

I would have personnaly think about the following :

nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
nSrcXOff = nSrcXOff - nGaussMatrixDim/2;
nSrcXOff2 = nSrcXOff + nGaussMatrixDim;

I'd be still interested to test it with sample images if you can upload some on a ftp/http server

I will set the ftp sever next week. Hopefully, it will work on Monday. I will let you know.

As far as downloading attachments in Trac is concerned, you can notice at the bottom of the page of the attachment a link called "Download in other formats: original format".

You can download a sample image and information about tests I made with this imagery. ftp access: photo.geo.fhm.edu user: gdal PW: gauss1

The download will be possible till Wednesday. Please read the readme.txt for more information.

comment:40 Changed 12 years ago by Even Rouault

I don't manage to reach that ftp server. "ping photo.geo.fhm.edu" returns unknown host. I also tried ftp.photo.geo.fhm.edu, etc, but no luck.

Note: if you want to pass some information only privately, you can email me at "even dot rouault at mines dash paris dot org"

comment:41 in reply to:  17 Changed 12 years ago by twinwizard

Replying to rouault:

The GDALRegenerateOverviewsMultiBand() should be supported too. But my previous comments must be addressed first.

Is there any urgent need to go ahead and support GDALRegenerateOverviewsMultiBand with Gaussian overview option?

comment:42 Changed 12 years ago by Even Rouault

I would not call that a "urgent" need but something interesting though. For that, in addition to overview.cpp, the relevant files to edit are frmt/gtiff/geotiff.cpp and frmt/gtiff/gt_overview.cpp as far as I remember. I think adding gaussian support should be rather trivial. But maybe we will observe some artifacts on tile boundaries due to the fact that we don't fetch data across tile. (We could probably already observe such artifacts on horizontal chuncks)

Side note. I've realized recently that there was a workaround to the implemention of gaussian overview. One way to do it is to set a VRT filter with the appropriate kernel values on the full resolution dataset (see http://trac.osgeo.org/gdal/browser/trunk/autotest/gdrivers/data/avfilt_nodata.vrt as an example of an average filter) and then compute nearest overview on that VRT. This should give more or less the same results.

Changed 12 years ago by twinwizard

Attachment: geotiff.cpp added

Changed 12 years ago by twinwizard

Attachment: gt_overview.cpp added

comment:43 in reply to:  42 Changed 12 years ago by twinwizard

Replying to rouault:

I would not call that a "urgent" need but something interesting though. For that, in addition to overview.cpp, the relevant files to edit are frmt/gtiff/geotiff.cpp and frmt/gtiff/gt_overview.cpp as far as I remember. I think adding gaussian support should be rather trivial. But maybe we will observe some artifacts on tile boundaries due to the fact that we don't fetch data across tile. (We could probably already observe such artifacts on horizontal chuncks)

I have added geotiff.cpp and gt_overview.cpp as attachments. They contain - as promised - only minor changes to apply for the Gaussian option. Also, I have successfully tested it with the image pyramid.tif with the options --config COMPRESS_OVERVIEW JPEG --config INTERLEAVE_OVERVIEW PIXEL.

I noticed that the help text of gdaladdo does noct contain the hint that Gaussian filtering is possible.

Side note. I've realized recently that there was a workaround to the implemention of gaussian overview. One way to do it is to set a VRT filter with the appropriate kernel values on the full resolution dataset (see http://trac.osgeo.org/gdal/browser/trunk/autotest/gdrivers/data/avfilt_nodata.vrt as an example of an average filter) and then compute nearest overview on that VRT. This should give more or less the same results.

I agree that it would be more or less the same. BTW, the kernel in the example is a box filter usually used for low pass filtering of an image.

comment:44 Changed 12 years ago by Even Rouault

Multiband case commited in r15414 (test in r15415) with adidtion of gauss parameter in gdaladdo text.

Note: See TracTickets for help on using tickets.