%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computes the integrated circular patterns in each colour 
% of an image stack
% April 2016
% Author: L. Muresan, lam94@cam.ac.uk
% Input: stack (First image CFP, second image GFP, third image RFP
% fourth image brightfield, fifth image DAPI)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters:
% r - Maximum radius for which profiles are computed - in pixels
r = 110
% w - width of circular bands
w =5

%% Select and read stack of 4 colours
sep = '/'
clear im
close all
[fn pt] = uigetfile('*.tif');
im = Read3d(strcat(pt, fn));
disp('Reading data...')

%%
% For the 5-image samples:
if size(im,3)<4
    error('Error: Less than four images in the stack!')
else
    CFP = im(:,:,1);
    GFP = im(:,:,2);
    RFP = im(:,:,3);
    BF = im(:,:,4);
    if size(im,3)==5
        DAPI = im(:,:,5);
    else
        DAPI = BF;
    end
end

figure; 
subplot(131);imagesc(CFP); title('CFP'); axis tight equal
subplot(132); imagesc(GFP); title('GFP'); axis tight equal
subplot(133); imagesc(RFP); title('RFP'); axis tight equal

%%
disp('Manual segmentation...')
figure;
BWx = roipoly((GFP/max(GFP(:)))+(RFP/max(RFP(:)))+(CFP/max(CFP(:))));
imshow(BWx, [])

%% Plot the result of the input perimeter

[y, x] = find(BWx>0);
ix = mean(x)
iy = mean(y)

h = figure; imshow(DAPI, []); hold on; plot(ix,iy, 'rx')
perBW = BWx-imerode(BWx, strel('disk',1));
[py px] = find(perBW>0);
plot(px,py, 'r.')
title('DAPI (BF) check of the perimeter')

print(h, strcat(pt,sep,'Fig',strtok(fn, '.'), '_perimeter.pdf'),'-dpdf');
imwrite(BWx, strcat(pt,sep,'Mask',strtok(fn, '.'), '.tif'), 'tif');

%% From center to periphery
disp('Computing profiles from centre to periphery...')
imgB=imgpolarcoordMod(DAPI, ix, iy);
imgG=imgpolarcoordMod(GFP,ix, iy);
imgR=imgpolarcoordMod(double(RFP), ix, iy);
imgC=imgpolarcoordMod(CFP, ix, iy);

figure; 
plot(sum(imgG,2),'g', 'LineWidth',2); xlabel('radius'); hold on; 
title('Profiles (Centre to periphery)')
plot(sum(imgR,2), 'r', 'LineWidth',2); 
plot(sum(imgC,2), 'c', 'LineWidth',2);
res = [(1:size(imgG,1))'-1 sum(imgG,2) sum(imgR,2) sum(imgC,2)];

%% From periphery to centre
disp('Computing profiles from periphery to centre...')
d = bwdist(~BWx);

c = 1;
clear tG tR tC dd
for dist = 0:w:max(d(:))-w
    dist;    
    id = find((d>=dist)&(d<dist+w));
    dd(c) = dist;    
    tR(c) = sum(RFP(id))/length(id);
    tG(c) = sum(GFP(id))/length(id);
    tC(c) = sum(CFP(id))/length(id);    
    c = c+1;
end

figure;
plot(fliplr(tG),'g', 'LineWidth',2); xlabel('radius'); hold on; 
title('Profiles (Periphery to Centre)')
plot(fliplr(tR), 'r', 'LineWidth',2); 
plot(fliplr(tC), 'c', 'LineWidth',2);
res2 = [dd', tG', tR', tC']; 

%% Saving results
disp('Saving results')
% pixel, G, R, C, CR, jR, jC 
csvwrite(strcat(pt,strtok(fn,'.'), 'CtoP.csv'), res)
csvwrite(strcat(pt,strtok(fn,'.'), 'PtoC.csv'), res2)
disp('Finished!')