[gdal-dev] GDALWarp API and paletted images
Thomas Sevaldrud
thomas at silentwings.no
Fri Jun 12 06:15:59 PDT 2015
I have those nodata-values there since the same code path is used for
elevation data and such, but it had no effect to remove it.
I have narrowed it down to the case where I use an overview as input. I
have some very large images with overviews, so I find the overview closest
matching to the output resolution of my image and make a VRT of the
overview level like this:
vrtDS = vrtDriver->Create("", ovrW, ovrH, _numBands, _dataType, NULL);
GDALSetProjection(vrtDS, _srcProjectionWKT.c_str());
geoTransform[XFM_PIXEL_SIZE_EW] = scaleX*basePixelSizeX;
geoTransform[XFM_PIXEL_SIZE_NS] = scaleY*basePixelSizeY;
vrtDS->SetGeoTransform(geoTransform);
for(int i = 1; i <= _numBands; i++)
{
GDALRasterBand* srcRootBand = _ds->GetRasterBand(i);
GDALRasterBand* srcBand = srcRootBand->GetOverview(overviewId);
VRTSourcedRasterBand* vrtBand =
(VRTSourcedRasterBand*)vrtDS->GetRasterBand(i);
vrtBand->AddSimpleSource(srcBand);
}
srcDS = vrtDS;
The strange thing is that this used to work when I had the same mistake as
you pointed out, where I added the bands twice, i.e the overview had 6
bands, but only three or four of them were used in the warp operation. When
fixing the code, I got this problem.
Is this the best way to do it btw, or is it a better way to make use of
overviews when cutting out and reprojecting parts of a big image?
Thanks again for all your help :-)
- Thomas
On Fri, Jun 12, 2015 at 1:35 PM, Even Rouault <even.rouault at spatialys.com>
wrote:
> Selon Thomas Sevaldrud <thomas at silentwings.no>:
>
> > Thanks again! The problem was in the warper setup, it used only the first
> > band of the image. This also revealed a bug in another part of my code
> > related to VRT's and Warp, since I had added double sets of bands there
> as
> > well. However, when fixing this, some old code stopped working :-)
> >
> > Basically, my problem is that I want to reproject an RGB image which may
> be
> > the VRT I made above, or another RGB(A) image, into an output RGBA. If I
> > don't have an input alpha channel, I still want the output to have alpha,
> > and I want any borders that result from the reprojection to have alpha =
> 0.
> > What is the correct way of setting up this? Currently I set up the bands
> > like this:
> >
> > warpOptions->nBandCount = srcNumBands;
> >
> > warpOptions->panSrcBands =
> > (int*)CPLMalloc(warpOptions->nBandCount*sizeof(int));
> > warpOptions->panDstBands =
> > (int*)CPLMalloc(warpOptions->nBandCount*sizeof(int));
> >
> > warpOptions->padfSrcNoDataReal =
> > (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));
> > warpOptions->padfSrcNoDataImag =
> > (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));
> >
> > warpOptions->padfDstNoDataReal =
> > (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));
> > warpOptions->padfDstNoDataImag =
> > (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));
> >
> > for (int j = 0; j < warpOptions->nBandCount; j++)
> > {
> > warpOptions->panSrcBands[j] = j+1;
> > warpOptions->panDstBands[j] = j+1;
> >
> > warpOptions->padfSrcNoDataReal[j] = _srcNodata;
> > warpOptions->padfSrcNoDataImag[j] = 0.0;
> >
> > warpOptions->padfDstNoDataReal[j] = _nodata;
> > warpOptions->padfDstNoDataImag[j] = 0.0;
> > }
> >
> > if (srcNumBands == 4)
> > warpOptions->nSrcAlphaBand = 4;
> >
> > if (dstNumBands == 4)
> > warpOptions->nDstAlphaBand = 4;
> >
> > char **papszWarpOptions = NULL;
> > papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST",
> "NO_DATA"
> > );
> >
> > In some cases this works, but in others I get only a black image after
> the
> > warp. So my question is basically if the above code should have worked.
> My
> > assumption here is that the GDAL Warper will regard all input pixels as
> > having full opacity when no alpha channel is specified, and that output
> > alpha is initialized to 0. Are these valid assumptions, or do I need to
> > create an alpha channel in the input image as well?
>
> I don't see anything fundamentally wrong in the above. If your source
> image is
> fully opaque, you don't need to set padfSrcNoDataReal and
> padfSrcNoDataImag. And
> if you have an output alpha channel, you don't need to set
> padfDstNoDataReal and
> padfDstNoDataImag too. In that cas you would just do papszWarpOptions =
> CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0" ); I'm not sure this
> will
> change the end result however.
>
> >
> > - Thomas
> >
> > On Thu, Jun 11, 2015 at 5:25 PM, Even Rouault <
> even.rouault at spatialys.com>
> > wrote:
> >
> > > Selon Thomas Sevaldrud <thomas at silentwings.no>:
> > >
> > > > Great, thanks!
> > > >
> > > > I tried this, but got only a red image as result, so I guess only the
> > > first
> > > > channel was used.
> > > >
> > > > This is the relevant code, where _ds is the input paletted data set
> > > >
> > > > vrtDS = vrtDriver->Create("", origW, origH, 3, GDT_Byte, NULL);
> > > > double geoTransform[6];
> > > > _ds->GetGeoTransform(geoTransform);
> > > > vrtDS->SetGeoTransform(geoTransform);
> > > > GDALSetProjection(vrtDS, _srcProjectionWKT.c_str());
> > > >
> > > > GDALRasterBand* srcBand = _ds->GetRasterBand(1);
> > > >
> > > > vrtDS->AddBand(GDT_Byte, NULL);
> > > > vrtDS->AddBand(GDT_Byte, NULL);
> > > > vrtDS->AddBand(GDT_Byte, NULL);
> > >
> > > You don't need the above 3 AddBand() since you called Create() with 3
> > > already
> > >
> > > > VRTSourcedRasterBand* vrtRBand =
> > > > (VRTSourcedRasterBand*)vrtDS->GetRasterBand(1);
> > > > VRTSourcedRasterBand* vrtGBand =
> > > > (VRTSourcedRasterBand*)vrtDS->GetRasterBand(2);
> > > > VRTSourcedRasterBand* vrtBBand =
> > > > (VRTSourcedRasterBand*)vrtDS->GetRasterBand(3);
> > > > vrtRBand->AddComplexSource(srcBand, -1, -1, -1, -1, -1, -1, -1, -1,
> 0.0,
> > > > 1.0, VRT_NODATA_UNSET, 1);
> > > > vrtGBand->AddComplexSource(srcBand, -1, -1, -1, -1, -1, -1, -1, -1,
> 0.0,
> > > > 1.0, VRT_NODATA_UNSET, 2);
> > > > vrtBBand->AddComplexSource(srcBand, -1, -1, -1, -1, -1, -1, -1, -1,
> 0.0,
> > > > 1.0, VRT_NODATA_UNSET, 3);
> > > >
> > > > Is there anything I have missed in the setup here?
> > >
> > > Looks good otherwise. Did you check that you configured properly the
> > > warper to
> > > use the 3 bands of the VRT ?
> > > >
> > > > - Thomas
> > > >
> > > >
> > > >
> > > > On Tue, Jun 9, 2015 at 1:44 PM, Even Rouault <
> even.rouault at spatialys.com
> > > >
> > > > wrote:
> > > >
> > > > > Thomas,
> > > > > >
> > > > > > I'm using the GDALWarp api from C++ to reproject and cut various
> > > input
> > > > > > images. In general this works very well for my purposes, except
> that
> > > for
> > > > > > paletted images I have to use NearestNeighbour resampling,
> > > > > >
> > > > > > I would like to use Bilinear or higher order resampling, and
> wonder
> > > if
> > > > > > there is any way to expand a paletted image to true RGB during
> the
> > > warp
> > > > > > process, so I can use these resampling methods.
> > > > >
> > > > > No, you must use an intermediate step before.
> > > > >
> > > > > >
> > > > > > Alternatively, is there a way to do this conversion (from code)
> > > before
> > > > > the
> > > > > > warp without having to run through the entire image in memory?
> > > > >
> > > > > You can use a in-memory VRT dataset to do on-the-fly expansion to
> RGB.
> > > > >
> > > > > With the AddComplexSource() method of VRTSourcedRasterBand class:
> > > > >
> > > > > CPLErr AddComplexSource( GDALRasterBand *poSrcBand,
> > > > > int nSrcXOff=-1, int
> nSrcYOff=-1,
> > > > > int nSrcXSize=-1, int
> > > nSrcYSize=-1,
> > > > > int nDstXOff=-1, int
> nDstYOff=-1,
> > > > > int nDstXSize=-1, int
> > > nDstYSize=-1,
> > > > > double dfScaleOff=0.0,
> > > > > double dfScaleRatio=1.0,
> > > > > double dfNoDataValue =
> > > > > VRT_NODATA_UNSET,
> > > > > int nColorTableComponent = 0);
> > > > >
> > > > > Set nColorTableComponent to 1,2,3 for respectively R,G,B. This is
> what
> > > > > gdal_translate will do if you use "gdal_translate -of VRT -expand
> rgb
> > > > > pct.tif
> > > > > out.vrt"
> > > > >
> > > > > Even
> > > > >
> > > > > --
> > > > > Spatialys - Geospatial professional services
> > > > > http://www.spatialys.com
> > > > >
> > > >
> > >
> > >
> > > --
> > > Spatialys - Geospatial professional services
> > > http://www.spatialys.com
> > >
> >
>
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150612/4543de2a/attachment-0001.html>
More information about the gdal-dev
mailing list