MATLAB处理自旋极化电子显微镜数据自动化（1）- imagestacks registration

1. 设置filter读取文件夹中的图像数据并且自动stack并对部分stack去除噪音背景—函数AutoStack()

function []= AutoStack()

bg = imread("bg.tif");

Folder = pwd; % Get current path
dataPath = fullfile(Folder,"All");
pattern_up = '_01';
pattern_down = '_02';
pattern_SPLEEM = '_03';

fileList = dir(dataPath);
up = {};
down = {};
SPLEEM = {};

for i = 1:length(fileList)
    fileName = fileList(i).name;

    % read LEEM_up
    if contains(fileName,pattern_up) && endsWith(fileName, '.tif')
        fullPath = fullfile(dataPath, fileName);
        up{end+1} = imread(fullPath);
    end

    % read LEEM_down
    if contains(fileName,pattern_down) && endsWith(fileName, '.tif')
        fullPath = fullfile(dataPath, fileName);
        down{end+1} = imread(fullPath);
    end

    % read SPLEEM
    if contains(fileName,pattern_SPLEEM) && endsWith(fileName, '.tif')
        fullPath = fullfile(dataPath, fileName);
        SPLEEM{end+1} = imread(fullPath);
    end

end

numFrames = length(up);
for i = 1:numFrames
    % Deducting hexagonal backgrounds
    up_hex = double(cell2mat(up(i)));
    down_hex = double(cell2mat(down(i)));
    SPLEEM_i = cell2mat(SPLEEM(i));
    up_i = up_hex./bg;
    up_i = uint16(up_i);
    down_i = down_hex./bg;
    down_i = uint16(down_i);
    imwrite(up_i,'Stack_LEEMup.tif',"WriteMode","append")
    imwrite(down_i,'Stack_LEEMdown.tif',"WriteMode","append")
    imwrite(SPLEEM_i,'Stack_SPLEEM.tif',"WriteMode","append")
end

end

2. 对imagestack按照指定区域特征registration—函数AutoDC()

function []= AutoDC()

%file = dir('*.tif');
% Specify the image stack file path
DCStackFile = "DCstack.tif";

% Get information about the image stack using imfinfo
imageInfo = imfinfo(DCStackFile);
numFrames = numel(imageInfo); % Number of frames in the image stack

% Preallocate a 3D array to store the image stack
imageStack1 = zeros(imageInfo(1).Height, imageInfo(1).Width, numFrames, 'uint16');

% Read each frame in the image stack
for i = 1:numFrames
    imageStack1(:, :, i) = imread(DCStackFile, i);
end

% Specify the image stack file path
SPLEEMStackFile = "Stack_SPLEEM.tif";
LEEMupStackFile = "Stack_LEEMup.tif";
LEEMdownStackFile = "Stack_LEEMdown.tif";

% Get information about the image stack using imfinfo
imageInfo = imfinfo(SPLEEMStackFile);
numFrames = numel(imageInfo); % Number of frames in the image stack

% Preallocate a 3D array to store the image stack
imageStack2 = zeros(imageInfo(1).Height, imageInfo(1).Width, numFrames, 'uint16');
imageStack3 = zeros(imageInfo(1).Height, imageInfo(1).Width, numFrames, 'uint16');
imageStack4 = zeros(imageInfo(1).Height, imageInfo(1).Width, numFrames, 'uint16');

% Read each frame in the image stack
for i = 1:numFrames
    imageStack2(:, :, i) = imread(SPLEEMStackFile, i);
    imageStack3(:, :, i) = imread(LEEMupStackFile, i);
    imageStack4(:, :, i) = imread(LEEMdownStackFile, i);
end


% Initialize the corrected image stack
correctedImageStack = zeros(size(imageStack1));
correctedSPLEEMStack = zeros(size(imageStack2));
correctedLEEMupStack = zeros(size(imageStack3));
correctedLEEMdownStack = zeros(size(imageStack4));
TransformationVector = zeros(1,2,numFrames);

% Loop over each frame in the image stack
for i = 1:numFrames
    % First frame is the reference frame
    if i == 1
        correctedImageStack(:,:,i) = imageStack1(:,:,i);
        referenceFrame = imageStack1(:,:,i);
    else
%         figure(1)
%         imshowpair(referenceFrame,imageStack(:,:,i))
        % Compute cross-correlation between current frame and reference frame
        crossCorr = normxcorr2(referenceFrame, imageStack1(:,:,i));
%         figure(2)
%         surf(crossCorr)
%         shading flat
        % Find the peak of the cross-correlation
        [~, peakLoc] = max(crossCorr(:));
        [ypeak, xpeak] = ind2sub(size(crossCorr), peakLoc);

        % Compute the drift vector
        yoffSet = ypeak-size(referenceFrame,1);
        xoffSet = xpeak-size(referenceFrame,2);
        driftVector = [xoffSet, yoffSet];
        TransformationVector(:,:,i) = driftVector;
        % Shift the current frame to correct for drift
        correctedImageStack(:,:,i) = imtranslate(imageStack1(:,:,i),-driftVector);
        correctedSPLEEMStack(:,:,i) = imtranslate(imageStack2(:,:,i),-driftVector);
        correctedLEEMupStack(:,:,i) = imtranslate(imageStack3(:,:,i),-driftVector);
        correctedLEEMdownStack(:,:,i) = imtranslate(imageStack4(:,:,i),-driftVector);
    end
end

correctedImageStack = uint16(correctedImageStack);
% Display the corrected image stack
% montage(correctedImageStack, 'DisplayRange', [])
for i = 1:numFrames
    imwrite(correctedImageStack(:,:,i),'DCstack_DC.tif',"WriteMode","append")
