基于B_spline 的非刚性形变

上图中的（f）图就使用基于B_spline的非刚性形变后的图像。刚性形变保持了像素的平直型，非刚性形变就会破坏像素的平直性，但是图像的信息又不会有太大的变化。

计算公式

T(x)=∑l=04∑m=04∑m=04Bl(u)Bm(v)Bn(w)ϕi+l,j+m,k+nT(x)=∑l=04∑m=04∑m=04Bl(u)Bm(v)Bn(w)ϕi+l,j+m,k+n

T(x) =\sum_{l=0}^{4} \sum_{m=0}^{4}\sum_{m=0}^{4} B_l(u)B_m(v)B_n(w) \phi_{i+l,j+m,k+n}

此公式是三维图像的坐标变化公式。如上图所示，将一副图像用一定间隔的网格分割成不同的区域。ϕi,j,kϕi,j,k\phi_{i,j,k} 是格点的坐标。原图像中的任意一个像素点的坐标变换到新图像中的坐标，又该点周围 4∗4∗44∗4∗44*4*4 的控制点来计算。
设网格点的间隔分别大小为：δx,δy,δzδx,δy,δz\delta_x,\delta_y, \delta_z
i,j,ki,j,ki,j,k 格点坐标：(int)(x/δx)(int)(x/δx)(int)(x/\delta_x)
u,v,wu,v,wu,v,w 分别问x,y,z方向上的相对位置： u=(x/δx)−(int)(x/δx)u=(x/δx)−(int)(x/δx)u =(x/\delta_x) - (int)(x/\delta_x)
B_spline 基函数：

B0(t)=(−t3+3t2−3t+1)/6B1(t)=(3t3−6t2+4)/6B2(t)=(−3t3+3t2+3t+1)/6B3(t)=(t3)/6B0(t)=(−t3+3t2−3t+1)/6B1(t)=(3t3−6t2+4)/6B2(t)=(−3t3+3t2+3t+1)/6B3(t)=(t3)/6

\begin{array}{lc}B_0(t) = (-t^3+3t^2-3t+1)/6 \\B_1(t) = (3t^3-6t^2+4)/6 \\B_2(t) = (-3t^3+3t^2+3t+1)/6 \\B_3(t) = (t^3)/6\end{array}

如果是二维图像，那么网格也变成二维，新坐标由领域的4*4个控制点来计算。
在实际计算新图像的时候，可以直接正向映射，把原图像中的每一像素点的坐标映射到新图像中。如果多个像素点映射到了同一个位置，那么就要考虑合并像素的问题，而新图像的有些位置没有像素映射过来。所有一般不采用正向映射，而采用反向映射。反向映射计算新图像中的每个像素点在原图像中的位置，然后赋予对应像素值，如果原位置不在整数上就通过插值计算。

C++ 使用opencv 实现代码


