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