A common 3 × 3 Template Function Based on gdal

Many template operators are often used for remote sensing image processing, such as smoothing and sharpening, Laplace operators, and souber operators. In fact, these algorithms are the same. They use a template window to move the image and then write the computing result into the image.

When viewing the gdal source code, there is a gdaldem tool, which has a template function similar to 3x3. I have transformed it to support any 3x3 Template operation.

/*** @ Brief 3*3 Template calculation processing function * @ Param hsrcband input image band * @ Param hdstband output image band * @ Param pfnalg algorithm callback function * @ Param pdata algorithm callback function parameter * @ Param pfnprogress progress bar callback function * @ Param pprogressdata progress bar callback function parameter */cplerr gdalgeneric3x3processing (gdalrasterbandh hsrcband, gdalrasterbandh hdstband, gdalgeneric3x3processingalg pfnalg, void * pdata, gdalprogressfunc pfnprogress, void * pprogressdata) {cplerr eerr; float * pafthreelinewin ;/* Three-line data storage space for input images */float * pafoutputbuf;/* one-line data storage space for output images */INT I, j; int bsrchasnodata, bdsthasnodata; float fsrcnodatavalue = 0.0, fdstnodatavalue = 0.0; int nxsize = Merge (hsrcband); int nysize = Merge (hsrcband); If (pfnprogress = NULL) pfnprogress = gdaldummyprogress; // initialize the progress bar counter if (! Pfnprogress (0.0, null, pprogressdata) {cplerror (ce_failure, cple_userinterrupt, "user terminated"); Return ce_failure;} // allocate memory space pafoutputbuf = (float *) cplmalloc (sizeof (float) * nxsize); pafthreelinewin = (float *) cplmalloc (3 * sizeof (float) * (nxsize + 1); fsrcnodatavalue = (float) gdalgetrasternodatavalue (hsrcband, & bsrchasnodata); fdstnodatavalue = (float) gdalgetrasternodatavalue (hdstband, & Bdsthasnodata); If (! Bdsthasnodata) fdstnodatavalue = 0.0; // move a 3x3 window pafwindow to traverse each pixel. // (the center pixel number is #4) /// 0 1 2 // 3 4 5 // 6 7 8 // pre-load the first two rows of data for (I = 0; I <2 & I <nysize; I ++) {gdalrasterio (hsrcband, gf_read, 0, I, nxsize, 1, pafthreelinewin + I * nxsize, nxsize, 1, gdt_float32, 0, 0 );} // does not include the boundary for (j = 0; j <nxsize; j ++) pafoutputbuf [J] = fdstnodatavalue; gdalrasterio (hdstband, gf_write, 0, 0, nxsize, 1, Pafoutputbuf, nxsize, 1, gdt_float32, 0, 0); If (nysize> 1) {gdalrasterio (hdstband, gf_write, 0, nysize-1, nxsize, 1, pafoutputbuf, nxsize, 1, gdt_float32, 0, 0);} int nline1off = 0 * nxsize; int nline2off = 1 * nxsize; int nline3off = 2 * nxsize; for (I = 1; I <nYSize-1; I ++) {// read the third row of Data eerr = gdalrasterio (hsrcband, gf_read, 0, I + 1, nxsize, 1, pafthreelinewin + nline3off, nxsize, 1, gdt_float 32, 0, 0); If (eerr! = Ce_none) goto end; // does not include boundary data pafoutputbuf [0] = fdstnodatavalue; If (nxsize> 1) pafoutputbuf [nxsize-1] = fdstnodatavalue; // use OpenMP to accelerate this sentence # pragma OMP parallel for (j = 1; j <nxsize-1; j ++) {float afwin [9]; afwin [0] = pafthreelinewin [nline1off + J-1]; afwin [1] = pafthreelinewin [nline1off + J]; afwin [2] = pafthreelinewin [nline1off + J + 1]; afwin [3] = pafthreelinewin [nline2off + J-1]; afwin [4] = pafthreelinewin [nline2off + J]; afwin [5] = pafthreelinewin [nline2off + J + 1]; afwin [6] = pafthreelinewin [nline3off + J-1]; afwin [7] = pafthreelinewin [nline3off + J]; afwin [8] = pafthreelinewin [nline3off + J + 1]; if (bsrchasnodata & (afwin [0] = fsrcnodatavalue | afwin [1] = fsrcnodatavalue | afwin [2] = fsrcnodatavalue | afwin [3] = Hangzhou | afwin [4] = fsrcnodatavalue | afwin [5] = Fsrcnodatavalue | afwin [6] = fsrcnodatavalue | afwin [7] = fsrcnodatavalue | afwin [8] = fsrcnodatavalue )) {// if one of the nine data items is nodata, set the current vertex to nodata pafoutputbuf [J] = fdstnodatavalue ;} else {// a qualified 3*3 window pafoutputbuf [J] = pfnalg (afwin, fdstnodatavalue, pdata) ;}// write a row of Data eerr = gdalrasterio (hdstband, gf_write, 0, I, nxsize, 1, pafoutputbuf, nxsize, 1, gdt_float32, 0, 0); If (eerr! = Ce_none) goto end; If (! Pfnprogress (1.0 * (I + 1)/nysize, null, pprogressdata) {cplerror (ce_failure, cple_userinterrupt, "user terminated"); eerr = ce_failure; goto end ;} int ntemp = nline1off; nline1off = nline2off; nline2off = nline3off; nline3off = ntemp;} pfnprogress (1.0, null, pprogressdata); eerr = ce_none; end: cplfree ); cplfree (pafthreelinewin); Return eerr ;}

To make this function generic, we have defined a function pointer to set the algorithm and a void parameter to set some parameters of the algorithm. The following code is used to prototype the callback function:

/*** @ Brief 3*3 Template algorithm callback function * @ Param pafwindow 3*3 window value array * @ Param fdstnodatavalue result nodata value * @ Param pdata algorithm callback function parameter *@ return 3*3 window calculation result value */typedef float (* gdalgeneric3x3processingalg) (float * pafwindow, float fdstnodatavalue, void * pdata );

The usage and parameter descriptions of functions are included in the Code. Here are two examples. First, do a Laplace operator. First, describe the Laplace operator. Laplacian Operator Formula: Delta ^ 2f (x, y) = f (x + 1, Y) + f (x-1, Y) + f (x, y + 1) + f (x, Y-1)-4f (x, y), as shown in figure 1 of the template:

Figure 1 Laplacian Operator

Now, we need to implement an algorithm function for template calculation and an algorithm parameter struct. The type of this algorithm function is like the callback function above, and the algorithm parameter struct is used to store template operators. The following is the code of the parameter struct of the algorithm, which is very simple. It is an array of double [9.

/*** @ Brief space filter algorithm Data Structure */typedef struct {/*! Template operator */Double dtemplate [9];} spacefilteralgdata;

Next, the implementation function of the template algorithm. We know that the template algorithm is to multiply each value of the template operator by the Image Element value and then sum it. The result is the image element value of the result image.

float SpaceFilterAlg (float* afWin, float fDstNoDataValue, void* pData){    SpaceFilterAlgData* psData = (SpaceFilterAlgData*)pData;    double fResult = 0.0;    for(int i=0; i<9; i++)        fResult += (psData->dTemplate[i]*afWin[i]);            return (float)fResult;}

The above Algorithm functions are also completed, and the following is the problem of starting to call. Before calling the API, We need to write a function that is used to create the parameter struct of the algorithm:

void* CreateSpaceFilterAlgData(double* adfTemplate){    SpaceFilterAlgData* pData =        (SpaceFilterAlgData*)CPLMalloc(sizeof(SpaceFilterAlgData));    if(adfTemplate != NULL)        memcpy(pData->dTemplate, adfTemplate, sizeof(double)*9);    else        memset(pData->dTemplate, 0, sizeof(double)*9);    return pData;}

With the above preparations, we will start to write a space filter function:

/*** @ Brief image space filtering, only applicable to 3*3 templates * @ Param pszsrcfile input file path * @ Param pszdstfile output file path * @ Param ptemplate 3*3 Template operator * @ Param edatatype output data type, the default value is gdt_unknown, which is consistent with the original image type * @ Param pszformat output file format. For details, refer to gdal support data type * @ Param pprocess progress bar pointer * @ return returned value, various error messages during calculation */INT imagespacefilter (const char * pszsrcfile, const char * pszdstfile, double * ptemplate, gdaldatatype edatatype, const char * pszformat, cprocessbase * PP Rocess) {If (pprocess! = NULL) {pprocess-> resetprocess (); pprocess-> setmessage ("executing space filter...") ;}if (pszsrcfile = NULL) {If (pprocess! = NULL) pprocess-> setmessage ("the source file path is empty! "); Return re_filenotexist;} If (pszdstfile = NULL) {If (pprocess! = NULL) pprocess-> setmessage ("the output file path is empty! "); Return region;} gdalallregister (); gdaldataseth hsrcdataset = NULL; gdaldataseth hdstdataset = NULL; export hsrcband = NULL; gdalrasterbandh hdstband = NULL; gdaldriverh hdriver = NULL; char ** papszcreateoptions = NULL; gdalprogressfunc pfnprogress = progress; // open the image hsrcdataset = gdalopen (pszsrcfile, ga_readonly); If (hsrcdataset = NULL) {If (pprocess! = NULL) pprocess-> setmessage ("the input file cannot be opened! "); Return re_filenotexist;} int nxsize = round (hsrcdataset); int nysize = round (hsrcdataset); int ndstbands = gdalgetrastercount (hsrcdataset); double adfgeotransform [6]; gdalgetgeotransform (hsrcdataset, adfgeotransform); hdriver = gdalgetdriverbyname (pszformat); If (hdriver = NULL) {If (pprocess! = NULL) pprocess-> setmessage ("You cannot create a file of the specified type. check whether this file type gdal supports creation! "); Gdalclose (hsrcdataset); Return re_filenotsupport;} hsrcband = gdalgetrasterband (hsrcdataset, 1); gdaldatatype edstdatatype = aggregate (hsrcband); If (edatatype! = Gdt_unknown) edstdatatype = edatatype; hdstdataset = gdalcreate (hdriver, temperature, nxsize, nysize, ndstbands, edstdatatype, papszcreateoptions); If (hdstdataset = NULL) {If (pprocess! = NULL) pprocess-> setmessage ("An error occurred while creating the output file. Check whether the file path is correct! "); Gdalclose (hsrcdataset); return failed; // image creation failed} failed (hdstdataset, adfgeotransform); gdalsetprojection (hdstdataset, gdalgetprojectionref (hsrcdataset); double success = 0; int bdsthasnodata = false; void * pdata = NULL; optional pfnalg = NULL; pdata = createspacefilteralgdata (ptemplate); pfnalg = spacefilteralg; For (INT I = 1; I <= ndstbands; I + +) {Hsrcband = gdalgetrasterband (hsrcdataset, I); hdstband = random (hdstdataset, I); If (bdsthasnodata) random (hdstband, random); If (random (hsrcband, hdstband, pfnalg, pdata, pfnprogress, pprocess )! = Ce_none) {If (pprocess! = NULL) {If (! Pprocess-> m_biscontinue) {gdalclose (hsrcdataset); gdalclose (hdstdataset); cplfree (pdata); csldestroy (papszcreateoptions); Return re_usercancel;} return re_failed ;}} gdalclose (hsrcdataset); gdalclose (hdstdataset); cplfree (pdata); csldestroy (papszcreateoptions); If (pprocess! = NULL) pprocess-> setmessage ("computing complete! "); Return re_success ;}

The descriptions of these functions are annotated in the code, basically all rely on gdal, there are no other libraries, there are those progress bars, what are the returned values, refer to my previous blog. At this point, we have completed all the functions of Space filtering, and we will call them below. Okay, we use the Laplace operator. The call code is as follows:

int LaplacianFilter(){    CConsoleProcess *pProgress = new CConsoleProcess();      double dTemplate[9] =     {        0, 1, 0,        1, -4, 1,        0, 1, 0    };        int iRev = ImageSpaceFilter("C:\\Test.img", "C:\\Lapiacian.tif", dTemplate, GDT_Unknown, "GTiff", pProgress);    if(iRev == RE_SUCCESS)        printf("Success!\n");    else        printf("Failed!\n");            delete pProgress;}

The original graph and result graph of program execution, as shown in figure 2 and Figure 3:

Figure 2 Original Image

Figure 3 Laplacian filtered image

Next, I will write an example using this 3 × 3 template to extract the slope and slope of DEM.