C++ GDAL——在原始影像中根据云检图生成抠掉云的原始影像

#include 
#include 
#include 
#include "shlwapi.h"
#include 
#include  
#include 
#pragma comment(lib,"imagehlp.lib")
#pragma comment(lib, "shlwapi.lib")
#include "CreatePyramid.h"
#include "gdal_priv.h"
#pragma warning(disable:4996)

#define RAM_SIZE  100

using namespace std;

int main(int argc, char* argv[])
{

    if (argc != 2)
    {
        cout << "argument is false" << endl;
        return 1;
    }
    string strError = argv[1]; int posError = strError.find_last_of("\\"); strError.erase(posError + 1); strError += "Error.log";
    std::ofstream finerror(strError, ios::app);
    FILE*pErrFile;
    if ((pErrFile = fopen(strError.c_str(), "r")) == NULL)
    {

        printf("Open Tsk File failed:\n%s\n", strError);

        return 1;

    }
    string strTsk = argv[1];
    std::ifstream fin(strTsk);
    if (!fin)
    {
        finerror << "ERROR: " << argv[1] << endl;
        std::cout << "Can not open ini file!  " << strTsk << std::endl;
        return 1;

    }


    string Argument; vector<string>strArg;
    while (getline(fin, Argument))
    {
        strArg.push_back(Argument);
    }
    string resultDom = strArg[2];
    int poResult = resultDom.find_last_of("\\");
    resultDom.erase(poResult+1);
    const char*resultDomPath = resultDom.c_str();
    if (!PathIsDirectory(resultDomPath))
    {
        MakeSureDirectoryPathExists((PCSTR)resultDomPath);
    }
    string TskLog = argv[1]; int posLog = TskLog.find_last_of("."); TskLog.erase(posLog); TskLog += ".log";
    ofstream LogWrite(TskLog);
    FILE*pLogFile;
    if ((pLogFile = fopen(TskLog.c_str(), "r")) == NULL)
    {

        printf("Open Tsk File failed:\n%s\n", TskLog);
        finerror << "ERROR: " << strArg [0]<< endl;
        return 1;

    }
    GDALDataset* poSrcDS;
    const char*strDom = strArg[0].c_str();
    GDALAllRegister();
    poSrcDS = (GDALDataset*)GDALOpen(strDom, GA_ReadOnly);
    if (poSrcDS == NULL)
    {
        LogWrite << strDom<<"  :no exit or errror" << endl;
        cout << "Can't Open the File and Image!" << endl;
        return EXIT_FAILURE;
    }
    LogWrite << strDom << "  Open succeed" << endl;
    int iBandCount = poSrcDS->GetRasterCount();
    int iSrcWidth = poSrcDS->GetRasterXSize();
    int iSrcHeight = poSrcDS->GetRasterYSize();

    const char* pszFormat = "GTiff";
    GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if (poDriver == NULL)
    {
        finerror << "ERROR: " << strArg[0] << endl;
        GDALClose((GDALDatasetH)poSrcDS);
        return 1;
    }
    double dGeoTransform[6];
    poSrcDS->GetGeoTransform(dGeoTransform);  
    const char* pszProj = poSrcDS->GetProjectionRef();


    GDALDataset* poSrcDS1;
    poSrcDS1 = (GDALDataset*)GDALOpen(strArg[1].c_str(), GA_ReadOnly);
    if (poSrcDS1 == NULL)
    {
        finerror << "ERROR: " << strArg[0] << endl;
        LogWrite << strArg[1]<< "  :no exit or errror" << endl;
        cout << "Can't Open the File and Image!" << endl;
        return EXIT_FAILURE;
    }
    LogWrite << strArg[1] << "  Open succeed" << endl;

    int iBandCount1 = poSrcDS1->GetRasterCount();
    int iSrcWidth1 = poSrcDS1->GetRasterXSize();
    int iSrcHeight1 = poSrcDS1->GetRasterYSize();


    const char* pszFormat1 = "GTiff";
    GDALDriver *poDriver1 = GetGDALDriverManager()->GetDriverByName(pszFormat1);
    if (poDriver1 == NULL)
    {
        finerror << "ERROR: " << strArg[0] << endl;
        GDALClose((GDALDatasetH)poSrcDS1);
        return 1;
    }
    double dGeoTransform1[6];
    poSrcDS1->GetGeoTransform(dGeoTransform1); 
    const char* pszProj1 = poSrcDS1->GetProjectionRef();


    GDALDataset* poDstDS;
    poDstDS = poDriver->Create(strArg[2].c_str(), iSrcWidth, iSrcHeight, iBandCount, GDT_UInt16, NULL);
    int nStepSize = (RAM_SIZE * 1024 * 1024) / (iSrcWidth * iBandCount);
    int nStepNum = iSrcHeight / nStepSize;
    if (iSrcHeight%nStepSize)nStepNum++;
    int pBandMap1[] = { 1 };
    int pBandMap[] = { 1, 2, 3, 4 };
    double nodata;
    nodata = poSrcDS->GetRasterBand(1)->GetNoDataValue();
    poDstDS->GetRasterBand(1)->SetNoDataValue(nodata);
    poDstDS->SetGeoTransform(dGeoTransform);
    poDstDS->SetProjection(pszProj);
    cout << "Nocld DOM Images are generating:" << strArg[2].c_str() << endl;
    for (int k = 0; k < nStepNum; k++)
    {
        cout << double(k*100 / nStepNum) << "%" << endl;
        int ybeg = max(0, min(nStepSize*k, iSrcHeight - 1));
        int yend = max(0, min(nStepSize*(k + 1), iSrcHeight));
        WORD *pImg = new WORD[(yend - ybeg)*iSrcWidth*iBandCount];       
        memset(pImg, 0, (yend - ybeg)*iSrcWidth*iBandCount*sizeof(WORD));
        BYTE *pImg1 = new BYTE[(yend - ybeg)*iSrcWidth*iBandCount1];
        memset(pImg1, 0, (yend - ybeg)*iSrcWidth*iBandCount1*sizeof(BYTE));
        poSrcDS->RasterIO(GF_Read, 0, ybeg, iSrcWidth, (yend - ybeg), pImg, iSrcWidth, (yend - ybeg), GDT_UInt16, iBandCount, pBandMap, 0, 0, 0);

        poSrcDS1->RasterIO(GF_Read, 0, ybeg, iSrcWidth1, (yend - ybeg), pImg1, iSrcWidth1, (yend - ybeg), GDT_Byte, iBandCount1, pBandMap1, 0, 0, 0);
        for (int j = 0; j < (yend - ybeg)*iSrcWidth*iBandCount1; ++j)
        {
            if (pImg1[j] == 255)
            {
                pImg[j] = 0;
                pImg[j + iSrcWidth*(yend - ybeg) * 1] = 0;
                pImg[j + iSrcWidth*(yend - ybeg) * 2] = 0;
                pImg[j + iSrcWidth*(yend - ybeg) * 3] = 0;
            }
        }
        poDstDS->RasterIO(GF_Write, 0, ybeg, iSrcWidth, (yend - ybeg), pImg, iSrcWidth, (yend - ybeg), GDT_UInt16, iBandCount, pBandMap, 0, 0, 0);
        delete[]pImg; pImg = NULL;
        delete[]pImg1; pImg1 = NULL;
    }
    cout << endl;
    GDALClose((GDALDatasetH)poSrcDS);
    GDALClose((GDALDatasetH)poSrcDS1); 
    GDALClose((GDALDatasetH)poDstDS);
    cout << "\nCreating Pyramids:" << endl;

    CConsoleProcess *pProgress = new CConsoleProcess();

    bool f = CreatePyramids(strArg[2].c_str(), pProgress);

    if (f == true)
    {
        cout << "***********************CreatePyramids Successed...***********************\n" << endl;
        LogWrite << "CreatePyramids Successed" << endl;
    }

    else
    {
        finerror << "ERROR: " << strArg[0] << endl;
        LogWrite << "CreatePyramids Failed" << endl;
        cout << "***********************CreatePyramids Failed...***********************\n" << endl;
    }
    delete pProgress;
    LogWrite.close();
    fclose(pLogFile);
    finerror.close();
    fclose(pErrFile);
    return 0;
}

/**
 * 在原始影像中根据云检图生成抠掉云的原始影像
 * 程序主要依赖的方法是gdal中rasterIO函数对tif影像的分块读写，使用性较强.
 */