<DIV>Hello again..</DIV> <DIV> </DIV> <DIV>I would like to repost my query to the mailing list regarding the thin plate</DIV> <DIV>splines... (I now have access to my code and the internet at the same time.)</DIV> <DIV> </DIV> <DIV>I was able to get reasonable results using the </DIV> <DIV>GDALGCPsToGeoTransform function to generate a geo transform, but I would</DIV> <DIV>like to have an exact fit at the GCPs. That drove me to delve further into the</DIV> <DIV>GDAL functionality where I ran across the warping functions.</DIV> <DIV> </DIV> <DIV>I tried to first use the gdal warper to be able to reproduce the results I got</DIV> <DIV>using the GDALGCPstoGeoTransform function with the warper. I was</DIV> <DIV>not able to do that, but here is the code that I wrote. My intention was to</DIV> <DIV>look at the TPS transform's ability to preserve the GCP values, but I could not</DIV> <DIV>get that to work, or even to reproduce the
results from GDALGCPsToGeoTransform</DIV> <DIV>using regular GCP warping.</DIV> <DIV> </DIV> <DIV>I apologize for putting this much code on the mailing but I guess it's the best way....</DIV> <DIV> </DIV> <DIV>NOTE: m_pRasterDataset is the dataset that will be used for rendering and is initially </DIV> <DIV> loaded as a JPEG dataset with a set of GCPs. (The gcpCount and gcp list</DIV> <DIV> are passed in the parameter list for clarity.)</DIV> <DIV> </DIV> <DIV> </DIV> <DIV>// This method transforms gcps and saves in the raster dataset.<BR>// 0 = use GDALGCPsToGeoTrans, 1 = use GCP warp, 2 = use TPS GCP warp<BR>bool CGDALGCPImageFile::transformGCPs(int gcpCount, const GDAL_GCP *pGCPList, int
iMethod)<BR>{<BR> // Ensure we have a good dataset and it hasn't already been loaded.<BR> bool bOK = m_pRasterDataset != NULL && !hasBeenLoaded();</DIV> <DIV> if(bOK)<BR> {<BR> // Get the regular transform.<BR> double adfGeoTransform[6];<BR> bOK = GDALGCPsToGeoTransform(gcpCount, pGCPList, adfGeoTransform, TRUE);<BR> if(iMethod == 0)<BR> {<BR> // Set it and return.<BR> m_pRasterDataset->SetGeoTransform(adfGeoTransform);<BR> }<BR> else<BR>
{<BR> // Set up the warping.<BR> char *pszSrcWKT = NULL, *pszDstWKT = NULL;<BR> OGRSpatialReference oSRS;<BR> oSRS.SetWellKnownGeogCS("WGS84");<BR> oSRS.exportToWkt(&pszSrcWKT);<BR> oSRS.exportToWkt(&pszDstWKT);</DIV> <DIV> </DIV> <DIV> double defGeoTrans[6] = {0., 1., 0., 0., 0., 1.};<BR> // Ensure we have set the geo transform as the default. </DIV> <DIV> CPLErr err = m_pRasterDataset->SetGeoTransform(defGeoTrans);</DIV>
<DIV> int i = 0;<BR> GDALWarpOptions *psWO = GDALCreateWarpOptions();</DIV> <DIV> psWO->eResampleAlg = GRA_NearestNeighbour;<BR> psWO->pfnProgress = GDALTProgressGCP;</DIV> <DIV> // Create the dest dataset.<BR> double adfDstGeoTransform[6] = {0., 1., 0., 0., 0., 1.};</DIV> <DIV> // Copy off the raster metadata into a new dataset.<BR> GDALDatasetH hSrcDS = CopyRasterWindow(m_pRasterDataset, 0, 0, <BR> m_pRasterDataset->GetRasterXSize(), m_pRasterDataset->GetRasterYSize());</DIV> <DIV> int nBands = GDALGetRasterCount(hSrcDS);<BR> vector<GDALColorInterp>
srcColors;<BR> for(int i = 1; i <= nBands; i++)<BR> {<BR> GDALRasterBand *poSrcBand = m_pRasterDataset->GetRasterBand(i);<BR> srcColors.push_back(poSrcBand->GetColorInterpretation());<BR> }</DIV> <DIV> if(bOK)<BR> {<BR> psWO->hSrcDS = m_pRasterDataset;<BR> psWO->hDstDS = NULL;</DIV> <DIV> psWO->nBandCount = GDALGetRasterCount(m_pRasterDataset);<BR> psWO->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);<BR> psWO->panDstBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);</DIV> <DIV> for(int i = 0; i
< psWO->nBandCount; i++)<BR> {<BR> psWO->panSrcBands[i] = i+1;<BR> psWO->panDstBands[i] = i+1;<BR> }</DIV> <DIV> // Use 0 for GCP regular warp and -1 for TPS warp.<BR> int iGCPWarp = (iMethod == 1) ? 0 : -1;<BR> psWO->pfnTransformer = GDALGenImgProjTransform;</DIV> <DIV> </DIV> <DIV> // Use NULLs as parameters to get from pixel/line to geocoordinates<BR> psWO->pTransformerArg = GDALCreateGenImgProjTransformer(psWO->hSrcDS, </DIV> <DIV> NULL, NULL, NULL, TRUE, -1.0,
iMethod);</DIV> <DIV> </DIV> <DIV> // Figure out the desired output bounds and the geo transform.<BR> int nPixels = 0, nLines = 0;<BR> CPLErr eErr = GDALSuggestedWarpOutput(m_pRasterDataset, psWO- </DIV> <DIV> >pfnTransformer, psWO->pTransformerArg, adfDstGeoTransform, &nPixels, </DIV> <DIV> &nLines);</DIV> <DIV> // Update the transformer to include the output geotransform. <BR> GDALSetGenImgProjTransformerDstGeoTransform(psWO->pTransformerArg, </DIV>
<DIV> adfDstGeoTransform);</DIV> <DIV> if(bOK)<BR> {<BR> int nBands = GDALGetRasterCount(m_pRasterDataset);<BR> psWO->dfWarpMemoryLimit = nPixels*nLines*4*8*nBands; <BR> psWO->nBandCount = nBands;<BR> psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));<BR> psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));</DIV> <DIV> for(int i = 0; i < psWO->nBandCount; i++)<BR> {<BR> psWO->panSrcBands[i] =
i+1;<BR> psWO->panDstBands[i] = i+1;<BR> }</DIV> <DIV> psWO->nDstAlphaBand = 0;</DIV> <DIV> // Create a warped VRT dataset.<BR> GDALDatasetH hDstDS = GDALCreateWarpedVRT(m_pRasterDataset, nPixels, </DIV> <DIV> nLines, adfDstGeoTransform, psWO);<BR> bOK = (hDstDS != NULL);</DIV> <DIV> if(bOK)<BR> {<BR> // Write out the projection definition. <BR> GDALSetProjection(hDstDS,
pszDstWKT);<BR> <BR> // Copy the color table, if required.<BR> GDALColorTableH hCT = GDALGetRasterColorTable(GDALGetRasterBand</DIV> <DIV> (m_pRasterDataset, 1));<BR> </DIV> <DIV> if(hCT != NULL)<BR> {<BR> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,1), hCT);<BR> }</DIV> <DIV> // Initialize and execute the warp operation.
<BR> GDALWarpOperation *oOperation = new GDALWarpOperation;<BR> psWO->hSrcDS = m_pRasterDataset;<BR> psWO->hDstDS = hDstDS; <BR> for(int i = 0; i < nBands; i++)<BR> {<BR> ((GDALRasterBand *) GDALGetRasterBand(hDstDS,i+1))->SetAccess</DIV> <DIV> (GA_Update); <BR> }<BR> <BR> eErr =
oOperation->Initialize(psWO);<BR> if (eErr == CE_None)<BR> {<BR> eErr = oOperation->ChunkAndWarpMulti(0, 0, GDALGetRasterXSize(hDstDS), </DIV> <DIV> GDALGetRasterYSize(hDstDS));<BR> } <BR> if(eErr == CE_None)<BR> {<BR> // Set the colors in the dest.<BR> GDALRasterBand *pBand =
NULL;<BR> GDALDataset *pDstDS = (GDALDataset *)hDstDS;</DIV> <DIV> for(int i = 1; i <= nBands; i++)<BR> {<BR> pBand = pDstDS->GetRasterBand(i);<BR> GDALColorInterp gColor = srcColors[i - 1];<BR> pBand->SetColorInterpretation(gColor);<BR> }</DIV> <DIV> // We can
overwrite the destination as the source dataset.<BR> delete m_pRasterDataset;<BR> m_pRasterDataset = (GDALDataset *)hDstDS;<BR> m_pRasterDataset->SetGeoTransform(adfDstGeoTransform);<BR> }<BR> }<BR> } <BR> }<BR> // Clean up.<BR>
GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);<BR> GDALDestroyWarpOptions(psWO);<BR> }<BR> }<BR> return bOK;<BR>}<BR></DIV> <DIV> </DIV> <DIV>I appreciate any guidance that can be provided!!</DIV> <DIV> </DIV> <DIV>Thanks,</DIV> <DIV> </DIV> <DIV>Chris Pozsgai</DIV><BR><BR><p> 
<hr size=1>We won't tell. Get more on <a href="http://us.rd.yahoo.com/evt=49980/*http://tv.yahoo.com/collections/265
">shows you hate to love</a><br>(and love to hate): <a href="http://us.rd.yahoo.com/evt=49980/*http://tv.yahoo.com/collections/265
">Yahoo! TV's Guilty Pleasures list.</a>