[gdal-dev] Boundary Artifact with ChunkAndWarpImage

Fri Apr 17 11:00:36 PDT 2015

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.

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