phase-mag-nonilm-1f.m

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

imgFile = 'bp.bmp';
magFile = 'di.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 = 500;
maxInt = 255/scale2;
magMat = magMat / scale2 * maxInt;

O1 = 2.5;
O2 = 2.5;
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;

Zfo = 1;
Zf = Zfo * O1 * O2;
ivo = 1;

k1 = (0 : (y2 - 1));
k2 = (0 : (x2 - 1));



py = exp(i * pi * (k1- y2/2).^2 * Zf * ivo ^2 / (y2 * y2)) .* sinc(Zf * ivo / (y2 * y2) * (k1- y2/2));
py = fft(py) .* exp(i * pi * k1);
px = py;
prgt = py.'* px;
iprgt = 1 ./ prgt;

T3 = 40;
ilmy = exp((  - 1)* (k1 - y2/2).^2 / (T3 * y2));
ilmx = exp((  - 1)* (k2 - y2/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;

%imgMatSup = double(abs(ilm) > 0.3);
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 = 100;
savestep = 1;

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);
%         midMat = midMat * sqrt(sum(sum(abs(midMat).^2)) /tot);
        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 = ifft2(fft2(midMat) .* prgt);
midMat1 = midMat;
phai = imag(log(midSp1 + 0.00000001));

h = 0.01;
count = 0;
for l = 1:20
%     midSp1 = fm .* midSp1 ./ (abs(midSp1) + 0.00000001);

    dfE = ifft2((midMat1 - ilm) .* support);    
    drSD = -2 * imag(fft2(dfE .* iprgt) .* midSp1);
%     dMidMat = i * h * ifft2(iprgt .* fft2(midSp1 .* drSD));    
    count
    R = sum(sum((midMat1 - fm) .* conj(midMat1 - fm) .* support))
    count = 0;
     
    for m = 1:200
        phai = phai - drSD * h;
        midSp1 = fm .* exp(i * phai);
        midMat1 = ifft2(fft2(midSp1) .* iprgt);
        sa(1) = real(sum(sum((midMat1 - ilm) .* conj((midMat1 - ilm)) .* support)));
        
        phai = phai - drSD * h;
        midSp1 = fm .* exp(i * phai);
        midMat1 = ifft2(fft2(midSp1) .* iprgt);
        sa(2) = real(sum(sum((midMat1 - ilm) .* conj((midMat1 - ilm)) .* support)));
        
        phai = phai - drSD * h;
        midSp1 = fm .* exp(i * phai);
        midMat1 = ifft2(fft2(midSp1) .* iprgt);
        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)))
            phai = phai + drSD * 3 * h;
            midSp1 = fm .* exp(i * phai);
            midMat1 = ifft2(fft2(midSp1) .* iprgt);
%             midSp1 = ifft2(fft2(midMat1) .* prgt);
            break;
        end
%         midSp1 = fm .* exp(i * phai);
%         midMat1 = ifft2(fft2(midSp1) .* iprgt);
    end       
end
        

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