[Gdal-dev] RE: GCPs and Thin Plate Splines - C++

[chris pozsgai]

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

   
  I appreciate any guidance that can be provided!!
   
  Thanks,
   
  Chris Pozsgai



 
---------------------------------
We won't tell. Get more on shows you hate to love
(and love to hate): Yahoo! TV's Guilty Pleasures list.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20070309/9281aa54/attachment.html