void ImageTransform::B_spline_form(Mat &srcimg, Mat &dstimg, Mat &offset)
{int delta_x = 32;//越大，支撑域越大，变形越温和int delta_y = 32;//y 方向的支撑域int grid_rows = (int)(srcimg.rows / delta_x) + 1 + 3;int grid_cols = (int)(srcimg.cols / delta_y) + 1 + 3;Mat noiseMat(grid_rows, grid_cols, CV_32FC2);default_random_engine rand_engine(static_cast<unsigned int>(std::time(NULL)));uniform_real_distribution<float> uniform_distribution(-20, 20);for (int row = 0; row < grid_rows; row++){for (int col = 0; col < grid_cols; col++){noiseMat.at<Vec2f>(row, col)[0] = uniform_distribution(rand_engine);noiseMat.at<Vec2f>(row, col)[1] = uniform_distribution(rand_engine);}}// B_spline 基函数auto B = [](int flag, float t)->double {if (flag == 0)return (1 - t*t*t + 3 * t*t - 3 * t) / 6.0;else if (flag == 1)return (4 + 3 * t*t*t - 6 * t*t) / 6.0;else if (flag == 2)return (1 - 3 * t*t*t + 3 * t*t + 3 * t) / 6.0;else if (flag == 3)return (t*t*t / 6.0);elsereturn 0.0;};// B_spline 变形 for (int x = 0; x < srcimg.rows; x++){for (int y = 0; y < srcimg.cols; y++){int i = (int)(x / delta_x)/* - 1*/;int j = (int)(y / delta_y)/* - 1*/;float u = (float)x / delta_x - i /*- 1*/;/************************/float v = (float)y / delta_y - j/* - 1*/;float pX[4], pY[4];for (int k = 0; k < 4; k++){pX[k] = (float)B(k, u);pY[k] = (float)B(k, v);}float Tx = 0;float Ty = 0;for (int m = 0; m < 4; m++){for (int n = 0; n < 4; n++){int control_point_x = i + m /*+1*/;int control_point_y = j + n /*+1*/;float temp = pY[n] * pX[m];Tx += temp*(noiseMat.at<Vec2f>(control_point_x, control_point_y)[0]);Ty += temp*(noiseMat.at<Vec2f>(control_point_x, control_point_y)[1]);}}offset.at<Vec2f>(x, y)[0] = Tx;offset.at<Vec2f>(x, y)[1] = Ty;}}//反向映射，双线性插值for (int row = 0; row < dstimg.rows; row++){for (int col = 0; col < dstimg.cols; col++){float src_x = row + offset.at<Vec2f>(row, col)[0];float src_y = col + offset.at<Vec2f>(row, col)[1];int x1 = int(src_x);int y1 = int(src_y);int x2 = x1 + 1;int y2 = y1 + 1;if (x1<0 || x1>(srcimg.rows - 2) || y1<0 || y1>(srcimg.cols - 2))//越界dstimg.at<Vec3b>(row, col) = Vec3b(0, 0, 0);else{Vec3b pointa = srcimg.at<Vec3b>(x1, y1);Vec3b pointb = srcimg.at<Vec3b>(x2, y1);Vec3b pointc = srcimg.at<Vec3b>(x1, y2);Vec3b pointd = srcimg.at<Vec3b>(x2, y2);uchar B = (uchar)((x2 - src_x)*(y2 - src_y)*pointa[0] - (x1 - src_x)*(y2 - src_y)*pointb[0] - (x2 - src_x)*(y1 - src_y)*pointc[0] + (x1 - src_x)*(y1 - src_y)*pointd[0]);uchar G = (uchar)((x2 - src_x)*(y2 - src_y)*pointa[1] - (x1 - src_x)*(y2 - src_y)*pointb[1] - (x2 - src_x)*(y1 - src_y)*pointc[1] + (x1 - src_x)*(y1 - src_y)*pointd[1]);uchar R = (uchar)((x2 - src_x)*(y2 - src_y)*pointa[2] - (x1 - src_x)*(y2 - src_y)*pointb[2] - (x2 - src_x)*(y1 - src_y)*pointc[2] + (x1 - src_x)*(y1 - src_y)*pointd[2]);dstimg.at<Vec3b>(row, col) = Vec3b(B, G, R);}}}}int main()
{Mat srcimg = imread("face.jpg");Mat dstimg(srcimg.rows, srcimg.cols, srcimg.type());Mat offset(srcimg.rows, srcimg.cols, CV_32FC2);ImageTransform::B_spline_form(srcimg, dstimg, offset);Mat diff = dstimg - srcimg;imshow("原图像", srcimg);imshow("B_spline_form", dstimg);imshow("diff", diff);//  imshow("offset", offset);waitKey(0);return 0;
}

运行结果

参考论文
参考文献
b-spline学习-系数计算及程序实践
双线性插值
图像配准之《常用图像变换模型》简述

基于B_spline 的非刚性形变相关推荐

《Master Opencv...读书笔记》非刚性人脸跟踪 II
上一篇博文中,我们了解了系统的功能和模块,明确了需要采集哪些类型的样本点及利用类的序列化的保存方式.这次将介绍几何约束模块,通过统计形态分析法(Statistical Shape Analysis, ...
软件工程 / 为什么基于接口而非实现编程？
基于接口而非实现编程(基于抽象而非实现编程)的目的是解耦. 这里面接口的含义可以理解为 dll 或者 so 文件对应的头文件中提供的函数列表,或者理解为C++中的抽象类. 该原则可以将接口和实现分离, ...
基于参考点的非支配遗传算法-NSGA-III（二）
上一篇我们讲了有关非支配遗传算法NSGA-III的非约束过程,接下来这一篇我们讲一下NSGA-III约束实现以及扩展自适应方法.同理,我们先列一下我们参考的博客.代码以及论文. 文章目录参考博客参 ...
非刚性配准（Non-rigid ICP ）
原文https://blog.csdn.net/linmingan/article/details/79270874?utm_medium=distribute.pc_relevant.none-ta ...
为什么基于接口而非实现编程？
如何解读原则中的"接口"二字? "基于接口而非实现编程"这条原则的英文描述是:"Program to an interface, not an imp ...
基于Android9的非root环境下frida-gadget持久化
基于Android9的非root环境下frida持久化博客: http://www.zhuoyue360.com 参考: 小肩膀安卓系统沙箱课程 https://bbs.pediy.com/thre ...
MATLAB 数学应用微分方程常微分方程求解非刚性ODE
本文介绍两个使用 ode45 来求解非刚性常微分方程的示例.MATLAB拥有三个非刚性 ODE 求解器. ode45 ode23 ode113 对于大多数非刚性问题,ode45 的性能最佳.但对于允许 ...
刚性微分方程与非刚性的大概区分（自看，求指点）
(这里的内容是作为自己的一个粗略的总结,不确定是否正确,希望有大佬能够更明确的指出其中的错误,作出指导) 对于刚性和非刚性微分方程的区分,可以简单的转变为在将原方程转换为常微分方程组后,进行一个简单的 ...
某系统采用基于优先权的非抢占式进程调度策略，完成一次进程调度和进程切换的系统时间开销为 1μs。
某系统采用基于优先权的非抢占式进程调度策略,完成一次进程调度和进程切换的系统时间开销为 1μs.在 T 时刻就绪队列中有 3 个进程 P1.P2 和 P3,其在就绪队列中的等待时间.需要的 CPU 时 ...
【JAVA程序设计】（C00073）基于SSH（非maven）便利店管理系统-有文档
@TOC 项目简介基于ssh框架非maven开发的便利店管理系统共分为三个角色:系统管理员.销售管理员角色包含以下功能: 系统管理.用户管理.商品管理.采购管理.库存管理.销售管理.财务管理(成本 ...

基于B_spline 的非刚性形变

基于B_spline 的非刚性形变相关推荐

最新文章

热门文章