Revision: 70132 http://sourceforge.net/p/brlcad/code/70132 Author: starseeker Date: 2017-08-25 21:11:09 +0000 (Fri, 25 Aug 2017) Log Message: ----------- Getting closer. For some reason GTiff works better than MEM as an intermediate for the projection, but with the gdal_translate steps added the hawaii example looks closer to the expected results.
Modified Paths: -------------- brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp Modified: brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp =================================================================== --- brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp 2017-08-25 15:56:53 UTC (rev 70131) +++ brlcad/trunk/src/libgcv/plugins/gdal/gdal.cpp 2017-08-25 21:11:09 UTC (rev 70132) @@ -49,6 +49,7 @@ #include <cpl_string.h> #include <cpl_multiproc.h> #include <ogr_spatialref.h> +#include "vrtdataset.h" #include "raytrace.h" #include "gcv/api.h" @@ -250,7 +251,7 @@ GDALAllRegister(); - state->hDataset = GDALOpenEx(source_path, GDAL_OF_READONLY | GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, NULL, NULL, NULL); + state->hDataset = GDALOpenEx(source_path, GDAL_OF_READONLY | GDAL_OF_VERBOSE_ERROR, NULL, NULL, NULL); if (!state->hDataset) { bu_log("GDAL Reader: error opening input file %s\n", source_path); BU_PUT(state, struct conversion_state); @@ -278,7 +279,7 @@ GDALDestroyGenImgProjTransformer(hTsfA); if (eErr != CE_None) return 0; - GDALDatasetH pjdata = GDALCreate(GDALGetDriverByName("MEM"), "", nPixels, nLines, GDALGetRasterCount(state->hDataset), + GDALDatasetH pjdata = GDALCreate(GDALGetDriverByName("GTiff"), "test_out.tiff", nPixels, nLines, GDALGetRasterCount(state->hDataset), GDALGetRasterDataType(GDALGetRasterBand(state->hDataset,1)), NULL); if (pjdata == NULL) return 0; GDALSetProjection(pjdata, dst_Wkt); @@ -310,16 +311,57 @@ //GDALDatasetH indata = state->hDataset; unsigned int xsize = GDALGetRasterXSize(indata); unsigned int ysize = GDALGetRasterYSize(indata); +#endif + + + double adfGeoTransform[6]; + VRTDataset *poVDS = (VRTDataset *)VRTCreate(GDALGetRasterXSize(indata), GDALGetRasterYSize(indata)); + poVDS->SetProjection(GDALGetProjectionRef(indata)); + GDALGetGeoTransform(indata, adfGeoTransform); + poVDS->SetGeoTransform(adfGeoTransform); + double adfWin[4]; + adfWin[0] = 0; + adfWin[1] = 0; + adfWin[2] = xsize; + adfWin[3] = ysize; + /* psOptions->adfSrcWin */ + GDALRasterBand *poSrcBand = ((GDALDataset *)indata)->GetRasterBand(1); + poVDS->AddBand(GDALGetRasterDataType(poSrcBand), NULL); + VRTSourcedRasterBand *poVRTBand = (VRTSourcedRasterBand *)poVDS->GetRasterBand(1); + double dfScale=1.0, dfOffset=0.0; + double dfScaleSrcMin = 0.0, dfScaleSrcMax = 0.0; + double dfScaleDstMin = 0.0, dfScaleDstMax = 255.999; + double adfCMinMax[2]; + GDALComputeRasterMinMax(poSrcBand, TRUE, adfCMinMax); + dfScaleSrcMin = adfCMinMax[0]; + dfScaleSrcMax = adfCMinMax[1]; + dfScale = (dfScaleDstMax - dfScaleDstMin) / (dfScaleSrcMax - dfScaleSrcMin); + dfOffset = -1 * dfScaleSrcMin * dfScale + dfScaleDstMin; + VRTSimpleSource* poSimpleSource; + VRTComplexSource* poSource = new VRTComplexSource(); + poSource->SetLinearScaling(dfOffset, dfScale); + poSimpleSource = poSource; + poSimpleSource->SetResampling(NULL); + poVRTBand->ConfigureSource( poSimpleSource, poSrcBand, FALSE, adfWin[0], adfWin[1], adfWin[2], adfWin[3], adfWin[0], adfWin[1], adfWin[2], adfWin[3]); + poVRTBand->AddSource(poSimpleSource); + GDALDatasetH hOutDS = GDALCreateCopy(GDALGetDriverByName("MEM"), "", (GDALDatasetH) poVDS, 0, NULL, GDALTermProgress, NULL); + GDALClose((GDALDatasetH)poVDS); + + indata = hOutDS; + +#if 1 unsigned short *uint16_array = (unsigned short *)bu_calloc(xsize*ysize, sizeof(unsigned short), "unsigned short array"); GDALRasterBandH band = GDALGetRasterBand(indata, 1); - int bmin, bmax; - double adfMinMax[2]; - adfMinMax[0] = GDALGetRasterMinimum(band, &bmin); - adfMinMax[1] = GDALGetRasterMaximum(band, &bmax); - if (!(bmin && bmax)) GDALComputeRasterMinMax(band, TRUE, adfMinMax); - bu_log("Min/Max: %f, %f\n", adfMinMax[0], adfMinMax[1]); + /* + int bmin, bmax; + double adfMinMax[2]; + adfMinMax[0] = GDALGetRasterMinimum(band, &bmin); + adfMinMax[1] = GDALGetRasterMaximum(band, &bmax); + if (!(bmin && bmax)) GDALComputeRasterMinMax(band, TRUE, adfMinMax); + bu_log("Min/Max: %f, %f\n", adfMinMax[0], adfMinMax[1]); + */ /* If we're going to DSP we need the unsigned short read. */ uint16_t *scanline = (uint16_t *)CPLMalloc(sizeof(uint16_t)*GDALGetRasterBandXSize(band)); @@ -330,7 +372,7 @@ /* This is the critical assignment point - if we get this * indexing wrong, data will not look right in dsp */ uint16_array[ysize*i+j] = scanline[j]; - bu_log("%d, %d: %d\n", i, j, scanline[j]); + //bu_log("%d, %d: %d\n", i, j, scanline[j]); } } } @@ -367,7 +409,7 @@ dsp->dsp_datasrc = RT_DSP_SRC_OBJ; dsp->dsp_xcnt = ysize; dsp->dsp_ycnt = xsize; - dsp->dsp_smooth = 0; + dsp->dsp_smooth = 1; dsp->dsp_cuttype = DSP_CUT_DIR_ADAPT; MAT_IDN(dsp->dsp_stom); dsp->dsp_stom[0] = dsp->dsp_stom[5] = 1000; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ------------------------------------------------------------------------------ Check out the vibrant tech community on one of the world's most engaging tech sites, Slashdot.org! http://sdm.link/slashdot _______________________________________________ BRL-CAD Source Commits mailing list brlcad-commits@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/brlcad-commits