Geometric Correction Code for gdal

A gdal rectification project rasterrectifier.rar was found on the Internet, with an error in it. The modified code is as follows:

A class is actually a function;

 

The Code is as follows:

//////////////////////////////////////// //////////////////////////////////////

# Include "./include/gdal. H"

Class crectifier
{
Public:
Crectifier ();
Crectifier (cstring filename );
Virtual ~ Crectifier ();

PRIVATE:
Int m_count; // number of control points
Gdal_gcp * m_gcps; // Control Point pair
Cstring m_scrfilename; // original file name
Cstring m_destfilename; // target file name

Public:
// Function: Add a Control Point
// Parameter: GCPs-Control Point pair; ncount-control point count
Void addgcps (double * GCPs, int ncount );
 
// Function: perform corrective actions
Bool transform ();

// Function: Set the source file to be retrieved.
Void setscrfilename (cstring filename );
// Function: Obtain the target file
Cstring getdestfilename ();
};

//////////////////////////////////////// ///////////////

# Include "stdafx. H"
# Include "rasterrectifier. H"
# Include "crectifier. H"

# Include "./include/cpl_conv.h"
# Include "./include/cpl_string.h"
# Include "./include/ogrsf_frmts.h"
# Include "./include/gdalwarper. H"

# Pragma comment (Lib, "lib // gdal_ I .lib ")

# Ifdef _ debug
# UNDEF this_file
Static char this_file [] =__ file __;
# Define new debug_new
# Endif

//////////////////////////////////////// //////////////////////////////
// Construction/destruction
//////////////////////////////////////// //////////////////////////////

Crectifier: crectifier ()
{
 
}

Crectifier: crectifier (cstring filename)
{
Setscrfilename (filename );
}

Crectifier ::~ Crectifier ()
{

}

Void crectifier: addgcps (double * GCPs, int ncount)
{
M_count = ncount;

M_gcps = new gdal_gcp [m_count];

For (INT I = 0; I <m_count; I ++)
{
Char pszid [5];
Sprintf (pszid, "GCP-% d", I );

M_gcps [I]. pszid = pszid;
M_gcps [I]. pszinfo = "";

M_gcps [I]. dfgcppixel = GCPs [4 * I];
M_gcps [I]. dfgcpline = GCPs [4 * I + 1];
M_gcps [I]. dfgcpx = GCPs [4 * I + 2];
M_gcps [I]. dfgcpy = GCPs [4 * I + 3];
M_gcps [I]. dfgcpz = 0;
}
}
 
Void crectifier: setscrfilename (cstring filename)
{
M_scrfilename = filename;
Int Index = m_scrfilename.reversefind ('.') + 1;
M_destfilename = m_scrfilename.left (INDEX) + "_ DEST ";
M_destfilename + = m_scrfilename.right (m_scrfilename.getlength ()-index );
}

Cstring crectifier: getdestfilename ()
{
Return m_destfilename;
}

