<div dir="ltr">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.<br><br>Best,<br>Mike</div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Apr 17, 2015 at 1:42 PM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@spatialys.com" target="_blank">even.rouault@spatialys.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">Le vendredi 17 avril 2015 20:00:36, Michael Aschenbeck a écrit :<br>
> Hey guys,<br>
><br>
> I realize this is borderline impossible to debug with what I've given you.<br>
> So I managed to reproduce it with the C++ code at the bottom of the email.<br>
> The image i used is:<br>
> <a href="https://www.dropbox.com/s/i9fxstyuf8eqe8t/testsingleband.tif?dl=0" target="_blank">https://www.dropbox.com/s/i9fxstyuf8eqe8t/testsingleband.tif?dl=0</a><br>
><br>
> Something that is a little weird is that I cropped this image so recreating<br>
> the problem is less clunky. When I did that, the artifacts were not<br>
> present!  So, that's the reason i'm including a 65mb image.  (Perhaps Matt<br>
> was on to something when he thought about gdal chunking bigger images.)<br>
>  Using this code and this image, you should see the first nonzero row is<br>
> too dark, and the second row is too bright.<br>
><br>
> Let me know if anyone can recreate the issue and/or has any further ideas.<br>
<br>
</span>I can now reproduce and finally found the reason. When you do :<br>
<span class=""><br>
psWarpOptions->padfSrcNoDataReal = (double<br>
*)CPLMalloc(psWarpOptions->nBandCount*sizeof(double));<br>
<br>
</span>nBandCount has not yet been set, so it is at its default value 0, and there's<br>
no nodata declared...<br>
<br>
If you do the following instead, it will work:<br>
<br>
psWarpOptions->nBandCount = GDALGetRasterCount(poDataset);<br>
psWarpOptions->padfSrcNoDataReal = (double<br>
*)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));<br>
psWarpOptions->padfSrcNoDataImag = (double<br>
*)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));<br>
<br>
(note that since you did CreateCopy() before, the original lines will remain<br>
above the area touched by your warped imagery. So you might want to erase<br>
everything with ((GDALDataset*)hDstDS)->GetRasterBand(1)->Fill(0); before)<br>
<div class="HOEnZb"><div class="h5"><br>
><br>
> Thanks again,<br>
> Mike<br>
><br>
><br>
><br>
> int pushUpToRight(void *hTransforArg, int bDstToSrc, int nPointCount,<br>
> double *x, double *y, double *z, int *panSuccess)<br>
> {<br>
> //First transform to georeferenced units:<br>
> GDALGenImgProjTransform(hTransforArg, 0, nPointCount, x, y, z, panSuccess);<br>
> //Now we are in georeferenced units<br>
><br>
>         //push up and to the right<br>
> for (int i = 0; i<nPointCount; i++)<br>
> {<br>
> if (panSuccess[i] != 0)<br>
> {<br>
> x[i] = x[i] + .01;<br>
> y[i] = y[i] + .651;<br>
> }<br>
> }<br>
>  double *geo = new double[6]; //hardcoded for this example<br>
> geo[0] = 451581.6; geo[1] = 0.4; geo[2] = 0.0; geo[3] =<br>
> 2420520.8000000003; geo[4]<br>
> = 0.0; geo[5] = -0.4;<br>
><br>
> //Now convert to pixels<br>
> for (int i = 0; i < nPointCount; i++)<br>
> {<br>
> x[i] = convertXToPixel(x[i], geo);<br>
> y[i] = convertYToPixel(y[i], geo);<br>
> }<br>
><br>
> return 1;<br>
> }<br>
> double convertXToPixel(double xGeo, double *adfGeoTransform) { return<br>
> ((xGeo - adfGeoTransform[0]) / adfGeoTransform[1]);}<br>
> double convertYToPixel(double yGeo, double *adfGeoTransform) { return<br>
> ((yGeo - adfGeoTransform[3]) / adfGeoTransform[5]);}<br>
><br>
> int main()<br>
> {<br>
> GDALDataset *poDataset;<br>
> GDALDriver *hDriver;<br>
> GDALDataset *hDstDS;<br>
><br>
> GDALAllRegister();<br>
><br>
> // Open the source file<br>
> poDataset = (GDALDataset *)GDALOpen("E:/Test/TestSingleBand.tif",<br>
> GA_ReadOnly);<br>
> CPLAssert(poDataset != NULL);<br>
><br>
> // Get output driver (GeoTIFF format)<br>
> hDriver = (GDALDriver*)GDALGetDriverByName("GTiff");<br>
> CPLAssert(hDriver != NULL);<br>
><br>
> // Create the output file.<br>
> hDstDS = hDriver->CreateCopy("E:/Test/out.tif",<br>
> poDataset,FALSE,NULL,NULL,NULL);<br>
> CPLAssert(hDstDS != NULL);<br>
><br>
> //Create the transformer<br>
> void *hTransformArgForward;<br>
> hTransformArgForward = GDALCreateGenImgProjTransformer2(poDataset, NULL,<br>
> NULL);<br>
><br>
>         //Create the warper and set various variables<br>
> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();<br>
><br>
> psWarpOptions->padfSrcNoDataReal = (double<br>
> *)CPLMalloc(psWarpOptions->nBandCount*sizeof(double));<br>
> for (int ii = 0; ii < psWarpOptions->nBandCount; ii++)<br>
> {<br>
> psWarpOptions->padfSrcNoDataReal[ii] = 0;<br>
> }<br>
> psWarpOptions->hSrcDS = poDataset;<br>
> psWarpOptions->hDstDS = hDstDS;<br>
> psWarpOptions->nBandCount = GDALGetRasterCount(poDataset);<br>
> psWarpOptions->pTransformerArg = hTransformArgForward;<br>
> psWarpOptions->panSrcBands = (int *)CPLMalloc(sizeof(int));<br>
> psWarpOptions->panSrcBands[0] = 1;<br>
> psWarpOptions->panDstBands = (int *)CPLMalloc(sizeof(int));<br>
> psWarpOptions->panDstBands[0] = 1;<br>
> psWarpOptions->eResampleAlg = GRA_Cubic;<br>
> psWarpOptions->pfnProgress = GDALTermProgress;<br>
><br>
> // Establish my transformer.<br>
>     GDALTransformerFunc myTransformer = &pushUpToRight;<br>
> psWarpOptions->pfnTransformer = myTransformer;<br>
><br>
> // Initialize and execute the warp operation.<br>
> GDALWarpOperation oOperation;<br>
> oOperation.Initialize(psWarpOptions);<br>
> oOperation.ChunkAndWarpImage(0, 0,<br>
> GDALGetRasterXSize(hDstDS),<br>
> GDALGetRasterYSize(hDstDS));<br>
><br>
>         //clean<br>
> GDALDestroyWarpOptions(psWarpOptions);<br>
> GDALClose(hDstDS);<br>
> GDALClose(poDataset);<br>
><br>
> return 0;<br>
> }<br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
> On Thu, Apr 16, 2015 at 1:10 PM, Michael Aschenbeck <<br>
><br>
> <a href="mailto:m.g.aschenbeck@gmail.com">m.g.aschenbeck@gmail.com</a>> wrote:<br>
> > Some good ideas, thanks for the suggestions.<br>
> ><br>
> > I added padfDstNoDataReal, as well as the imaginary values.  I also set<br>
> > INIT_DEST to NODATA, although from what I read that was more of an option<br>
> > that affected speed.  Unfortunately my output was identical to the<br>
> > previous run and the artifacts remain present.<br>
> ><br>
> > On Thu, Apr 16, 2015 at 11:59 AM, Matt Hanson <<a href="mailto:matt.a.hanson@gmail.com">matt.a.hanson@gmail.com</a>><br>
> ><br>
> > wrote:<br>
> >> Mike,<br>
> >> Are you setting padfDstNoDataReal as well?     Because GDAL performs the<br>
> >> warp in chunks I'm wondering if previously written data, that should be<br>
> >> NoData, is not recognized as such with later chunks.   However, I'm<br>
> >> still not clear on the exact details of how that works internally, Even<br>
> >> would have a better idea if this even makes sense.<br>
> >><br>
> >> The code I'm using to set WarpOptions for doing warping on multiple<br>
> >> input images to a single outout looks like the below.   Also, are you<br>
> >> setting INIT_DEST option to "NODATA" (it's one of the key-value options<br>
> >> that are assigned to psWarpOptions->papszWarpOptions).<br>
> >><br>
> >> Hope this helps!<br>
> >><br>
> >> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();<br>
> >> psWarpOptions->hDstDS = imgout.GetGDALDataset();<br>
> >> psWarpOptions->nBandCount = imgout.NumBands();<br>
> >> psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) *<br>
> >> psWarpOptions->nBandCount ); psWarpOptions->panDstBands = (int *)<br>
> >> CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );<br>
> >> psWarpOptions->padfSrcNoDataReal = (double *) CPLMalloc(sizeof(double) *<br>
> >> psWarpOptions->nBandCount ); psWarpOptions->padfSrcNoDataImag = (double<br>
> >> *) CPLMalloc(sizeof(double) * psWarpOptions->nBandCount );<br>
> >> psWarpOptions->padfDstNoDataReal = (double *) CPLMalloc(sizeof(double) *<br>
> >> psWarpOptions->nBandCount ); psWarpOptions->padfDstNoDataImag = (double<br>
> >> *) CPLMalloc(sizeof(double) * psWarpOptions->nBandCount ); for (unsigned<br>
> >> int b=0;b<imgout.NumBands();b++) { psWarpOptions->panSrcBands[b] = b+1;<br>
> >> psWarpOptions->panDstBands[b] = b+1; psWarpOptions->padfSrcNoDataReal[b]<br>
> >> = images[0][b].NoDataValue(); psWarpOptions->padfDstNoDataReal[b] =<br>
> >> imgout[b].NoDataValue(); psWarpOptions->padfSrcNoDataImag[b] = 0.0;<br>
> >> psWarpOptions->padfDstNoDataImag[b] = 0.0; }<br>
> >><br>
> >> On Thu, Apr 16, 2015 at 1:42 PM, Michael Aschenbeck <<br>
> >><br>
> >> <a href="mailto:m.g.aschenbeck@gmail.com">m.g.aschenbeck@gmail.com</a>> wrote:<br>
> >>> Thanks for your reply.<br>
> >>><br>
> >>> I'm using gdal 1.10.1.  I'm using cubic resampling, but i have noticed<br>
> >>> the dark line artifact with bilinear as well.  I unfortunately do not<br>
> >>> reproduce the problem with a basic gdalwarp with cubic resampling.<br>
> >>><br>
> >>> I have written this code to do some registration. It's hard to tell in<br>
> >>> the example, but my transformation actually pushed pixels up a bit. The<br>
> >>> transformation I use handles interior pixels nicely and extends beyond<br>
> >>> the image, so there is nothing 'weird' going on with that near the<br>
> >>> boundary that I can think of.<br>
> >>><br>
> >>> I was hoping there is a warp option i have not yet thought of.<br>
> >>><br>
> >>> Thanks again,<br>
> >>> Mike<br>
> >>><br>
> >>><br>
> >>><br>
> >>> On Thu, Apr 16, 2015 at 11:10 AM, Even Rouault <<br>
> >>><br>
> >>> <a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>> wrote:<br>
> >>>> Le jeudi 16 avril 2015 19:01:57, Michael Aschenbeck a écrit :<br>
> >>>> > Hello,<br>
> >>>> ><br>
> >>>> > I'm using ChunkAndWarpMulti to warp an image.  The warping is<br>
> >>>> > working nicely, however, at the boundary I seem to be getting some<br>
> >>>> > artifacts.<br>
> >>>> ><br>
> >>>> > The first artifact i see is a DARK boundary of pixels in some<br>
> >>>><br>
> >>>> location.  My<br>
> >>>><br>
> >>>> > guess is that the interpolator is interpolating with blackfill<br>
> >>>><br>
> >>>> (intensity<br>
> >>>><br>
> >>>> > zero pixels).  Note that I am using the following setup:<br>
> >>>> > psWarpOptions->padfSrcNoDataReal = (double *)<br>
> >>>> > CPLMalloc(psWarpOptions->nBandCount*sizeof(double));<br>
> >>>> > for (int ii = 0; ii < psWarpOptions->nBandCount; ii++)<br>
> >>>> > {<br>
> >>>> > psWarpOptions->padfSrcNoDataReal[ii] = 0;<br>
> >>>> > }<br>
> >>>> > which i thought was supposed to treat zeros as nodata.  It doesn't<br>
> >>>><br>
> >>>> seem to<br>
> >>>><br>
> >>>> > be doing what I think it should.<br>
> >>>> ><br>
> >>>> > In some cases, I am also noticing a BRIGHT band strip of pixels<br>
> >>>><br>
> >>>> adjacent to<br>
> >>>><br>
> >>>> > the dark boundary pixels.  I don't have any thoughts on where this<br>
> >>>><br>
> >>>> artifact<br>
> >>>><br>
> >>>> > is coming from.<br>
> >>>> ><br>
> >>>> > Below you can find links to an example.  The orange you see is just<br>
> >>>><br>
> >>>> the<br>
> >>>><br>
> >>>> > background color of my viewer.  Zero pixels are set to transparent<br>
> >>>> > so<br>
> >>>><br>
> >>>> you<br>
> >>>><br>
> >>>> > can see the dark boundary artifact.  (Note that all of the orange<br>
> >>>><br>
> >>>> region is<br>
> >>>><br>
> >>>> > covered with zero-intensity pixels, so we haven't reached the<br>
> >>>><br>
> >>>> boundary of<br>
> >>>><br>
> >>>> > the file, just the boundary of the non-zero pixels.)  Sorry if<br>
> >>>> > that's confusing.<br>
> >>>><br>
> >>>> > The before image:<br>
> >>>> <a href="https://www.dropbox.com/s/3kmyi8wu0qybsq9/before.JPG?dl=0" target="_blank">https://www.dropbox.com/s/3kmyi8wu0qybsq9/before.JPG?dl=0</a><br>
> >>>><br>
> >>>> > The after image with the artifacts:<br>
> >>>> > <a href="https://www.dropbox.com/s/fo0m8q95b26s61m/after.JPG?dl=0" target="_blank">https://www.dropbox.com/s/fo0m8q95b26s61m/after.JPG?dl=0</a><br>
> >>>><br>
> >>>> Mike,<br>
> >>>><br>
> >>>> Which gdal version ? Which resampling method ?<br>
> >>>> Do you reproduce with gdalwarp ?<br>
> >>>> If so, providing input file + full gdalwarp command line would help.<br>
> >>>><br>
> >>>> Even<br>
> >>>><br>
> >>>> --<br>
> >>>> Spatialys - Geospatial professional services<br>
> >>>> <a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a><br>
> >>><br>
> >>> _______________________________________________<br>
> >>> gdal-dev mailing list<br>
> >>> <a href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
> >>> <a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
<br>
--<br>
Spatialys - Geospatial professional services<br>
<a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a><br>
</div></div></blockquote></div><br></div>