![4fef8a9da0ece80107a0813cafe65f3b.png](https://i-blog.csdnimg.cn/blog_migrate/a10aca86db6ea754f8dc1782134ac1cd.jpeg)
close all;
clear;
clc;
I=imread('sx1.jpg');
I=rgb2gray(I);
imshow(I);
I=double(I);
[m,n]=size(I);
I_Ex=zeros(m,n);
C=zeros(m,1);
R=floor(m/2-1);
% R=299;
k=40/151; % 27/151;
a=0:(1/R):1;
b=a.^3.5 *k; %2
c=1.-b;
for i=1:R
C(R+1-i)=c(i)*sqrt(R^2-(i)^2)/R;
C(R+1+i)=c(i)*sqrt(R^2-(i)^2)/R;
% C(R+1-i)=sqrt(R^2-(i)^2)/R
% C(R+1+i)=sqrt(R^2-(i)^2)/R;
end
C(R+1)=1;
yc=floor(m/2);
xc=floor(n/2); % 231;
count=0;
for i=2:m-10
for j=2:n-10
%%%经纬方法,存在失真
xp= (j-xc)*C(i)+xc;
yp= i;
xs=floor(xp);
ys=floor(yp);
a0=xs+1-xp;
b0=xp-xs;
a=ys+1-yp;
b=yp-ys;
p1=I(ys,xs);
% p2=I(ys,xs+1);
% p3=I(ys+1,xs);
% p4=I(ys+1,xs+1);
%
% p5=p1*a0+p2*b0;
% p6=p3*a0+p4*b0;
I_Ex(i,j)=p1;%p5*a+p6*b;
if p1>100
count=count+1;
end
end
end
I_Ex=uint8(I_Ex);
figure;
imshow(I_Ex);
% imwrite(I_Ex,'D:b.bmp')
count
m*n