SingleMoleculeAnalysisScripts/virtualAlignmentFunction.m at main · neuhoffmj/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
    function offset = virtualAlignmentFunction(fn, numChannels, useChannels, channelFrames)
    filename = fn;
    contrastBoost = 0.6;

    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);

    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

    offset = [xoffset, yoffset];
    figure;

    imagesc(cor(img_height-15:img_height+15, img_width/numChannels-15:img_width/numChannels+15))
    colorbar
    colormap turbo
    hold on
    viscircles([16 16], 0.8, 'Color','k');
    title(sprintf('Channels %d, %d', useChannels(1), useChannels(2)))
    % 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
end