function F = improcMULT2(I,numxs,stage,dark,LOthreshold,plotting,bfrac,dfrac,filename,stret,rim)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ImProcMULT2.m
%
% Same as imProcMULT.m, but with plotting functions mostly disabled. 
% ImProcMULT2 analyses each cross-section separately rather than in bulk 
% because of evolving gray scale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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)
bright=find(I>250);
I(bright)=1;
K=I;
%improve contrast stretch
if stret == 1
   I=stretch(I);
end
%smooth (adaptive filtering)
J=wiener2(I,[10 10]); %20 is good
% 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 %dark = 2 is 'both' selected in gui
            vesDNmaxD(m)=0;
            vesDNmaxB(m)=0
            break
        else
            vesDNmax(m)=0;
	    vesDNmin=5;      
        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=5;
        end
    elseif  isempty(borders) == 1
            vesDNmax(m)=0;
            vesDNmin=5;
    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=1;%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 %ch
     	  ymax=vesDNmaxD(m)*ones(1,length(x));
	  else %ch
		  ymax=vesDNmax(m)*ones(1,length(x)); %ch
	  end
     ymin=vesDNmin*ones(1,length(x));
%     plot(x,ymax,'r',x,ymin,'g')
%     legend('data (windowed)','amgydale DN upperbound', 'amygdale DN 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) & (J<vesDNmax));
    textproc=['processing: ' filename];
    disp(sprintf(textproc));
    %disp(sprintf('amygdale DN cutoff: %0.3g',vesDNmax));
elseif dark == 0
    if (length(vesDNmax) < 3) & (stage == 1)
    	display(sprintf('CAUTION! You should increase your Number of Cross Sections'))
	length(vesDNmax)
	keyboard
    end  
    vesDNmax=mean(vesDNmax)+bfrac*std(vesDNmax);
    amyg = find(J>vesDNmax);
    textproc=['processing: ' filename];
    disp(sprintf(textproc));
    %disp(sprintf('amygdale DN cutoff: %0.3g',vesDNmax));
else
    vesDNmaxD=mean(vesDNmaxD)-dfrac*std(vesDNmaxD);
    amyg = find((J>vesDNmin) & (J<vesDNmaxD));
    vesDNmaxB=mean(vesDNmaxB)+bfrac*std(vesDNmaxB);
    amyg2 = find(J>vesDNmaxB);
    textproc=['processing: ' filename];
    disp(sprintf(textproc));
end

%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);
%puncturing, to avoid getting pancakes
F=imfill(closebw,'holes');
Fsmall= F(475:525,475:525);
if (length(find(Fsmall == 1))==prod(size(Fsmall))) & (sum(sum(Fsmall)) > 1);
   %pancaked!
   closebw(250,1:end)=0;
   disp(sprintf('punctured'))
   F=imfill(closebw,'holes');
end
%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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (stage == 2)
    figure,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('Finished Product')
	final=Jw;
    final(find(perim==1))= 256;
    subplot(236),imshow(final);
	title('Wiener image w/ amygdale contour')
end