Bool crectifier: Transform ()
{
Gdaldataseth hsrcds, hdstds;
Gdaldriverh hdriver;
Gdaldatatype EDT;

Try
{
Gdalallregister ();

// Open the source file.

Hsrcds = gdalopen (m_scrfilename, ga_readonly );
Cplassert (hsrcds! = NULL );
Gdalsetgcps (hsrcds, m_count, m_gcps, null );

// Create output with same datatype as first input band.

EDT = gdalgetrasterdatatype (gdalgetrasterband (hsrcds, 1 ));

// Get output driver (geotiff format)

Hdriver = gdalgetdriverbyname ("gtiff ");
Cplassert (hdriver! = NULL );

// Get source coordinate system.

Const char * pszsrcwkt = NULL;
//
// Pszsrcwkt = gdalgetprojectionref (hsrcds );
// Cplassert (pszsrcwkt! = NULL & strlen (pszsrcwkt)> 0 );
//
/// Setup output coordinate system that is input by user.
//
Char * pszdstwkt = NULL;
//
// Otargetsrs. exporttowkt (& pszdstwkt );
// Cplassert (pszdstwkt! = NULL & strlen (pszdstwkt)> 0 );

// Create a transformer that maps from source pixel/line coordinates
// To destination georeferenced coordinates (not destination
// Pixel line). We do that by omitting the destination Dataset
// Handle (setting it to null ).

Void * htransformarg;

Htransformarg =
Gdalcreategenimgprojtransformer (hsrcds, pszsrcwkt, null, pszdstwkt,
False, 0, 1 );
Cplassert (htransformarg! = NULL );

// Get approximate output georeferenced bounds and resolution for file.

Double adfdstgeotransform [6];
Int npixels = 0, nlines = 0;
Cplerr eerr;

Eerr = gdalsuggestedwarpoutput (hsrcds,
Gdalgenimgprojtransform, htransformarg,
Adfdstgeotransform, & npixels, & nlines );
Cplassert (eerr = ce_none );

Gdaldestroygenimgprojtransformer (htransformarg );

// Obtain the number of bands

Int nbandcount = gdalgetrastercount (hsrcds );

// Create the output file.

Hdstds = gdalcreate (hdriver, m_destfilename, npixels, nlines,
Nbandcount, EDT, null );

Cplassert (hdstds! = NULL );

// Write out the projection definition.

Gdalsetprojection (hdstds, pszdstwkt );

Gdalsetgeotransform (hsrcds, adfdstgeotransform );
Gdalsetgeotransform (hdstds, adfdstgeotransform );

// Obtain the color table in band-by-band Mode
Gdalcolortableh HCT;

For (INT I = 1; I <= nbandcount; I ++) // the band starts from 1.
{
HCT = gdalgetrastercolortable (gdalgetrasterband (hsrcds, I ));
If (Hct! = NULL)
Gdalsetrastercolortable (gdalgetrasterband (hdstds, I), HCT );
}

// Setup warp options.

Gdalwarpoptions * pswarpoptions = gdalcreatewarpoptions ();

Pswarpoptions-> hsrcds = hsrcds;
Pswarpoptions-> hdstds = hdstds;

Char ** papszwarpoptions = NULL;
Papszwarpoptions = cslsetnamevalue (papszwarpoptions,
"Tfw", "yes ");

Pswarpoptions-> papszwarpoptions = papszwarpoptions;

Pswarpoptions-> nbandcount = nbandcount;

// Apply for memory

Pswarpoptions-> pansrcbands =
(Int *) cplmalloc (sizeof (INT) * pswarpoptions-> nbandcount );
Pswarpoptions-> pandstbands =
(Int *) cplmalloc (sizeof (INT) * pswarpoptions-> nbandcount );

For (Int J = 0; j <nbandcount; j ++)
{
Pswarpoptions-> pansrcbands [J] = J + 1;
Pswarpoptions-> pandstbands [J] = J + 1;
}

Pswarpoptions-> pfnprogress = gdaltermprogress;

// Establish Reprojection transformer.

Pswarpoptions-> ptransformerarg =
Gdalcreategenimgprojtransformer (hsrcds,
Gdalgetprojectionref (hsrcds ),
Hdstds,
Gdalgetprojectionref (hdstds ),
False, 0.0, 1 );
Pswarpoptions-> pfntransformer = gdalgenimgprojtransform;

// Initialize and execute the warp operation.

Gdalwarpoperation ooperation;

Ooperation. initialize (pswarpoptions );
Ooperation. chunkandwarpimage (0, 0,
Gdalgetrasterxsize (hdstds ),
Gdalgetrasterysize (hdstds ));

Csldestroy (papszwarpoptions );

Gdaldestroygenimgprojtransformer (pswarpoptions-> ptransformerarg );
Gdaldestroywarpoptions (pswarpoptions );

Gdalclose (hdstds );
Gdalclose (hsrcds );

Return true;
}
Catch (...)
{
Return false;
}
}

 

//////////////////////////////////////// //////////////////////////////////////// ////////

 

The above code is accurate after testing.

 

Generally, the following two errors may occur in online code:

1. After correction, the image is flipped.

The reason is: The gdalsetgeotransform () function sets the number of times the radiology transformation coefficient is used;

2. After correction, the image pixels are all 0.

Cause: No projection information is initialized when the result affects initialization.

 

Of course, we recommend that you encapsulate the above Code as a gdal command line application.

 

The correction function is implemented in gdalwarp.exe in the gdalapplication. However, this program is mainly used for re-projection, and the calling parameters are somewhat confusing.