phase-mag-nonilm.m

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

imgFile = 'sc.bmp';
magFile = 'de.bmp';


name0 = imgFile(1:(length(imgFile)-4));
imgMat = double(imread(imgFile));
magMat = double(imread(magFile));
dimen1 = double(size(imgMat, 1));
dimen2 = double(size(imgMat, 2));
dimen3 = double(size(imgMat, 3));
imgMat = imgMat(:,:,1);

scale1 = 250;
scale2 = 200;
maxInt = 255/scale2;
magMat = magMat / scale2 * maxInt;

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

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


py = exp( -i * k1.^2 / (T1 * y2));
px = exp( -i * k2.^2 / (T2 * x2));
prgt = py.'* px;
iprgt = conj(prgt);
T3 = 40;
ilmy = exp((  - 1)* (k1 - y2/2).^2 / (T3 * y2));
ilmx = exp((  - 1)* (k2 - x2/2).^2 / (T3 * x2));
ilm  = ilmy.' * ilmx;
ilm  = ilm .* (abs(ilm) > 0.3);

ilmup = ilm(1:M, 1:x2);
ilmlt = ilm((M+1):y1, 1:N);
ilmrt = ilm((M+1):y1, (x1+1):x2);
ilmdn = ilm((y1+1):y2, 1:x2);
ilmmd = 1 ./ (ilm((M+1):y1, (N+1):x1) + 0.00000001);
iilm = 1 ./ (ilm + 0.00000001);

imgMatSup = [imgMat; zeros(M, dimen2)];
imgMatSup = [zeros(M, dimen2); imgMatSup];
imgMatSup = [zeros(y2, N), imgMatSup];
imgMatSup = [imgMatSup, zeros(y2, N)];

support = double((imgMatSup == 0));

magMatSup = [magMat; maxInt * ones(M, dimen2)];
magMatSup = [maxInt * ones(M, dimen2); magMatSup];
magMatSup = [maxInt * ones(y2, N), magMatSup];
magMatSup = [magMatSup, maxInt * ones(y2, N)];

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

spec = ifft2(fft2(imgMatSup) .* prgt);
fm = abs(spec);
tot = sum(sum(fm .^ 2));
midMat = zeros(y2, x2);
reform = exp(i * pi/2);

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

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

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

beta = 0.9;
steps = 50;
savestep = 1;
R = 1:steps;
midSp = ifft2(fft2(midMat) .* prgt);

for k = 1:steps
    for l = 1:savestep

        midSp = fm .* midSp ./ (abs(midSp) + 0.00000001);
        imgMatSup = ifft2(fft2(midSp) .* iprgt);

        midMat(1:M, 1:x2) = midMat(1:M, 1:x2) - beta * imgMatSup(1:M, 1:x2) + maxInt * ilmup * beta;
        midMat((M+1):y1, 1:N) = midMat((M+1):y1, 1:N) - beta * imgMatSup((M+1):y1, 1:N) + maxInt * ilmlt * beta;
        midMat((M+1):y1, (x1+1):x2) = midMat((M+1):y1, (x1+1):x2) - beta * imgMatSup((M+1):y1, (x1+1):x2) + maxInt * ilmrt * beta;
        midMat((y1+1):y2, 1:x2) = midMat((y1+1):y2, 1:x2) - beta * imgMatSup((y1+1):y2, 1:x2) + maxInt * ilmdn * beta;
        midMat((M+1):y1, (N+1):x1) =  imgMatSup((M+1):y1, (N+1):x1);

        midSp = ifft2(fft2(midMat) .* prgt);
        %imshow(uint8(imag(log(0.0000001 + midMat .* iilm)) * scale1 / pi));
        imshow(uint8(abs(midMat) * scale2 / maxInt))
        %imshow(uint8(imag(log(0.0000001 + midMat((M+1):y1, (N+1):x1) .* ilmmd)) * scale1 / pi));
%         imshow(uint8(abs(midMat((M+1):y1, (N+1):x1)) .* ilmmd * scale2 / maxInt));
    end
    R(k) = sum(sum(abs(abs(midSp) - fm)));
%     name = strcat('.\\sp-dm-b08-s26\\p-', num2str(k));
%     name = strcat(name, '.jpg');
%     name1 = strcat('.\\sp-dm-b08-s26\\m-', 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


midSp1 = midSp;
h = 0.001;
count = 0;
for l = 1:5
    midSp1 = fm .* midSp1 ./ (abs(midSp1) + 0.00000001);
    midMat1 = ifft2(fft2(midSp1) .* iprgt);
    dfE = ifft2((midMat1 - ilm) .* support);    
    drSD = -2 * imag(fft2(dfE .* iprgt) .* midSp1);
    dMidMat = i * h * ifft2(iprgt .* fft2(midSp1 .* drSD));    
    count
    R(l) = sum(sum((abs(midSp1) - fm) .^ 2));
    count = 0;
    %phai = imag(log(midSp1);
    while (1)        
        midMat1 = midMat1 - dMidMat;
        sa(1) = real(sum(sum((midMat1 - ilm) .* conj((midMat1 - ilm)) .* support)));
        midMat1 = midMat1 - dMidMat;
        sa(2) = real(sum(sum((midMat1 - ilm) .* conj((midMat1 - ilm)) .* support)));
        midMat1 = midMat1 - dMidMat;
        sa(3) = real(sum(sum((midMat1 - ilm) .* conj((midMat1 - ilm)) .* support)));
        %imshow(uint8(abs(midMat1) * scale2 / maxInt));
        imshow(uint8(abs(midMat((M+1):y1, (N+1):x1)) .* ilmmd * scale2 / maxInt));
        
        count = count + 1;
        if ((sa(2) > sa(1)) && (sa(3) > sa(2)))
            midMat1 = midMat1 + 3 * dMidMat;
            midSp1 = ifft2(fft2(midMat1) .* prgt);
            break;
        end         
    end    
end
        

R = R/sum(sum(fm));
save('.\\sp-dm-b08-s26\\R', 'R', '-mat');