三维薄板样条，用于三维模型变形（c++）

刚接触变形算法，其实不想碰三维变形算法，菜鸡一个，但是导师一直让试试，那就试试，事实证明薄板样条并不特别适用于古建筑点云变形，第一次写文章只为记录一下我曾经学过，不然按照本人的记忆力，基本都忘完。

薄板样条一般用于图像匹配，二维图像处理无敌，三维真的不是那么太强，起码基本只能二维变形，将三维视为二维，忽略第三维信息，用二维图像提取轮廓，在用相似三角形进行变形，着实凑巧，这里说的是基于轮廓线与薄板样条的三维自动化建模方法毕竟论文种可以完全避免同时三个维度的变形，咱也不知道说的对不对。经过这些天的死磕研究，自己写的三维薄板样条的代码总是并不是特别好的变形。简单的圆柱体、长方体、正方体。提供控制点和形变坐标就可以变形，然而使用一些点云的时候，最后通过薄板样条计算出的点各式各样，五花八门。

失败的就不放了，丢人现眼。本来不想碰拉普拉斯的，看来还得去自己编。

static double tps_base_func(double r)//用于计算基函数的r
{if (r == 0.0)return 0.0;elsereturn r*r * log(r*r);
}
static void calc_tps()
{// You We need at least 3 points to define a planeif (control_points.size() < 3)return;unsigned p = control_points.size();// Allocate the matrix and vectormatrix<double> mtx_l(p + 4, p + 4);//L矩阵  matrix<double> mtx_v(p + 4, 3);//包含了v 和  omatrix<double> mtx_orig_k(p, p);//K矩阵matrix<double>w(p, 4);double a = 0.0;for (unsigned i = 0; i < p; ++i){for (unsigned j = i + 1; j < p; ++j){mtx_l(i, j) = mtx_l(j, i) =mtx_orig_k(i, j) = mtx_orig_k(j, i) = tps_base_func((control_points[i] - control_points[j]).len()); //tps_base_func(elen);//???//Uij传给上下对角//a += elen * 2; // same for upper & lower tri}}a /= (double)(p*p);//a=a/p^2// Fill the rest of Lfor (unsigned i = 0; i < p; ++i){// diagonal: reqularization parameters (lambda * a^2)对角线重新量化参数mtx_l(i, i) = mtx_orig_k(i, i) = regularization * (a*a);//变形程度 目前为0// P (p x 3, upper right) 已改 4*P 增加Ymtx_l(i, p + 0) = 1.0;mtx_l(i, p + 1) = control_points[i].x;mtx_l(i, p + 2) = control_points[i].y;mtx_l(i, p + 3) = control_points[i].z;//fill P up right and down leftmtx_l(p + 0, i) = 1.0;mtx_l(p + 1, i) = control_points[i].x;mtx_l(p + 2, i) = control_points[i].y;mtx_l(p + 3, i) = control_points[i].z;}// O (3 x 3, lower right)  改4*4for (unsigned i = p; i < p + 4; ++i){for (unsigned j = p; j < p + 4; ++j){mtx_l(i, j) = 0.0;}}// Fill the right hand vector V   maybe is matrixfor (unsigned j = 0; j < 3; ++j){for (int i = 0; i < 4; i++){mtx_v(p + i, j) = 0;}}for (size_t i = 0; i < TransformationPoint.size(); i++){mtx_v(i, 0) = TransformationPoint[i].x;mtx_v(i, 1) = TransformationPoint[i].y;mtx_v(i, 2) = TransformationPoint[i].z;}cout << "--------------------L矩阵-----------------" << endl;for (unsigned i = 0; i < p+4; i++){for (unsigned j = 0; j < p+4; j++){cout << mtx_l(i, j) << "\t";}cout << endl;}cout << " ----------------v矩阵-----------------------" << endl;for (int i = 0; i < p+4; i++){for (int j = 0; j < 3; j++){cout << mtx_v(i, j) << "\t";}cout << endl;}if (0!= LU_Solve(mtx_l, mtx_v)&&1== LU_Solve(mtx_l, mtx_v))//0为成功  1为奇异矩阵 2 行数不对等{puts("Singular matrix! Aborting.");cout << "计算不成功,奇异矩阵" << endl;exit(1);}cout << "---------------------------系数----------------------------------" << endl;for (int i = 0; i < p + 4; i++){for (int j = 0; j < 3; j++){cout << mtx_v(i, j) << "\t";}cout << endl;}//WxU  WyU  WzU权重分配float WxU = 0;float WyU = 0;float WzU = 0;for (int i = 0; i < EnteredPoint.size(); i++){double x= mtx_v(p + 0, 0) + mtx_v(p + 1, 0)*EnteredPoint[i].x+ mtx_v(p + 2, 0)*EnteredPoint[i].y + mtx_v(p + 3, 0)*EnteredPoint[i].z;double y= mtx_v(p + 0, 1) + mtx_v(p + 1, 1)*EnteredPoint[i].x+ mtx_v(p + 2, 1)*EnteredPoint[i].y + mtx_v(p + 3, 1)*EnteredPoint[i].z;double z= mtx_v(p + 0, 2) + mtx_v(p + 1, 2)*EnteredPoint[i].x+ mtx_v(p + 2, 2)*EnteredPoint[i].y + mtx_v(p + 3, 2)*EnteredPoint[i].z;for (int j = 0; j < p; j++)//权重系数很小{WxU += mtx_v(j, 0) * tps_base_func((EnteredPoint[i] - control_points[j]).len());WyU += mtx_v(j, 1) * tps_base_func((EnteredPoint[i] - control_points[j]).len());WzU += mtx_v(j, 2) * tps_base_func((EnteredPoint[i] - control_points[j]).len());}//cout <<"v    "<< x+WxU << "   " << y + WyU << " " << z + WzU << "" << endl;//cout <<  WxU << " " <<  WyU << "  " <<  WzU << endl;std::cout << "计算第" << i << "次" << "   " << "还剩" << EnteredPoint.size() - i << "待计算" << std::endl;cout <<" v  "<<x+ WxU << " " <<y+ WyU << "  " << z+WzU << endl;AfterDeformPoint.push_back(Vec(x + WxU, y + WyU, z + WzU));}std::cout << "计算完成" << endl;// Calc bending energyfor (int i = 0; i < p; ++i){w(i, 0) = mtx_v(i, 0);//    弯曲能量部分待研究w(i, 1) = mtx_v(i, 1);w(i, 2) = mtx_v(i, 2);}matrix<double> be = prod(prod<matrix<double> >(trans(w), mtx_orig_k), w);bending_energy = be(0, 0);cout << "弯曲能量为=  " << bending_energy << endl;}
void ReadControlPoint()
{ifstream infile("E:\\ExperimentalData\\expData\\test\\gd.txt");float a, b, c;while (infile >> a >> b >> c){control_points.push_back(Vec(a, b, c));}}
void ReadTransformationPoint()
{ifstream tranfile("E:\\ExperimentalData\\expData\\test\\mv.txt");//C:\\Users\\Administrator\\Desktop//E:\\ExperimentalData\\expData\\test\\zhenshi700.txtfloat a, b, c,N1,N2,N3,N4,N5,N6;while (tranfile>>a>>b>>c)//tranfile >> a>> b >> c>>N1>>N2>>N3>>N4>>N5>>N6{TransformationPoint.push_back(Vec(a, b, c));}
}
void reaeAllPoints()
{//c_obj 定义一个网格对象int mask;     //定义mask//注意your filePath不能有中文路径vcg::tri::io::ImporterOBJ<MyMesh>::Open(m_obj, "E:\\ExperimentalData\\expData\\test\\524.obj", mask);//为每个点计算法线并归一化vcg::tri::RequirePerVertexNormal(m_obj);vcg::tri::UpdateNormal<MyMesh>::PerVertexNormalized(m_obj);std::vector <MyVertex> &vs = m_obj.vert;std::vector<MyEdge> &es = m_obj.edge;std::vector<MyFace> &fs = m_obj.face;for (auto i = 0; i < m_obj.vert.size(); ++i){MyVertex v = vs[i];EnteredPoint.push_back(Vec(v.P().X(), v.P().Y(), v.P().Z()));}
}
void OutputPoints()
{MyVertex v;for (int i = 0; i < EnteredPoint.size(); i++){v.P().X() = AfterDeformPoint[i].x;v.P().Y() = AfterDeformPoint[i].y;v.P().Z() = AfterDeformPoint[i].z;c_obj.vert.push_back(v);//cout << "x=  "<< AfterDeformPoint[i].x << "y=  " << AfterDeformPoint[i].y << "z= " << AfterDeformPoint[i].z << endl;}/*for (size_t i = 0; i < c_obj.vert.size(); i++){cout <<"v"<<"  "<< c_obj.vert[i].P().X() << "  " << c_obj.vert[i].P().Y() << "   " << c_obj.vert[i].P().Z()<< endl;}*/c_obj.edge = m_obj.edge;c_obj.face = m_obj.face;//vcg::tri::BuildMeshFromCoordVectorIndexVector(c_obj, c_obj.vert, c_obj.face);//vcg::tri::Octahedron(c_obj);vcg::tri::io::ExporterOBJ<MyMesh>::Save(c_obj, "C:\\Users\\Administrator\\Desktop\\c_obj.obj", 0);
}
void outDataToTxt()
{ofstream location_out;///*/注意路径前面不能有！空格否则会读不进去坐标//location_out.open("E:\\ExperimentalData\\expData\\test\\location_out.txt", std::ios::out | std::ios::app);  //以写入和在文件末尾添加的方式打开.txt文件，没有的话就创建该文件。if (!location_out.is_open())cout << "open file error" << endl;for (int i = 0; i < EnteredPoint.size(); i++){location_out << "v " << AfterDeformPoint[i].x << " " << AfterDeformPoint[i].y << " " << AfterDeformPoint[i].z << endl;}location_out.close();}int main()
{reaeAllPoints();ReadControlPoint();ReadTransformationPoint();calc_tps();outDataToTxt();system("pause");return 0;}

代码水平有限，也是自己摸索。真实点云作变形，计算出的模型顶点坐标会是特别大的数简直离谱。懂的大佬有没有知道古建筑三维变形最好用什么变形算法~难道只有拉普拉斯吗？

三维薄板样条，用于三维模型变形（c++）相关推荐

TPS(薄板样条) 2D 插值
如何测量2个2D形状之间的相似性?2D形状的形变使这个问题变得非常复杂. 很多情况下,我们无法准确的描述这种形变,使用2D插值是一个很好的办法. TPS(薄板样条)插值是常用的2D插值方法.它的物理意 ...
2022CVPR：一种用于空间变形的文本注意网络鲁棒场景文本图像超分辨率
TATT:A Text Attention Network for Spatial Deformation Robust Scene Text Image Super-resolution 单位:香港 ...
图像配准系列之“Sift特征点+薄板样条变换+FFD变换”配准方法
上篇文章中我们讲了"Sift+TPS"的配准方法: 图像配准系列之"Sift特征点+薄板样条变换"配准方法我们知道,TPS薄板样条变换(简称TPS变换)与FF ...
Thin Plate Spline TPS薄板样条变换基础理解
什么是图像扭曲问题? 给定控制点和相应位移点稀疏对应集,我们需要找到一个映射,且两点之间的尽可能平滑. 一维空间举例,绿色为对应集,需找到蓝色曲线映射,满足形变后控制点重合且之间连线平滑 ...
三维匹配_基于三维模型的目标识别和分割在杂乱的场景中的应用
作者:仲夏夜之星来源:3D视觉工坊公众号链接: 基于三维模型的目标识别和分割在杂乱的场景中的应用在杂波和遮挡情况下,对自由形式物体的识别及分割是一项具有挑战性的任务.本文提出了一种新的基于三维模 ...
osg 三维gis开发_OSG三维模型初探
最近在研究OSG开发,准备用OSG+OSGEARTH开发一套三维地形GIS系统,目前研究在VS2008下把OSG-2.8.3(Debug和Release)和OSEARTH-2.0.0(Release, ...
Thin Plate Spline薄板样条
自变量 x = [ x 1 , x 2 ] \mathbf{x}=[x_1,x_2] x=[x1,x2] 是2维空间中的一个点,函数值 y = [ y 1 , y 2 ] \mathbf{y}=[ ...
基于MATLAB三维重构开题报告,三维模型开题报告范文
1. 毕业设计开题报告立体依据怎么写为了点明论文的研究对象.研究内容.研究目的,对总标题加以补充.解说,有的论文还可以加副标题.非凡是一些商榷性的论文,一般都有一个副标题,如在总标题下方,添上&qu ...
三维激光扫描后处理软件_三维激光扫描应用于工程测量
三维扫描技术是继GIS技术后的又一项新测绘技术,通过对物体进行快速.非接触式的激光扫描,从而得到物体的高精度三维点云数据,生成物体的三维虚拟模型. 在使用全站仪和GPS这样的传统测量仪器所获得的测量数 ...

三维薄板样条，用于三维模型变形（c++）

三维薄板样条，用于三维模型变形（c++）相关推荐

最新文章

热门文章