diff --git a/demos/3D/demo3D61/comparemodels.m b/demos/3D/demo3D61/comparemodels.m
new file mode 100644
index 0000000..2f925a6
--- /dev/null
+++ b/demos/3D/demo3D61/comparemodels.m
@@ -0,0 +1,3 @@
+options = struct('verbose',true,'includenuclear',true, ...
+    'includecell',true,'includeprot',true,'hd_threshold',10);
+slml2report('images1to5/model.mat','images11to15/model.mat',options);
diff --git a/demos/3D/demo3D61/demo3D61.m b/demos/3D/demo3D61/demo3D61.m
index 1aa1f70..7981b9a 100644
--- a/demos/3D/demo3D61/demo3D61.m
+++ b/demos/3D/demo3D61/demo3D61.m
@@ -59,6 +59,8 @@
     options.labels{length(options.labels)+1} = 'LAMP2';
     options.masks{i} = [directory filesep 'LAM_cell' num2str(i) '_mask_t1.tif'];
 end
+disp(dna_paths)
+pause
 
 %set the dimensionality
 dimensionality = '3D';
diff --git a/img2shapespace.m b/img2shapespace.m
new file mode 100644
index 0000000..5ee9440
--- /dev/null
+++ b/img2shapespace.m
@@ -0,0 +1,178 @@
+function answer = img2shapespace(varargin)
+% img2shapespace.m
+%
+% Train 3D generative SPHARM-RPDM cell shape model and creates a report
+% using the trained model.
+%
+% Input
+% -----
+% * a directory of raw or synthetic cell shape images (cell shapes file
+%       format 'LAM_cell[1-9]_ch1_t1.tif') or OME.TIFFs
+% * an options structure to override default options
+%       see http://www.cellorganizer.org/docs/2.9.0/chapters/cellorganizer_options.html#img2slml
+%       important options include
+%       .model.name for name of the model itself
+%       .model.filename for name of the file in which to store the model
+%       .model.resolution for setting image/model pixel size
+%       .downsampling for setting downsampling before model building
+%       .shape_evolution for setting whether report includes example
+%           shape evolution
+%
+% Output
+% ------
+% * a valid SLML model file
+% * a report with an embedded shape space
+
+% Soham Chakraborti (sohamc@andrew.cmu.edu)
+% R.F. Murphy (murphy@cmu.edu)
+%
+% Copyright (C) 2020 Murphy Lab
+% Computational Biology Department
+% School of Computer Science
+% Carnegie Mellon University
+%
+% 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.cbd.cmu.edu or
+% send email to murphy@cmu.edu
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if isdeployed()
+    %when method is deployed
+    if length(varargin) == 1
+        text_file = varargin{1};
+    else
+        error('Deployed function takes only 1 argument. Exiting method.');
+        return
+    end
+    
+    [filepath, name, ext] = fileparts(text_file);
+    
+    if ~exist(text_file, 'file')
+        warning('Input file does not exist. Exiting method.');
+        return
+    end
+    
+    disp(['Attempting to read input file ' text_file]);
+    fid = fopen(text_file, 'r' );
+    
+    disp('Evaluating lines from input file');
+    while ~feof(fid)
+        line = fgets(fid);
+        disp(line);
+        try
+            eval(line);
+        catch err
+            disp('Unable to parse line');
+            getReport(err)
+            return
+        end
+    end
+    fclose(fid);
+    
+    if ~exist('options', 'var')
+        options = {};
+    end
+else
+    if nargin == 2
+        dnaImagesDirectoryPath = varargin{1};
+        cellImagesDirectoryPath = varargin{2};
+        options = struct([]);
+    else
+        dnaImagesDirectoryPath = varargin{1};
+        cellImagesDirectoryPath = varargin{2};
+        options = varargin{3};
+    end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% img2slml %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+defaultoptions = struct();
+defaultoptions.verbose = false;
+defaultoptions.debug = false;
+defaultoptions.display = false;
+defaultoptions.model.id = uuidgen();
+%is this needed?  what does it do?
+defaultoptions.train = struct( 'flag', 'cell' );
+defaultoptions.cell.class = 'cell_membrane';
+defaultoptions.cell.type = 'spharm_rpdm';
+defaultoptions.labels = 'unique';
+
+options = ml_initparam( options, defaultoptions);
+
+if ~isfield(options,'model.filename')
+    options.model.filename = [options.model.name '.mat'];
+end    
+
+if ~isfield(options,'model.resolution')
+    disp ('No resolution specified for the input images;')
+    disp ('Using default of [0.1, 0.1, 0.1]');
+    % this is the resolution of the input image
+    options.model.resolution = [0.1, 0.1, 0.1];
+end
+
+if ~isfield(options,'downsampling')
+    disp ('No downsampling specified prior to model creation;')
+    disp ('Using default of no downsampling');
+    options.downsampling = [1, 1, 1];
+end
+
+% postprocess of parameterization: alignment
+options.model.spharm_rpdm.postprocess = ~false;
+% alignment method: 'major_axis' or 'foe'
+options.model.spharm_rpdm.alignment_method = 'major_axis';
+% plane of rotation: 'xy', 'xz', 'yz' or 'xyz'
+options.model.spharm_rpdm.rotation_plane = 'xy';
+
+% degree of the descriptor
+options.model.spharm_rpdm.maxDeg = 31;
+
+% latent dimension for the model
+options.model.spharm_rpdm.latent_dim = 15;
+options.model.spharm_rpdm.segmindnaImagesDirectoryPathfraction = 0.1;
+
+if exist(dnaImagesDirectoryPath,'var') && ~isempty(dnaImagesDirectoryPath)
+    options.model.name = 'SPHARM-RPDM-NucCellShapeModel';
+    options.dnaImagesDirectoryPathleus.class = 'nuc_membrane';
+    options.dnaImagesDirectoryPathleus.type = 'spharm_rpdm';
+    options.model.spharm_rpdm.components = {'nuc', 'cell'};
+else
+    options.model.spharm_rpdm.components = {'cell'};
+    options.model.name = 'SPHARM-RPDM-CellShapeModel';
+end
+
+tic; answer = img2slml( '3D', dnaImagesDirectoryPath, ...
+    cellImagesDirectoryPath, [], options ); toc,
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if exist( [pwd filesep options.model.filename] )
+    answer = slml2info( {[pwd filesep options.model.filename]}, options );
+    
+    if exist( ['./' options.model.name '-report'])
+        rmdir(['./' options.model.name '-report'],'s')
+    end
+    
+    try
+        movefile( './report' , ['./' options.model.name '-report']);
+        web(['./' options.model.name '-report' filesep 'index.html'])
+    catch
+        disp('Couldn''t rename report folder and open report in browser');
+    end
+else
+    disp ('Failed to find the model file. Exiting application.' );
+    answer = false;
+end
+end%img2shapespace
diff --git a/spharm_rpdm_v2.m b/spharm_rpdm_v2.m
index 2e2adfa..4242a80 100644
--- a/spharm_rpdm_v2.m
+++ b/spharm_rpdm_v2.m
@@ -1,594 +1,594 @@
-function answer = spharm_rpdm_v2(spharm_obj_files,options)
-%call to spharm_rpdm
-% 1/26/2021 R.F. Murphy - save indices of objects that were successfully parameterized
-% 2/8/2021 R.F. Murphy - set a default hd_thresh
-
-%IMG2SLML2
-answer = false;
-
-if ~exist( [ pwd filesep 'log'], 'dir' )
-    mkdir( 'log' );
-end
-c=clock;
-logfile = '';
-
-for i=1:1:length(c)
-    logfile = ['',logfile,num2str(c(i))]; %#ok<AGROW>
-end
-
-logfile = [ pwd filesep 'log' filesep logfile, '.log' ];
-diary( logfile )
-
-disp('Running img2slml')
-disp('Check number of input arguments')
-dnaImagesDirectoryPath = spharm_obj_files;
-cellImagesDirectoryPath = spharm_obj_files;
-
- disp('Checking and getting default parameters')
- options = get_cellorganizer_default_parameters( 'training', options );
- 
- [spharm_obj_files, spharm_obj_files, ...
-  ~,labels] = check_images_directory_paths( ...
-        dnaImagesDirectoryPath, cellImagesDirectoryPath, ...
-        [], options );
- options = clean_up_training_input_arguments( '3D', options );
- 
- if options.verbose
-    fprintf( 1, '%s\n', 'Checking the existence of temporary folder' );
-end
- 
- % %TRAIN THE GENERATIVE MODEL
-disp(' '); print_large_title('Training generative model' );
-options.dimensionality = '3D';
-
-%IMG2SLML2
-
-
-%IMG2MODEL2
-
-model = [];
-
-options = ml_initparam( options, struct( 'display', false ,...
-    'debug', false, ...
-    'verbose', false, ...
-    'paramdir', [pwd filesep 'param'], ...
-    'tempparent', [pwd filesep 'temp'],...
-    'hd_thresh', 20));
-
-check_required_options(options, {'dimensionality'})
-check_required_options(options.model, {'resolution'})
-
-disp('Setting up data');
-
-%Could be cleaned up
-[spharm_obj_files, spharm_obj_files,~, options] = setup_data(spharm_obj_files, ...
-    spharm_obj_files, {}, options);
-%Could be cleaned up 
-
-%Could be cleaned up
-disp('Setting up model options');
-options = setup_model_options(spharm_obj_files, spharm_obj_files, {}, options);
-%Could be cleaned up
-
-
-paramfiles = cell(size(spharm_obj_files));
-isdone = false(size(spharm_obj_files));
-
-
-disp(' '); print_large_title('Processing images');
-for i = 1:length(spharm_obj_files)
-
-    paramfiles{i} = [options.paramdir filesep 'param' num2str(i) '.mat'];
-    
-    tmpfile = [paramfiles{i} '.tmp'];
-    if exist(tmpfile, 'file')
-        continue
-    elseif exist(paramfiles{i}, 'file')
-        isdone(i) = true;
-        continue
-    end
-    
-    system(['touch "' tmpfile '"']);
-    [spharm_obj,options.dna_image_path] = readfileifnonblank(spharm_obj_files,i);
-    [immask,options.crop_image_path] = readfileifnonblank(options.masks,i);
-
-    
-    if ~isa(spharm_obj, 'uint8' )
-        spharm_obj = uint8(spharm_obj);
-    end
-    if ~isa(immask, 'uint8' )
-        immask = uint8(immask);
-    end
-    
-    savedir = [options.paramdir filesep 'param' num2str(i)];
-
-    try
-        [cell_params] = img2param(spharm_obj, spharm_obj, spharm_obj, immask, savedir, options);
-        %Serena 03/21 - eliminate NaN values
-        if ~isempty(cell_params)
-            if (sum(isnan(cell_params.cell.fvec(:)))>0 || sum(isnan(cell_params.cell.vertices(:)))>0 || sum(isnan(cell_params.cell.faces(:)))>0 || sum(isnan(cell_params.cell.sph_verts(:)))>0)
-                isdone(i)=false;
-            else
-                save(paramfiles{i}, '-struct', 'cell_params')
-                isdone(i) = true;
-            end
-            delete(tmpfile);
-        end
-    catch the_error
-        warning(['Unable to extract parameters for cell ' num2str(i) ...
-            ': it will be ignored.']);
-        getReport( the_error )
-        disp( 'Check the images exist or that you are using the correct options.' );
-    end
-end
-goodparamfiles = paramfiles(isdone);
-
-model=param2model(goodparamfiles,options);
-
-
-if isempty( model )
-    warning( ['Method img2model returned an empty model. Exiting method.']);
-    answer = false;
-    return
-end
-
-%
-disp('in spharm_rpdm_v2: save indices of good objects to the model file');
-%
-model.cellShapeModel.parameterization_successful = isdone;
-
-%PARSE GENERATIVE MODEL INTO SLML INSTANCE
-disp(' '); print_large_title('Parse and clean generative model' );
-
-disp(upper('Adding parameters to model structure'));
-model = parse_and_clean_generative_model( model, options );
-    
-disp(upper('Adding documentation to model structure'));
-model = add_documentation_to_model( model, options );
-    
-disp(upper('Adding parameters to model structure'));
-model = parse_and_clean_generative_model( model, options );
-    
-disp(' '); print_simple_title('Clean up workspace and environment');
-model = clean_up_and_wrap_up_model( model, options );
-answer = clean_up( options );
-
-diary off
-
-
-end
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                             HELPER METHODS                              %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-function check_required_options(options, required_fields)
-ismissing = false(size(required_fields));
-for j = 1:length(required_fields)
-    if ~isfield(options, required_fields{j})
-        ismissing(j) = true;
-    end
-end
-if any(ismissing)
-    error(['Missing options fields: ' strjoin(required_fields(ismissing), ', ')])
-end
-end
-
-function [dna_images_list, cell_images_list, protein_images_list, options] = setup_data(dna_images, cell_images, prot_images, options)
-
-options = ml_initparam(options, struct('masks', []));
-
-if ~exist(options.paramdir, 'dir')
-    mkdir(options.paramdir);
-end
-
-disp(' '); print_simple_title('Creating list of nuclear membrane images');
-dna_images_list = {};
-label = -1;
-temp_labels = {};
-if isa(dna_images, 'cell' ) && ...
-        ~any(cellfun( @(x)(isa(x,'function_handle')), dna_images ))
-    for i=1:1:length((dna_images))
-        dataset = dna_images{i};
-        temp = ml_ls( dataset );
-        if length(temp) ~= 1
-            label = label+1;
-        end
-        if isempty(dna_images_list)
-            disp('Adding first dataset to list');
-            dna_images_list = temp;
-            for j=1:1:length(dna_images_list)
-                disp(['Adding file ' dna_images_list{j}]);
-                temp_labels{length(temp_labels)+1} = num2str(label);
-            end
-        else
-            disp('Adding another dataset to list')
-            for j=1:1:length(temp)
-                dna_images_list{length(dna_images_list)+1} = ...
-                    temp{j};
-                disp(['Adding file ' temp{j}]);
-                temp_labels{length(temp_labels)+1} = num2str(label);
-            end
-        end
-    end
-    if isempty(options.labels)
-        options.labels = temp_labels;
-    end
-    clear temp_labels;
-else
-    if iscell( dna_images ) && ... 
-        all(cellfun( @(x)(isa(x,'function_handle')), dna_images ))
-        dataset = dna_images;
-        label = label + 1;
-        disp('Adding datasets to list')
-        for j=1:1:length(dataset)
-            dna_images_list{length(dna_images_list)+1} = ...
-                dataset{j};
-            disp(['Adding file function handle ' num2str(j) ' to list']);
-            temp_labels{length(temp_labels)+1} = num2str(label);
-        end
-        options.labels = temp_labels;
-        clear dataset
-    else
-        dataset = dna_images;
-        temp = ml_ls( dataset );
-        label = label + 1;
-        disp('Adding datasets to list')
-        for j=1:1:length(temp)
-            dna_images_list{length(dna_images_list)+1} = ...
-                temp{j};
-            disp(['Adding file ' temp{j}]);
-            temp_labels{length(temp_labels)+1} = num2str(label);
-        end
-        options.labels = temp_labels;
-        clear temp_labels;
-    end
-end
-
-if isempty(dna_images_list)
-    disp('List is empty or no nuclear membrane images found');
-end
-
-disp(' '); print_simple_title('Creating list of cell membrane images');
-cell_images_list = {};
-if isa(cell_images, 'cell' ) && ...
-        ~any(cellfun( @(x)(isa(x,'function_handle')), cell_images ))
-    for l=1:1:length((cell_images))
-        dataset = cell_images{l};
-        temp = ml_ls( dataset );
-        
-        if isempty(cell_images_list)
-            disp('Adding first dataset to list');
-            cell_images_list = temp;
-            for j=1:1:length(cell_images_list)
-                disp(['Adding file ' cell_images_list{j}]);
-            end
-        else
-            disp('Adding another dataset to list')
-            for j=1:1:length(temp)
-                cell_images_list{length(cell_images_list)+1} = ...
-                    temp{j};
-                disp(['Adding file ' temp{j}]);
-            end
-        end
-    end
-else
-    if iscell( cell_images ) && ...
-            all(cellfun( @(x)(isa(x,'function_handle')), cell_images ))
-        dataset = cell_images;
-        label = label + 1;
-        disp('Adding datasets to list')
-        for j=1:1:length(dataset)
-            cell_images_list{length(cell_images_list)+1} = ...
-                dataset{j};
-            disp(['Adding file function handle ' num2str(j) ' to list']);
-        end
-        clear dataset
-    else
-        dataset = cell_images;
-        temp = ml_ls( dataset );
-        disp('Adding datasets to list')
-        for j=1:1:length(temp)
-            cell_images_list{length(cell_images_list)+1} = ...
-                temp{j};
-            disp(['Adding file ' temp{j}]);
-        end
-    end
-end
-
-if isempty(cell_images_list)
-    disp('List is empty or no cell membrane images found');
-    cell_images_list = cell(size(dna_images_list));
-end
-
-disp(' '); print_simple_title('Creating list of protein pattern images');
-protein_images_list = {};
-if isa(prot_images, 'cell' ) && ...
-        ~any(cellfun( @(x)(isa(x,'function_handle')), prot_images ))
-    for i=1:1:length((prot_images))
-        dataset = prot_images{i};
-        temp = ml_ls( dataset );
-        
-        if isempty(protein_images_list)
-            disp('Adding first dataset to list');
-            protein_images_list = temp;
-            for j=1:1:length(protein_images_list)
-                disp(['Adding file ' protein_images_list{j}]);
-            end
-        else
-            disp('Adding another dataset to list')
-            for j=1:1:length(temp)
-                protein_images_list{length(protein_images_list)+1} = ...
-                    temp{j};
-                disp(['Adding file ' temp{j}]);
-            end
-        end
-    end
-else
-    if iscell( prot_images ) && ... 
-            all(cellfun( @(x)(isa(x,'function_handle')), prot_images ))
-        dataset = prot_images;
-        label = label + 1;
-        disp('Adding datasets to list')
-        for j=1:1:length(dataset)
-            protein_images_list{length(protein_images_list)+1} = ...
-                dataset{j};
-            disp(['Adding file function handle ' num2str(j) ' to list']);
-            temp_labels{length(temp_labels)+1} = num2str(label);
-        end
-        clear dataset
-    else
-        dataset = prot_images;
-        temp = ml_ls( dataset );
-        disp('Adding datasets to list')
-        for j=1:1:length(temp)
-            protein_images_list{length(protein_images_list)+1} = ...
-                temp{j};
-            disp(['Adding file ' temp{j}]);
-        end
-    end
-end
-
-if isempty(protein_images_list)
-    disp('List is empty or no protein images found');
-    protein_images_list = cell(size(dna_images_list));
-end
-
-if ischar(options.masks)
-    options.masks = ml_ls(options.masks);
-end
-if isempty(options.masks)
-    options.masks = cell(size(dna_images_list));
-end
-
-if isempty(cell_images_list) && ~isempty(options.protein.type) && ...
-        strcmp(options.protein.type, 'standardized_map_half-ellipsoid')
-    disp('Setting standardized_map_half-ellipsoid model' );
-    [cell_images_list, options] = tcell_setup_options(options);
-end
-
-disp(' ' ); disp('Saving dataset and label information')
-options.dataset.nuclear_membrane_images = dna_images_list;
-options.dataset.cell_membrane_images = cell_images_list;
-options.dataset.protein_images = protein_images_list;
-options.dataset.labels = options.labels;
-if isfield( options, 'labels' )
-    options = rmfield( options, 'labels' );
-end
-
-model.dataset = options.dataset;
-
-
-end
-
-function options = setup_model_options(dna_images, cell_images, prot_images, options)
-%this function sets the default model options for cellorganizer
-component_struct = struct('type', '', ...
-    'name', '', ...
-    'id', '');
-
-options = ml_initparam(options, struct('nucleus', []));
-% xruan 01/05/2016 change ml_initparam(options, component_struct); to ml_initparam(options.nucleus, component_struct);
-options.nucleus = ml_initparam(options.nucleus, component_struct);
-options.documentation.numimgs = numel(dna_images);
-
-if strcmpi(options.dimensionality, '2D')
-    if isempty(options.nucleus.type)
-        options.nucleus.type = 'medial axis';
-    end
-    
-    if ~all(cellfun(@isempty, cell_images)) && isempty(options.cell.type)
-        options.cell.type = 'ratio';
-    end
-    
-    if ~all(cellfun(@isempty, prot_images)) && isempty(options.protein.type)
-        options.protein.type = 'vesicle';
-    end
-    
-elseif strcmpi(options.dimensionality, '3D')
-    if isempty(options.nucleus.type)
-        options.nucleus.type = 'cylindrical_surface';
-    end
-    
-    if ~all(cellfun(@isempty, cell_images)) && isempty(options.cell.type)
-        options.cell.type = 'ratio';
-    end
-    
-    if ~all(cellfun(@isempty, prot_images)) && isempty(options.protein.type)
-        options.protein.type = 'vesicle';
-    end
-else
-    error('Unsupported dimensionality. Exiting method.')
-end
-end
-
-function [img,filename] = readfileifnonblank(files,i)
-if ~isempty(files)
-    filename = files{i};
-    
-    if ~isempty( filename )
-        if ~strcmpi(class(filename), 'function_handle')
-            disp(['Reading file ' filename] );
-        end
-    end
-    img = ml_readimage(filename);
-else
-    filename = [];
-    img = [];
-end
-end
-
-function filename = copypathifnonblank(files,i)
-if ~isempty(files)
-    filename = files{i};
-else
-    filename = [];
-end
-end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                             HELPER METHODS                              %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-%%%%%%%%%%%%%%%%%%%%%%%%% HELPER METHOD - CLEAN UP %%%%%%%%%%%%%%%%%%%%%%%%
-function answer = clean_up( param )
-%removes temp folder if debug set to true
-if ~param.debug
-    disp( 'Removing temporary folder' );
-    if exist( [ pwd filesep 'temp' ] )
-        rmdir( [ pwd filesep 'temp' ], 's'  );
-    end
-end
-
-disp('Checking if model file exists on disk')
-if exist( [param.model.filename(1:end-3) 'mat'] )
-    answer = true;
-else
-    answer = false;
-end
-
-if isfield( param.protein, 'class' ) && ...
-        strcmpi( param.protein.class, 'standardized_voxels' ) && ...
-        exist( [param.model.filename(1:end-3) 'mat'] )
-    delete( [param.model.filename(1:end-3) 'mat'] )
-end
-end%clean_up
-
-function [dnaImagesDirectoryPath, cellImagesDirectoryPath, ...
-    proteinImagesDirectoryPath, labels] = check_images_directory_paths( ...
-    dnaImagesDirectoryPath, cellImagesDirectoryPath, ...
-    proteinImagesDirectoryPath, param )
-
-if ismember( param.train.flag, {'nuclear'} )
-    if strcmpi( param.cell.class, 'framework' ) && ...
-            strcmpi( param.cell.type, 'pca' ) && ...
-            strcmpi( param.cell.class, 'framework' ) && ...
-            strcmpi( param.cell.type, 'pca' )
-        proteinImagesDirectoryPath = {};
-    else
-        cellImagesDirectoryPath = {};
-        proteinImagesDirectoryPath = {};
-    end
-elseif ismember( param.train.flag, {'cell'} )
-    if strcmpi( param.cell.class, 'framework' ) && ...
-            strcmpi( param.cell.type, 'pca' ) && ...
-            strcmpi( param.cell.class, 'framework' ) && ...
-            strcmpi( param.cell.type, 'pca' )
-        proteinImagesDirectoryPath = {};
-    end
-    
-    if strcmpi( param.cell.class, 'cell_membrane' ) && ...
-            strcmpi( param.cell.type, 'spharm_rpdm' )
-        proteinImagesDirectoryPath = {};
-        dnaImagesDirectoryPath = cellImagesDirectoryPath;
-    end
-    
-    if strcmpi( param.cell.class, 'cell_membrane' ) && ...
-            strcmpi( param.cell.type, 'ratio' )
-        warning('Unable to train a cell membrane ratio model without a nuclear membrane');
-        dnaImagesDirectoryPath = {};
-        cellImagesDirectoryPath = {};
-        proteinImagesDirectoryPath = {};
-    end
-elseif ismember( param.train.flag, {'framework'} )
-    proteinImagesDirectoryPath = {};
-elseif ismember( param.train.flag, {'protein'} )
-    if isempty( proteinImagesDirectoryPath )
-        warning('Unable to train a protein shape model without images');
-        dnaImagesDirectoryPath = {};
-        cellImagesDirectoryPath = {};
-        proteinImagesDirectoryPath = {};
-    end
-else %param.train.flag == 'all'
-    if isempty( proteinImagesDirectoryPath )
-        warning('Unable to train a protein shape model without images');
-        dnaImagesDirectoryPath = {};
-        cellImagesDirectoryPath = {};
-        proteinImagesDirectoryPath = {};
-        labels = {};
-    end
-end
-
-disp('Checking if using multiple datasets');
-check_images_directory = @(x)( (isempty(x)) || (~isempty(x) && all(cellfun(@iscell,x))) );
-if iscell( dnaImagesDirectoryPath ) && ...
-        ( check_images_directory(dnaImagesDirectoryPath) && ...
-        check_images_directory(cellImagesDirectoryPath) && ...
-        check_images_directory(proteinImagesDirectoryPath) )
-    disp('Multiple datasets found');
-    disp('Checking consistency accross datasets')
-    if numel(unique(nonzeros([length(dnaImagesDirectoryPath), ...
-            length(cellImagesDirectoryPath), ...
-            length(proteinImagesDirectoryPath)]))) == 1
-        disp('All nonempty datasets have the same length')
-    end
-    
-    labels = {};
-elseif ~isempty(dnaImagesDirectoryPath)
-    disp('Only one dataset found');
-    labels = {};
-else
-    dnaImagesDirectoryPath = {};
-    cellImagesDirectoryPath = {};
-    proteinImagesDirectoryPath = {};
-    labels = {};
-end
-end%check_images_directory_paths
-
-function imgDir = parse_ometiff_deployed(arr)
-file_array = {};
-ch_num = {};
-time_num = {};
-
-disp('Parsing cell array')
-for i = 1:(length(arr))
-    split      = strsplit(arr{i},':');
-    file_array = [file_array ml_ls(split{1})];
-    
-    %Check if given array has delimiter
-    if length(strsplit(arr{1},':')) > 1
-        ch_num     = [ch_num split{2}];
-        if length(strsplit(arr{1},':')) > 1
-            time_num   = [time_num split{3}];
-        end
-    end
-end
-
-%If no delimiter, no need to get ometiff func. handles
-if isempty(ch_num)
-    imgDir = file_array;
-    return
-end
-
-imgDir = {};
-channel_num = str2num(ch_num{1});
-for i = 1:length(file_array)
-    disp(['Parsing image ' file_array{i} ]);
-    temp = get_list_of_function_handles_from_ometiff( [file_array{i}], channel_num);
-    for j=1:1:length(temp)
-        imgDir{length(imgDir)+1} = temp{j};
-    end
-end
+function answer = spharm_rpdm_v2(spharm_obj_files,options)
+%call to spharm_rpdm
+% 1/26/2021 R.F. Murphy - save indices of objects that were successfully parameterized
+% 2/8/2021 R.F. Murphy - set a default hd_thresh
+
+%IMG2SLML2
+answer = false;
+
+if ~exist( [ pwd filesep 'log'], 'dir' )
+    mkdir( 'log' );
+end
+c=clock;
+logfile = '';
+
+for i=1:1:length(c)
+    logfile = ['',logfile,num2str(c(i))]; %#ok<AGROW>
+end
+
+logfile = [ pwd filesep 'log' filesep logfile, '.log' ];
+diary( logfile )
+
+disp('Running img2slml')
+disp('Check number of input arguments')
+dnaImagesDirectoryPath = spharm_obj_files;
+cellImagesDirectoryPath = spharm_obj_files;
+
+ disp('Checking and getting default parameters')
+ options = get_cellorganizer_default_parameters( 'training', options );
+ 
+ [spharm_obj_files, spharm_obj_files, ...
+  ~,labels] = check_images_directory_paths( ...
+        dnaImagesDirectoryPath, cellImagesDirectoryPath, ...
+        [], options );
+ options = clean_up_training_input_arguments( '3D', options );
+ 
+ if options.verbose
+    fprintf( 1, '%s\n', 'Checking the existence of temporary folder' );
+end
+ 
+ % %TRAIN THE GENERATIVE MODEL
+disp(' '); print_large_title('Training generative model' );
+options.dimensionality = '3D';
+
+%IMG2SLML2
+
+
+%IMG2MODEL2
+
+model = [];
+
+options = ml_initparam( options, struct( 'display', false ,...
+    'debug', false, ...
+    'verbose', false, ...
+    'paramdir', [pwd filesep 'param'], ...
+    'tempparent', [pwd filesep 'temp'],...
+    'hd_thresh', 20));
+
+check_required_options(options, {'dimensionality'})
+check_required_options(options.model, {'resolution'})
+
+disp('Setting up data');
+
+%Could be cleaned up
+[spharm_obj_files, spharm_obj_files,~, options] = setup_data(spharm_obj_files, ...
+    spharm_obj_files, {}, options);
+%Could be cleaned up 
+
+%Could be cleaned up
+disp('Setting up model options');
+options = setup_model_options(spharm_obj_files, spharm_obj_files, {}, options);
+%Could be cleaned up
+
+
+paramfiles = cell(size(spharm_obj_files));
+isdone = false(size(spharm_obj_files));
+
+
+disp(' '); print_large_title('Processing images');
+for i = 1:length(spharm_obj_files)
+
+    paramfiles{i} = [options.paramdir filesep 'param' num2str(i) '.mat'];
+    
+    tmpfile = [paramfiles{i} '.tmp'];
+    if exist(tmpfile, 'file')
+        continue
+    elseif exist(paramfiles{i}, 'file')
+        isdone(i) = true;
+        continue
+    end
+    
+    system(['touch "' tmpfile '"']);
+    [spharm_obj,options.dna_image_path] = readfileifnonblank(spharm_obj_files,i);
+    [immask,options.crop_image_path] = readfileifnonblank(options.masks,i);
+
+    
+    if ~isa(spharm_obj, 'uint8' )
+        spharm_obj = uint8(spharm_obj);
+    end
+    if ~isa(immask, 'uint8' )
+        immask = uint8(immask);
+    end
+    
+    savedir = [options.paramdir filesep 'param' num2str(i)];
+
+    try
+        [cell_params] = img2param(spharm_obj, spharm_obj, spharm_obj, immask, savedir, options);
+        %Serena 03/21 - eliminate NaN values
+        if ~isempty(cell_params)
+            if (sum(isnan(cell_params.cell.fvec(:)))>0 || sum(isnan(cell_params.cell.vertices(:)))>0 || sum(isnan(cell_params.cell.faces(:)))>0 || sum(isnan(cell_params.cell.sph_verts(:)))>0)
+                isdone(i)=false;
+            else
+                save(paramfiles{i}, '-struct', 'cell_params')
+                isdone(i) = true;
+            end
+            delete(tmpfile);
+        end
+    catch the_error
+        warning(['Unable to extract parameters for cell ' num2str(i) ...
+            ': it will be ignored.']);
+        getReport( the_error )
+        disp( 'Check the images exist or that you are using the correct options.' );
+    end
+end
+goodparamfiles = paramfiles(isdone);
+
+model=param2model(goodparamfiles,options);
+
+
+if isempty( model )
+    warning( ['Method img2model returned an empty model. Exiting method.']);
+    answer = false;
+    return
+end
+
+%
+disp('in spharm_rpdm_v2: save indices of good objects to the model file');
+%
+model.cellShapeModel.parameterization_successful = isdone;
+
+%PARSE GENERATIVE MODEL INTO SLML INSTANCE
+disp(' '); print_large_title('Parse and clean generative model' );
+
+disp(upper('Adding parameters to model structure'));
+model = parse_and_clean_generative_model( model, options );
+    
+disp(upper('Adding documentation to model structure'));
+model = add_documentation_to_model( model, options );
+    
+disp(upper('Adding parameters to model structure'));
+model = parse_and_clean_generative_model( model, options );
+    
+disp(' '); print_simple_title('Clean up workspace and environment');
+model = clean_up_and_wrap_up_model( model, options );
+answer = clean_up( options );
+
+diary off
+
+
+end
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                             HELPER METHODS                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function check_required_options(options, required_fields)
+ismissing = false(size(required_fields));
+for j = 1:length(required_fields)
+    if ~isfield(options, required_fields{j})
+        ismissing(j) = true;
+    end
+end
+if any(ismissing)
+    error(['Missing options fields: ' strjoin(required_fields(ismissing), ', ')])
+end
+end
+
+function [dna_images_list, cell_images_list, protein_images_list, options] = setup_data(dna_images, cell_images, prot_images, options)
+
+options = ml_initparam(options, struct('masks', []));
+
+if ~exist(options.paramdir, 'dir')
+    mkdir(options.paramdir);
+end
+
+disp(' '); print_simple_title('Creating list of nuclear membrane images');
+dna_images_list = {};
+label = -1;
+temp_labels = {};
+if isa(dna_images, 'cell' ) && ...
+        ~any(cellfun( @(x)(isa(x,'function_handle')), dna_images ))
+    for i=1:1:length((dna_images))
+        dataset = dna_images{i};
+        temp = ml_ls( dataset );
+        if length(temp) ~= 1
+            label = label+1;
+        end
+        if isempty(dna_images_list)
+            disp('Adding first dataset to list');
+            dna_images_list = temp;
+            for j=1:1:length(dna_images_list)
+                disp(['Adding file ' dna_images_list{j}]);
+                temp_labels{length(temp_labels)+1} = num2str(label);
+            end
+        else
+            disp('Adding another dataset to list')
+            for j=1:1:length(temp)
+                dna_images_list{length(dna_images_list)+1} = ...
+                    temp{j};
+                disp(['Adding file ' temp{j}]);
+                temp_labels{length(temp_labels)+1} = num2str(label);
+            end
+        end
+    end
+    if isempty(options.labels)
+        options.labels = temp_labels;
+    end
+    clear temp_labels;
+else
+    if iscell( dna_images ) && ... 
+        all(cellfun( @(x)(isa(x,'function_handle')), dna_images ))
+        dataset = dna_images;
+        label = label + 1;
+        disp('Adding datasets to list')
+        for j=1:1:length(dataset)
+            dna_images_list{length(dna_images_list)+1} = ...
+                dataset{j};
+            disp(['Adding file function handle ' num2str(j) ' to list']);
+            temp_labels{length(temp_labels)+1} = num2str(label);
+        end
+        options.labels = temp_labels;
+        clear dataset
+    else
+        dataset = dna_images;
+        temp = ml_ls( dataset );
+        label = label + 1;
+        disp('Adding datasets to list')
+        for j=1:1:length(temp)
+            dna_images_list{length(dna_images_list)+1} = ...
+                temp{j};
+            disp(['Adding file ' temp{j}]);
+            temp_labels{length(temp_labels)+1} = num2str(label);
+        end
+        options.labels = temp_labels;
+        clear temp_labels;
+    end
+end
+
+if isempty(dna_images_list)
+    disp('List is empty or no nuclear membrane images found');
+end
+
+disp(' '); print_simple_title('Creating list of cell membrane images');
+cell_images_list = {};
+if isa(cell_images, 'cell' ) && ...
+        ~any(cellfun( @(x)(isa(x,'function_handle')), cell_images ))
+    for l=1:1:length((cell_images))
+        dataset = cell_images{l};
+        temp = ml_ls( dataset );
+        
+        if isempty(cell_images_list)
+            disp('Adding first dataset to list');
+            cell_images_list = temp;
+            for j=1:1:length(cell_images_list)
+                disp(['Adding file ' cell_images_list{j}]);
+            end
+        else
+            disp('Adding another dataset to list')
+            for j=1:1:length(temp)
+                cell_images_list{length(cell_images_list)+1} = ...
+                    temp{j};
+                disp(['Adding file ' temp{j}]);
+            end
+        end
+    end
+else
+    if iscell( cell_images ) && ...
+            all(cellfun( @(x)(isa(x,'function_handle')), cell_images ))
+        dataset = cell_images;
+        label = label + 1;
+        disp('Adding datasets to list')
+        for j=1:1:length(dataset)
+            cell_images_list{length(cell_images_list)+1} = ...
+                dataset{j};
+            disp(['Adding file function handle ' num2str(j) ' to list']);
+        end
+        clear dataset
+    else
+        dataset = cell_images;
+        temp = ml_ls( dataset );
+        disp('Adding datasets to list')
+        for j=1:1:length(temp)
+            cell_images_list{length(cell_images_list)+1} = ...
+                temp{j};
+            disp(['Adding file ' temp{j}]);
+        end
+    end
+end
+
+if isempty(cell_images_list)
+    disp('List is empty or no cell membrane images found');
+    cell_images_list = cell(size(dna_images_list));
+end
+
+disp(' '); print_simple_title('Creating list of protein pattern images');
+protein_images_list = {};
+if isa(prot_images, 'cell' ) && ...
+        ~any(cellfun( @(x)(isa(x,'function_handle')), prot_images ))
+    for i=1:1:length((prot_images))
+        dataset = prot_images{i};
+        temp = ml_ls( dataset );
+        
+        if isempty(protein_images_list)
+            disp('Adding first dataset to list');
+            protein_images_list = temp;
+            for j=1:1:length(protein_images_list)
+                disp(['Adding file ' protein_images_list{j}]);
+            end
+        else
+            disp('Adding another dataset to list')
+            for j=1:1:length(temp)
+                protein_images_list{length(protein_images_list)+1} = ...
+                    temp{j};
+                disp(['Adding file ' temp{j}]);
+            end
+        end
+    end
+else
+    if iscell( prot_images ) && ... 
+            all(cellfun( @(x)(isa(x,'function_handle')), prot_images ))
+        dataset = prot_images;
+        label = label + 1;
+        disp('Adding datasets to list')
+        for j=1:1:length(dataset)
+            protein_images_list{length(protein_images_list)+1} = ...
+                dataset{j};
+            disp(['Adding file function handle ' num2str(j) ' to list']);
+            temp_labels{length(temp_labels)+1} = num2str(label);
+        end
+        clear dataset
+    else
+        dataset = prot_images;
+        temp = ml_ls( dataset );
+        disp('Adding datasets to list')
+        for j=1:1:length(temp)
+            protein_images_list{length(protein_images_list)+1} = ...
+                temp{j};
+            disp(['Adding file ' temp{j}]);
+        end
+    end
+end
+
+if isempty(protein_images_list)
+    disp('List is empty or no protein images found');
+    protein_images_list = cell(size(dna_images_list));
+end
+
+if ischar(options.masks)
+    options.masks = ml_ls(options.masks);
+end
+if isempty(options.masks)
+    options.masks = cell(size(dna_images_list));
+end
+
+if isempty(cell_images_list) && ~isempty(options.protein.type) && ...
+        strcmp(options.protein.type, 'standardized_map_half-ellipsoid')
+    disp('Setting standardized_map_half-ellipsoid model' );
+    [cell_images_list, options] = tcell_setup_options(options);
+end
+
+disp(' ' ); disp('Saving dataset and label information')
+options.dataset.nuclear_membrane_images = dna_images_list;
+options.dataset.cell_membrane_images = cell_images_list;
+options.dataset.protein_images = protein_images_list;
+options.dataset.labels = options.labels;
+if isfield( options, 'labels' )
+    options = rmfield( options, 'labels' );
+end
+
+model.dataset = options.dataset;
+
+
+end
+
+function options = setup_model_options(dna_images, cell_images, prot_images, options)
+%this function sets the default model options for cellorganizer
+component_struct = struct('type', '', ...
+    'name', '', ...
+    'id', '');
+
+options = ml_initparam(options, struct('nucleus', []));
+% xruan 01/05/2016 change ml_initparam(options, component_struct); to ml_initparam(options.nucleus, component_struct);
+options.nucleus = ml_initparam(options.nucleus, component_struct);
+options.documentation.numimgs = numel(dna_images);
+
+if strcmpi(options.dimensionality, '2D')
+    if isempty(options.nucleus.type)
+        options.nucleus.type = 'medial axis';
+    end
+    
+    if ~all(cellfun(@isempty, cell_images)) && isempty(options.cell.type)
+        options.cell.type = 'ratio';
+    end
+    
+    if ~all(cellfun(@isempty, prot_images)) && isempty(options.protein.type)
+        options.protein.type = 'vesicle';
+    end
+    
+elseif strcmpi(options.dimensionality, '3D')
+    if isempty(options.nucleus.type)
+        options.nucleus.type = 'cylindrical_surface';
+    end
+    
+    if ~all(cellfun(@isempty, cell_images)) && isempty(options.cell.type)
+        options.cell.type = 'ratio';
+    end
+    
+    if ~all(cellfun(@isempty, prot_images)) && isempty(options.protein.type)
+        options.protein.type = 'vesicle';
+    end
+else
+    error('Unsupported dimensionality. Exiting method.')
+end
+end
+
+function [img,filename] = readfileifnonblank(files,i)
+if ~isempty(files)
+    filename = files{i};
+    
+    if ~isempty( filename )
+        if ~strcmpi(class(filename), 'function_handle')
+            disp(['Reading file ' filename] );
+        end
+    end
+    img = ml_readimage(filename);
+else
+    filename = [];
+    img = [];
+end
+end
+
+function filename = copypathifnonblank(files,i)
+if ~isempty(files)
+    filename = files{i};
+else
+    filename = [];
+end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                             HELPER METHODS                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%% HELPER METHOD - CLEAN UP %%%%%%%%%%%%%%%%%%%%%%%%
+function answer = clean_up( param )
+%removes temp folder if debug set to true
+if ~param.debug
+    disp( 'Removing temporary folder' );
+    if exist( [ pwd filesep 'temp' ] )
+        rmdir( [ pwd filesep 'temp' ], 's'  );
+    end
+end
+
+disp('Checking if model file exists on disk')
+if exist( [param.model.filename(1:end-3) 'mat'] )
+    answer = true;
+else
+    answer = false;
+end
+
+if isfield( param.protein, 'class' ) && ...
+        strcmpi( param.protein.class, 'standardized_voxels' ) && ...
+        exist( [param.model.filename(1:end-3) 'mat'] )
+    delete( [param.model.filename(1:end-3) 'mat'] )
+end
+end%clean_up
+
+function [dnaImagesDirectoryPath, cellImagesDirectoryPath, ...
+    proteinImagesDirectoryPath, labels] = check_images_directory_paths( ...
+    dnaImagesDirectoryPath, cellImagesDirectoryPath, ...
+    proteinImagesDirectoryPath, param )
+
+if ismember( param.train.flag, {'nuclear'} )
+    if strcmpi( param.cell.class, 'framework' ) && ...
+            strcmpi( param.cell.type, 'pca' ) && ...
+            strcmpi( param.cell.class, 'framework' ) && ...
+            strcmpi( param.cell.type, 'pca' )
+        proteinImagesDirectoryPath = {};
+    else
+        cellImagesDirectoryPath = {};
+        proteinImagesDirectoryPath = {};
+    end
+elseif ismember( param.train.flag, {'cell'} )
+    if strcmpi( param.cell.class, 'framework' ) && ...
+            strcmpi( param.cell.type, 'pca' ) && ...
+            strcmpi( param.cell.class, 'framework' ) && ...
+            strcmpi( param.cell.type, 'pca' )
+        proteinImagesDirectoryPath = {};
+    end
+    
+    if strcmpi( param.cell.class, 'cell_membrane' ) && ...
+            strcmpi( param.cell.type, 'spharm_rpdm' )
+        proteinImagesDirectoryPath = {};
+        dnaImagesDirectoryPath = cellImagesDirectoryPath;
+    end
+    
+    if strcmpi( param.cell.class, 'cell_membrane' ) && ...
+            strcmpi( param.cell.type, 'ratio' )
+        warning('Unable to train a cell membrane ratio model without a nuclear membrane');
+        dnaImagesDirectoryPath = {};
+        cellImagesDirectoryPath = {};
+        proteinImagesDirectoryPath = {};
+    end
+elseif ismember( param.train.flag, {'framework'} )
+    proteinImagesDirectoryPath = {};
+elseif ismember( param.train.flag, {'protein'} )
+    if isempty( proteinImagesDirectoryPath )
+        warning('Unable to train a protein shape model without images');
+        dnaImagesDirectoryPath = {};
+        cellImagesDirectoryPath = {};
+        proteinImagesDirectoryPath = {};
+    end
+else %param.train.flag == 'all'
+    if isempty( proteinImagesDirectoryPath )
+        warning('Unable to train a protein shape model without images');
+        dnaImagesDirectoryPath = {};
+        cellImagesDirectoryPath = {};
+        proteinImagesDirectoryPath = {};
+        labels = {};
+    end
+end
+
+disp('Checking if using multiple datasets');
+check_images_directory = @(x)( (isempty(x)) || (~isempty(x) && all(cellfun(@iscell,x))) );
+if iscell( dnaImagesDirectoryPath ) && ...
+        ( check_images_directory(dnaImagesDirectoryPath) && ...
+        check_images_directory(cellImagesDirectoryPath) && ...
+        check_images_directory(proteinImagesDirectoryPath) )
+    disp('Multiple datasets found');
+    disp('Checking consistency accross datasets')
+    if numel(unique(nonzeros([length(dnaImagesDirectoryPath), ...
+            length(cellImagesDirectoryPath), ...
+            length(proteinImagesDirectoryPath)]))) == 1
+        disp('All nonempty datasets have the same length')
+    end
+    
+    labels = {};
+elseif ~isempty(dnaImagesDirectoryPath)
+    disp('Only one dataset found');
+    labels = {};
+else
+    dnaImagesDirectoryPath = {};
+    cellImagesDirectoryPath = {};
+    proteinImagesDirectoryPath = {};
+    labels = {};
+end
+end%check_images_directory_paths
+
+function imgDir = parse_ometiff_deployed(arr)
+file_array = {};
+ch_num = {};
+time_num = {};
+
+disp('Parsing cell array')
+for i = 1:(length(arr))
+    split      = strsplit(arr{i},':');
+    file_array = [file_array ml_ls(split{1})];
+    
+    %Check if given array has delimiter
+    if length(strsplit(arr{1},':')) > 1
+        ch_num     = [ch_num split{2}];
+        if length(strsplit(arr{1},':')) > 1
+            time_num   = [time_num split{3}];
+        end
+    end
+end
+
+%If no delimiter, no need to get ometiff func. handles
+if isempty(ch_num)
+    imgDir = file_array;
+    return
+end
+
+imgDir = {};
+channel_num = str2num(ch_num{1});
+for i = 1:length(file_array)
+    disp(['Parsing image ' file_array{i} ]);
+    temp = get_list_of_function_handles_from_ometiff( [file_array{i}], channel_num);
+    for j=1:1:length(temp)
+        imgDir{length(imgDir)+1} = temp{j};
+    end
+end
 end
