phase-mag-b-sup.m

相位衬度成像中由振幅得到相位的恢复算法

imgFile = 'schrodinger.bmp';
magFile = 'de Broglie.bmp';
%imgFile = '1.jpg';
maxInt = 1;
name0 = imgFile(1:(length(imgFile)-4));
imgMat = double(imread(imgFile));
magMat = double(imread(magFile));
dimen1 = size(imgMat, 1);
dimen2 = size(imgMat, 2);
dimen3 = size(imgMat, 3);
imgMat = 255 - imgMat(:,:,1);
magMat = magMat / 255 * 1;


O1 = 2.8;
O2 = 2.8;
M = int32(dimen1 * (O1 - 1)/2);
N = int32(dimen2 * (O2 - 1)/2);
y1 = M + dimen1;
y2 = y1 + M;
x1 = N + dimen2;
x2 = x1 + N;
whiteSup = 0;

T1 = 40;
k1 = 0 : (y2 - 1);
T2 = x2/y2 * T1;
k2 = 0 : (x2 - 1);

% omiga1 = sqrt(2 * pi / T1 / M / z / lambda);
% omiga2 = omiga1 * sqrt( T1 * M ) / sqrt(T2 * N);
% res1 = 1 / (omiga1 * M);
% res2 = 1 / (omiga2 * N);
%n = sqrt(z * lambda / res1 / res1);

% py = exp( i * double(k1.^2 / (T1 * y2)));
% px = exp( i * double(k2.^2 / (T2 * x2)));
py = exp( -i * double(k1.^2 / (T1 * y2)));
px = exp( -i * double(k2.^2 / (T2 * x2)));
prgt = py.'* px;
%prgt = fft2(prgt1);
iprgt = conj(prgt);
%iprgt = fft2(iprgt);


imgMatSup = [imgMat; whiteSup * ones(M, dimen2)];
imgMatSup = [whiteSup * ones(M, dimen2); imgMatSup];
imgMatSup = [whiteSup * ones(y2, N), imgMatSup];
imgMatSup = [imgMatSup, whiteSup * ones(y2, N)];

magMatSup = [magMat; zeros(M, dimen2)];
magMatSup = [zeros(M, dimen2); magMatSup];
magMatSup = [zeros(y2, N), magMatSup];
magMatSup = [magMatSup, zeros(y2, N)];

imgMatSup = exp(i * imgMatSup/255 * pi);
imgMatSup = imgMatSup .* magMatSup;

spec = ifft2(fft2(imgMatSup) .* prgt);
fm = abs(spec);

midMat = zeros(y2, x2);
reform = exp(i * pi/2);

initPhSup = [rand(dimen1, dimen2) * 3; zeros(M, dimen2)];
initPhSup = [zeros(M, dimen2); initPhSup];
initPhSup = [zeros(y2, N), initPhSup];
initPhSup = [initPhSup, zeros(y2, N)];
midSp = ifft2(fft2(exp(i * initPhSup )) .* prgt);

midSp = fm .* midSp ./ (abs(midSp) + 0.00001);

imgMatSup = ifft2(fft2(midSp) .* iprgt);
midMat(1:M, 1:x2) = 1;
midMat(M:y1, 1:N) = 1;
midMat(M:y1, x1:x2) = 1;
midMat(y1:y2, 1:x2) = 1;
midMat(M:y1, N:x1) = imgMatSup(M:y1, N:x1); % .* (imgMatSup(M:y1, N:x1) > 0);

beta = 0.8;
steps = 30;
savestep = 30;

for k = 1:steps
    for l = 1:savestep
        
        midSp = ifft2(fft2(midMat) .* prgt);       
        midSp = fm .* midSp ./ (abs(midSp) + 0.00001);
        imgMatSup = ifft2(fft2(midSp) .* iprgt);
%         aver = (sum(sum(imgMatSup(1:M, 1:x2))) + ...
%                 sum(sum(imgMatSup((M+1):y1, (x1+1):x2))) + ...
%                 sum(sum(imgMatSup((M+1):y1, 1:N))) + ...               
%                 sum(sum(imgMatSup((y1+1):y2, 1:x2))) ) /...
%                 (dimen1 * dimen2 * (O1 - 1) * (O2 - 1));
%         imgMatSup = imgMatSup / aver;
        midMat(1:M, 1:x2) = midMat(1:M, 1:x2) - beta * imgMatSup(1:M, 1:x2) ;
        midMat((M+1):y1, 1:N) = midMat((M+1):y1, 1:N) - beta * imgMatSup((M+1):y1, 1:N) ;
        midMat((M+1):y1, (x1+1):x2) = midMat((M+1):y1, (x1+1):x2) - beta * imgMatSup((M+1):y1, (x1+1):x2) ;
        midMat((y1+1):y2, 1:x2) = midMat((y1+1):y2, 1:x2) - beta * imgMatSup((y1+1):y2, 1:x2) ;
%         midMat((M+1):y1, (N+1):x1) = (midMat((M+1):y1, (N+1):x1) - beta * imgMatSup((M+1):y1, (N+1):x1)) .* ...
%                              ((imgMatSup((M+1):y1, (N+1):x1)) < 0) + ...
%                              imgMatSup((M+1):y1, (N+1):x1) .* (imgMatSup((M+1):y1, (N+1):x1) > 0);
%       tmp = abs(imgMatSup((M+1):y1, (N+1):x1));
%       midMat((M+1):y1, (N+1):x1) =  (tmp >  0) .* imgMatSup((M+1):y1, (N+1):x1) + ...
%                                      (tmp == 0) .* imgMatSup((M+1):y1, (N+1):x1) * reform;
%        tmp = (tmp > 0) .* tmp + (tmp == 0); 
%       midMat((M+1):y1, (N+1):x1) =  imgMatSup((M+1):y1, (N+1):x1) ./ (tmp + 0.00001);
       %midMat((M+1):y1, (N+1):x1) =  imgMatSup((M+1):y1, (N+1):x1);
       %imshow(uint8(imag(log(0.00001 + midMat)) * 255 / pi));
       midMat((M+1):y1, (N+1):x1) =  imgMatSup((M+1):y1, (N+1):x1);
       %imshow(uint8(255 - imag(log(0.00001 + midMat((M+1):y1, (N+1):x1))) * 255 / pi));
       imshow(uint8(abs(midMat((M+1):y1, (N+1):x1)) * 255 / maxInt));
    end
    name = strcat(name0, num2str(k));
    name = strcat(name, '.jpg');
    name1 = strcat(name0, 'int');
    name1 = strcat(name1, num2str(k));
    name1 = strcat(name1, '.jpg');
    imwrite(uint8(255 - imag(log(0.00001 + midMat((M+1):y1, (N+1):x1))) * 255 / pi), name,'jpg');
    imwrite(uint8(abs(midMat((M+1):y1, (N+1):x1)) * 255 / maxInt), name1,'jpg');
end