function F = improcMULT(I,numxs,stage,dark,LOthreshold,plotting,bfrac,dfrac,stret,rim)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ImProcMULT.m
%
% Builds on ImProcGRID. ImProcGRID works great on single images, but
% doesn't do so well processing a bunch (because of evolving gray scales).
% ImProcMULT analyses each cross-section separately rather than in bulk as
% in ImProcGRID. This allows removal of the large and very dark edges,
% which screw the analyses by ravaging the mean.
%
% NOTE: At this stage, the code only detects vesicles *darker* than matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Improve images by adaptive smoothing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%location of cross-section pixels
x=repmat([0 1000],1,numxs); %the images are 1000 pixels wide (and tall)
bux=x; %backup x
int=floor(1000/(numxs+1)); % n+1 intervals yield n cross-sections across
% the image and one is at the bottom (draw it)
y=linspace(1,numxs+1,numxs+1);
y=y(1:end-1); %removes the cross-section that is at the bottom of the image
y=sort([(int*y) (int*y)]); %corresponding y-coordinates
%remove scale (very bright pixels: DN=255)
%I=imread('05017-12-2_rec0075.bmp');
%K=imread('05017-12-2_rec0075raw.bmp');
bright=find(I>250);
I(bright)=1;
K=I;
% figure(1), subplot(211), imshow(K)
% line(x,y,'Color', 'r')
% figure(2), subplot(211), imshow(I)
% line(x,y,'Color', 'r')
%improve contrast stretch
I=stretch(I);
%smooth (adaptive filtering)
J=wiener2(I,[20 20]);
figure(3), subplot(211), imshow(J)
line(x,y,'Color', 'r')
Jw=J; %save in the background...
%plot how good the adaptive filtering is
datai=improfile(I,x,y);
datak=improfile(K,x,y);
figure(1), subplot(212), plot(datak)
figure(2), subplot(212), plot(datak)
hold on
plot(datai,'r');
figure(3), subplot(212), plot(datai)
hold on
dataj=improfile(J,x,y);
plot(dataj,'r')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEPARATE cross-sections for analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Note: For each number of cross sections n, there are 2n-1 cross sections
%taken in total because there is a diagonal one taken between horizontal
%ones.
XStot=2*numxs-1; %total number of cross-sections including the
%diagonal ones
figc=XStot+3; %aside: this is a figure counter to ensure no overlap
data=reshape(dataj,1000,XStot); %convert data vector into a data matrix.
%each column contains one cross-section
%find edges of image using seismic picker
baddata=find(isnan(data)); % looks for NaN values in data
data(baddata)=0; %reset NaN from data matrix
meandata=mean(data); %finds the mean of each column. Needed coz I have to
% add meandata at the end
data=data-repmat(meandata,1000,1); %the mean is removed here
%coz the picker works better
%that way.
[pTime] = pickerRatioMULT(data, 3, 5, 2, stage);
%find amydgale borders
%narrow the search window
x=zeros(2,1); %matrix of zeros
for i=1:length(pTime)
if length(pTime{i})==0
i=i+1;
end
begend=[pTime{i}(1); pTime{i}(end)];
x=[x,begend];
end
x(:,1)=[]; %remove that matrix of zeros
begend=x;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ANALYZE each cross-section individually (since they are of different
%length now)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
matdata=data; %rename data matrix
vesDNmax=[];
vesDNmaxD=[];
vesDNmaxB=[];
for m=1:min(size(data))
data=matdata(begend(1,m):begend(2,m),m);
%add former mean, and subtract new mean
data=data+meandata(m); %the old mean is added here
newmean=mean(data);
data=data-newmean; %remove new mean for seismic picker
[pTime] = pickerRatioSING(data, 3, 5, 2, stage);
%plot(data)
if pTime == 0 %checks if amygdales are indeed detected
if dark == 2
vesDNmaxD(m)=0;
vesDNmaxB(m)=0
break
else
vesDNmax(m)=0;
vesDNmin=20;
end
end
borders = find(diff(pTime)>10); %finds where large intervals between
%borders are.
if max(pTime)-min(pTime) < 10 %checks if only one peak is present
if dark == 2
vesDNmaxD(m)=0;
vesDNmaxB(m)=0
else
vesDNmax(m)=0;
vesDNmin=20;
end
elseif isempty(borders) == 1
vesDNmax(m)=0;
vesDNmin=20;
else
pT = pTime(borders); %find the corresponding distance across
pT1 = pTime(borders);
pT2 = pTime(borders+1);
pT = [pT1;pT2];
pT = sort(pT); %pT is the coordinates of the auto-detected
%"amygdale boundaries"
%set conditions for vesicles DARKER THAN MATRIX: The mean of the data
%between borders must be less than the average of the borders. However,
%ignore data that is too broad (> 200 px) and too dark (<30), as it doesn't
%constitute empty vesicles, but the outside dark area.
flak=[];
%now find vesicles darker than matrix
if dark == 1
vesDNmax=darker(vesDNmax, m, pT, data, flak, newmean);
elseif dark == 0
vesDNmax=brighter(vesDNmax, m, pT, data, flak, newmean);
else
vesDNmaxD=darker(vesDNmaxD, m, pT, data, flak, newmean);
vesDNmaxB=brighter(vesDNmaxB, m, pT, data, flak, newmean);
end
end %check one peak
vesDNmin=20;%min(vesDN)-1+meandata;
%plot those values
%first remean the data
data=data+newmean;
subplot(211),plot(data)
xlabel('pixel across');
ylabel('DN');
title('Original data (if red line = 0, XS not used)')
hold on;
x=linspace(1,length(data),length(data));
if dark == 2
ymax=vesDNmaxD(m)*ones(1,length(x));
else
ymax=vesDNmax(m)*ones(1,length(x));
end
ymin=vesDNmin*ones(1,length(x));
plot(x,ymax,'r',x,ymin,'g')
legend('data (windowed)','amgydale DN upperbound','Dark amydale lowerbound')
end %end of m=1:min(size(data))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Apply results to entire image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%remove all "0" values in vesDNmax
if dark == 2
nodata=find(vesDNmaxD==0);
vesDNmaxD(nodata)=[];
nodata=find(vesDNmaxB==0);
vesDNmaxB(nodata)=[];
else
nodata=find(vesDNmax==0);
vesDNmax(nodata)=[];
end
%identify vesDNmax as the mean of those found & identify such pixels in
%image
if dark == 1
vesDNmax=mean(vesDNmax)-dfrac*std(vesDNmax);
amyg = find((J>vesDNmin) & (JvesDNmax);
else
vesDNmaxD=mean(vesDNmaxD)-dfrac*std(vesDNmaxD);
amyg = find((J>vesDNmin) & (JvesDNmaxB);
end
disp(sprintf('amygdale DN cutoff: %0.3g', vesDNmax));
%identify such pixels in image
J(amyg)=256;
if dark == 2
J(amyg2)=256;
end
Wp=J; %save the image in the background..
figure(figc+1), subplot(211), imshow(J)
x=bux;
line(x,y,'Color', 'r')
subplot(212), plot(dataj)
hold on
dataj=improfile(J,x,y);
plot(dataj,'r')
%%%%%%%%%%%%%%%%%%%%%%
% BINARIZE
%%%%%%%%%%%%%%%%%%%%%%
graypix=find(Wp<255);
Wp(graypix)=0;
J=im2bw(Wp); %binarize
Wnb=J; %save in the background...
if (dark == 1) | (dark == 2) %boundary only present when detecting dark amygdales
Wnb=roiset(J,rim); %loads function rim - eliminates the white boundary
J=Wnb;
end
figure(figc+2),subplot(211), imshow(J)
line(x,y,'Color', 'r')
subplot(212), plot(dataj)
hold on
dataj=improfile(J,x,y);
plot(improfile(J,x,y),'r')
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PUTZ
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%morphological opening: removes "snowflakes" having a radius less than 2
%pixels by opening it with a disk-shaped structuring element having a 2
%pixel radius.
se = strel('disk',2);
F = imopen(J,se);
figure(figc+3), subplot(211),imshow(F)
line(x,y,'Color', 'r')
subplot(212), plot(dataj)
hold on
dataf=improfile(F,x,y);
plot(dataf,'r')
%morphological closing: fills open amygdales by first filling the gaps
%between them and smoothing their outer edges: http://tinyurl.com/imclose
se=strel('disk',5);
closebw=imclose(F,se);
F=imfill(closebw,'holes');
%F=rounded(F,LOthreshold,plotting); %detects roundedness of amydales to remove
%non-amydgale detection.
perim=bwperim(F); %find the perimeter of amygdales
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Compare all steps of image processing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(figc+4),subplot(231),imshow(K);
title('original image')
subplot(232),imshow(Jw);
title('Wiener smoothed')
subplot(233),imshow(Wp);
title('Wiener image w/ detected amygdales')
subplot(234),imshow(Wnb);
title('Rim boundary removed')
subplot(235),imshow(F);
title('Final product')
final=Jw;
final(find(perim==1))= 256;
subplot(236),imshow(final);
title('Wiener image w/ amygdale contour')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Compare with automatic edge detector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
amygedge(F,Jw,final,LOthreshold,plotting,rim)