散斑照明下的盲卷积问题求解（仿真代码）

chunqiulei

设备描述：光从光纤出射经准直，照射到毛玻璃后，形成散斑光，在经过样品，在经过毛玻璃，最后的形成图案经相机采集。
目标：恢复样品的形貌特征。
问题：未知光照射样品，求解样品特征。
仿真代码如下：
%% 浑浊层分辨率增强成像
%  上述为仿真的名字： 其意思 光通过一个散射介质后照射到 物体上， 经物镜 收集成像在CCD 上


%================================================================
% 第一步我们定义了激光器的波长、 成像系统的有效像元尺寸 由散射介质造成的NA  
% 散斑NA的估计值 散斑的变化量 沿着一个方向最大的变化值
%  waveLength = 0.6328e-6;    % 采用He-Ne 激光器照明
%  
%  pixesize = 54e-6;          % 在成像系统中有效的像素尺寸
%  
%  diffuserNA = 0.5e-3;       % Low-pass filter of the diffuser
%  illuminationNA = 3e-3;     % 散射斑的数值孔径
%  
%  % 定义移动距离
%  shiftstepsize = 0.75;
%  
%  % 迭代次数
%  niterative = 30;
%  
%  % size of image (assumed to be square)
%  InputImgSize = 135;
%  
%  nPattern = (spiralradius*2+1)^2;   % number of illumination patterns
clear all
clc

%  定义一个初始图像
biaozhun = double(imread('分辨率板_1.bmp'));
biaozhun1 = imcrop(biaozhun,[256-67,256-67,135,135]);

%% Step 1: Set Parameters
WaveLength = 0.632e-6; % Wave length(red)
PixelSize = 54e-6; % Effective pixel size of imaging system
DiffuserNA = 0.5 * 1e-3; % Low-pass filter of the diffuser
IlluminationNA = 3 * 1e-3; % Numerical aperture of speckle patterns
ShiftStepSize = 0.75; % Step size of pattern's shift
SpiralRadius = 20; % Pattern number = (20*2+1)^2
nIterative = 30;
InputImgSize = 135; % Size of image (assumed to be square)
nPattern = (SpiralRadius*2+1)^2; % Number of illumination patterns
InputImgSizeHalf = floor(InputImgSize / 2);
InputImgCenter = ceil(InputImgSize / 2);
MaxShiftInPixel = round(ShiftStepSize * SpiralRadius) * 2;

%% Step 2: Input object image preparation
InputImg = biaozhun1;
InputImg = InputImg./max(max(InputImg)); % Image intensity normalization
InputImgFT = fftshift(fft2(InputImg));
InputImgSizeX = size(InputImg, 1);
InputImgSizeY = size(InputImg, 2);
% Show prepared input object image and its Fourier transform
figure;
subplot(1,2,1); imshow(abs(InputImg),[]); title('Ground Truth');
subplot(1,2,2); imshow(log(abs(InputImgFT) + 1),[]); title('Ground Truth FT');

%% Step 3: Optical transfer function (OTF) and target object image generation
CTF = ones(InputImgSizeX,InputImgSizeY);
CutoffFreqX = DiffuserNA*(1/WaveLength)*InputImgSizeX*PixelSize;
CutoffFreqY = DiffuserNA*(1/WaveLength)*InputImgSizeY*PixelSize;
[GridX, GridY] = meshgrid(1:InputImgSizeX,1:InputImgSizeY);
CTF = CTF.*(((GridY- (InputImgSizeY+1)/2)/CutoffFreqY).^2+((GridX-(InputImgSizeX+1)/2)/CutoffFreqX).^2<=1);
CTFSignificantPix = numel(find(abs(CTF)>eps(class(CTF))));
ifftscale = numel(CTF)/CTFSignificantPix;
aPSF = fftshift(ifft2(ifftshift(CTF)));
iPSF = ifftscale*abs(aPSF).^2;
OTF = fftshift(fft2(ifftshift(iPSF)));
LRTempTargetImgFT=fftshift(fft2(imresize(InputImg,1))).* OTF; % Low resolution target image FT
InputImgLR=abs(ifft2(ifftshift(LRTempTargetImgFT))); % Low resolution target image
figure; subplot(1,2,1); imshow(InputImgLR,[]); title('Low-resolution input image');
subplot(1,2,2); imshow(log(abs(LRTempTargetImgFT) + 1),[]); title('Low-resolution input image FT');

