遥感图像几何校正 matlab,算法仿真《基于畸变率的图像几何校正》matlab7.0

% 试验结果参考 相册中的图片

%只是用于镜头视角

%function DistortionCorrect();

clear all;

clc;

tic

pi=3.141592653589793;

I0=imread('pic.bmp');

I1=rgb2gray(I0);

[height0,width0]=size(I1);

height=2*height0;

width =2*width0;

I2=zeros(height,width);

a1=height/4;

a2=height*3/4;

b1=width/4;

b2=width*3/4;

for i=1:height

for j=1:width

if i<=a1 || i>a2 || j<=b1 || j>b2;

I2(i,j)=255;

else

I2(i,j)=I1(i-a1,j-b1);

end

end

end

Inew=DistortionRate(I2);

% I2=uint8(I2);

Inew=uint8(Inew);

%figure,imshow(I2);

figure,imshow(Inew);

toc

time=toc

%%====================================================

function In=DistortionRate(G);

pi=3.141592653589793;

[height,width]=size(G);

In=zeros(height,width);

%===============

% r_m=sqrt((height^2+width^2)/4);

r_m =sqrt((height/2)^2 + (width/2)^2);

theta_m = 6.63*pi/18.0;%(66.3/90)*(pi/2);

h = r_m/tan(theta_m);

Xc=height/2;

Yc=width/2;

% Xc=488+480;

% Yc=670+640;

for i=1:height

for j=1:width

r=sqrt(((i-Xc)^2+(j-Yc)^2)/1.0);

theta=atan(r/h);

a=tan(theta);

if abs(a)>0.00000001;

D=(theta-a)/a;

x1=ceil((1+D)*(i-Xc));

y1=ceil((1+D)*(j-Yc));

x2=floor((1+D)*(i-Xc));

y2=floor((1+D)*(j-Yc));

if x1+Xc

In(i,j)=(G(x1+Xc,y1+Yc) + G(x2+Xc,y1+Yc) + G(x1+Xc,y2+Yc) + G(x2+Xc,y2+Yc)) / 4;

end

else

x=i;

y=j;

if x+Xc

In(i,j)=G(x+Xc,y+Yc);

end

end

%         if x+Xc

%             In(i,j)=G(x+Xc,y+Yc);

%        end

end

end

%%==================================================