matlab練習程序（圖像扭曲算法）

方法就是先對圖像按照cellsize設置網格，一般是16*16或32*32。

然后對每個網格做投影變換，最后把所有格子拼起來就行了。

單個格子類似下圖：

在實際編程的時候這里沒有采用常見的反變換采樣法，而是采用了正向變換的方式直接處理，投影公式見這里。

正向變換后得到待采樣點集，再對點集重新進行一次柵格化，就能得到沒有空洞的圖像了。

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);

原圖：

扭曲結果：

扭曲后點集（就是在該點集上柵格化）：

圖像直接投到網格點雲上：

參考：http://www.360doc.com/content/06/0823/16/3500_188533.shtml