基于B_spline 的非刚性形变

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

• 计算公式

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

此公式是三维图像的坐标变化公式。如上图所示，将一副图像用一定间隔的网格分割成不同的区域。 ϕi,j,k ϕ i , j , k 是格点的坐标。原图像中的任意一个像素点的坐标变换到新图像中的坐标，又该点周围 4∗4∗4 4 ∗ 4 ∗ 4 的控制点来计算。

设网格点的间隔分别大小为： δx,δy,δz δ x , δ y , δ z

i,j,k i , j , k 格点坐标： (int)(x/δx) ( i n t ) ( x / δ x )

u,v,w u , v , w 分别问x,y,z方向上的相对位置： u=(x/δx)−(int)(x/δx) u = ( x / δ x ) − ( i n t ) ( x / δ x )

B_spline 基函数：

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

如果是二维图像，那么网格也变成二维，新坐标由领域的4*4个控制点来计算。

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

• C++ 使用opencv 实现代码

void ImageTransform::B_spline_form(Mat &srcimg, Mat &dstimg, Mat &offset)
{
    int delta_x = ;//越大，支撑域越大，变形越温和
    int delta_y = ;//y 方向的支撑域
    int grid_rows = (int)(srcimg.rows / delta_x) +  + ;
    int grid_cols = (int)(srcimg.cols / delta_y) +  + ;
    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(-, );

    for (int row = ; row < grid_rows; row++)
    {
        for (int col = ; col < grid_cols; col++)
        {
            noiseMat.at<Vec2f>(row, col)[] = uniform_distribution(rand_engine);
            noiseMat.at<Vec2f>(row, col)[] = uniform_distribution(rand_engine);
        }
    }

    // B_spline 基函数
    auto B = [](int flag, float t)->double {
        if (flag == )
            return ( - t*t*t +  * t*t -  * t) / ;
        else if (flag == )
            return ( +  * t*t*t -  * t*t) / ;
        else if (flag == )
            return ( -  * t*t*t +  * t*t +  * t) / ;
        else if (flag == )
            return (t*t*t / );
        else
            return ;
    };

    // B_spline 变形 
    for (int x = ; x < srcimg.rows; x++)
    {
        for (int y = ; 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[], pY[];
            for (int k = ; k < ; k++)
            {
                pX[k] = (float)B(k, u);
                pY[k] = (float)B(k, v);
            }

            float Tx = ;
            float Ty = ;
            for (int m = ; m < ; m++)
            {
                for (int n = ; n < ; 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)[]);
                    Ty += temp*(noiseMat.at<Vec2f>(control_point_x, control_point_y)[]);
                }
            }
            offset.at<Vec2f>(x, y)[] = Tx;
            offset.at<Vec2f>(x, y)[] = Ty;
        }
    }

    //反向映射，双线性插值
    for (int row = ; row < dstimg.rows; row++)
    {
        for (int col = ; col < dstimg.cols; col++)
        {
            float src_x = row + offset.at<Vec2f>(row, col)[];
            float src_y = col + offset.at<Vec2f>(row, col)[];
            int x1 = int(src_x);
            int y1 = int(src_y);
            int x2 = x1 + ;
            int y2 = y1 + ;
            if (x1< || x1>(srcimg.rows - ) || y1< || y1>(srcimg.cols - ))//越界
                dstimg.at<Vec3b>(row, col) = Vec3b(, , );
            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[] - (x1 - src_x)*(y2 - src_y)*pointb[] - (x2 - src_x)*(y1 - src_y)*pointc[] + (x1 - src_x)*(y1 - src_y)*pointd[]);
                uchar G = (uchar)((x2 - src_x)*(y2 - src_y)*pointa[] - (x1 - src_x)*(y2 - src_y)*pointb[] - (x2 - src_x)*(y1 - src_y)*pointc[] + (x1 - src_x)*(y1 - src_y)*pointd[]);
                uchar R = (uchar)((x2 - src_x)*(y2 - src_y)*pointa[] - (x1 - src_x)*(y2 - src_y)*pointb[] - (x2 - src_x)*(y1 - src_y)*pointc[] + (x1 - src_x)*(y1 - src_y)*pointd[]);

                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();
    return ;
}
           

运行结果

• 参考论文

  

  

• 参考文献

  b-spline学习-系数计算及程序实践

  双线性插值

  图像配准之《常用图像变换模型》简述