\ No newline at end of file
diff --git a/utilities/.DS_Store b/utilities/.DS_Store
deleted file mode 100644
index a9c590d..0000000
Binary files a/utilities/.DS_Store and /dev/null differ
diff --git a/utilities/._consolidateshapespacemodels.m b/utilities/._consolidateshapespacemodels.m
index 8d84442..0323dc5 100644
Binary files a/utilities/._consolidateshapespacemodels.m and b/utilities/._consolidateshapespacemodels.m differ
diff --git a/utilities/._model2info_v2.m b/utilities/._model2info_v2.m
index 384a82a..28f4060 100644
Binary files a/utilities/._model2info_v2.m and b/utilities/._model2info_v2.m differ
diff --git a/utilities/._models2report_v2.m b/utilities/._models2report_v2.m
index 8d84442..d34713a 100644
Binary files a/utilities/._models2report_v2.m and b/utilities/._models2report_v2.m differ
diff --git a/utilities/3D/.DS_Store b/utilities/3D/.DS_Store
deleted file mode 100644
index dddfd22..0000000
Binary files a/utilities/3D/.DS_Store and /dev/null differ
diff --git a/utilities/3D/._spharm_obj_percell_3D.m b/utilities/3D/._spharm_obj_percell_3D.m
index 1e08a93..23aafd9 100644
Binary files a/utilities/3D/._spharm_obj_percell_3D.m and b/utilities/3D/._spharm_obj_percell_3D.m differ
diff --git a/utilities/3D/spharm-ppm/._spharm_obj_model.m b/utilities/3D/spharm-ppm/._spharm_obj_model.m
index c541516..a2b5025 100644
Binary files a/utilities/3D/spharm-ppm/._spharm_obj_model.m and b/utilities/3D/spharm-ppm/._spharm_obj_model.m differ
diff --git a/utilities/3D/spharm-ppm/spharm_obj_model.m b/utilities/3D/spharm-ppm/spharm_obj_model.m
index 99f4c47..72ea172 100644
--- a/utilities/3D/spharm-ppm/spharm_obj_model.m
+++ b/utilities/3D/spharm-ppm/spharm_obj_model.m
@@ -2,6 +2,7 @@
 % 1/26/2021 R.F. Murphy - save normalized nuclear distances into model
 % (only for objects that were successfully parameterized)
 % 3/2/2021 R.F. Murphy - remove unused ppm code
