function find_cell_codes( dna_image_file, ...
    cell_image_file, ...
    prot_image_file, ...
    mask_image_file, savepath, param,currfile )
%FIND_CELL_CODES Calculates cell codes

% Tao Peng
%
% Copyright (C) 2011-2013 Murphy Lab
% Lane Center for Computational Biology
% School of Computer Science
% Carnegie Mellon University
%
% March 12, 2012 Added another preprocessing routine. This routine mimics
%                the preprocessing from A. Shariff, G. K. Rohde and R. F. Murphy (2010) A
%                Generative Model of Microtubule Distributions, and Indirect Estimation of
%                its Parameters from Fluorescence Microscopy Images. Cytometry Part A 77A:457-466.
% March 19, 2012 I. Cao-Berg Added Active Countour 3D segmentation
% March 22, 2012 I. Cao-Berg Added param.downsample so that user can select a
%                downsampling scale. The default is [1 1 1] which is
%                the one in the original generative model paper
% March 23, 2012 I. Cao-Berg Fixed a bug when getting cellnum from impath
% March 27, 2012 I. Cao-Berg When calling ml_parsecell instead of sending a slice
%                of the nucleus, now we are sending a z-projection of the nucleus
% April 5, 2012 I. Cao-Berg Added new preprocessing method
% April 9, 2012 R.F. Murphy Save cell-nucleus height ratio and nuc bottomslice
% June 19, 2012 I. Cao-Berg Changed index name in the preprocessing for-loops that had the same
%                           index name as the main loop
% July 26, 2012 I. Cao-Berg If the number of slices in the nuclear image that have
%                           fluorescence is less than 4, then ignore such image
% May 15, 2013 I. Cao-Berg Updated method to support wildcards
% June 4, 2013 D. Sullivan Added check if masks were specified to avoid
%                          index out of bounds code when they're not
%
%%
% June 7-13 2013 D. Sullivan changed to per-cell type computations. Major
%                            changes include passing single files instead
%                            of directories and passing in a tempdirectory
%                            structure as part of param.
%%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published
% by the Free Software Foundation; either version 2 of the License,
% or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA.
%
% For additional information visit http://murphylab.web.cmu.edu or
% send email to murphy@cmu.edu

if nargin == 2
    param = [];
end

try
    downsample = param.downsample;
catch
    downsample = [5 5 1];
end

%%%%celldirlist = ml_dir([impath '/cell*']);
% preprocessingFolder = [ pwd filesep 'temp' filesep 'preprocessing' ];
% if ~exist( preprocessingFolder )
%     mkdir( preprocessingFolder );
% end
% 
% dna_image_files = ml_ls( dnaImagesDirectoryPath );
% cell_image_files = ml_ls( cellImagesDirectoryPath );
% prot_image_files = ml_ls( protImagesDirectoryPath );
% mask_image_files = ml_ls( maskImagesDirectoryPath );
% celldirlist = ml_dir([preprocessingFolder '/cell*.mat']);

%D. Sullivan 6/6/13, no longer an issue with refactored code
%D. Sullivan 6/4/13 - added check if masks were specified
%     if isempty(mask_image_file)
%         no masks specified ignore
%         masks_specified = 0;
%     elseif length(mask_image_files)~=length(cell_image_files)
%         %masks specified, but don't have proper number. No current way to
%         %register which image gets a mask and which doesn't so ignore
%         warning('Number of masks specified ~= number of cells specified, ignoring masks');
%         masks_specified = 0;
%     else
%         masks_specified = 1;
%     end
% for i = 1:length(cell_image_files)
%     param.cellnum = i;
%     dna_image_file = dna_image_files{i};
%     cell_image_file = cell_image_files{i};
%     prot_image_file = prot_image_files{i};
%     %D. Sullivan 6/4/13 - added check if masks were specified
%     if masks_specified==1
%         %mask properly specified, read mask in
%         mask_image_file = mask_image_files{i};
%     else
%         mask_image_file = [];
%     end
    
    temporaryResultsFile = [savepath filesep 'cellcodes_' num2str(currfile) '.mat'];
    nucbody = {}; cellbody = {};
    disp( ['Image:' num2str(currfile)]);
    if ~exist( temporaryResultsFile, 'file' )
%         try
            %icaoberg april 5, 2012
            %check if preprocessing results exist
            %%%%        cellfile = [ preprocessingFolder filesep 'cell' celldirlist{i}(5:end) '.mat' ];
            cellfile = [ param.preprocessingFolder filesep 'cell' num2str(currfile) '.mat' ]; %%%%
            if ~exist( cellfile, 'file' )
                [segdna,segcell] = seg_cell_and_nucleus( ...
                    dna_image_file, ...
                    cell_image_file, ...
                    prot_image_file, ...
                    mask_image_file, ...
                    downsample, param.display, param );
            else
                disp('Preprocessed results found.');
                load( cellfile );
            end
            
            %icaoberg 26/7/2012
            nnzpixel = squeeze(sum(sum(segdna)));
            
            if length(find(nnzpixel)>0) < 4
                disp(['Ignoring image ' num2str(currfile)]);
                %D. Sullivan 6/6/13 want to return instead of continue now
                %that this code does one cell at a time.
