clear all
	
%% Plot the example images

SourceFileCell = {...
	'./ExtractedStacks/Cond_25/Image_2.mat',...
	'./ExtractedStacks/Cond_19/Image_2.mat',...
	'./ExtractedStacks/Cond_15/Image_1.mat',...
	'./ExtractedStacks/Cond_17/Image_1.mat',...
	'./ExtractedStacks/Cond_21/Image_6.mat',...
	};

Ser2P_blurRange = 0.0001; % micrometers
Ser5P_blurRange = 0.0001; % micrometers

SourceFileCell = {...
	'./ExtractedStacks/Cond_25/Image_2.mat',...
	'./ExtractedStacks/Cond_19/Image_2.mat',...
	'./ExtractedStacks/Cond_23/Image_5.mat',...
	'./ExtractedStacks/Cond_15/Image_12.mat',...
    './ExtractedStacks/Cond_17/Image_5.mat',...
    './ExtractedStacks/Cond_21/Image_6.mat',...
	};

zoomLocations = {...
    [23.2,18],...
    [9.8,7.63],...
    [6.2,13.4],...
    [12.5,5.8],...
    [16.5,12.1],...
    [14.6,3.6]};
zoomRange = 3.0; % in microns


Ser2P_blurRange = 0.0001; % micrometers
Ser5P_blurRange = 0.0001; % micrometers

% imgInds = {[1,7],[2,8]}; %[1,3]
imgRanges = {...
	[22,22+32,35,35+32],...
	[8.5,8.5+32,42,42+32],...
	[21,21+32,37,37+32],...
    [13,13+32,18,18+32],...
	[20,20+32,49,49+32],...
	[26,26+32,19,19+32]};
zCoordinate = [19,5,10,18,10,9];
rotate_flag = [false,false,false,false,false,false];


plotTitles = {'Control','3 h RHB','6 h RHB','12 h RHB','24 h RHB','48 h RHB'};

scaleBar = 10.0; % in microns

Ser5P_lim = [30,900];
Ser2P_lim = [120,1500];


numPlots = numel(imgRanges);

figure(1)
clf