+% 4/24/2021 R.F. Murphy - merge in Serena's changes
 
     options_spharm = options.options_spharm;
     spharm_obj_dir = [pwd filesep 'spharm_input'];
@@ -32,6 +33,7 @@
 anglestheta = [];
 anglesphi = [];
 mappos_x = [];
+distcodes=[];
 paramfiles = ml_ls([paramdir filesep '*.mat']);
 for i=1:length(paramfiles)
     pp = load(paramfiles{i});
@@ -40,10 +42,14 @@
     normdists = [normdists; pp.prot.seg.normdists];
     anglestheta = [anglestheta; pp.prot.seg.angles.theta];
     anglesphi = [anglesphi; pp.prot.seg.angles.phi];
-    mappos_x = [mappos_x; pp.prot.seg.mappos_x];    
+    mappos_x = [mappos_x; pp.prot.seg.mappos_x];
+    distcodes=[distcodes;pp.prot.seg.distcodes];
 end
 spatial.normdists = normdists(isdone)';
 spatial.anglestheta = anglestheta(isdone)';
 spatial.anglesphi = anglesphi(isdone)';
 spatial.mappos_x = mappos_x(isdone,:);
+spatial.distcodes=distcodes(isdone,:);
+beta=ml_logreg(spatial.mappos_x(:,2:6),spatial.distcodes(:,3));
+spatial.beta=beta;
 end
