hio.m

相位衬度成像中恢复算法经典程序

imgFile1 = 'bp.jpg';
imgFile = 'di.jpg';
scale = 255;
scale2 = 250;
imgMat = double(imread(imgFile)) / scale;
phMat = double(imread(imgFile1));
dimen1 = size(imgMat, 1);
dimen2 = size(imgMat, 2);
dimen3 = size(imgMat, 3);
imgMat = imgMat(:,:,1);
imgMat = imgMat .* exp(i * phMat(:,:,1)/scale2 * pi);


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

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

midMat = zeros(y2, x2);
initPhSup = [rand(dimen1, dimen2) * 100; zeros(M, dimen2)];
initPhSup = [zeros(M, dimen2); initPhSup];
initPhSup = [zeros(y2, N), initPhSup];
initPhSup = [initPhSup, zeros(y2, N)];

% initial estimate
midSp = fft2(initPhSup);
midSp = fm .* midSp ./ (abs(midSp) + 0.00001);
imgMatSup = ifft2(midSp);%, 'symmetric');
midMat(1:M, 1:x2) = 0;
midMat(M:y1, 1:N) = 0;
midMat(M:y1, x1:x2) = 0;
midMat(y1:y2, 1:x2) = 0;

midMat(M:y1, N:x1) = imgMatSup(M:y1, N:x1);

beta = 0.8;
steps = 1;
savestep = 700;
gn = 0;
for k = 1:savestep 
    for l = 1:steps
        imgMatSup = ifft2(midSp);%, 'symmetric');
%         R1(k) = (sum(sum(abs(midMat)))) / (sum(sum(abs(midMat((M+1):y1, (N+1):x1)))) - 1);
%         R3(k) = (sum(sum(abs(imgMatSup)))) / (sum(sum(abs(imgMatSup((M+1):y1, (N+1):x1)))) - 1);
        R1(k) = ((sum(sum(abs(midMat((M+1):y1, (N+1):x1)))) / sum(sum(abs(midMat)))));
        R3(k) = ((sum(sum(abs(imgMatSup((M+1):y1, (N+1):x1)))) / sum(sum(abs(imgMatSup)))));
%         beta = sum(R3(1:k))/sum(R1(1:k)) * 0.3;

        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) = imgMatSup((M+1):y1, (N+1):x1);       
        
        midSp = fft2(midMat);
        R2(k) = sum(sum((abs(midSp) - fm) .^ 2));
        midSp = fm .* midSp ./ (abs(midSp) + 0.00001);
        imshow(uint8(abs(scale * midMat)));
         
        %imshow(uint8(abs(scale * midMat((M+1):y1, (N+1):x1))));
        %imshow(uint8(80-imag(log(midMat((M+1):y1, (N+1):x1))) * scale2 / pi));
    end
end
R1 = R1 / (O1 * O2 - 1);
R3 = R3 / (O1 * O2 - 1);
R2 = R2 / sum(sum(fm .^ 2)); 

A = midMat ./ (imgMatSup + 0.000000001);
surf(imag(log(A)))
surf(abs(A))