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 <[email protected]> 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 < > > > > [email protected]> 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 <[email protected] > > > > > > > > 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 < > > >> > > >> [email protected]> 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 < > > >>> > > >>> [email protected]> 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 > > >>> [email protected] > > >>> http://lists.osgeo.org/mailman/listinfo/gdal-dev > > -- > Spatialys - Geospatial professional services > http://www.spatialys.com >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