\ No newline at end of file
diff --git a/utilities/3D/spharm_obj_percell_3D.m b/utilities/3D/spharm_obj_percell_3D.m
index 17b8ceb..c4189d8 100644
--- a/utilities/3D/spharm_obj_percell_3D.m
+++ b/utilities/3D/spharm_obj_percell_3D.m
@@ -6,6 +6,8 @@
 % 1/27/2021 R.F. Murphy - also save normdists and angles; don't suppress
 %                           objects with zero distcodes
 % 2/1/2021 R.F. Murphy - mask protein image before object finding
+% 4/24/2021 R.F. Murphy - merge in Serena's fixes for cell and nuclear
+% segmentation and normdists
 
 min_obj_size = options.min_obj_size;
 max_obj_size = options.max_obj_size;
@@ -58,8 +60,19 @@
 options = ml_initparam(options, struct('use_geodesic_distance', false));
 
 
-nucEdge = bwperim(seg.nuc,18);
-cellEdge = bwperim(seg.cell,18);
+%Serena 03/21
+se = strel('line', 3,5);
+seg.cell=imfill(seg.cell,'holes');
+cellerode = imerode(seg.cell, se);
+cellEdge = logical(abs(imsubtract(cellerode, seg.cell)));
+cellEdge=imfill(cellEdge,'holes');
+
+seg.nuc=imfill(seg.nuc,'holes');
+nucerode=imerode(seg.nuc,se);
+nucEdge= logical(abs(imsubtract(nucerode, seg.nuc)));
+
+%nucEdge = bwperim(seg.nuc,18);
+%cellEdge = bwperim(seg.cell,18);
 obj_center_image = zeros(size(cellEdge));
             
 for p=1:size(centers)