%                 continue
                return
            end
            
            
            %icaoberg april 5, 2012
            stacknumber = 15;
            [segdna, segcell] = ml_rescaleImage2Cell( segdna, segcell, stacknumber );
            
            %fills the empty slices of DNA with the closest
            z = find(sum(sum(segdna,1),2));
            segdna(:,:,1:min(z)-1)=repmat(segdna(:,:,min(z)),[1,1,min(z)-1]);
            segdna(:,:,max(z)+1:end)=repmat(segdna(:,:,max(z)),[1,1,size(segdna,3)-max(z)]);
            
            %needed for model of nuclear position
            cellnucheightratio = stacknumber/length(z);
            nucbottomslice = (min(z)-1)/stacknumber;
            
            %preprocess DNA image
            se = strel('disk',3);
            crossectionArea = [];
            %icaoberg june 19, 2012 changed index from i to j
            for j=1:1:size(segdna,3)
                mask = segdna(:,:,j);
                mask = imclose(mask,se);
                mask = imfill(mask,'holes');
                obj = ml_findmainobj2d_bw(mask);
                crossectionArea(j) = nnz(mask);
                nucbody{j}=obj;
            end
            [maxArea,equatorZ] = max(crossectionArea);
            
            %preprocess cell image
            %icaoberg june 19, 2012 changed index from i to j
            for j=1:1:size(segcell,3)
                mask = segcell(:,:,j);
                mask = imclose(mask,se);
                mask = imfill(mask,'holes');
                if ~all(mask(:) == 0)
                    obj = ml_findmainobj2d_bw(mask);
                    cellbody{j}=obj;
                end
                
            end
            
            cellcodes = {};
            for k = 1:1:length(cellbody)
                disp(['Calculating cell code on stack number: ' num2str(k) ]);
                
                method = 0; 
                if method == 0
                    %calculate projection here
                    %let nucbody become a z-projection and calculate distance according to the projection
                    mask = sum( segdna, 3 );
                    mask(find( mask ~= 0)) = 255;
                    mask = logical( mask );
                    obj = ml_findmainobj2d_bw(mask);
                    crossectionArea = nnz(mask);
                    nucbody = obj;
                    
                    if ~isempty(nucbody)
                        cellcodes{k} = ml_parsecell({},cellbody{k},nucbody,...
                            1,size(segdna(:,:,k)),...
                            {'da','nucarea','nuccenter','nucmangle','nuchitpts',...
                            'nuccontour','nucellhitpts','nucdist','nucelldist','nucecc',...
                            'cellarea','cellcenter','cellmangle','cellcontour',...
                            'cellhitpts','celldist','cellecc','nucedge','celledge'},0);
                    end
                elseif method == 1
                    %calculate projection here
                    %let nucbody become a z-projection and calculate distance according to the projection
%                     mask = sum( segdna, 3 );
%                     mask(find( mask ~= 0)) = 255;
%                     mask = logical( mask );
%                     obj = ml_findmainobj2d_bw(mask);
%                     crossectionArea = nnz(mask);
                    %nucbody = obj;
                    nuclearbody = nucbody{k};
                    
                    if ~isempty(nucbody)
                        cellcodes{k} = ml_parsecell({},cellbody{k},nuclearbody,...
                            1,size(segdna(:,:,k)),...
                            {'da','nucarea','nuccenter','nucmangle','nuchitpts',...
                            'nuccontour','nucellhitpts','nucdist','nucelldist','nucecc',...
                            'cellarea','cellcenter','cellmangle','cellcontour',...
                            'cellhitpts','celldist','cellecc','nucedge','celledge'},0);
                    end
                elseif method == 2
                                        %calculate projection here
                    %let nucbody become a z-projection and calculate distance according to the projection
                    mask = sum( segdna, 3 );
                    mask(find( mask ~= 0)) = 255;
                    mask = logical( mask );
                    obj = ml_findmainobj2d_bw(mask);
                    crossectionArea = nnz(mask);
                    %nucbody = obj;
                    nuclearbody = nucbody{k};
                    
                    if ~isempty(nucbody)
                        cellcodes{k} = ml_parsecell({},cellbody{k},nuclearbody,...
                            1,size(segdna(:,:,k)),...
                            {'da','nucarea','nuccenter','nucmangle','nuchitpts',...
                            'nuccontour','nucellhitpts','nucdist','nucelldist','nucecc',...
                            'cellarea','cellcenter','cellmangle','cellcontour',...
                            'cellhitpts','celldist','cellecc','nucedge','celledge'},0);
                    end
                end
            end
            
            %icaoberg april 5, 2012
            save(temporaryResultsFile, 'cellcodes', 'equatorZ', 'cellnucheightratio', 'nucbottomslice');
%             param = rmfield( param, 'cellnum' );
%         catch
%             disp( ['Ignoring image ' num2str(i)] );
%         end
    end
% end
end
