SingleMoleculeAnalysisScripts/virtualAlignment.m at main · NilsSomm/SingleMoleculeAnalysisScripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
clear
close all

filename = ('20240729_2.5nMPho4dN7pmcbf1smFRETNucS_50mw5hz_2bin_bgsub.tif');
% filename = 'robinmovie.tif';

numChannels = 2;
% First channel is the find maxima channel
useChannels = [1, 2];
% Frames to average over for each channel
channelFrames = [[1 10]; [24 30]];

% useChannels = [1, 1];
% channelFrames = [[1 10]; [1 10]];

testAlignment = 0;
testhshift = 1;
testwshift = 1; %shifts channel 2 up and to the left
contrastBoost = 0.8;

info = imfinfo(filename);
img_height = info(1,1).Height;
img_width = info(1,1).Width;

channelWidth = round(img_width/numChannels);

channel1avg = imgAvg(filename, info, useChannels(1), channelWidth, channelFrames(1,:), img_height, img_width);
channel2avg = imgAvg(filename, info, useChannels(2), channelWidth, channelFrames(2,:), img_height, img_width);

channelX = channel2avg;

%%Test misalignment
if testAlignment
    for i = 1:img_height
        for j = 1:channelWidth
            if ((i+testhshift)>0)&&((i+testhshift)<=img_height)...
                    &&((j+testwshift)>0)&&((j+testwshift)<=channelWidth)
                channelX(i, j) = channel2avg(i+testhshift, j+testwshift);
            end
        end
    end
    channel2avg = channelX;
end

compositeImage = zeros(size(channel1avg, 1),size(channel1avg, 2),3);
channel1Scaled = rescale(channel1avg);
channel2Scaled = rescale(channel2avg);
channel1Contrast = channel1Scaled * 1/contrastBoost;
channel1Contrast(channel1Contrast>1) = 1;
channel2Contrast = channel2Scaled * 1/contrastBoost;
channel2Contrast(channel2Contrast>1) = 1;

cor = xcorr2(channel1Scaled, channel2Scaled);
[val,ind] = max(cor(:));
[ij,ji] = ind2sub(size(cor),ind);

xoffset = -(ji - channelWidth) % Positive means Channel 2 (Cy3) is shifted to the right
yoffset = -(ij - img_height) % Positive means Channel 2 (Cy3) is shifted down

figure;
imagesc(cor(img_height-15:img_height+15, img_width/2-15:img_width/2+15))
colorbar
colormap turbo
% set(gca,'ColorScale','log')

compositeImage(:, :, 1) = channel1Contrast;
% compositeImage(:, :, 2) = ones(size(channel1avg, 1),size(channel1avg, 2));
compositeImage(:, :, 2) = channel2Contrast;
figure;
subplot(1,2,1)
imshow(compositeImage)
title('Original Alignment')
subplot(1,2,2)
channelShift = channel2Contrast;


for i = 1:img_height
    for j = 1:channelWidth
        if ((i+yoffset)>0)&&((i+yoffset)<=img_height)...
                &&((j+xoffset)>0)&&((j+xoffset)<=channelWidth)
            channelShift(i, j) = channel2Contrast(i+yoffset, j+xoffset);
        end
    end
end

alignedComposite = compositeImage;
alignedComposite(:, :, 2) = channelShift;
imshow(alignedComposite)
title('Adjusted Alignment')

function avgimg = imgAvg(filename, info, channel, channelWidth, boundaryFrames, img_height, img_width)
    firstframe = boundaryFrames(1);
    lastframe = boundaryFrames(2);
    avgimg = zeros(img_height, channelWidth, 'uint16');
    imgstk = zeros(img_height, img_width, lastframe-firstframe+1, 'uint16');
    for k = firstframe:lastframe
        imgstk(:,:,k+1-firstframe) = imread(filename, k, 'Info', info);
    end
    % create the merged image from the average of the image stack
    for i = 1:img_height
        for j = (channel-1)*channelWidth+1:(channel-1)*channelWidth+channelWidth
            avgimg(i,j-(channel-1)*channelWidth) = mean(imgstk(i,j,:)); %imgmrg(i,j) = mean(imgstk(i,j,:));
        end
    end
end