[gdal-dev] Boundary Artifact with ChunkAndWarpImage

[Even Rouault]

Le vendredi 17 avril 2015 20:00:36, Michael Aschenbeck a écrit :
> Hey guys,
> 
> I realize this is borderline impossible to debug with what I've given you.
> So I managed to reproduce it with the C++ code at the bottom of the email.
> The image i used is:
> https://www.dropbox.com/s/i9fxstyuf8eqe8t/testsingleband.tif?dl=0
> 
> Something that is a little weird is that I cropped this image so recreating
> the problem is less clunky. When I did that, the artifacts were not
> present!  So, that's the reason i'm including a 65mb image.  (Perhaps Matt
> was on to something when he thought about gdal chunking bigger images.)
>  Using this code and this image, you should see the first nonzero row is
> too dark, and the second row is too bright.
> 
> Let me know if anyone can recreate the issue and/or has any further ideas.

I can now reproduce and finally found the reason. When you do :

psWarpOptions->padfSrcNoDataReal = (double
*)CPLMalloc(psWarpOptions->nBandCount*sizeof(double));

nBandCount has not yet been set, so it is at its default value 0, and there's 
no nodata declared...

If you do the following instead, it will work:

psWarpOptions->nBandCount = GDALGetRasterCount(poDataset);
psWarpOptions->padfSrcNoDataReal = (double
*)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));
psWarpOptions->padfSrcNoDataImag = (double
*)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));

(note that since you did CreateCopy() before, the original lines will remain 
above the area touched by your warped imagery. So you might want to erase 
everything with ((GDALDataset*)hDstDS)->GetRasterBand(1)->Fill(0); before)

