下载文档原格式
实验二图像几何变换与插值一、实验目的巩固图像处理编程的步骤格式，理解数据插值及图像几何变换的原理，掌握图像几何变换的实现方法。
二、实验内容1、理解反向变换的实现思路2、图像缩放及插值Matlab取整命令：floor, ceil, round分别实验最近邻插值和双线性插值f=imread('lena.bmp');beishu=0.5;[row,col]=size(f);r1=round(row*beishu);c1=round(col*beishu);b=zeros(r1,c1);for i=1:r1for j=1:c1i1=round(i/beishu);j1=round(j/beishu);if i1<1i1=1;endif j1<1j1=1;endb(i,j)=f(i1,j1);endendb=uint8(b);figure;imshow(f);figureimshow(b);3、图像旋转及插值以图像中心为轴实现任意角度（逆时针为正）的图像旋转，分别实验两种插值算法f=imread('lena.bmp');B=imrotate(f,45,'nearest','crop');C=imrotate(f,45,'bilinear','crop');figure;subplot(121);imshow(f);title('原图像');subplot(122);imshow(B);title('最近邻插值');figure;subplot(121);imshow(f);title('原图像');subplot(122);imshow(C);title('双线性插值');原图像最近邻插值处理原图像双线性插值处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matlab练习程序（图像放大/缩小，最邻近插值）ccccl;w=0.6; %放大或缩小的宽度h=1.4; %放大或缩小的高度img=imread('Corner.png');imshow(img);[m n]=size(img);imgn=zeros(h*m,w*n);rot=[h 0 0;0 w 0;0 0 1]; %变换矩阵x=h*u,y=w*vinv_rot=inv(rot);for x=1:h*mfor y=1:w*npix=[x y 1]*inv_rot;imgn(x,y)=img(round(pix(1)),round(pix(2)));endendfigure,imshow(uint8(imgn)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matlab练习程序（图像放大/缩小，双线性插值）cccm=1.8; %放大或缩小的高度n=2.3; %放大或缩小的宽度img=imread('lena.jpg');imshow(img);[h w]=size(img);imgn=zeros(h*m,w*n);rot=[m 0 0;0 n 0;0 0 1]; %变换矩阵for i=1:h*mfor j=1:w*npix=[i j 1]/rot;float_Y=pix(1)-floor(pix(1));float_X=pix(2)-floor(pix(2));if pix(1) < 1%边界处理pix(1) = 1;endif pix(1) > hpix(1) = h;endif pix(2) < 1pix(2) =1;endif pix(2) > wpix(2) =w;endpix_up_left=[floor(pix(1)) floor(pix(2))];%四个相邻的点pix_up_right=[floor(pix(1)) ceil(pix(2))];pix_down_left=[ceil(pix(1)) floor(pix(2))];pix_down_right=[ceil(pix(1)) ceil(pix(2))];value_up_left=(1-float_X)*(1-float_Y);%计算临近四个点的权重value_up_right=float_X*(1-float_Y);value_down_left=(1-float_X)*float_Y;value_down_right=float_X*float_Y;%按权重进行双线性插值imgn(i,j)=value_up_left*img(pix_up_left(1),pix_up_left(2))+ ...value_up_right*img(pix_up_right(1),pix_up_right(2))+ ...value_down_left*img(pix_down_left(1),pix_down_left(2))+ ...value_down_right*img(pix_down_right(1),pix_down_right(2));endendfigure,imshow(uint8(imgn)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matlab练习程序（图像旋转，最邻近插值）ccccl;H=1; %索引pix中第一个元素，即高度W=2; %索引pix中第二个元素，即宽度jiaodu=45; %要旋转的角度，旋转方向为顺时针img=imread('Corner.png'); %这里v为原图像的高度，u为原图像的宽度imshow(img); %这里y为变换后图像的高度，x为变换后图像的宽度[v u]=size(img);theta=jiaodu/180*pi;rot=[cos(theta) -sin(theta) 0;sin(theta) cos(theta) 0;0 0 1];inv_rot=inv(rot);pix1=[1 1 1]*rot; %变换后图像左上点的坐标pix2=[1 u 1]*rot; %变换后图像右上点的坐标pix3=[v 1 1]*rot; %变换后图像左下点的坐标pix4=[v u 1]*rot; %变换后图像右下点的坐标height=round(max([abs(pix1(H)-pix4(H))+0.5 abs(pix2(H)-pix3(H))+0.5])); %变换后图像的高度width=round(max([abs(pix1(W)-pix4(W))+0.5 abs(pix2(W)-pix3(W))+0.5])); %变换后图像的宽度imgn=zeros(height,width);delta_y=abs(min([pix1(H) pix2(H) pix3(H) pix4(H)])); %取得y方向的负轴超出的偏移量delta_x=abs(min([pix1(W) pix2(W) pix3(W) pix4(W)])); %取得x方向的负轴超出的偏移量for y=1-delta_y:height-delta_yfor x=1-delta_x:width-delta_xpix=[y x 1]*inv_rot; %用变换后图像的点的坐标去寻找原图像点的坐标，%否则有些变换后的图像的像素点无法完全填充if pix(H)>=0.5 && pix(W)>=0.5 && pix(H)<=v && pix(W)<=uimgn(y+delta_y,x+delta_x)=img(round(pix(H)),round(pix(W)));endendendfigure,imshow(uint8(imgn))%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% matlab练习程序（图像旋转，双线性插值）cccjiaodu=45; %要旋转的角度，旋转方向为顺时针img=imread('lena.jpg'); %这里v为原图像的高度，u为原图像的宽度imshow(img); %这里y为变换后图像的高度，x为变换后图像的宽度[h w]=size(img);theta=jiaodu/180*pi;rot=[cos(theta) -sin(theta) 0;sin(theta) cos(theta) 0;0 0 1];pix1=[1 1 1]*rot; %变换后图像左上点的坐标pix2=[1 w 1]*rot; %变换后图像右上点的坐标pix3=[h 1 1]*rot; %变换后图像左下点的坐标pix4=[h w 1]*rot; %变换后图像右下点的坐标height=round(max([abs(pix1(1)-pix4(1))+0.5 abs(pix2(1)-pix3(1))+0.5])); %变换后图像的高度width=round(max([abs(pix1(2)-pix4(2))+0.5 abs(pix2(2)-pix3(2))+0.5])); %变换后图像的宽度imgn=zeros(height,width);delta_y=abs(min([pix1(1) pix2(1) pix3(1) pix4(1)])); %取得y方向的负轴超出的偏移量delta_x=abs(min([pix1(2) pix2(2) pix3(2) pix4(2)])); %取得x方向的负轴超出的偏移量for i=1-delta_y:height-delta_yfor j=1-delta_x:width-delta_xpix=[i j 1]/rot; %用变换后图像的点的坐标去寻找原图像点的坐标，%否则有些变换后的图像的像素点无法完全填充float_Y=pix(1)-floor(pix(1));float_X=pix(2)-floor(pix(2));if pix(1)>=1 && pix(2)>=1 && pix(1) <= h && pix(2) <= wpix_up_left=[floor(pix(1)) floor(pix(2))]; %四个相邻的点pix_up_right=[floor(pix(1)) ceil(pix(2))];pix_down_left=[ceil(pix(1)) floor(pix(2))];pix_down_right=[ceil(pix(1)) ceil(pix(2))];value_up_left=(1-float_X)*(1-float_Y); %计算临近四个点的权重value_up_right=float_X*(1-float_Y);value_down_left=(1-float_X)*float_Y;value_down_right=float_X*float_Y;imgn(i+delta_y,j+delta_x)=value_up_left*img(pix_up_left(1),pix_up_left(2))+ ...value_up_right*img(pix_up_right(1),pix_up_right(2))+ ...value_down_left*img(pix_down_left(1),pix_down_left(2))+ ...value_down_right*img(pix_down_right(1),pix_down_right(2));endendendfigure,imshow(uint8(imgn)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%实验一 Matlab图像基本操作一、实验目的熟悉利用Matlab进行图像处理的基本操作，了解图像数据的存储形式及进行图像处理编程的步骤方法。
Matlab中插值函数汇总和使用说明MATLAB中的插值函数为interp1，其调用格式为： yi= interp1(x,y,xi,'method')其中x，y为插值点，yi为在被插值点xi处的插值结果；x,y为向量， 'method'表示采用的插值方法，MATLAB提供的插值方法有几种： 'method'是最邻近插值， 'linear'线性插值； 'spline'三次样条插值； 'cubic'立方插值．缺省时表示线性插值注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。
例如：在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为12，9，9，10，18 ，24，28，27，25，20，18，15，13，推测中午12点（即13点）时的温度．x=0:2:24;y=[12 9 9 10 18 24 28 27 25 20 18 15 13];a=13;y1=interp1(x,y,a,'spline')结果为： 27.8725若要得到一天24小时的温度曲线，则：xi=0:1/3600:24;yi=interp1(x,y,xi, 'spline');plot(x,y,'o' ,xi,yi)命令1 interp1功能一维数据插值（表格查找）。
该命令对数据点之间计算内插值。
它找出一元函数f(x)在中间点的数值。
其中函数f(x)由所给数据决定。
x：原始数据点Y：原始数据点xi：插值点Yi：插值点格式(1)yi = interp1(x,Y,xi)返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。
参量x 指定数据Y 的点。
若Y 为一矩阵，则按Y 的每列计算。
yi是阶数为length(xi)*size(Y,2)的输出矩阵。
Matlab中插值函数汇总和使用说明以下内容的参考地址：/blog/article.asp?id=28904已知向量x,y,通过x,y求出向量xi的插值的值1、分段线性插值yi=interp1(x,y,xi)或yi=interp1(x,y,'linear')2、多项式插值多项式为y=p(1)*x+p(2)*x^2+...+p(n)*x^n,则可以用p=polyfit(x,y,n)，（n是多项式的最高次数）求出系数向量p,然后用y=polyval(p,xi)求出向量xi的函数值y3、三次样条插值yi=interp1(x,y,xi,'spline')或yi=spline(x,y,xi)或pp=interp1(x,y,'spline','pp')，然后pp=spline(x,y)-->yi=ppval(pp,xi)4、分段三次埃尔米特yi=interp1(x,y,xi,'pchip')或yi=pchip(x,y,xi)或pp=interp1(x,y,'pchip','pp') pp=pchip(x,y)-->yi=pppval(pp.xi)以下内容的参考地址：/s/blog_6056e8800100k3xf.html命令1 interp1功能一维数据插值（表格查找）。
该命令对数据点之间计算内插值。
它找出一元函数f(x)在中间点的数值。
其中函数f(x)由所给数据决定。
x：原始数据点Y：原始数据点xi：插值点Yi：插值点格式(1)yi = interp1(x,Y,xi)返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。
参量x 指定数据Y 的点。
若Y 为一矩阵，则按Y 的每列计算。
yi 是阶数为length(xi)*size(Y,2)的输出矩阵。
(2)yi = interp1(Y,xi)假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。
magnify matlab源程序以下是一个简单的 MATLAB 源程序，用于实现图片的放大功能： % 以下代码实现了图片的放大功能% 输入的变量 im 是一个三维 RGB 图像矩阵，factor 是放大因子 % 输出的变量 im_scaled 是一个放大后的三维 RGB 图像矩阵function im_scaled = magnify(im, factor)% 读取输入的图片大小[height, width, ~] = size(im);% 通过放大因子计算新的图片大小new_height = round(height * factor);new_width = round(width * factor);% 通过双线性插值法计算新的像素值，并把结果存储到输出矩阵中 for i = 1:new_heightfor j = 1:new_widthy = y_ratio * (i-1) + 1;x = x_ratio * (j-1) + 1;y1 = floor(y);y2 = ceil(y);x1 = floor(x);x2 = ceil(x);if x1 == x2x2 = x1 + 1;endif y1 == y2y2 = y1 + 1;endf11 = double(im(y1, x1, :));f12 = double(im(y2, x1, :));f21 = double(im(y1, x2, :));f22 = double(im(y2, x2, :));dx = x - x1;dy = y - y1;im_scaled(i,j,:) = uint8((1-dx)*(1-dy)*f11 + dx*(1-dy)*f21 +(1-dx)*dy*f12 + dx*dy*f22);endend% 显示原始图片和放大后的图片subplot(1, 2, 1);imshow(im);title('原始图片');subplot(1, 2, 2);imshow(im_scaled);title('放大后的图片');end该程序的功能是将输入的 RGB 图片矩阵进行放大，输出结果是放大后的图片矩阵。
% Define a function for shrinking and zooming images by bilinear interpolation% I is the original image and N is the factor for shrinking firstly and zooming later function interpolation_shrink_zoom_image(I, N% First to shrink the original image% To measure the size of original image I[P,Q]=size(I;% Rebuild the original image I to I_TEMP for bilinear interpolationI_TEMP=zeros(P+2,Q+2; I_TEMP(2:P+1,2:Q+1=I;I_TEMP(1,2:Q+1=I(1,:; I_TEMP(P+2,2:Q+1=I(P,:;I_TEMP(2:P+1,1=I(:,1; I_TEMP(2:P+1,Q+2=I(:,Q;I_TEMP(1,1=I(1,1; I_TEMP(1,P+2=I(1,P;I_TEMP(P+2,1=I(P,1; I_TEMP(P+2,Q+2=I(P,Q;I_TEMP=double(I_TEMP;% To fix the size of shrunken image% R is the number of rows and C is the number of columnsR=floor(P/N; C=floor(Q/N;% Define the shrunken imageI_S=zeros(R,C;% For bilinear interpolation f(x,y in 4 points f(x1,y1,f(x2,y1,f(x1,y2,f(x2,y2 % f(x,y = f(x1,y1*(x2-x*(y2-y/((x2-x1*(y2-y1% + f(x2,y1*(x-x1*(y2-y/((x2-x1*(y2-y1% + f(x1,y2*(x2-x*(y-y1/((x2-x1*(y2-y1% + f(x2,y2*(x-x1*(y-y1/((x2-x1*(y2-y1% For pixel matrix, x2=x1+1, y2=y1+1, x=x1+a, y=y1+b% Let x1=u, y1=v, then x2=u+1, y2=v+1, x=u+a, y=v+b% Then f(u+a, v+b = f(u,v*(1-a*(1-b + f(u+1,v*a*(1-b% + f(u,v+1*(1-a*b + f(u+1,v+1*a*b% Method: from original image I get suitable gray value by surrounding 4 pixels,% then fill the gray value into certain location in shrunken image I_S.for i=1:R;for j=1:C;% Step1: get x,y,u,v,a,b. The size of I to size of I_S ratio = N,% (x,y is pixel "in" I and (i,j is pixel in I_S, then x/i=y/j=N.x=i*N; y=j*N; u=floor(x; v=floor(y; a=x-u; b=y-v;% Step2: considering the situation on bound of I, using I_TEMP to get suitable value of f(x,y by bilinear interpolationI_S(i,j = I_TEMP(u+1,v+1*(1-a*(1-b + I_TEMP(u+2,v+1*a*(1-b +I_TEMP(u+1,v+2*(1-a*b + I_TEMP(u+2,v+2*a*b;end ;end ;% Display the original image I and the shrunken image I_SI_S=uint8(I_S;figure; imshow(I,[]; title(['Original image: ',num2str(P,'*',num2str(Q]; figure; imshow(I_S,[]; title(['Interpolation shrinking image:' ,num2str(R,'*',num2str(C];% Second to zoom the shrunken image% R*C is the size of shrunken image as shown before% Rebuild the shrunken image I_S to I_TEMP for bilinear interpolation I_TEMP=zeros(R+2,C+2; I_TEMP(2:R+1,2:C+1=I_S;I_TEMP(1,2:C+1=I_S(1,:; I_TEMP(R+2,2:C+1=I_S(R,:;I_TEMP(2:R+1,1=I_S(:,1; I_TEMP(2:R+1,C+2=I_S(:,C;I_TEMP(1,1=I_S(1,1; I_TEMP(1,R+2=I_S(1,R;I_TEMP(R+2,1=I_S(R,1; I_TEMP(R+2,R+2=I_S(R,C;I_TEMP=double(I_TEMP;% Define the zoomed imageI_Z=zeros(R*N,C*N;% Zooming the shrunken image by bilinear interpolationfor i=1:R*N;for j=1:C*N;x=i/N; y=j/N; u=floor(x; v=floor(y; a=x-u; b=y-v;% Step2: considering the situation on bound of I, using I_TEMP to get suitable value of f(x,y by bilinear interpolationI_Z(i,j = I_TEMP(u+1,v+1*(1-a*(1-b + I_TEMP(u+2,v+1*a*(1-b +I_TEMP(u+1,v+2*(1-a*b + I_TEMP(u+2,v+2*a*b;end ;end ;% Display the original image I and the shrunken image I_SI_Z=uint8(I_Z;figure; imshow(I,[]; title(['Original image: ',num2str(P,'*',num2str(Q]; figure; imshow(I_Z,[]; title(['Interpolation zooming image:' ,num2str(R*N,'*',num2str(C*N];end。
插值方法晚上做一个曲线拟合，结果才开始用最小二乘法拟合时，拟合出来的东西太难看了！于是尝试用其他方法。
经过一番按图索骥，终于发现做曲线拟合的话，采用插值法是比较理想的方法。
尤其是样条插值，插完后线条十分光滑。
方法付后，最关键的问题是求解时要积分，放这里想要的时候就可以直接过来拿，不用死去搜索啦。
呵呵插值方法的Matlab实现一维数据插值MATLAB中用函数interp1来拟合一维数据，语法是YI = INTERP1(X,Y,XI,方法)其中（X，Y）是已给的数据点，XI 是插值点，其中方法主要有'linear' -线性插值，默认'pchip' -逐段三次Hermite插值'spline' -逐段三次样条函数插值其中最后一种插值的曲线比较平滑例：x=0:.12:1; x1=0:.02:1;%（其中x=0:.12:1表示显示的插值点，x1=0:.02:1表示插值的步长）y=(x.^2-3*x+5).*exp(-5*x).*sin(x);plot(x,y,'o'); hold on;y1=interp1(x,y,x1,'spline');plot(x1,y1,':')如果要根据样本点求函数的定积分，而函数又是比较光滑的，则可以用样条函数进行插值后再积分，在MATLAB中可以编写如下程序：function y=quadspln(x0,y0,a,b)f=inline('interp1(x0,y0,x,''spline'')','x','x0','y0');y=quadl(f,a,b,1e-8,[],x0,y0);现求sin(x)在区间[0，pi]上的定积分，只取5点x0=[0,0.4,1,2,pi];y0=sin(x0);I=quadspln(x0,y0,0,pi)结果得到的值为2.01905，精确值为2求一段matlab插值程序悬赏分：20 - 解决时间：2009-12-26 19:57已知5个数据点：x=[0.25 0.5 0.75 1] y=[0 0.3104 0.6177 0.7886 1] ，求一段matlab插值程序，求过这5个数据点的插值多项式，并在x-y坐标中画出y=f（x）图形，并且求出f （x）与x轴围成图形的面积（积分），不胜感激！使用Lagrange 插值多项式的方法：首先把下面的代码复制到M文件中，保存成lagranfunction [C,L]=lagran(X,Y)% input - X is a vector that contains a list of abscissas% - Y is a vector that contains a list of ordinates% output - C is a matrix that contains the coefficients of the lagrange interpolatory polynomial%- L is a matrix that contains the lagrange coefficients polynomialw=length(X);n=w-1;L=zeros(w,w);for k=1:n+1V=1;for j=1:n+1if k~=jV=conv(V,poly(X(j)))/(X(k)-X(j));endendL(k,:)=V;endC=Y*L;然后在命令窗口中输入以下内容：x=[0 0.25 0.5 0.75 1];y=[0 0.3104 0.6177 0.7886 1];lagran(x,y)ans =3.3088 -6.3851 3.3164 0.7599 0得到的数据就是多项式各项的系数，注意最后一个是常数项，即x^0，所以表达式为:f=3.3088*x.^4-6.3851*x.^3+3.3164*x.^2 +0.7599*x求面积就是积分求解>> f=@(x)3.3088*x.^4-6.3851*x.^3+3.3164*x.^2 +0.7599*x;>> quad(f,0,1)ans =0.5509这些点肯定是通过这个多项式的！MATLAB插值与拟合§1曲线拟合实例：温度曲线问题气象部门观测到一天某些时刻的温度变化数据为：试描绘出温度变化曲线。
使⽤双线性插值法放⼤图像（matlab实现）双线性插值的概念及公式可以参考百度，这⾥仅对算法原理进⾏简单的说明：双线性插值计算公式：f(i+u,j+v) = (1-u)(1-v)f(i,j)+u(1-v)f(i+1,j)+(1-u)vf(i,j+1)+ uvf(i+1,j+1)这个公式表明了如何利⽤矩阵中的四个像素值计算新的像素值，这些新的像素值就组成了放⼤后的图像。
下图是如何将3x3的图像放⼤为4x4的图像：原图像表⽰为3x3的矩阵（像素值处在⿊线的交叉点上），如何计算4x4矩阵的值呢？（像素值处在红⾊虚线交叉点及红线与⿊线的交点上）⽐如新图像B的第⼀列与原图像A的第⼀列的对应关系是：B(1,1) = A(1,1)B(1,2) = A(1,1.66667)B(1,3) = A(1,2.33334)B(1,4) = A(1,3.00001)应⽤上⾯的公式，实际上就是⽤A的含有⼩数点的位置的像素值来计算B的像素值，但含有⼩数点位的像素是不存在的，这⾥称为虚拟位置。
⽤原图像A的值就能计算出放⼤后B的值,是不是很神奇？实际上可以这样认为：双线性插值就是把放⼤后的图像再压缩到原来图像的尺⼨⼤⼩，计算原图像中虚拟的像素值，等同于计算放⼤后图像的像素值，对于本例来说，B图像的步长相当于A图像步长的(3-1)/(4-1)=0.66667倍。
下⾯我们就可以利⽤这个⽐率来对应B中像素位置与A中虚拟像素位置的关系。
B(1,1) = A(1,1) (1-1)*0.66667+1=1B(1,2) = A(1,1.66667) (2-1)*0.66667+1=1.66667B(1,3) = A(1,2.33334) (3-1)*0.66667+1=2.33334B(1,4) = A(1,3.00001) (4-1)*0.66667+1=3.00001根据上⾯的对应关系，我们就可以⽤代码实现了。
现在还有⼀个问题：我们计算虚拟像素值是需要周围四个原像素值，⽐如上列中的（下图中红圈圈住的部分）A(1,3) = (1-0)(1-0)A(1,3) + (1-0)0A(1,4) + 0(1-0)A(2,3) + 00A(2,4)显然这⾥的A(1,4)和A(2,4)是⽆法索引到得，因为原图像是3x3的矩阵。
