图像接缝裁剪（seam carving）算法实现-SIGGRAPH 2007

seam carving是SIGGRAPH 2007数字图形学年会上，以色列两位教授提出的算法，用于实现“内容保留”的图像伸缩。 出自论文《Seam Carving for Content-Aware Image Resizing》，作者的个人主页有对该算法的描述。

常规的图像伸缩基于采样和插值，图像内容会产生拉伸或压缩，比如图像中有人脸等内容，通常产生糟糕的效果。seam carving的基本思想是：图像中不同的位置“能量”不同，有重要内容的地方能量大，反之，可有可无的位置能量就小。这个能量函数可以是梯度，这有一定道理，梯度大的地方亮度变化剧烈、纹理丰富，一般是重要内容；梯度小的位置亮度变化小，比如蓝天，删除一条缝隙对图像影响较小。

考虑横向或者纵向寻找一条接缝，怎么样寻找这样一条能量最小的接缝？可以通过动态规划的方法实现。看了wiki百科：seam carving和这篇博客Eric Yuan's Blog，你会发现代码实现还是很简单的。

放上代码(Matlab):

main：

%%%%%%%%%%%%%%%%%%使用动态规划的方法实现seam carving%%   2016.12.09 南京%%%%%%%%%%%%%%%%%clc;clear;Im=imread('ori.jpg');gauss=fspecial('gaussian',[3,3],2);Blur=imfilter(Im,gauss,'same');[m,n,c]=size(Im);%灰度图Gray=rgb2gray(Im);%求梯度图hy=fspecial('sobel');hx=hy';Iy=imfilter(double(Gray),hy,'replicate');Ix=imfilter(double(Gray),hx,'replicate');Gradient=sqrt(Ix.^2+Iy.^2);%归一化max1=max(max(Gradient)');Gradient=Gradient./max1;%能量值图Energy=zeros(m,n);%路径图Path=zeros(m,n);tmp=0;for i=1:m    for j=1:n        if(i==1)            Energy(i,j)=Gradient(i,j);            Path(i,j)=0;        else            if(j==1)                tmp=which_min2(Energy(i-1,j),Energy(i-1,j+1));                Energy(i,j)=Gradient(i,j)+Energy(i-1,j+tmp);                Path(i,j)=tmp;            elseif(j==n)                tmp=which_min2(Energy(i-1,j-1),Energy(i-1,j));                Energy(i,j)=Gradient(i,j)+Energy(i-1,j-1+tmp);                Path(i,j)=tmp-1;            else                tmp=which_min3(Energy(i-1,j-1),Energy(i-1,j),Energy(i-1,j+1));                Energy(i,j)=Gradient(i,j)+Energy(i-1,j-1+tmp);                Path(i,j)=tmp-1;            end        end    endend%归一化%归一化max2=max(max(Energy)');Energy=Energy./max2;%能量最小路径的最后一行的纵坐标lastCol=find(Energy(m,:)==min(Energy(m,:)));col=lastCol(1);%描画出分割线Line=Im;for i=m:-1:1    Line(i,col,:)=[0,255,0];    col=col+Path(i,col);end%消除路径上的点Im=Im(:);for i=m:-1:1    Im(1*i+col)=[];    Im(2*i+col)=[];    Im(3*i+col)=[];    col=col+Path(i,col);endIm=reshape(Im,m,n-1,3);%Gradient是梯度图figure,imshow(Gradient);title('Gradient Image');%Energy是能量图figure,imshow(Energy);title('Energy Image');%Line是标注了分割线的图figure,imshow(Line);title('Image with Seam');%在Im上将分割线切掉figure,imshow(Im);title('after Cut Seam');

which_min2函数：

function tmp=which_min(x,y)if(min([x,y])==x)    tmp=0;else    tmp=1;end

which_min3函数：

function tmp=which_min3(x,y,z)if(min([x,y,z])==y)    tmp=1;elseif(min([x,y,z])==x)    tmp=0;    else    tmp=2;end

效果：

原图

能量（梯度）图

累加能量图

将要裁剪的接缝

删除接缝处的像素之后（因为仅删除一条线，和原图相比，变化甚微）

在网上搜集的C++版（使用OpenCV）

#include "opencv2/opencv.hpp"#include "opencv2/highgui/highgui.hpp"#include <stdio.h>#define NUM 220void calculateEnergy(cv::Mat& srcMat,cv::Mat& dstMat,cv::Mat& traceMat){srcMat.copyTo(dstMat);   //不用“=”，防止两个矩阵指向的都是同一个矩阵，现在只需要传里面的数值for (int i = 1;i < srcMat.rows;i++)  //从第2行开始计算{//第一列if (dstMat.at<float>(i-1,0) <= dstMat.at<float>(i-1,1)){dstMat.at<float>(i,0) = srcMat.at<float>(i,0) + dstMat.at<float>(i-1,0);traceMat.at<float>(i,0) = 1; //traceMat记录当前位置的上一行应取那个位置，上左为0，上中1，上右为2}else{dstMat.at<float>(i,0) = srcMat.at<float>(i,0) + dstMat.at<float>(i-1,1);traceMat.at<float>(i,0) = 2;}//中间列for (int j = 1;j < srcMat.cols-1;j++){float k[3];k[0] = dstMat.at<float>(i-1,j-1);k[1] = dstMat.at<float>(i-1,j);k[2] = dstMat.at<float>(i-1,j+1);int index = 0;if (k[1] < k[0])index = 1;if (k[2] < k[index])index = 2; dstMat.at<float>(i,j) = srcMat.at<float>(i,j) + dstMat.at<float>(i-1,j-1+index);traceMat.at<float>(i,j) = index;}//最后一列if (dstMat.at<float>(i-1,srcMat.cols-1) <= dstMat.at<float>(i-1,srcMat.cols-2)){dstMat.at<float>(i,srcMat.cols-1) = srcMat.at<float>(i,srcMat.cols-1) + dstMat.at<float>(i-1,srcMat.cols-1);traceMat.at<float>(i,srcMat.cols-1) = 1; }else{dstMat.at<float>(i,srcMat.cols-1) = srcMat.at<float>(i,srcMat.cols-1) + dstMat.at<float>(i-1,srcMat.cols-2);traceMat.at<float>(i,srcMat.cols-1) = 0;}}}// 找出最小能量线void getMinEnergyTrace(const cv::Mat& energyMat,const cv::Mat& traceMat,cv::Mat& minTrace){int row = energyMat.rows - 1;// 取的是energyMat最后一行的数据，所以行标是rows-1int index = 0;// 保存的是最小那条轨迹的最下面点在图像中的列标// 获得index，即最后那行最小值的位置for (int i = 1;i < energyMat.cols;i++){if (energyMat.at<float>(row,i) < energyMat.at<float>(row,index)){index = i;} // end if} // end for i = ...// 以下根据traceMat，得到minTrace，minTrace是多行一列矩阵{minTrace.at<float>(row,0) = index;int tmpIndex = index;for (int i = row;i > 0;i--){int temp = traceMat.at<float>(i,tmpIndex);// 当前位置traceMat所存的值if (temp == 0) // 往左走{tmpIndex = tmpIndex - 1;}else if (temp == 2) // 往右走{tmpIndex = tmpIndex + 1;} // 如果temp = 1，则往正上走，tmpIndex不需要做修改minTrace.at<float>(i-1,0) = tmpIndex;}}}// 删掉一列void delOneCol(cv::Mat& srcMat,cv::Mat& dstMat,cv::Mat& minTrace,cv::Mat& beDeletedLine){for (int i = 0;i < dstMat.rows;i++){int k = minTrace.at<float>(i,0);for (int j = 0;j < k;j++){dstMat.at<cv::Vec3b>(i,j)[0] = srcMat.at<cv::Vec3b>(i,j)[0];dstMat.at<cv::Vec3b>(i,j)[1] = srcMat.at<cv::Vec3b>(i,j)[1];dstMat.at<cv::Vec3b>(i,j)[2] = srcMat.at<cv::Vec3b>(i,j)[2];}for (int j = k;j < dstMat.cols-1;j++){if (j == dstMat.cols-1){int a = 1;}dstMat.at<cv::Vec3b>(i,j)[0] = srcMat.at<cv::Vec3b>(i,j+1)[0];dstMat.at<cv::Vec3b>(i,j)[1] = srcMat.at<cv::Vec3b>(i,j+1)[1];dstMat.at<cv::Vec3b>(i,j)[2] = srcMat.at<cv::Vec3b>(i,j+1)[2];}{beDeletedLine.at<cv::Vec3b>(i,0)[0] = srcMat.at<cv::Vec3b>(i,k)[0];beDeletedLine.at<cv::Vec3b>(i,0)[1] = srcMat.at<cv::Vec3b>(i,k)[1];beDeletedLine.at<cv::Vec3b>(i,0)[2] = srcMat.at<cv::Vec3b>(i,k)[2];}}}void run(cv::Mat& image,cv::Mat& outImage,cv::Mat& outMinTrace,cv::Mat& outDeletedLine){cv::Mat image_gray(image.rows,image.cols,CV_8U,cv::Scalar(0));cv::cvtColor(image,image_gray,CV_BGR2GRAY); //彩色图像转换为灰度图像cv::Mat gradiant_H(image.rows,image.cols,CV_32F,cv::Scalar(0));//水平梯度矩阵cv::Mat gradiant_V(image.rows,image.cols,CV_32F,cv::Scalar(0));//垂直梯度矩阵cv::Mat kernel_H = (cv::Mat_<float>(3,3) << 0, 0, 0, 0, 1, -1, 0, 0, 0); //求水平梯度所使用的卷积核（赋初始值）cv::Mat kernel_V = (cv::Mat_<float>(3,3) << 0, 0, 0, 0, 1, 0, 0, -1, 0); //求垂直梯度所使用的卷积核（赋初始值）cv::filter2D(image_gray,gradiant_H,gradiant_H.depth(),kernel_H);cv::filter2D(image_gray,gradiant_V,gradiant_V.depth(),kernel_V);cv::Mat gradMag_mat(image.rows,image.rows,CV_32F,cv::Scalar(0));cv::add(cv::abs(gradiant_H),cv::abs(gradiant_V),gradMag_mat);//水平与垂直滤波结果的绝对值相加，可以得到近似梯度大小////如果要显示梯度大小这个图，因为gradMag_mat深度是CV_32F，所以需要先转换为CV_8U//cv::Mat testMat;//gradMag_mat.convertTo(testMat,CV_8U,1,0);//cv::imshow("Image Show Window2",testMat);//计算能量线cv::Mat energyMat(image.rows,image.cols,CV_32F,cv::Scalar(0));//累计能量矩阵cv::Mat traceMat(image.rows,image.cols,CV_32F,cv::Scalar(0));//能量最小轨迹矩阵calculateEnergy(gradMag_mat,energyMat,traceMat); //找出最小能量线cv::Mat minTrace(image.rows,1,CV_32F,cv::Scalar(0));//能量最小轨迹矩阵中的最小的一条的轨迹getMinEnergyTrace(energyMat,traceMat,minTrace);//显示最小能量线cv::Mat tmpImage(image.rows,image.cols,image.type());image.copyTo(tmpImage);for (int i = 0;i < image.rows;i++){int k = minTrace.at<float>(i,0);tmpImage.at<cv::Vec3b>(i,k)[0] = 0;tmpImage.at<cv::Vec3b>(i,k)[1] = 0;tmpImage.at<cv::Vec3b>(i,k)[2] = 255;}cv::imshow("Image Show Window (A)",tmpImage);//删除一列cv::Mat image2(image.rows,image.cols-1,image.type());cv::Mat beDeletedLine(image.rows,1,CV_8UC3);//记录被删掉的那一列的值delOneCol(image,image2,minTrace,beDeletedLine);cv::imshow("Image Show Window",image2);image2.copyTo(outImage);minTrace.copyTo(outMinTrace);beDeletedLine.copyTo(outDeletedLine);}void recoverOneLine(cv::Mat& inImage,cv::Mat&inTrace,cv::Mat& inDeletedLine,cv::Mat& outImage){cv::Mat recorvedImage(inImage.rows,inImage.cols+1,CV_8UC3);for (int i = 0; i < inImage.rows; i++){int k = inTrace.at<float>(i);for (int j = 0; j < k; j++){recorvedImage.at<cv::Vec3b>(i,j)[0] = inImage.at<cv::Vec3b>(i,j)[0];recorvedImage.at<cv::Vec3b>(i,j)[1] = inImage.at<cv::Vec3b>(i,j)[1];recorvedImage.at<cv::Vec3b>(i,j)[2] = inImage.at<cv::Vec3b>(i,j)[2];}recorvedImage.at<cv::Vec3b>(i,k)[0] = inDeletedLine.at<cv::Vec3b>(i,0)[0];recorvedImage.at<cv::Vec3b>(i,k)[1] = inDeletedLine.at<cv::Vec3b>(i,0)[1];recorvedImage.at<cv::Vec3b>(i,k)[2] = inDeletedLine.at<cv::Vec3b>(i,0)[2];for (int j = k + 1;j < inImage.cols + 1; j++){recorvedImage.at<cv::Vec3b>(i,j)[0] = inImage.at<cv::Vec3b>(i,j-1)[0];recorvedImage.at<cv::Vec3b>(i,j)[1] = inImage.at<cv::Vec3b>(i,j-1)[1];recorvedImage.at<cv::Vec3b>(i,j)[2] = inImage.at<cv::Vec3b>(i,j-1)[2];}}//显示恢复的轨迹cv::Mat tmpImage(recorvedImage.rows,recorvedImage.cols,recorvedImage.type());recorvedImage.copyTo(tmpImage);for (int i = 0;i < tmpImage.rows;i++){int k = inTrace.at<float>(i,0);tmpImage.at<cv::Vec3b>(i,k)[0] = 0;tmpImage.at<cv::Vec3b>(i,k)[1] = 255;tmpImage.at<cv::Vec3b>(i,k)[2] = 0;}cv::imshow("Image Show Window (B)",tmpImage);recorvedImage.copyTo(outImage);}int main(int argc,char* argv){cv::Mat image = cv::imread("1.jpg");cv::namedWindow("Original Image");cv::imshow("Original Image",image);cv::Mat tmpMat;image.copyTo(tmpMat);cv::Mat traces[NUM];cv::Mat deletedLines[NUM];cv::Mat outImage;cv::waitKey(2000);for (int i = 0;i < NUM;i++){run(tmpMat,outImage,traces[i],deletedLines[i]);tmpMat = outImage;cv::waitKey(50);}cv::Mat tmpMat2;outImage.copyTo(tmpMat2);for (int i = 0; i < NUM; i++){recoverOneLine(tmpMat2,traces[NUM-i-1],deletedLines[NUM-i-1],outImage);tmpMat2 = outImage;cv::waitKey(50);}cv::waitKey(115000);return 0;}