念念不忘,必有回响-1
条评论废了很大的功夫，从什么都不知道到完成光束法平差以及密集匹配，把POS数据添加上得到与POS数据精度较为一直的坐标已经快要耗尽自己的全部精力了，不过走完这一步一切都变得简单了。经过一段时间的努力终于把结合控制点的平差方法给添加上去了。
想想是不是有必要讲讲将控制点加入进行光束法平差的方法，实际上通过ceres solver我构建了两类方程，第一类为加密点方程，在此方程中同时调整相机参数和三维点的坐标，另一类为控制点误差方程，在控制点误差方程中仅仅调整相机参数，在此两类方程的约束下进行光束法平差，实际上的代码为：
struct BundlerReprojectionError
{
BundlerReprojectionError(double observed_x,double observed_y,double flen)
: _observed_x(observed_x), _observed_y(observed_y) ,_focalLength(flen){}
//************************************
// 函数名: operator()
// 功能描述: 重载括号操作符计算残差，在这里我们并不对相机内参进行优化，假设相机内参定标得到
// 作者: W.W.Frank
// 返回值: bool
// 输入参数： const T* const camera,const T* const point3D 相机参数（0，1，2为角元素，3，4，5为线元素，6为焦距），得到的三维点坐标
// 输出参数： T* residuals 得到的残差
//************************************
template<class T>
bool operator()(const T* const camera,const T* const point3D,T* residuals) const
{
T p[3], t[3];
t[0] = point3D[0]-camera[3];
t[1] = point3D[1]-camera[4];
t[2] = point3D[2]-camera[5];
T dPhi = camera[0];
T dOmega = camera[1];
T dKappa = camera[2];
T a1 = cos(dPhi)*cos(dKappa) - sin(dPhi)*sin(dOmega)*sin(dKappa);
T a2 = -cos(dPhi)*sin(dKappa) - sin(dPhi)*sin(dOmega)*cos(dKappa);
T a3 = -sin(dPhi)*cos(dOmega);
T b1 = cos(dOmega)*sin(dKappa);
T b2 = cos(dOmega)*cos(dKappa);
T b3 = -sin(dOmega);
T c1 = sin(dPhi)*cos(dKappa) + cos(dPhi)*sin(dOmega)*sin(dKappa);
T c2 = -sin(dPhi)*sin(dKappa) + cos(dPhi)*sin(dOmega)*cos(dKappa);
T c3 = cos(dPhi)*cos(dOmega);
p[0]=a1*t[0]+a2*t[1]+a3*t[2];
p[1]=b1*t[0]+b2*t[1]+b3*t[2];
p[2]=c1*t[0]+c2*t[1]+c3*t[2];
T xp = - p[0] / p[2];
T yp = - p[1] / p[2];
T predicted_x = _focalLength * xp;
T predicted_y = _focalLength * yp;
//误差项的估计
residuals[0] = predicted_x - T(_observed_x);
residuals[1] = predicted_y - T(_observed_y);
return true;
}
//************************************
// 函数名: Create()
// 功能描述: create 函数，为了隐藏该结构体的构造函数
// 作者: W.W.Frank
// 返回值: bool
// 输入参数： const double observed_x, const double observed_y 观测的二维点的坐标
// 输出参数： /
//************************************
static ceres::CostFunction* Create(const double observed_x, const double observed_y,const double flen)
{
return (new ceres::AutoDiffCostFunction<BundlerReprojectionError,2,6,3>(
new BundlerReprojectionError(observed_x,observed_y,flen)));
}
double _observed_x;
double _observed_y;
double _focalLength;
};
//控制点误差结构体，控制点只调整
struct BundlerReprojectionErrorConst
{
BundlerReprojectionErrorConst(double observed_x,double observed_y,double flen,GeoInfo_Point3f pnd3D)
: _observed_x(observed_x), _observed_y(observed_y) ,_focalLength(flen),_Xpnt(pnd3D.x),_Ypnt(pnd3D.y),_Zpnt(pnd3D.z){}
//************************************
// 函数名: operator()
// 功能描述: 重载括号操作符计算残差，在这里我们并不对相机内参进行优化，假设相机内参定标得到
// 作者: W.W.Frank
// 返回值: bool
// 输入参数： const T* const camera ；0，1，2为角元素，3，4，5为线元素，6为焦距），得到的三维点坐标（控制点坐标进行调整）
// 输出参数： T* residuals 得到的残差
//************************************
template<class T>
bool operator()(const T* const camera,T* residuals) const
{
T p[3], t[3];
t[0] = _Xpnt-camera[3];
t[1] = _Ypnt-camera[4];
t[2] = _Zpnt-camera[5];
T dPhi = camera[0];
T dOmega = camera[1];
T dKappa = camera[2];
T a1 = cos(dPhi)*cos(dKappa) - sin(dPhi)*sin(dOmega)*sin(dKappa);
T a2 = -cos(dPhi)*sin(dKappa) - sin(dPhi)*sin(dOmega)*cos(dKappa);
T a3 = -sin(dPhi)*cos(dOmega);
T b1 = cos(dOmega)*sin(dKappa);
T b2 = cos(dOmega)*cos(dKappa);
T b3 = -sin(dOmega);
T c1 = sin(dPhi)*cos(dKappa) + cos(dPhi)*sin(dOmega)*sin(dKappa);
T c2 = -sin(dPhi)*sin(dKappa) + cos(dPhi)*sin(dOmega)*cos(dKappa);
T c3 = cos(dPhi)*cos(dOmega);
p[0]=a1*t[0]+a2*t[1]+a3*t[2];
p[1]=b1*t[0]+b2*t[1]+b3*t[2];
p[2]=c1*t[0]+c2*t[1]+c3*t[2];
T xp = - p[0] / p[2];
T yp = - p[1] / p[2];
T predicted_x = _focalLength * xp;
T predicted_y = _focalLength * yp;
//误差项的估计
residuals[0] = predicted_x - T(_observed_x);
residuals[1] = predicted_y - T(_observed_y);
return true;
}
//************************************
// 函数名: Create()
// 功能描述: create 函数，为了隐藏该结构体的构造函数
// 作者: W.W.Frank
// 返回值: bool
// 输入参数： const double observed_x, const double observed_y 观测的二维点的坐标
// 输出参数： /
//************************************
static ceres::CostFunction* Create(const double observed_x, const double observed_y,const double flen,const GeoInfo_Point3f pnt)
{
return (new ceres::AutoDiffCostFunction<BundlerReprojectionErrorConst,2,6>(
new BundlerReprojectionErrorConst(observed_x,observed_y,flen,pnt)));
}
double _observed_x;
double _observed_y;
double _focalLength;
double _Xpnt;
double _Ypnt;
double _Zpnt;
};
通过构造的控制点误差方程直接添加到光束法平差的方程中就可以自动的进行迭代求解，虽然这样可以将控制点方程添加进行求解，然而在实际求解的过程中求解的精度与控制点的密度和控制点的分布具有较大的关系，如果控制点分布较差，且控制点较少，则求解精度难以得到保证，采用三幅影像5个控制点能够取得较高的平面精度。
实际上文章取了这么一个文艺的名字不仅仅是想要说说我控制点怎么添加上去的，主要是讲讲做这么一件事情的心路历程，研究生读下来好像做了很多东西，一开始做无人机，后来做航空摄影测量，做高光谱图像处理，做点云处理，然后到毕设做回高光谱图像处理，兜兜转转这么一个大圈，发现最后果然还是应了自己的专业：摄影测量，虽然在读硕士，可是很多时候也不是想做什么就能够做什么的，毕竟有很多方向不太能够出成果，如果进入这个方向做了三年没有取得什么研究成果那岂不是太可惜了，然后选了一个能够出成果的方向好好做，可是却不是自己想要做的。很多时候都是这样，为了成果不停的去换方向，然后什么都没有做好。我觉得做任何一个东西，研究也好软件也罢，如果不坚持恐怕什么都做不出来吧。
讲真话，搞了三年所谓的科研之后我感觉自己有些憎恨科研了，我不知道我的研究到底有什么意义，所做的工作不过是在很成熟的商业软件下做的并不怎么样的工作罢了，可是毕竟这是我的成果，如果一切都这么简单那我们做的又还有什么意义呢，起码我还不知道有国内做的很好的商也软件吧。不过这也成为了马太效应产生的原因吧，既然人家都做好了，我们何必要自己做，买现成的不就好了么，自己做的又没有人家的好…..好吧可是虽然不一定有商业软件做的好，可是毕竟我做出来了，坚持终于还是有点回报的吧，虽然只是一点小小的涟漪。