無(wú)人機(jī)影像 + DEM 計(jì)算四個(gè)角點(diǎn)坐標(biāo)（剛性變換）

1. 像空間坐標(biāo)（x,y,-f）

2. 像空間坐標(biāo)畸變糾正 deltax,deltay
   

3. 已知(x,y)，求解(X,Y, Z)或者(Lat,Lon)
   
   這里的Z是DEM上獲取的坐標(biāo)和Zs為相機(jī)坐標(biāo)的高程，如果均為已知的情況下，則可以求解(X,Y)，這里的(X,Y,Z)為地固地心坐標(biāo)，單位為米。平地的情況只需要獲取行高即可求解(X,Y)，接著使用proj庫(kù)將地固地心坐標(biāo)轉(zhuǎn)化為經(jīng)緯度坐標(biāo)即可。

4. 地理配準(zhǔn)
   這里直接采用**gdal_translate和gdal_wrap**，gdal_translate轉(zhuǎn)換過(guò)程如下，大概就是將jpg進(jìn)行地理配準(zhǔn)。請(qǐng)注意，GDAL的影像起點(diǎn)是左上角，但是我們的相機(jī)模型是左下角，所以需要變換Y軸。轉(zhuǎn)換之后的tif只有專(zhuān)業(yè)的QGIS之類(lèi)的軟件才能讀取。
   

   下面是QGIS讀取的效果。但是為了geoserver能夠識(shí)別，還需要轉(zhuǎn)換，這時(shí)候需要gdal_wrap,這個(gè)時(shí)候就很關(guān)鍵，我們需要設(shè)置其為透明。

   gdalwarp -r cubic -ovr AUTO -dstalpha D:\code\roadProj\public\demo\DSC00002\test.tif D:\code\roadProj\public\demo\DSC00002\test3_geotiff.tif
   

   

   gdal_translate.exe -of GTiff -gcp 0 5304 102.1265090139 29.6453703982 -gcp 7952 5304 102.1164515460 29.6474820822 -gcp 7952 0 102.1131839750 29.6445496193 -gcp 0 0 102.1233217949 29.6424804051 -ovr AUTO -co GCPs_Creation=YES D:\code\roadProj\public\demo\DSC00002\DSC00002.jpg D:\code\roadProj\public\demo\DSC00002\test.tif
   

   5. geoserver發(fā)布，具體操作比較簡(jiǎn)單

代碼邏輯

1. 下面是求解影像四個(gè)角點(diǎn)經(jīng)緯度的簡(jiǎn)單思路，主要還是共線(xiàn)方程，代碼中的1000還是得根據(jù)距離地面的高度，即需要DEM的高程值才能求解得到較為精確的精度。

2. 具體實(shí)現(xiàn)分為相機(jī)模型（固定的參數(shù)不部分），大疆無(wú)人機(jī)是WGS84橢球，EPSG:4978是地心地固的轉(zhuǎn)換參數(shù)。

