From fbb28dcc251a75de2052e9e630b6b03b4b671a5a Mon Sep 17 00:00:00 2001 From: Ted Zhang Date: Wed, 28 Apr 2021 12:51:50 -0400 Subject: [PATCH] cellorganizer v2.9.1 --- demos/3D/demo3D61/comparemodels.m | 3 + demos/3D/demo3D61/demo3D61.m | 2 + img2shapespace.m | 178 +++ spharm_rpdm_v2.m | 1186 ++++++++--------- utilities/.DS_Store | Bin 30724 -> 0 bytes utilities/._consolidateshapespacemodels.m | Bin 212 -> 312 bytes utilities/._model2info_v2.m | Bin 268 -> 268 bytes utilities/._models2report_v2.m | Bin 212 -> 310 bytes utilities/3D/.DS_Store | Bin 14340 -> 0 bytes utilities/3D/._spharm_obj_percell_3D.m | Bin 268 -> 268 bytes utilities/3D/spharm-ppm/._spharm_obj_model.m | Bin 786 -> 786 bytes utilities/3D/spharm-ppm/spharm_obj_model.m | 8 +- utilities/3D/spharm_obj_percell_3D.m | 28 +- .../spharm_rpdm/SPHARM-RPDM-package/.DS_Store | Bin 6148 -> 0 bytes utilities/consolidateshapespacemodels.m | 70 +- utilities/model2info_v2.m | 5 + utilities/models2report_v2.m | 12 +- utilities/opentstool/.DS_Store | Bin 6148 -> 0 bytes utilities/show_spatial_distribution.m | 49 + utilities/total_variation_joint_pca.m | 37 + 20 files changed, 942 insertions(+), 636 deletions(-) create mode 100644 demos/3D/demo3D61/comparemodels.m create mode 100644 img2shapespace.m delete mode 100644 utilities/.DS_Store delete mode 100644 utilities/3D/.DS_Store delete mode 100644 utilities/3D/spharm_rpdm/SPHARM-RPDM-package/.DS_Store delete mode 100644 utilities/opentstool/.DS_Store create mode 100644 utilities/show_spatial_distribution.m create mode 100644 utilities/total_variation_joint_pca.m 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 -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 +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 a9c590d48662565404e6b979701cce01c4329174..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 30724 zcmeHQ3vgUlc|K<)OS{%eTG>wQ2sn%FvcVNIOkSnbIdSg+L#)G?UO4+R`b36sC_b?M!KAC>h!gls?k$J7@2` z=fC&dE4h?51#Pn0$Nm5Rod0?L|37E1l~Q{5!t{_*b4sZ^o_sA8mC6;A8o=u|Jc9?X z^O<ip{N)o&b|sh=KePBa%9yC>%F=cmIH&6&y8`y=&w zHGDrYc~7N$dc3r5PrcThsMN~y2Dvg_2FM$ZRch1q)4S^nwdsbTGiER3^M!n2ymZ<5 z^V>%T2W}i1y7{5Of%8KeUw_rmZ3j-wo;`Q|Ef0O)qnIro8dhl@ zy{*^a<^6a=jti=yX4SAdrRMQ&T5VOc2ded{7QZm$_RHW41239~LGHlcUf|RLT~<}# zR#Z)$#OHCfh))R$>Q+#Bk3#_-qPj8O1n{PTHw(-LB+ePs)T!C}^xoF^aJCCu4t(Ot zUBIB@wF9mSsHwS`p|Pp@ES(h(OtUr9j*>G@4=9%5&~>~vEp2rKpX#8|sGf2uSUM|z zS0u0^f$u>Hw4!9b>wA#*6_!^dup)tM2{1kk;K|p~*6RiPYC;G8FoooGrjX`=_ee_2 zcLS~0(bnq)`)W{uV^HyGQFW`t&=wuv4-g!$qpjBq_SK>*ZPAr(jjUT0iZ-JIZ@_>n zJ!)T9{;o)1MFL3)?8jGm)EBl|7)K?7Dp1x!ZMP}a@lL10NWjagqVIJ+nzQx2PJN@I z`X)BEtb@|wWG%FfoUH~QbEtA`|6wbtvx+(;d0}0b;Kjv>*MUDAgYRk8Yzt}%)zvV@ zn#X?=>TI$CZ2y-SZf^{C5p^1=Q1;&Xv9st%U1CoGYZlyc{IUo8Lukpk*xeZWEVMhU zCLy7+sX@!WarOEbW+}v+1AJY*sWs~g-q(^N>mQ2YM(ZG%qd_}Uvs3$SH%27iaaNfSD#Cb8%- zLi)uJi+oTD)J0RxNb79-<>lfDkL&E7gY~D5k>Yt7@NOp}^uSbj;feVHdSU^*RUl7x zNoLVyNywFcC?wR^!XOSUJJQKGzS--c0MA@Dg%K|W_@|(-CSwi<)h={V5`3`B)A^&0 zm*zo(@qj%Y)7Oz2R`knEc%fOf>&GF@1kDtUdf>d8Ouz{Z!oBdgyv0Wiw#tY~t3L*g zV6%{A+4$lDt~AvF(-EUJ^QtEfnuy}^yc%%l@U%UtNxYs>!%|KriY$pEWVQCj=P_?D zZC74DaZTk{@3p(47XL^s?qwZ7(rG`wb^}-2`a&NlaaMYq9 zISOM-y9%#I(g$9UC4|=nM>%7udGud1)=?9S-FdDv@D0dd3iew@7uH=%E3kL(x@Kat zOQ%nr2IM5x!||9c2OZ&UQer0NT_gn|OLcSvGFisBDNPu#`^X`BldT>mHk9Rq+t)?_uAtlAaSOw}{l2ld z$8Hk04FpVoDcF=RtnMmy7mwa^^w{y@cxlZd_S+B5H7oU6O`}=t6CkaoGatI-p15?y;?n7Z&XO;-s8>knI^9# zE*-zGQJ-(_F~p8nK;J@@78>Qd%OF}=Ebj&BMoAZbt!MqlE3Uf!rtLdly=VWy!_V}V z`qr#nS1jI#9sgRha8CjxAOimC0)PNNr-a+?OjA3I=P# zW&Qof7bYujUZ^zBmNw_}=G#q|_ZQzVjlu4$Pm|ho>7#iIf92)<{rB_Ap7}EVPPY)Q zzNWu;qyer<6D`0S28+cb58_BdsnHs0+fcFpXmw(;TrEu(>^QI3QS3jqXa!rkqxI!0 zZ!Pv8#}550SZVa>!PeKW8SU?XfFmEAXfzMam1~&QAGU_yi}~KmYrCeWF{i>aYf62+ z{RjH*$d~>+e_vx}V!nLnG>%II-yh27jUMqmsog$ZY1YGwW4u_Z)D|l^l>wF6!x@Fm zfthU{z1F@p3~wIYWs%u%u0EokP#;sDQlC{{LiO;x`nvjt z`j+}P_3wI>?$IUPr`PD0>1*}%`UX9yU!iyD-FlBcs9&$g^bvhjKcJ^{S)bIW^qgMQ zXY^TpPXC~OhkiuAQ~#KLuYO#AK>vdNp#F9JyZV#*_w|$dGy3!T3;J38oc^-@Tm5(X zYx?x|guee6h!R$t9PrAE4{p{Hk4{m07b4EYY=ZqhsL;#;A=K#jNoe_|nTZZYsV! zaJ9G#9)&lQ&*W*+rA`ZdWpPU$-cVX)KVDP3E)PE_b=Z%C8@J}+d8M5BSRC3oQh9-_*aYe@8#9|5SfY|7FZjX_(GK z{}nzMLv>ftfBx$Wzj*i?PwtAYMfrLBn!Ad&J#T7aLvL|xzE-W*PJsrUWGy;a)2LP# zkl{Af-6qrJnVm)4Q%hQo}|H>^&sAFX&(kZc4zK0pyf=5EZ4B-8l4EF_mQkpWP0Dze(Un=hqsb zMk?&Ivd3pFU;?vC>Ef~fAt(=kvbzVt(`-R+pYMaZ%D(57d5p%6J9{RHnYW`GAn5m? z>tyNMj(IT0gP38=23UhHAJ?@9=#Xh==2+4;!g_QGoXB56Kd4eOaa@~G>Rmv?qD&c6 zKlD^~?QADDl5`|zSGNi2ve%$K=;ZTjT9j>%q6t}#sr}}h;GHHi%p^6$N?y>&hctK+ z@?@gtiK;mpUTWetrnw24Qu-GHNqXR1iKoL!DWAZQ=owSWbLd^O{-+!h@?$?*@X`w$ z+NIS;VB_rfc)Z%gG*|jP0-Ldz-PkjjLto@UH`ZX946!G;3P%ko*01ztEn@ z&e)kXtBh;2U7^osVR3zNvATg(H+F2-PoA4^Wl_?m7Wcm=z=DG2K|`%EdxUu z;svsoK|f}a9Waaqgk?`4=kh5`VI7Q?>0kxWT;RF9kDG-c*Z`kT+fEg9R=1Xp;0TZB z*bA4hrE#P#zOWOVdudB8sHn*qN|G@_GQns}7SYpMsR&;L^3Ys0yA0tl*{X zo3D@4l01arOVEhfP7t~TjF@wKiF57>)B9nB2&-f7vMy%d?O;)ZI*$FHa)D1zup-Ii3OZv?)rM>~P~h!5x=Uf!PcCEG{OH zfyIT2TP16wI5pWl2~3_UII&4i0xJtj5Jt?H&6_kqE6a23mPE- z4~FEr^9gY6fs&ktM+Cf(E_Ydwf+ID6fz>(_w`#Bl9)96Q={TP8=$5b(vpoKsvpjD0 z(5WG!+O2w)kh07?X%btOm6Rf_*pq&n$m@>M=NW~C6bo&C0*_KgyA(PM1;%M?^-A5D zSeptfwAT4BJN|-}tFIsVz*Ml65`&2}l?Eep)cZ;5=xQ_PCS(;vC9W9sOCDG`XLx5* zz&+eE4yjP;6X297mfg@!=+VStmuD+sjun-q9M>Ms^fa1S<_HOw!r8hMc)-Yec+f;e zcLouDbRcUS5sx$Q#X3eu7oT2Dj(4lOP#C)Ikl<%Eb-ziy7o6N*Y z7q{5L9fRFZVg*QP^)QtcGw#7W51OMo<7%H7w+RZg=Y|@RHfrLFJ7y-XuvxY50oYO| zylcaN4c5&ia3Y&tE{?Eqbr@(S?%U(I{EMGeg8N}(<&()(iqS{>$TM#mSPis6yf(r( zBH)W5ynB=_=^cn30B%(c_X3MLY+?j^tP2vBQeqaA z4KF06zw&#KRDxzqT7iXo;Dp3d;e^z3@1i9a&5E|fHIYwA1r>Bjxxwz4Y zJnzf>9M3$B^SD^{ck~&%ImSF)jLvY5i&zq4XT`V_3i-vLg=?zcdLwK%29e^J2`H~y z#-M~#s^{ob5x?f7TZ&z~{N@e6foCOjXuQdVw|cN;0ZTe|g%eoxg0_v07j$D3LJuSG zF=?=852rksSi%Hmhryfy7Ls)tD(M61uu>Om{sz(!v{Ez?gYoJ?z(_3&Pz#J#Ho)8R zlf7*=P?(}7WW^;^aIq%D4w6xx1a=eEBn>RBg2slBMgu7Z?*YMxRYeS1R@p=E6ZC>P zaYybX{wo`9n4YsJ;rtSFqW61gun23Kv&r*j&oGKS2^#cgFO4~jI00T*UfI;2XPrqe zNdpUzpus3a?eNqu3$u*z>F)x!uZ3%iv0%qAivrdWLE#cFh2Wfa7A34}f>NnPsR>!w zQZ|zp$i-?3@n)u+3+hxY3HP2wLcegncoxm~? z!NpqL!kSUHo1H*U0n0UBNTeskv3Ugj*cN4p@k|78#|v;7nY{lsT4bop1eza^P7@;6z@wTpVG&*S&yg!OQK5lX=*KXZ2^!AoyMjcoSjF zi3wVi{ry=3^5eFGw?)f%)xN=rR|QbvO2gH381~grC#Ad@v6jLRQdXz$1vehLu58ta@uA5+cj+;sz0= z0$LdnMc{U}#>M@blItzwiGU7f*#Lc>=ExX{F+t#N!bm}5YS_ddfA(fgY+JJSH_vDBjqh${hzMc~TL1DMhX6)D zaV5Qku8p6#vg;Q3OL>Wld8xqUnit%+jQJzi#US@eXprl$xW8J5Lo4=6GvhPvxB^zj ztV{m_5)s$kXoYJ*6DI}}#i=wH8KDLMhDs)qJZnBz0tt?dt-D-Yn}Pj5V7Vq*xZ}A# z5}$2yyo<136FuCK)?{+#z26?Q5MdA2Dj5;&wy6D$s3)Qt4|PYUX6>(;lRBWMNmyHN znOG28FR{X|TIZ~=9VS9Fv8Eum5mwjp)X&2j#C6x|w*SDXApeNLL{^suBdl%{jSido zjA!($)H^9KK8m?tf(XHLn^uZ<;u9syyo1>VKQj;YWQU9vf>OQ(%U%Unl-OBLX64!N zLUQVgzXwMNny&P4J*X`_aFGPF;Di*_Hh`k@6dcnqHM`V3qm$qw1}%8ly8f$4dNB@y zgA*}2vHB=Aot@%}T($T$QqtAWetOwCh=X7fYx&)Bgvj*!CJKQ<02 zcxk(`c4p? z(hBZa<6S`nc9+0T$#YFSez&uh^mwpRdi|{QR=K#{0yoM)5_2`N+r@M_WtWi^>bX#A zHEACvl3xV2uoN`}>q34MY05L6w3+Y_*OFMSKP(qZ*pc0A4m8K@%+!vGSqWSzP;*mi z2oaZClkrxwx-BB;YsPthT$>5O>d-#a`CRP~BII3e}nM;vY zWGPM|K4w_FiES=sF~p2)&zF}mzPNaOUcAUc67n;i>QbQW?VJyzSF#IxlYu0XxXARx z=$y69zNfXBBIc=R@LKW54)iD&9b=`Zy0JnEQn$>KGPkCidD z=UzATO@V7AVsSYl{bJguWqM`|(;hhS+}T@dnVz{^Ord4fb-U4W>Nu{J@w=uV6QCy2 zAN)OV#$CyUR@DIMTngzjYFJM1HZdhwH+j&B?^8OLj3I*elCcWjh246S{S{htKb=?| z&q#7{_)>y{U&s8H95Y3Z%fO3##*zP$cv8kEDC5j~xJ!Ft7A9npi7Vyt*i%c5C1nw2 zmziV4va`SrGjXI0OohVx2HK9M{ncA)1YndG65##vCsCWm_OsbH_rpo4T>S@-N=QRW zAoMEdx>xJrb;;Q_YPT!Q7RC^(tnXm=SnAcRBko-{Q?SI`^{qdhPG52o>N16HiWXC# z1(&%2a#{3qsYYUaa@BMBOvRI&{^ypgF;aLlt@z)dfY z@2Qylm?^mO3*n?(dcFd=gx4g$B3_v>q@~;IVosXIZ2B;XT3*g?rT%~MMI$S>-ef-^ zvl`+mv{L_@NEPcINUL6{|8W))RfecRS)Z)b|2#|{>^Io^AFv|vth6ijzo~=qZ%geB z@4ddfYTjKXmWLtMkd^wM`y#dr5xWEZOQ`-YnAJS7gh4{b!H;R7)@Gzby@LYMU zNMJ<*-=z}J?R&@e?gn!<^-priqq_rmp1_k2!d3P=N4lRT))7t zbZeCV^?v|p7Y0wQ=iT)Va6!*16<(}%xt01Kc?M3#$@*}m{zp%T`wr{|g*_c%IcfmL O7S--C`F%(0|NjD%7?0%u diff --git a/utilities/._consolidateshapespacemodels.m b/utilities/._consolidateshapespacemodels.m index 8d844422df04854bba487efa48b29295fb665e6d..0323dc52c37c2c10c9dd28c60a4f90046fa6dbd2 100644 GIT binary patch literal 312 zcmZus%}T>S5T2E)ZH|JTbMTrtY`0-GUSgYwU?{XHdXq7(5^Q3cY#+$C@u56)vndGD zftl}z`344#FHYeIAf?=$FWqN%U!4ad(;R@E7&MwmgOEm9V)_vVojO!emMwbUCVmXB z|Lt|xJ6`}?=3k{xb-h-8sWiQjrkj_~boQWCxoTCD(MN03zh>S7od4}_w$xCii%P+3 zy)Igf_|2Z~8>K3T4;`BV#(K0L?Wu*c%@v{?hXIO!ay|){^iE9?qPUAAejRxj^VkzO SLr$rp diff --git a/utilities/._model2info_v2.m b/utilities/._model2info_v2.m index 384a82aab43ab43da5735f65a48d03e5c5d3f7ba..28f40602ec4be40fc051ef5cba1c5d100c03cf05 100644 GIT binary patch delta 22 ccmeBS>S3C2fQP%LC4m757N6zfn|Q_n08ZBj2><{9 delta 22 ccmeBS>S3C2fahADQvw4JxR$nWupNWQ29*!X)K_D@x+!i3C3Cu+hKz$4# zaSsLt{^a~zy~KiooK(Hs#N-@^nm<6wYGQ+g5Hmxo@1*;_&i->|t+~1V*~^6!CK@VF zbP%>G)Hg7&Ft#=`urM(pFGAbJY delta 129 zcmdnSbcJz(3Kzp71_lNuAigs3qMD9lNJtQn4FXqycm@#LfH+8i5y%FMPn;uRQ>brX zU}0=+nwFesnwV^09h{h!Sd?jPZ0hFhVr1Z?YvyL|q-$c~W})lkYH6wKXyNK);cV#W J;^gSU001&-8VCRY diff --git a/utilities/3D/.DS_Store b/utilities/3D/.DS_Store deleted file mode 100644 index dddfd225ff88fcdc97f53c22c9698f3b25b006af..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 14340 zcmeHN4RBml6+Y)BY2V9knwQd_WSe%kq?AHa$R@9f8C4Ud;yvaeV?eLAh7)7m6y26GMonX^ZzvHR ziN&n}y=EvJ9U6&-(?)wV$_i6jwKbGVM5EghmYFsa@m>9C<3Ku7kt3^rx0OhxJ9%XT zx;|nV+YPGA95uSAMyhHiYi(hY9x>zldx_9Epw?icrcW28 z&-X4;^@yfdi0Fmhi&TA+rcV>m#`;TC{iLcZJn9YBw`h7;)yqV*v3{kd>L?I!xmTj^P13gg4_S@D{ukKZUpBQT#mKf%oEZJb_=vui!)Yb^IoN3s2&s_yc?l zPvH;oaeNwIz*q2(qMr&|vFu@|j`ZNGi6OC4*+Ha)DQn65{NaUWfJx`lerTpu&d#VX zKc!lr3oe0Jcho`qj0mWQP+XWJxTYsSqq=Dg=M%o3>w`7O=yR9ZrallM)RlrBINf4KYxV zQ+j|Er!ZuP6CxDkfPx&*ZW(Qag1mL23-IEAArY4TOB^V1ppXOeD5;3w^I#kPgh!-Mwsjzl7w{SFWAB%m>% z&+1IX(_u4iq&Tq|F{sGZeWd0m#yS!s@rcFi=@A(kwWAH_wd3Q>!Bvfc^uCjqYmIga zrLb&Ch5<5+^Hh0yclP%6?{dXe=K|14_|XN%Y1%_G=GHG{xUKE-&MP->dyauQk7xFr z4A_!j7KdHCiRtk)Lu@IHK_L!qI8|+}Z=R~^4gDj7=Jg|H`k<(*)*m1^b<%3_DF=#XLXG}v$nf-ld35EXN80kSf-o59eV$F2JR@99s!Ub`XTzOdxVA z?!asCS{%SIfye|V@p`oIAb#8ike|al@h&`uU$SB2HwZXBhEL&-@hA8^{)PbK@9_8d z2EL8&;Jb=PnWfBD<|r2_7b`8w8l_v=E(Vu{le;s6wbOS{D9cyQ5OPlu3EI%TVdJLO z53B^8BLj$MfnOLp@+C4Osrq!b(p7GlSmlGv+BD&bK#IH{y=@8MS{$+k$VI~k> zxXIt_(aT8WG#+jBFZ1ZAib$JhwMQ>!CJ>%?xhLe&ftx^hywktVuQO3114w3`&$6f| zUD8j&v+zq=(y!B!{wJ2vLY|8OT+A&Ymy<1|9lLNHUWME6BiO?i@&N9|K^(#;#%VEM zk0~6-!*~Nu;7#~xJc6G^wxEyU-L#l>{E3gH1lGYl({4WK4(?t$n`Jm&GlTNvdOj>@^9pys(qMpgx?NS6|(~ zRbR79ji3CFYHV?-@y)y%=Q{crCr%AbeZ26KM{P;sW`a@-o1>HPP9bZTqoz4i ztjUFI%=+c1OQls{&@7rXwkZ25)Y+7Pq_Y*2mo&GhGTZL|N^`3l9}hM)HwJ<$gG_S^ zHn+Ak23lKM3sg6zn$c5Kb^9V^7K~I@;GqA^DYr7RG{5)y@$NIf;L==BS1gf}>B-8; zk`2PNSWs70NnwAn(T?Yk9O)bokHxyO=VLH&Rg!Wki9D&UY7Yt z*tuiT8*|z1>C|`lWQtHZc0$n0h^k0Qc|Qc0XzLcfj-sV=yW-hPJN8BhP$ytEjgpN_ zc#zV3jBLzj|wdc@XPoW|ct8Ub4?Og~7w>r0N2SCjWHOj?pfKc*Qn z4_=aM6{$jMrXM+|Oa4C>KcuWDOzJD78ItD`y+31fsdMf#Z))5zHz}H>c;TEdt>5A2 zoKC*cZdwQH;40{Wb$KhF$&Bn3_0oJZJyK}4BQQvDT`S|%x}v24mlNTfniy>g&EX>e0t{Y8x+nQ z!>ux{k7%99Ei>!Rp2VqzF*Xa-%e>nJZ2ddxoCe!T!rJ*AgnSuem{e}Rnp$VK_x4l$ zX+9%tZ>BbA^%eAju6jEBcZ-03Q(U~sa951{`c+ms<7A3 z86RK1yfM(kA>WENf_$r1t}K9jchd5wrzqsRi!xo=%z_JlRF4r0M=3xMzY*ZH-X-KyeH?vS>tNsc5&NeJp``q zoAYLQy^@@u4soFa8TU)FV_H>fI4oxi~6bvP%;Xe`Q~zrg2{oi{eY z;tPFV?z*8Hb)w+XMz5Ek6+gR%7`)1|*X z{Mz3~$7ZR7Ui_ynHDC-GT)34e_*;~9Jn|Ac?W zxA1)h6jaVr>XrG*0wt(4E31_*Ws9<1=@FyKdT#cJc;;W}&3u)@$vIFVvYXibzjej$ z|4&{?DYwLd5(hri9KiCfo~{lWLdL0>KN^J+0`%^smp!c9ga|na{jnnap$~QyZV45&&)^6a?!%FyH+?kiLb_xGGQMapwg8Q;%knS3C2fM;!8O9BHBocd=kKJknL09(}vLjV8( delta 22 ccmeBS>S3C2fXD2wd;$XyG}s!*O+4cO08jD;(*OVf diff --git a/utilities/3D/spharm-ppm/._spharm_obj_model.m b/utilities/3D/spharm-ppm/._spharm_obj_model.m index c5415162253c3c002418a3ad7106fc8d1e29c715..a2b5025c4d7e12ec1f599b2d64979cbd181cfa72 100644 GIT binary patch delta 25 fcmbQlHi>Nm4Nm4XIfHU`ArZBU5)>dBq^d;bQM(+E91?CBVL z+M|rdumqGvgO=OlaVnFk9DhB}7%}g3{D}LA6^Dy+@czTS^NrL0BjcpQLQdP?@F{ea zT)IUTv=dHCkW zPVa~Oc^2ZA&|mNNwX7+{r-+i4>ueYY+o xC+=DcS5T0!dwy4mPc#;Lrg2kJsSmHqtJoE)5ZBZ~`3;xh!?!Jd_;l+3G;>kzx zo86(!Hi;(@k(n~{ZDwb7CtpIQOGIY;I2#cSiKq`{>rDlw( z3YycBju8z5$0X@(?)Zv#?2D^-EI(*)3j%ul^ zVZD~x8X}|mPPjZ3Wk4BF27a6Y^lX;?c0jetfHI&AtQg?$gM~6Cf@MJe>45P^0ALSp zFBtP#f^)23B3K5*48(*~U`RD`#4sTpek*Z_U>Pu^laa%Rk)4g4P>kOl>sudACJCrk z8Bhjd2D;&L$o>Dk{rn#%>6J2|4E!qwOfQ|JW89MT*4oW+ueH!?C=17x0aqy)_*RTq dZpHghFYsF)029G7AUqKLBj9OJqYV5i17E{CWLy9M 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