> 
> Thanks again,
> Mike
> 
> 
> 
> int pushUpToRight(void *hTransforArg, int bDstToSrc, int nPointCount,
> double *x, double *y, double *z, int *panSuccess)
> {
> //First transform to georeferenced units:
> GDALGenImgProjTransform(hTransforArg, 0, nPointCount, x, y, z, panSuccess);
> //Now we are in georeferenced units
> 
>         //push up and to the right
> for (int i = 0; i<nPointCount; i++)
> {
> if (panSuccess[i] != 0)
> {
> x[i] = x[i] + .01;
> y[i] = y[i] + .651;
> }
> }
>  double *geo = new double[6]; //hardcoded for this example
> geo[0] = 451581.6; geo[1] = 0.4; geo[2] = 0.0; geo[3] =
> 2420520.8000000003; geo[4]
> = 0.0; geo[5] = -0.4;
> 
> //Now convert to pixels
> for (int i = 0; i < nPointCount; i++)
> {
> x[i] = convertXToPixel(x[i], geo);
> y[i] = convertYToPixel(y[i], geo);
> }
> 
> return 1;
> }
> double convertXToPixel(double xGeo, double *adfGeoTransform) { return
> ((xGeo - adfGeoTransform[0]) / adfGeoTransform[1]);}
> double convertYToPixel(double yGeo, double *adfGeoTransform) { return
> ((yGeo - adfGeoTransform[3]) / adfGeoTransform[5]);}
> 
> int main()
> {
> GDALDataset *poDataset;
> GDALDriver *hDriver;
> GDALDataset *hDstDS;
> 
> GDALAllRegister();
> 
> // Open the source file
> poDataset = (GDALDataset *)GDALOpen("E:/Test/TestSingleBand.tif",
> GA_ReadOnly);
> CPLAssert(poDataset != NULL);
> 
> // Get output driver (GeoTIFF format)
> hDriver = (GDALDriver*)GDALGetDriverByName("GTiff");
> CPLAssert(hDriver != NULL);
> 
> // Create the output file.
> hDstDS = hDriver->CreateCopy("E:/Test/out.tif",
> poDataset,FALSE,NULL,NULL,NULL);
> CPLAssert(hDstDS != NULL);
> 
> //Create the transformer
> void *hTransformArgForward;
> hTransformArgForward = GDALCreateGenImgProjTransformer2(poDataset, NULL,
> NULL);
> 
>         //Create the warper and set various variables
> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
> 
> psWarpOptions->padfSrcNoDataReal = (double
> *)CPLMalloc(psWarpOptions->nBandCount*sizeof(double));
> for (int ii = 0; ii < psWarpOptions->nBandCount; ii++)
> {
> psWarpOptions->padfSrcNoDataReal[ii] = 0;
> }
> psWarpOptions->hSrcDS = poDataset;
> psWarpOptions->hDstDS = hDstDS;
> psWarpOptions->nBandCount = GDALGetRasterCount(poDataset);
> psWarpOptions->pTransformerArg = hTransformArgForward;
> psWarpOptions->panSrcBands = (int *)CPLMalloc(sizeof(int));
> psWarpOptions->panSrcBands[0] = 1;
> psWarpOptions->panDstBands = (int *)CPLMalloc(sizeof(int));
> psWarpOptions->panDstBands[0] = 1;
> psWarpOptions->eResampleAlg = GRA_Cubic;
> psWarpOptions->pfnProgress = GDALTermProgress;
> 
> // Establish my transformer.
>     GDALTransformerFunc myTransformer = &pushUpToRight;
> psWarpOptions->pfnTransformer = myTransformer;
> 
> // Initialize and execute the warp operation.
> GDALWarpOperation oOperation;
> oOperation.Initialize(psWarpOptions);
> oOperation.ChunkAndWarpImage(0, 0,
> GDALGetRasterXSize(hDstDS),
> GDALGetRasterYSize(hDstDS));
> 
>         //clean
> GDALDestroyWarpOptions(psWarpOptions);
> GDALClose(hDstDS);
> GDALClose(poDataset);
> 
> return 0;
> }
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> On Thu, Apr 16, 2015 at 1:10 PM, Michael Aschenbeck <
> 
> m.g.aschenbeck at gmail.com> wrote:
> > Some good ideas, thanks for the suggestions.
> > 
> > I added padfDstNoDataReal, as well as the imaginary values.  I also set
> > INIT_DEST to NODATA, although from what I read that was more of an option
> > that affected speed.  Unfortunately my output was identical to the
> > previous run and the artifacts remain present.
> > 
> > On Thu, Apr 16, 2015 at 11:59 AM, Matt Hanson <matt.a.hanson at gmail.com>
> > 
> > wrote:
> >> Mike,
> >> Are you setting padfDstNoDataReal as well?     Because GDAL performs the
> >> warp in chunks I'm wondering if previously written data, that should be
> >> NoData, is not recognized as such with later chunks.   However, I'm
> >> still not clear on the exact details of how that works internally, Even
> >> would have a better idea if this even makes sense.
> >> 
> >> The code I'm using to set WarpOptions for doing warping on multiple
> >> input images to a single outout looks like the below.   Also, are you
> >> setting INIT_DEST option to "NODATA" (it's one of the key-value options
> >> that are assigned to psWarpOptions->papszWarpOptions).
> >> 
> >> Hope this helps!
> >> 
> >> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
> >> psWarpOptions->hDstDS = imgout.GetGDALDataset();
> >> psWarpOptions->nBandCount = imgout.NumBands();
> >> psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) *
> >> psWarpOptions->nBandCount ); psWarpOptions->panDstBands = (int *)
> >> CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
> >> psWarpOptions->padfSrcNoDataReal = (double *) CPLMalloc(sizeof(double) *
> >> psWarpOptions->nBandCount ); psWarpOptions->padfSrcNoDataImag = (double
> >> *) CPLMalloc(sizeof(double) * psWarpOptions->nBandCount );
> >> psWarpOptions->padfDstNoDataReal = (double *) CPLMalloc(sizeof(double) *
> >> psWarpOptions->nBandCount ); psWarpOptions->padfDstNoDataImag = (double
> >> *) CPLMalloc(sizeof(double) * psWarpOptions->nBandCount ); for (unsigned
> >> int b=0;b<imgout.NumBands();b++) { psWarpOptions->panSrcBands[b] = b+1;
> >> psWarpOptions->panDstBands[b] = b+1; psWarpOptions->padfSrcNoDataReal[b]
> >> = images[0][b].NoDataValue(); psWarpOptions->padfDstNoDataReal[b] =
> >> imgout[b].NoDataValue(); psWarpOptions->padfSrcNoDataImag[b] = 0.0;
> >> psWarpOptions->padfDstNoDataImag[b] = 0.0; }
> >> 
> >> On Thu, Apr 16, 2015 at 1:42 PM, Michael Aschenbeck <
> >> 
> >> m.g.aschenbeck at gmail.com> wrote:
> >>> Thanks for your reply.
> >>> 
> >>> I'm using gdal 1.10.1.  I'm using cubic resampling, but i have noticed
> >>> the dark line artifact with bilinear as well.  I unfortunately do not
> >>> reproduce the problem with a basic gdalwarp with cubic resampling.
> >>> 
> >>> I have written this code to do some registration. It's hard to tell in
> >>> the example, but my transformation actually pushed pixels up a bit. The
> >>> transformation I use handles interior pixels nicely and extends beyond
> >>> the image, so there is nothing 'weird' going on with that near the
> >>> boundary that I can think of.
> >>> 
> >>> I was hoping there is a warp option i have not yet thought of.
> >>> 
> >>> Thanks again,
> >>> Mike
> >>> 
> >>> 
> >>> 
> >>> On Thu, Apr 16, 2015 at 11:10 AM, Even Rouault <
> >>> 
> >>> even.rouault at spatialys.com> wrote:
> >>>> Le jeudi 16 avril 2015 19:01:57, Michael Aschenbeck a écrit :
> >>>> > Hello,
> >>>> > 
> >>>> > I'm using ChunkAndWarpMulti to warp an image.  The warping is
> >>>> > working nicely, however, at the boundary I seem to be getting some
> >>>> > artifacts.
> >>>> > 
> >>>> > The first artifact i see is a DARK boundary of pixels in some
> >>>> 
> >>>> location.  My
> >>>> 
> >>>> > guess is that the interpolator is interpolating with blackfill
> >>>> 
> >>>> (intensity
> >>>> 
> >>>> > zero pixels).  Note that I am using the following setup:
> >>>> > psWarpOptions->padfSrcNoDataReal = (double *)
> >>>> > CPLMalloc(psWarpOptions->nBandCount*sizeof(double));
> >>>> > for (int ii = 0; ii < psWarpOptions->nBandCount; ii++)
> >>>> > {
> >>>> > psWarpOptions->padfSrcNoDataReal[ii] = 0;
> >>>> > }
> >>>> > which i thought was supposed to treat zeros as nodata.  It doesn't
> >>>> 
> >>>> seem to
> >>>> 
> >>>> > be doing what I think it should.
> >>>> > 
> >>>> > In some cases, I am also noticing a BRIGHT band strip of pixels
> >>>> 
> >>>> adjacent to
> >>>> 
> >>>> > the dark boundary pixels.  I don't have any thoughts on where this
> >>>> 
> >>>> artifact
> >>>> 
> >>>> > is coming from.
> >>>> > 
> >>>> > Below you can find links to an example.  The orange you see is just
> >>>> 
> >>>> the
> >>>> 
> >>>> > background color of my viewer.  Zero pixels are set to transparent
> >>>> > so
> >>>> 
> >>>> you
> >>>> 
> >>>> > can see the dark boundary artifact.  (Note that all of the orange
> >>>> 
> >>>> region is
> >>>> 
> >>>> > covered with zero-intensity pixels, so we haven't reached the
> >>>> 
> >>>> boundary of
> >>>> 
> >>>> > the file, just the boundary of the non-zero pixels.)  Sorry if
> >>>> > that's confusing.
> >>>> 
> >>>> > The before image:
> >>>> https://www.dropbox.com/s/3kmyi8wu0qybsq9/before.JPG?dl=0
> >>>> 
> >>>> > The after image with the artifacts:
> >>>> > https://www.dropbox.com/s/fo0m8q95b26s61m/after.JPG?dl=0
> >>>> 
> >>>> Mike,
> >>>> 
> >>>> Which gdal version ? Which resampling method ?
> >>>> Do you reproduce with gdalwarp ?
> >>>> If so, providing input file + full gdalwarp command line would help.
> >>>> 
> >>>> Even
> >>>> 
> >>>> --
> >>>> Spatialys - Geospatial professional services
> >>>> http://www.spatialys.com
> >>> 
> >>> _______________________________________________
> >>> gdal-dev mailing list
> >>> gdal-dev at lists.osgeo.org
> >>> http://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com