struct CameraModel {double f = 7538.508;  // 像素為單位double u0 = 3982.417; // 像素為單位double v0 = 2671.637;double pixelSize = 4.5e-6; // 米為單位double k1 = 2.470920e-9;   // 徑向畸變系數(shù)double k2 = -2.767172e-16;double k3 = 2.479935e-23;double k4 = -6.583598e-31;double p1 = 1.388595e-8; // 偏心畸變系數(shù)double p2 = 1.781812e-7;double alpha = -4.697031e-4; // CCD非正方形比例系數(shù)double beta = -1.300023e-4;double width = 7952; // 影像的高度和寬度double height = 5304;
};``````cpp
///
/// \brief The ComputeBoundingBox class
/// 計(jì)算影像的包圍盒的經(jīng)緯度坐標(biāo) + 高程，然后貼地
///
class ComputeBoundingBox {
public:ComputeBoundingBox() {transTool.init("EPSG:4326","EPSG:4978"); // 橢球坐標(biāo)->地心坐標(biāo) XYZ}////// \brief resetR/// \param phi 俯仰角/// \param omega 橫滾角/// \param kappa 旋轉(zhuǎn)角///void resetR(double phi, double omega, double kappa) {this->phi = degreeToRadian(phi);this->omega = degreeToRadian(omega);this->kappa = degreeToRadian(kappa);}////// \brief resetCamera/// \param lat 維度/// \param lon 經(jīng)度/// \param height 高程///void resetCamera(const QString &imgNumber, double lat, double lon, double height) {this->imgNumber = imgNumber;this->cameraLat = lat;this->cameraLon = lon;this->height = height;}////// \brief compute 計(jì)算相機(jī)位置 + 旋轉(zhuǎn)矩陣///void compute();private:////// \brief computeLatlon 根據(jù)像點(diǎn)坐標(biāo)計(jì)算經(jīng)緯度坐標(biāo)(共線(xiàn)方程的邏輯)/// @param vector2d uv x軸，y軸坐標(biāo)///Eigen::Vector2d computeLatlon(Eigen::Vector2d &uv);////// \brief degreeToRadian 度轉(zhuǎn)弧度/// \param degree/// \return///double degreeToRadian(double degree) { return degree * M_PI / 180.0; }//    double computeDeltaX();//    double computeDeltaY();Eigen::Vector3d cameraGeo;Eigen::Matrix3d rMatrix;Transform transTool;CameraModel intrinsic; // 內(nèi)參數(shù)矩陣QString imgNumber;     // 影像編號(hào)double cameraLat;      // 相機(jī)外參數(shù)矩陣double cameraLon;double height;double phi;double omega;double kappa;
};

Eigen::Vector2d ComputeBoundingBox::computeLatlon(Eigen::Vector2d &uv) {double u = uv.x();double v = uv.y();Eigen::Vector3d cameraSpace(u, v, -this->intrinsic.f); // 像空間坐標(biāo)double r = qSqrt(pow(u - this->intrinsic.u0, 2) + pow(v - this->intrinsic.v0, 2));// (x-x0) * (k1*r^2 + k2*r^4 + k3*r^6 + k4*r8) + p1*(r^2 + 2*(x-x)^2) + 2p2*(x-x0)(y-y0) + alpha*(x-x0) +// beta*(y-y0)double deltaX = (u - this->intrinsic.u0) * (this->intrinsic.k1 * pow(r, 2) + this->intrinsic.k2 * pow(r, 4) +this->intrinsic.k3 * pow(r, 6) + this->intrinsic.k4 * pow(r, 8)) +this->intrinsic.p1 * (pow(r, 2) + 2 * pow((u - this->intrinsic.u0), 2)) +2 * this->intrinsic.p2 * (u - this->intrinsic.u0) * (v - this->intrinsic.v0) +this->intrinsic.alpha * (u - this->intrinsic.u0) + this->intrinsic.beta * (v - this->intrinsic.v0);double deltaY = (v - this->intrinsic.v0) * (this->intrinsic.k1 * pow(r, 2) + this->intrinsic.k2 * pow(r, 4) +this->intrinsic.k3 * pow(r, 6) + this->intrinsic.k4 * pow(r, 8)) +this->intrinsic.p2 * (pow(r, 2) + 2 * pow((v - this->intrinsic.v0), 2)) +2 * this->intrinsic.p1 * (u - this->intrinsic.u0) * (v - this->intrinsic.v0);Eigen::Vector3d cameraOffset(deltaX - this->intrinsic.u0, deltaY - this->intrinsic.v0,0);                              // 像點(diǎn)坐標(biāo)偏移Eigen::Vector3d cameraSpaceTrue = cameraSpace + cameraOffset; // 實(shí)際的像點(diǎn)坐標(biāo)[最后一位該如何求解]Eigen::Vector3d worldCoordBa =this->rMatrix * cameraSpaceTrue * this->intrinsic.pixelSize; // (xBa, yBa, zBa) pixelSize這個(gè)參數(shù)多余worldCoordBa =Eigen::Vector3d(worldCoordBa.x() / worldCoordBa.z() * 1000, worldCoordBa.y() / worldCoordBa.z() * 1000, 0);//    qDebug() << worldCoordBa.x() << " " << worldCoordBa.y() << " " << worldCoordBa.z();Eigen::Vector3d worldCoord = worldCoordBa + this->cameraGeo; // 真正的坐標(biāo)//    std::cout << worldCoord.x() << " " << worldCoord.y() << " " << worldCoord.z();//    PJ_COORD latlonh = proj_coord(cameraLon, cameraLat, height, 0);PJ_COORD geoxyz = proj_coord(worldCoord.x(), worldCoord.y(), worldCoord.z(), 0); // 地心坐標(biāo)PJ_COORD latlon = this->transTool.backward(geoxyz);//    std::cout << "lat:" << latlon.lpz.lam << " ,lon:" << latlon.lpz.phi << ", " << latlon.lpz.z;return Eigen::Vector2d(latlon.lp.phi, latlon.lp.lam); // lat and lon
}

void ComputeBoundingBox::compute() {PJ_COORD latlonh = proj_coord(cameraLon, cameraLat, height, 0);PJ_COORD geoxyz = transTool.forward(latlonh);this->cameraGeo = Eigen::Vector3d(geoxyz.xyz.x, geoxyz.xyz.y, geoxyz.xyz.z); // 地心坐標(biāo)// 計(jì)算繞X軸旋轉(zhuǎn)的旋轉(zhuǎn)矩陣 Rx(φ)Eigen::Matrix3d Rx;Rx << 1, 0, 0, 0, std::cos(phi), -std::sin(phi), 0, std::sin(phi), std::cos(phi);// 計(jì)算繞Y軸旋轉(zhuǎn)的旋轉(zhuǎn)矩陣 Ry(ω)Eigen::Matrix3d Ry;Ry << std::cos(omega), 0, std::sin(omega), 0, 1, 0, -std::sin(omega), 0, std::cos(omega);// 計(jì)算繞Z軸旋轉(zhuǎn)的旋轉(zhuǎn)矩陣 Rz(κ)Eigen::Matrix3d Rz;Rz << std::cos(kappa), -std::sin(kappa), 0, std::sin(kappa), std::cos(kappa), 0, 0, 0, 1;// 計(jì)算總的旋轉(zhuǎn)矩陣 R_total = Rz(κ) * Ry(ω) * Rx(φ)this->rMatrix = Rz * Ry * Rx;// 像素點(diǎn)畸變糾正Eigen::Vector2d lb(0, 0);                                          // 左下角為起點(diǎn)Eigen::Vector2d rb(this->intrinsic.width, 0);                      // 右下角坐標(biāo)Eigen::Vector2d rt(this->intrinsic.width, this->intrinsic.height); // 右上角坐標(biāo)Eigen::Vector2d lt(0, this->intrinsic.height);                     // 左上角成果std::vector<Eigen::Vector2d> latlonVec = {lb, rb, rt, lt};qDebug() << "####################" << this->imgNumber << "########################";for (Eigen::Vector2d &pixelCoord : latlonVec) {pixelCoord = this->computeLatlon(pixelCoord);qDebug() << "lat: " << QString::number(pixelCoord.x(), 'f', 10)<< ",lon:" << QString::number(pixelCoord.y(), 'f', 10);}
}

效果圖