www.pudn.com > dstile-0.2.rar > CSTrans.cpp

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
using namespace std;
#include "exceptions.h"
#include "giskit.h"
#include "CSTrans.h"

CSTrans::CSTrans() :
    m_srcGT(),
    m_dstGT(),
    m_rta(0),
    m_ata(0)
{
}

CSTrans::~CSTrans() {
    if (m_ata) GDALDestroyApproxTransformer(m_ata);
    if (m_rta) GDALDestroyReprojectionTransformer(m_rta);
}

void CSTrans::Create(const string& srcProj, const string& dstProj, double maxError) {
    m_rta = GDALCreateReprojectionTransformer(ProjToSRS(srcProj).c_str(), ProjToSRS(dstProj).c_str());
    if (maxError >= 0.001) m_ata = GDALCreateApproxTransformer(GDALReprojectionTransform, m_rta, maxError);
}

void CSTrans::Destroy() {
    if (m_ata) {
	GDALDestroyApproxTransformer(m_ata);
	m_ata = 0;
    }
    if (m_rta) {
	GDALDestroyReprojectionTransformer(m_rta);
	m_rta = 0;
    }
}

void CSTrans::SetSrcGT(const GeoT& gt) {
    m_srcGT = gt;
    GDALInvGeoTransform(const_cast(gt.m_mat), m_srcIGT.m_mat);
}

void CSTrans::SetDstGT(const GeoT& gt) {
    m_dstGT = gt;
    GDALInvGeoTransform(const_cast(gt.m_mat), m_dstIGT.m_mat);
}

int CSTrans::STransform(void *arg, int dts, int n, double *x, double *y, double *z, int *ps) {
    double tx, ty;
    int i;
    CSTrans *t = reinterpret_cast(arg);
    double srcMat[6];
    double dstMat[6];
	
    if (dts) {
        for (i = 0; i < 6; ++i) srcMat[i] = t->m_dstGT.m_mat[i];
        for (i = 0; i < 6; ++i) dstMat[i] = t->m_srcIGT.m_mat[i];
    } else {
        for (i = 0; i < 6; ++i) srcMat[i] = t->m_srcGT.m_mat[i];
        for (i = 0; i < 6; ++i) dstMat[i] = t->m_dstIGT.m_mat[i];
    }
    
    for (i = 0; i < n; ++i) {
        tx = x[i]; ty = y[i];
        x[i] = srcMat[0] + tx * srcMat[1] + ty * srcMat[2];
        y[i] = srcMat[3] + tx * srcMat[4] + ty * srcMat[5];
    }

    if (t->m_ata) GDALApproxTransform(t->m_ata, dts, n, x, y, z, ps);
    else GDALReprojectionTransform(t->m_rta, dts, n, x, y, z, ps);

    for (i = 0; i < n; ++i) {
        tx = x[i]; ty = y[i];
        x[i] = dstMat[0] + tx * dstMat[1] + ty * dstMat[2];
        y[i] = dstMat[3] + tx * dstMat[4] + ty * dstMat[5];
    }
}