[gdal-dev] Boundary Artifact with ChunkAndWarpImage

Fri Apr 17 13:43:18 PDT 2015

That did the trick!  Sorry to burden the mailing list with a stupid error
like this.  However, I really appreciate you taking the time to help.

Best,
Mike

On Fri, Apr 17, 2015 at 1:42 PM, Even Rouault <even.rouault at spatialys.com>
wrote:

> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150417/dbfd9b87/attachment-0001.html>