@@ -79,7 +92,16 @@
 %distcodes(nullidx,:)=[];
 %angles.theta(nullidx,:)=[];
 %angles.phi(nullidx,:)=[];
-normdists = distcodes(:,1)./sum(abs(distcodes(:,1:2)),2);
+%normdists = distcodes(:,1)./sum(abs(distcodes(:,1:2)),2);
+
+for i=1:length(distcodes)
+    %Serena 03/21
+    if and(distcodes(i,1)==0,distcodes(i,2)==0)
+        normdists(i,1)=0;
+        continue;
+    end
+    normdists(i,1)=distcodes(i,1)/sum(abs(distcodes(i,1:2)),2);
+end
 
 disp('Mapping positions')
 x = ml_mappos([normdists angles.theta angles.phi]);
diff --git a/utilities/3D/spharm_rpdm/SPHARM-RPDM-package/.DS_Store b/utilities/3D/spharm_rpdm/SPHARM-RPDM-package/.DS_Store
deleted file mode 100644
index 396a50c..0000000
Binary files a/utilities/3D/spharm_rpdm/SPHARM-RPDM-package/.DS_Store and /dev/null differ
diff --git a/utilities/consolidateshapespacemodels.m b/utilities/consolidateshapespacemodels.m
index 0eda6af..07115fc 100644
--- a/utilities/consolidateshapespacemodels.m
+++ b/utilities/consolidateshapespacemodels.m
@@ -1,35 +1,35 @@
-function consolidated_model = consolidateshapespacemodels(model1,model2)
-% input models should be the cellShapeModels substruct
-
-X1=model1.X;
-coeff1=model1.coeff;
-score1=model1.train_score; 
-
-X2=model2.X;
-coeff2=model2.coeff;
-score2=model2.train_score; 
-
-consolidateX=[X1;X2];
-
-%default value in cellorganizer
-latent_dim=100;
-
-[scales,coeff,score,latent,tsquared,explained,mu,train_score_consolidate,train_explained,train_coeff_consolidate] = CalculatePCA(consolidateX, latent_dim);
-spharm1=model1.all_spharm_descriptors;
-spharm2=model2.all_spharm_descriptors;
-
-%Create new model with required features 
-consolidated_model.coeff=train_coeff_consolidate;
-consolidated_model.X=consolidateX;
-consolidated_model.train_score=train_score_consolidate;
-consolidated_model.numimgs1=model1.numimgs;
-consolidated_model.numimgs2=model2.numimgs;
-consolidated_model.numimgs=model1.numimgs+model2.numimgs;
-consolidated_model.all_spharm_descriptors=cat(3,spharm1,spharm2);
-consolidated_model.hausdorff_distances=horzcat(model1.hausdorff_distances,model2.hausdorff_distances);
-if ismember(model1.components,model2.components)
-    consolidated_model.components=model2.components;
-else
-    consolidated_model.components={};
-    consolidated_model.components={model1.components,model2.components};
-end
+function consolidated_model = consolidateshapespacemodels(model1,model2)
+% input models should be the cellShapeModels substruct
+
+X1=model1.X;
+coeff1=model1.coeff;
+score1=model1.train_score; 
+
+X2=model2.X;
+coeff2=model2.coeff;
+score2=model2.train_score; 
+
+consolidateX=[X1;X2];
+
+%default value in cellorganizer
+latent_dim=100;
+
+[scales,coeff,score,latent,tsquared,explained,mu,train_score_consolidate,train_explained,train_coeff_consolidate] = CalculatePCA(consolidateX, latent_dim);
+spharm1=model1.all_spharm_descriptors;
+spharm2=model2.all_spharm_descriptors;
+
+%Create new model with required features 
+consolidated_model.coeff=train_coeff_consolidate;
+consolidated_model.X=consolidateX;
+consolidated_model.train_score=train_score_consolidate;
+consolidated_model.numimgs1=model1.numimgs;
+consolidated_model.numimgs2=model2.numimgs;
+consolidated_model.numimgs=model1.numimgs+model2.numimgs;
+consolidated_model.all_spharm_descriptors=cat(3,spharm1,spharm2);
+consolidated_model.hausdorff_distances=horzcat(model1.hausdorff_distances,model2.hausdorff_distances);
+if ismember(model1.components,model2.components)
+    consolidated_model.components=model2.components;
+else
+    consolidated_model.components={};
+    consolidated_model.components={model1.components,model2.components};
+end
diff --git a/utilities/model2info_v2.m b/utilities/model2info_v2.m
index 00f43b5..233f4fa 100644
--- a/utilities/model2info_v2.m
+++ b/utilities/model2info_v2.m
@@ -33,6 +33,7 @@ function model2info_v2( model, fileID, options )
 % 2020/10/16 R. F. Murphy added printout of Hausdorff distances for
 %                       SPHARM-RPDM models
 % 2021/03/05 R.F. Murphy modified to use new name for hd values