for pp = 1:numPlots
	
	pp
	
	thisFilePath = SourceFileCell{pp};
	loadStruct = load(thisFilePath,...
		'imgStack','imgSize','pixelSize','zStepSize');
	imgStack = loadStruct.imgStack;
	imgSize = loadStruct.imgSize;
	pixelSize = loadStruct.pixelSize;
	zStepSize = loadStruct.zStepSize;

	thisImg = imgStack;
	thisPixelSize = pixelSize;
	
	% --- Removal of DNA background
	Ser2P_img = double(thisImg{2}(:,:,zCoordinate(pp)));
	Ser5P_img = double(thisImg{1}(:,:,zCoordinate(pp)));
	
	if rotate_flag(pp)
		Ser2P_img = Ser2P_img';
		Ser5P_img = Ser5P_img';
	end

	Ser5P_img = imgaussfilt(Ser5P_img,...
		Ser5P_blurRange./thisPixelSize);
	Ser2P_img = imgaussfilt(Ser2P_img,...
		Ser2P_blurRange./thisPixelSize);
	
	thisImgRange = round(imgRanges{pp}./thisPixelSize)+1;
	Ser2P_img = Ser2P_img(...
		thisImgRange(1):thisImgRange(2),...
		thisImgRange(3):thisImgRange(4));
	Ser5P_img = Ser5P_img(...
		thisImgRange(1):thisImgRange(2),...
		thisImgRange(3):thisImgRange(4));
	thisSize = size(Ser5P_img);
	
	Ser5P_img = (Ser5P_img-Ser5P_lim(1))./(diff(Ser5P_lim));
	Ser2P_img = (Ser2P_img-Ser2P_lim(1))./(diff(Ser2P_lim));
	
	
	subplot(numPlots,3,(pp-1).*3+1)
	imagesc([0,thisSize(2)].*thisPixelSize,...
		[0,thisSize(1)].*thisPixelSize,...
		Ser5P_img,[0,1])
	axis equal tight
	set(gca,'XTick',[],'YTick',[])
	set(gca,'YDir','normal')
	hold on
	plot([1,1+scaleBar],thisSize(1).*thisPixelSize-[1,1],...
		'w-','LineWidth',3) % Plot scale bar
	textObject = ...
		text(1+0.5.*scaleBar,thisSize(2).*thisPixelSize-3,...
		sprintf('%d \\mum',scaleBar));
	set(textObject,'Color',[1,1,1],'HorizontalAlignment','center',...
		'FontSize',12,'FontName','Helvetica Neue')
	
    if pp == 1
        title('Pol II Ser5P','FontWeight','normal')
    end
    ylabel(plotTitles{pp})
    
    colormap(inferno)
    
    hold on
    thisZoomLoc = zoomLocations{pp};
    plot(thisZoomLoc(1)+zoomRange.*[-0.5,+0.5,+0.5,-0.5,-0.5],...
        thisZoomLoc(2)+zoomRange.*[-0.5,-0.5,+0.5,+0.5,-0.5],...
        'w-','LineWidth',1)

	
	subplot(numPlots,3,(pp-1).*3+2)
	imagesc([0,thisSize(2)].*thisPixelSize,...
		[0,thisSize(1)].*thisPixelSize,...
		Ser2P_img,[0,1])
	axis equal tight
	set(gca,'XTick',[],'YTick',[])
	set(gca,'YDir','normal')
    if pp == 1
        title('Pol II Ser2P','FontWeight','normal')
    end
	
    hold on
    thisZoomLoc = zoomLocations{pp};
    plot(thisZoomLoc(1)+zoomRange.*[-0.5,+0.5,+0.5,-0.5,-0.5],...
        thisZoomLoc(2)+zoomRange.*[-0.5,-0.5,+0.5,+0.5,-0.5],...
        'w-','LineWidth',1)
    
	% --- Ser5P / Ser2P
	% Prepare the color channels for RGB plotting
    
    thisZoomLoc = zoomLocations{pp};
    thisImgRange = ...
        [thisZoomLoc(2)+zoomRange.*[-0.5,+0.5],...
        thisZoomLoc(1)+zoomRange.*[-0.5,+0.5]];
    thisImgRange = round(thisImgRange./thisPixelSize)+1;
    
	Ser2P_img = Ser2P_img(...
		thisImgRange(1):thisImgRange(2),...
		thisImgRange(3):thisImgRange(4));
	Ser5P_img = Ser5P_img(...
		thisImgRange(1):thisImgRange(2),...
		thisImgRange(3):thisImgRange(4));
	thisSize = size(Ser5P_img);
    
    
	this_magenta_plot = Ser5P_img;
	this_green_plot = Ser2P_img;
	thisSize = size(Ser5P_img);
	
    this_magenta_plot = this_magenta_plot-min(this_magenta_plot(:));
    this_magenta_plot = this_magenta_plot./max(this_magenta_plot(:));
%     this_green_plot = this_green_plot-min(this_green_plot(:));
%     this_green_plot = this_green_plot./max(this_green_plot(:));

    
	subplot(numPlots,3,(pp-1).*3+3)
	
	redChannel = this_magenta_plot;
	blueChannel = this_magenta_plot;
	greenChannel = this_green_plot;
	rgb_img = zeros(thisSize(1),thisSize(2),3);
	rgb_img(:,:,1) = redChannel;
	rgb_img(:,:,2) = greenChannel;
	rgb_img(:,:,3) = blueChannel;
	
	image([0,thisSize(2)].*thisPixelSize,...
		[0,thisSize(1)].*thisPixelSize,...
		rgb_img)
	set(gca,'XTick',[],'YTick',[])
	set(gca,'YDir','normal')
    if pp == 1
        title('Detail merge','FontWeight','normal')
    end
	axis equal tight
	
%     thisZoomLoc = zoomLocations{pp};
%     set(gca,'XLim',thisZoomLoc(1)+zoomRange.*[-0.5,+0.5],...
%         'YLim',thisZoomLoc(2)+zoomRange.*[-0.5,+0.5])
    
end