Enhances spots and creates gimplike screen version

screens.m

clear all;

base=5;

if(base==1)
b3=imread('MAX_base_1_ch3.tif');
b4=imread('MAX_base_1_ch4.tif');
b5=imread('MAX_base_1_ch5.tif');
b6=imread('MAX_base_1_ch6.tif');
finalname='b1-color.tif';
end;

if(base==2)
b3=imread('MAX_base_2_ch3.tif');
b4=imread('MAX_base_2_ch4.tif');
b5=imread('MAX_base_2_ch5.tif');
b6=imread('MAX_base_2_ch6.tif');
finalname='b2-color.tif';
end;

if(base==3)
b3=imread('MAX_base_3_ch3.tif');
b4=imread('MAX_base_3_ch4.tif');
b5=imread('MAX_base_3_ch5.tif');
b6=imread('MAX_base_3_ch6.tif');
finalname='b3-color.tif';
end;

if(base==4)
b3=imread('MAX_base_4_ch3.tif');
b4=imread('MAX_base_4_ch4.tif');
b5=imread('MAX_base_4_ch5.tif');
b6=imread('MAX_base_4_ch6.tif');
finalname='b4-color.tif';
end;

if(base==5)
b3=imread('MAX_base_5_ch3.tif');
b4=imread('MAX_base_5_ch4.tif');
b5=imread('MAX_base_5_ch5.tif');
b6=imread('MAX_base_5_ch6.tif');
finalname='b5-color.tif';
end;

[r,c]=size(b3(:,:));
b3th=imtophat(b3,strel('disk',6)); b4th=imtophat(b4,strel('disk',6));
b5th=imtophat(b5,strel('disk',6)); b6th=imtophat(b6,strel('disk',6));

clear b3 b4 b5 b6;

t3=(prctile(b3th(:),100)-prctile(b3th(:),99.5))/2; t4=(prctile(b4th(:),100)-prctile(b4th(:),99.5))/2;
t5=(prctile(b5th(:),100)-prctile(b5th(:),99.5))/2; t6=(prctile(b6th(:),100)-prctile(b6th(:),99.5))/2;

ts = zeros(4,1,'double');
ts(1)=t3; ts(2)=t4; ts(3)=t5; ts(4)=t6;

clear t3 t4 t5 t6;

tsmax=max(ts);
m=mean(ts);
st=std(ts);

for i=1:4
  if(ts(i)>(m+st*0.75)) ts(i)=round((m)); end;
end

bools=[];
%maxbinval=61680;%
%minbinval=4112;%
maxbinval=240;
minbinval=16;
mediumval=round((maxbinval- minbinval)/2);
bmax=[];
bmax(1)=max(b3th(:)); bmax(2)=max(b4th(:)); bmax(3)=max(b5th(:)); bmax(4)=max(b6th(:));
bmin=[];
bmin(1)=min(b3th(:)); bmin(2)=min(b4th(:)); bmin(3)=min(b5th(:)); bmin(4)=min(b6th(:));


gb1=zeros(size(b3th(:,:)),'uint8');
gb2=zeros(size(b4th(:,:)),'uint8');
gb3=zeros(size(b5th(:,:)),'uint8');
gb4=zeros(size(b6th(:,:)),'uint8');

b3th=(b3th-bmin(1))*( 255/ (bmax(1)-bmin(1)));
b4th=(b4th-bmin(2))*( 255/ (bmax(2)-bmin(2)));
b5th=(b5th-bmin(3))*( 255/ (bmax(3)-bmin(3)));
b6th=(b6th-bmin(4))*( 255/ (bmax(4)-bmin(4)));

gb1(:,:)=uint8(b3th); gb1(:,:,3)=uint8(b3th);  gb1(:,:,2)=uint8(b3th);
gb2(:,:)=uint8(b4th); gb2(:,:,3)=uint8(b4th);  gb2(:,:,2)=uint8(b4th);
gb3(:,:)=uint8(b5th); gb3(:,:,3)=uint8(b5th);  gb3(:,:,2)=uint8(b5th);
gb4(:,:)=uint8(b6th); gb4(:,:,3)=uint8(b6th);  gb4(:,:,2)=uint8(b6th);

gbts = zeros(4,1,'double');
gbts(1)=(prctile(gb1(:),100)-prctile(gb1(:),99.5))/2;
gbts(2)=(prctile(gb2(:),100)-prctile(gb2(:),99.5))/2;
gbts(3)=(prctile(gb3(:),100)-prctile(gb3(:),99.5))/2;
gbts(4)=(prctile(gb4(:),100)-prctile(gb4(:),99.5))/2;

newgb1=gb1; newgb2=gb2; newgb3=gb3; newgb4=gb4;

% indices=find(gb1>gbts(1));
% for i=1:length(indices)
%   gbindex=gb1(indices(i));
%   newgb1(indices(i)) = gbindex+round(abs( gbindex-gbts(1) )/2);
% end
% 
% indices=find(gb2>gbts(2));
% for i=1:length(indices)
%   gbindex=gb2(indices(i));
%   newgb2(indices(i)) = gbindex+round(abs( gbindex-gbts(2) )/2);
% end
% 
% indices=find(gb3>gbts(3));
% for i=1:length(indices)
%   gbindex=gb3(indices(i));
%   newgb3(indices(i)) = gbindex+round(abs( gbindex-gbts(3) )/2);
% end
% 
% indices=find(gb4>gbts(4));
% for i=1:length(indices)
%   gbindex=gb4(indices(i));
%   newgb4(indices(i)) = gbindex+round(abs( gbindex-gbts(4) )/2);
% end

clear gb1 gb2 gb3 gb4;

newgb1=uint16(newgb1); newgb2=uint16(newgb2);
newgb3=uint16(newgb3); newgb4=uint16(newgb4);

mgreen=newgb1; mmagenta=newgb2; morange=newgb3; mcyan=newgb4;
mgreen(:,:,1)=103;   mgreen(:,:,2)=255;  mgreen(:,:,3)=0;
mmagenta(:,:,1)=255; mmagenta(:,:,2)=0;  mmagenta(:,:,3)=255;
morange(:,:,1)=255;  morange(:,:,2)=97;  morange(:,:,3)=9;
mcyan(:,:,1)=0;      mcyan(:,:,2)=255;   mcyan(:,:,3)=255;


multip1=((newgb1).*(mgreen))./255;   %multip1=uint8(multip1); %imshow(multip1); figure;
multip2=((newgb2).*(mmagenta))./255; %multip2=uint8(multip2); %imshow(multip2); figure;
multip3=((newgb3).*(morange))./255;  %multip3=uint8(multip3); %imshow(multip3); figure;
multip4=((newgb4).*(mcyan))./255;    %multip4=uint8(multip4); %imshow(multip4);

multall=255-(255-multip1).*(255-multip2)./255;
multall=255-(255-multall).*(255-multip3)./255;
multall=255-(255-multall).*(255-multip4)./255;
multall=uint8(multall); imshow(multall);

imwrite(multall,finalname);
%imagesc(gb1); figure;imagesc(gb2); figure;imagesc(gb3); figure;imagesc(gb4);%