matlab练习程序（图像扭曲算法）

clear all;close all;clc;

img=imread('lena.jpg');
imshow(img)

[h,w]=size(img);

cellsize = 32;                              %网格大小，不同的大小会产生不同尺度的噪声
G = rand(2,h/cellsize+1,w/cellsize+1)-0.5;  %每个网格顶点的随机方向向量
G = G*8;                                    %方向向量乘上距离，扭曲大小

pc = zeros(h*w,3);                          %扭曲后点集
num=1;
for y=1:cellsize:h-cellsize+1
    for x=1:cellsize:w-cellsize+1
        
        src = [x y;x y+cellsize-1; x+cellsize-1 y+cellsize-1;x+cellsize-1 y];           %扭曲前一个格子的四个坐标点
        
        indx = floor(x / cellsize) + 1;
        indy = floor(y / cellsize) + 1;
        
        dst = src + [G(:,indx,indy)';G(:,indx,indy+1)'; G(:,indx+1,indy+1)';G(:,indx+1,indy)']; %扭曲后一个格子的四个坐标点
        
        X = zeros(8,8);
        Y = zeros(8,1);                        
        X(1:4,1) = src(:,1);                        %计算投影矩阵
        X(1:4,2) = src(:,2);
        X(1:4,3) = 1;
        
        X(1:4,7) = -src(:,1).*dst(:,1);
        X(1:4,8) = -src(:,2).*dst(:,1);
        
        X(5:8,4) = src(:,1);
        X(5:8,5) = src(:,2);
        X(5:8,6) = 1;
        
        X(5:8,7) = -src(:,1).*dst(:,2);
        X(5:8,8) = -src(:,2).*dst(:,2);
        
        Y(1:4) = dst(:,1);
        Y(5:8) = dst(:,2);   
        A = inv(X)*Y;
        
        %得到扭曲后待采样点集
        for yy=y:y+cellsize-1
            for xx=x:x+cellsize-1          
                x1 = (A(1)*xx+A(2)*yy+A(3)) / (A(7)*xx+A(8)*yy+1);
                y1 = (A(4)*xx+A(5)*yy+A(6)) / (A(7)*xx+A(8)*yy+1);     
                pc(num,:) = [x1 y1 double(img(xx,yy))];
                num = num+1;
            end
        end     
    end
end

%点集投到一个较大的栅格(3*3)中，提升后续栅格化速度
maxx = max(pc(:,1));
minx = min(pc(:,1));
maxy = max(pc(:,2));
miny = min(pc(:,2));
gridx = floor((maxx - minx)/3 )+ 1;
gridy = floor((maxy - miny)/3 )+ 1;
grid = cell(gridy,gridx);                  
for i=1:length(pc)
    indx = floor((pc(i,1) - minx)/3) + 1;
    indy = floor((pc(i,2) - miny)/3) + 1;    
    grid{indy,indx} = [grid{indy,indx};pc(i,:)];
end

%对点集按(1*1)重新栅格化得到扭曲后图像
imgre = zeros(h,w);
for y=1:h
    for x=1:w
        indx = floor((x - minx)/3) + 1;
        indy = floor((y - miny)/3) + 1; 

        pctmp = [];        %在(3*3)*9格子中找最邻近点
        for dy = indy-1:indy+1
            for dx = indx-1:indx+1
                if dx>=1 && dy>=1 && dx<=gridx && dy<=gridy
                    pctmp = [pctmp;grid{dy,dx}];
                end
            end
        end
            
        if isempty(pctmp)==false
            d = [x y] - pctmp(:,1:2);
            [~,ind] = min(sqrt(d(:,1).^2+d(:,2).^2));
            imgre(y,x) = pctmp(ind,3);
        end
    end
end

figure;
plot(pc(:,1),pc(:,2),'.')
axis equal;

figure;
imshow(imgre',[])
imwrite(uint8(imgre'),'imgre.jpg');

%投到点云上看看
figure;
I=uint8(pc(:,3));   
x=pc(:,1);
y=pc(:,2);
z=zeros(h*w,1);

pcshow([x y z],[I I I]);    %如果是彩色图：pcshow([x y z],I);