%% Step 4: Mother Speckle pattern generation
% generate random speckle pattern
MotherSpeckleSizeX = InputImgSizeX + MaxShiftInPixel;
MotherSpeckleSizeY = InputImgSizeY + MaxShiftInPixel;
randomAmplitude = imnoise(ones(MotherSpeckleSizeX,MotherSpeckleSizeY),'speckle',0.5);     % 散斑图样的复振幅模
randomPhase = imnoise(ones(MotherSpeckleSizeX,MotherSpeckleSizeY),'speckle',0.5);         % 散斑图案的相位
randomPhase = 0.5*pi*randomPhase./max(max(randomPhase));     % 对散斑相位进行归一化
randomSpeckle = randomAmplitude.*exp(1j.*randomPhase);       % 散斑的分布
randomSpeckleFT = fftshift(fft2(randomSpeckle));             % 散斑复振幅的归一化
% speckle pattern lowpass filter
specklePatternCTF = ones(MotherSpeckleSizeX,MotherSpeckleSizeY);
CutoffFreqX = IlluminationNA * (1/WaveLength)*MotherSpeckleSizeX* PixelSize;
CutoffFreqY = IlluminationNA * (1/WaveLength)*MotherSpeckleSizeY* PixelSize;
[GridX, GridY] = meshgrid(1:MotherSpeckleSizeX,1:MotherSpeckleSizeY);
specklePatternCTF = specklePatternCTF.*(((GridY-(MotherSpeckleSizeY+1)/2) /CutoffFreqY).^2+((GridX-(MotherSpeckleSizeX+1)/2)/CutoffFreqX).^2<=1);
% lowpassed speckle intensity
aMotherSpeckle = ifft2(ifftshift(specklePatternCTF.*randomSpeckleFT));
iMotherSpeckle = (abs(aMotherSpeckle)).^2;
MotherSpeckle = iMotherSpeckle./max(max(iMotherSpeckle));

%% Step 5: Shifted speckle pattern generation
LocationX = zeros(1, nPattern);
LocationY = zeros(1, nPattern);
ShiftY = zeros(1, nPattern);
ShiftX = zeros(1, nPattern);
SpiralPath = spiral(2*SpiralRadius+1);      
SpiralPath = flipud(SpiralPath);           % 在这里矩阵翻转，主要是由于矩阵的
for iShift = 1:nPattern
[iRow, jCol] = find(SpiralPath == iShift);
LocationX(iShift) = iRow;
LocationY(iShift) = jCol;
end;
ShiftedSpecklePatterns = zeros(size(InputImg,1), size(InputImg,2), nPattern);
for iPattern = 1:nPattern
    ShiftedMotherSpeckle = subpixelshift(MotherSpeckle, ...
    LocationX(1,iPattern) * ShiftStepSize * PixelSize,...
    LocationY(1,iPattern) * ShiftStepSize * PixelSize,...
    PixelSize); % Generate speckle pattern sequence
    ShiftedSpecklePatterns(:,:,iPattern) = ShiftedMotherSpeckle(MotherSpeckleSizeX/2-InputImgSizeX/2:MotherSpeckleSizeX/2+InputImgSizeX/2-1,MotherSpeckleSizeX/2-InputImgSizeX/2:MotherSpeckleSizeX/2+InputImgSizeX/2-1);
    if iPattern < (nPattern)
        ShiftX(iPattern+1)=(LocationX(1,iPattern+1)-LocationX(1,iPattern)).* ShiftStepSize;
        ShiftY(iPattern+1) = (LocationY(1,iPattern+1)-LocationY(1,iPattern)).* ShiftStepSize;
    end
     disp(iPattern);
end

%% Step 6: Low-resolution target image generation
TargetImgs = zeros(size(ShiftedSpecklePatterns));
nTargetImgs = nPattern;% Total number of target low-resolution images
for iTargetImg = 1:nTargetImgs
    TargetImgs(:,:,iTargetImg) = abs(ifft2(ifftshift(fftshift(fft2(InputImg.*ShiftedSpecklePatterns(:,:,iTargetImg))).*OTF)));
    disp(iTargetImg);
end;