+% 2021/04/24 R.F. Murphy merged in Serena's show_spatial_distribution
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODEL.NAME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 header2html( fileID, 'Model name');
@@ -225,6 +226,10 @@ function model2info_v2( model, fileID, options )
     close all
     
     movefile( '*.png', './html' );
+    
+    %Spatial Distribution
+    show_spatial_distribution(model.spatial,fileID)
+
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%% IS TCell MODEL? %%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/utilities/models2report_v2.m b/utilities/models2report_v2.m
index 4c49463..b60f2e6 100644
--- a/utilities/models2report_v2.m
+++ b/utilities/models2report_v2.m
@@ -320,13 +320,17 @@ function models2report_v2(models, param, classlabels, fileID)
 %SPHARM_OBJ_MODEL COMPARISON
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     if strcmpi(models{1}.proteinModel.type,'spharm_obj')
-%         header2html(fileID, 'Spharm Model Comparison'); 
-%         total_variation= ppm_total_variation(models);
-%         text2html(fileID, ['Total_Variation=''%s'';\n', total_variation]); 
         
         models_ShapeModel{1} = models{1}.proteinModel.spharm_obj_model.cellShapeModel;
         models_ShapeModel{2} = models{2}.proteinModel.spharm_obj_model.cellShapeModel;