end

% % mean specific slice: 1 to 175
% correctedSPLEEMStack = correctedSPLEEMStack(:,:,1:175);
% correctedLEEMupStack = correctedLEEMupStack(:,:,1:175);
% correctedLEEMdownStack = correctedLEEMdownStack(:,:,1:175);

SPLEEM_DC = mean(correctedSPLEEMStack,3); 
SPLEEM_DC = uint16(SPLEEM_DC);

LEEMup_DC = mean(correctedLEEMupStack,3);
LEEMup_DC = uint16(LEEMup_DC);

LEEMdown_DC = mean(correctedLEEMdownStack,3);
LEEMdown_DC = uint16(LEEMdown_DC);

LEEM_DC = LEEMup_DC + LEEMdown_DC;

% imwrite cannot write 32-bit data
imwrite(SPLEEM_DC,'DC_SPLEEM.tif')
imwrite(LEEMup_DC,'DC_LEEMup.tif')
imwrite(LEEMdown_DC,'DC_LEEMdown.tif')
imwrite(LEEM_DC,'DC_LEEM.tif')

% t = Tiff('SPLEEM_DC.tiff','w');
% 
% % Setup tags
% % Lots of info here:
% % http://www.mathworks.com/help/matlab/ref/tiffclass.html
% tagstruct.ImageLength     = size(SPLEEM_DC,1);
% tagstruct.ImageWidth      = size(SPLEEM_DC,2);
% tagstruct.Photometric     = Tiff.Photometric.MinIsBlack;
% tagstruct.BitsPerSample   = 32;
% tagstruct.SamplesPerPixel = 1;
% tagstruct.RowsPerStrip    = 16;
% tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
% tagstruct.Software        = 'MATLAB';
% t.setTag(tagstruct)
% t.write(SPLEEM_DC);
% t.close();

end

3 对不同自旋极化条件下图像间手动Drift-Correction并同时扣除不均匀光照背景—函数DC_3image()

function []= DC_3image()

% Based on imageJ positions
p_ip1 = [662 900]; % Location of the reference in ip1 -the smaller one
p_ip2 = [663 898]; % Location of the reference in ip2
p_pp = [663 905]; % Location of the reference in pp

Folder = pwd; % Get current path
ip1_LEEM = imread(Folder+"\LEEM_DC\ip90_LEEM_DC.tif");
ip2_LEEM = imread(Folder+"\LEEM_DC\ip180_LEEM_DC.tif");
pp_LEEM  = imread(Folder+"\LEEM_DC\pp_LEEM_DC.tif");
ip1_SPLEEM = imread(Folder+"\SPLEEM_DC\ip90_SPLEEM_DC.tif");
ip2_SPLEEM = imread(Folder+"\SPLEEM_DC\ip180_SPLEEM_DC.tif");
pp_SPLEEM  = imread(Folder+"\SPLEEM_DC\pp_SPLEEM_DC.tif");

% Drift correction
% imtranslate(image,[x y])
ip1_SPLEEM = imtranslate(ip1_SPLEEM,p_pp-p_ip1);
ip2_SPLEEM = imtranslate(ip2_SPLEEM,p_pp-p_ip2);
ip1_LEEM = imtranslate(ip1_LEEM,p_pp-p_ip1);
ip2_LEEM = imtranslate(ip2_LEEM,p_pp-p_ip2);

% crop image
[ip1_LEEM_crop,rectout] = imcrop(ip1_LEEM);
ip2_LEEM_crop = imcrop(ip2_LEEM,rectout);
pp_LEEM_crop = imcrop(pp_LEEM,rectout);
ip1_SPLEEM_crop = imcrop(ip1_SPLEEM,rectout);
ip2_SPLEEM_crop = imcrop(ip2_SPLEEM,rectout);
pp_SPLEEM_crop = imcrop(pp_SPLEEM,rectout);

% Deducting inhomogeneous backgrounds
f = 200;
ip1_LEEM_bg = imgaussfilt(ip1_LEEM_crop,f);
ip1_LEEM_dbg = double(ip1_LEEM_crop)./double(ip1_LEEM_bg);
ip1_LEEM = mat2gray(ip1_LEEM_dbg);
ip2_LEEM_bg = imgaussfilt(ip2_LEEM_crop,f);
ip2_LEEM_dbg = double(ip2_LEEM_crop)./double(ip2_LEEM_bg);
ip2_LEEM = mat2gray(ip2_LEEM_dbg);
pp_LEEM_bg = imgaussfilt(pp_LEEM_crop,f);
pp_LEEM_dbg = double(pp_LEEM_crop)./double(pp_LEEM_bg);
pp_LEEM = mat2gray(pp_LEEM_dbg);

% convert to uint16
ip1_LEEM = im2uint16(ip1_LEEM);
ip2_LEEM = im2uint16(ip2_LEEM);
pp_LEEM = im2uint16(pp_LEEM);

% save DC result
imwrite(ip1_LEEM,'Stack_LEEM.tif',"WriteMode","append")
imwrite(ip2_LEEM,'Stack_LEEM.tif',"WriteMode","append")
imwrite(pp_LEEM,'Stack_LEEM.tif',"WriteMode","append")
imwrite(ip1_SPLEEM_crop,'Stack_SPLEEM.tif',"WriteMode","append")
imwrite(ip2_SPLEEM_crop,'Stack_SPLEEM.tif',"WriteMode","append")
imwrite(pp_SPLEEM_crop,'Stack_SPLEEM.tif',"WriteMode","append")

end