%% Step 7: Iterative reconstruction
ImgRecovered = mean(TargetImgs,3); % Initial guess of the object,这里描述的是将采集到的图进行平均
MotherSpeckleRecovered = ones(InputImgSizeX + MaxShiftInPixel, InputImgSizeY + MaxShiftInPixel); % Initial guess of the speckle pattern
for iterative = 1:nIterative
    for iShift = 1:nPattern
        display([iterative iShift])
        MotherSpeckleRecovered = subpixelshift(MotherSpeckleRecovered,ShiftX(1,iShift)*PixelSize, ShiftY(1,iShift)*PixelSize, PixelSize);
        MotherSpeckleRecoveredCropped =MotherSpeckleRecovered(round(MotherSpeckleSizeX/2-InputImgSizeX/2):round(MotherSpeckleSizeX/2+InputImgSizeX/2-1),round(MotherSpeckleSizeY/2 -InputImgSizeY/2):round(MotherSpeckleSizeY/2+InputImgSizeY/2-1));
        TempTargetImg = ImgRecovered .* MotherSpeckleRecoveredCropped;
        TempTargetImgCopy = TempTargetImg;  % Use it to recover Itn
        TempTargetImgFT = fftshift(fft2(TempTargetImg));  % 2D Fourier transform
        LRTempTargetImgFT = OTF .* TempTargetImgFT;  % Lowpass the image
        LRTempTargetImg = ifft2(ifftshift(OTF .* TempTargetImgFT)); % Inverse FT
        LRTempTargetImg_AmpUpdated = TargetImgs(:,:,iShift);
        LRTempTargetImg_AmpUpdated_FT = fftshift(fft2(LRTempTargetImg_AmpUpdated));
        TempTargetImgFT =TempTargetImgFT+conj(OTF)./(max(max((abs(OTF)).^2))).*(LRTempTargetImg_AmpUpdated_FT- LRTempTargetImgFT);
        % Update the target image
        if (iterative > 2)
            OTF = OTF + conj(TempTargetImgFT)./(max(max((abs(TempTargetImgFT)).^2))).*(LRTempTargetImg_AmpUpdated_FT-LRTempTargetImgFT);
        end
        TempTargetImg = ifft2(ifftshift(TempTargetImgFT));
        ImgRecovered = ImgRecovered+MotherSpeckleRecoveredCropped.*(TempTargetImg-TempTargetImgCopy)./(max(max(MotherSpeckleRecoveredCropped))).^2;
        % Update the object
        ImgRecovered = abs(ImgRecovered);
        MotherSpeckleRecoveredCropped = MotherSpeckleRecoveredCropped+ImgRecovered.*(TempTargetImg-TempTargetImgCopy)./(max(max(ImgRecovered))).^2;
        MotherSpeckleRecoveredCropped = abs(MotherSpeckleRecoveredCropped);
        MotherSpeckleRecovered(round(MotherSpeckleSizeX/2 - InputImgSizeX/2):round (MotherSpeckleSizeX/2+InputImgSizeX/2-1),round(MotherSpeckleSizeY/2-InputImgSizeY/2):round(MotherSpeckleSizeY/2+InputImgSizeY/2-1))= abs(MotherSpeckleRecoveredCropped);
    end
    MotherSpeckleRecovered = subpixelshift(MotherSpeckleRecovered,-sum(ShiftX)*PixelSize, -sum(ShiftY)*PixelSize, PixelSize);
    ImgRecoveredFT = fftshift(fft2(ImgRecovered));
    figure; subplot(1,2,1); imshow(abs(ImgRecovered),[]); title('Recovered image');
    subplot(1,2,2); imshow(log(abs(ImgRecoveredFT)+1),[]); title('Recovered FT');
    pause(0.5)
    subplot(1,2,2); imshow(log(abs(ImgRecoveredFT)+1),[]); title('Recovered FT');
    pause(0.5)
end




function output_image=subpixelshift(input_image,xshift,yshift,spsize)
[m,n,num]=size(input_image);
[FX,FY]=meshgrid(1:m,1:n);
    for i=1:num
        Hs=exp(-1j*2*pi.*(FX.*(xshift(1,i)/spsize)/m+FY.*(yshift(1,i)/spsize)/n));
        output_image(:,:,i)=abs(ifft2(ifftshift(fftshift(fft2(input_image(:,:,i))).*Hs)));
    end
end