-        compare_shape_space_spharm_obj(models_ShapeModel,fileID,param);
+        consolidated_model=compare_shape_space_spharm_obj(models_ShapeModel,fileID,param);
+        header2html(fileID, 'Spharm Model Comparison');
+        total_variation=consolidated_model.total_variation;
+        text2html(fileID, ['Total_Variation=''%s'';\n', total_variation]);
+        %spatial model comparision
+        models_SpatialModel{1}=models{1}.proteinModel.spharm_obj_model.spatial;
+        models_SpatialModel{2}=models{2}.proteinModel.spharm_obj_model.spatial;
+        compare_spatial_model_spharm_obj(models_SpatialModel,fileID,param);
     end
 end
 
diff --git a/utilities/opentstool/.DS_Store b/utilities/opentstool/.DS_Store
deleted file mode 100644
index 5a5e04f..0000000
Binary files a/utilities/opentstool/.DS_Store and /dev/null differ
diff --git a/utilities/show_spatial_distribution.m b/utilities/show_spatial_distribution.m
new file mode 100644
index 0000000..6c92854
--- /dev/null
+++ b/utilities/show_spatial_distribution.m
@@ -0,0 +1,49 @@
+function [ output_args ] = show_spatial_distribution( model,fileID)
+%Author : Serena Abraham
+
+
+
+    f = figure('visible','off');
+    spatial=model;
+    [x1,y1,z1]=sph2cart(spatial.anglestheta,spatial.anglesphi,spatial.normdists);
+    [X,Y,Z] = sphere;
+    scatter3(X(:),Y(:),Z(:),'filled','MarkerFaceColor',[0.7 0.7 0.7]);
+    hold on;
+    scatter3(x1,y1,z1,'filled','b');
+    legend('cell boundary','model');
+    title('Spatial Distribution of objects');
+    hold off;
+    saveas( f, 'show_spatial_distribution.png', 'png');
+    I = imread( 'show_spatial_distribution.png' );
+    I = imresize( I, 0.50 );
+    imwrite( I, 'show_spatial_distribution_thumbnail.png' );
+    img2html(fileID,'show_spatial_distribution.png','show_spatial_distribution_thumbnail.png','Spatial Distribution across both models');
+    
+    
+    
+    %spatial positions fit
+    %Generate all possible angles
+    angle_range=linspace(0,6.28319,36);
+    three_angle_range=nchoosek(angle_range,2);
+    normdist_range=linspace(0.1,1,10);
+    append_range=normdist_range(1).*ones(length(three_angle_range),1);
+    final_range=horzcat(append_range,three_angle_range);
+    for i=2:length(normdist_range)
+        append_range=normdist_range(i).*ones(length(three_angle_range),1);
+        final_range=vertcat(final_range,horzcat(append_range,three_angle_range));
+    end
+    f1 = figure('visible','on');
+    final_range2=ml_mappos(final_range);
+    y_a=1./(1+exp(-final_range2(:,2:6)*spatial.beta));
+    [x_1,y_1,z_1]=sph2cart(final_range(:,2),final_range(:,3),final_range(:,1));
+    scatter3(x_1,y_1,z_1,'filled','CData',y_a);
+    title('Spatial Distribution of Model  fitted to a Logistic Regression Model');
+    colorbar;
+    saveas( f1, 'spatial_model_fit.png', 'png');
+    I = imread( 'spatial_model_fit.png' );
+    I = imresize( I, 0.50 );
+    imwrite( I, 'spatial_model_fit_thumbnail.png' );
+    img2html(fileID,'spatial_model_fit.png','spatial_model_fit.png','Spatial Distribution fitted to a Logistic Regression Model');
+
+
+end
\ No newline at end of file
diff --git a/utilities/total_variation_joint_pca.m b/utilities/total_variation_joint_pca.m
new file mode 100644
index 0000000..93bc52e
--- /dev/null
+++ b/utilities/total_variation_joint_pca.m
@@ -0,0 +1,37 @@
+function [vol_norm]= total_variation_joint_pca(X1,X2,coeff)
+
+score1=X1*coeff;
+score2=X2*coeff;
+vol_norm=0;
+if isequal(score1,score2)
+    vol_norm=1;
+else 
+ [m1,s1]=normfit(score1(:,1:2));
+ [m2,s2]=normfit(score2(:,1:2));
+ x1=mvnrnd(m1,s1,100);
+ x2=mvnrnd(m2,s2,100);
+ y1=mvnpdf(x1,m1);
+ y2=mvnpdf(x2,m2);
+ xi = linspace(min(min(x1(:,1),x2(:,1))),max(max(x1(:,1),x2(:,1))),100);
+ yi = linspace(min(min(x1(:,2),x2(:,2))),max(max(x1(:,2),x2(:,2))),100);
+ dx = xi(2) - xi(1);
+ dy = yi(2) - yi(1);
+
+
+ [X, Y] = meshgrid(xi,yi);
+ Z1=griddata(x1(:,1),x1(:,2),y1,X,Y);
+ dv = Z1(~isnan(Z1))*dx*dy;
+ vol_1=sum(dv);
+ Z1_norm=Z1/vol_1;
+
+ Z2=griddata(x2(:,1),x2(:,2),y2,X,Y);
+ dv = Z2(~isnan(Z2))*dx*dy;
+ vol_2=sum(dv);
+ Z2_norm=Z2/vol_2;
+
+ Z_shared=min(Z1_norm,Z2_norm);
+ dv = Z_shared(~isnan(Z_shared))*dx*dy;
+
+ vol_norm=sum(dv);
+end
+end
\ No newline at end of file