model2diffeomorphicInstance

MODEL2DIFFEORMORPHICINSTANCE Helper method that synthesizes a framework
from a diffeomorphic model

0001 function [nucimg, cellimg] = model2diffeomorphicInstance( model, param )
0002 %MODEL2DIFFEORMORPHICINSTANCE Helper method that synthesizes a framework
0003 %from a diffeomorphic model
0004 
0005 %icaoberg 10/1/2012
0006 %
0007 %Edited:
0008 %
0009 %6/25/13 D. Sullivan - moved this to its own file so it can be called
0010 %externally from model2framework.
0011 %6/25/13 D. Sullivan - added support for param.randomwalk, a flag that
0012 %indicates that a random walk is desired.
0013 %9/6/13 added support for the chunk_start method.
0014 %9/22/13 grj - Added missing chunk_finish ^^
0015 %9/22/13 grj - Moved framefolder variable construction outside of
0016 %              conditional part of code
0017 % Copyright (C) 2007-2013  Murphy Lab
0018 % Carnegie Mellon University
0019 %
0020 % This program is free software; you can redistribute it and/or modify
0021 % it under the terms of the GNU General Public License as published
0022 % by the Free Software Foundation; either version 2 of the License,
0023 % or (at your option) any later version.
0024 %
0025 % This program is distributed in the hope that it will be useful, but
0026 % WITHOUT ANY WARRANTY; without even the implied warranty of
0027 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0028 % General Public License for more details.
0029 %
0030 % You should have received a copy of the GNU General Public License
0031 % along with this program; if not, write to the Free Software
0032 % Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
0033 % 02110-1301, USA.
0034 %
0035 % For additional information visit http://murphylab.web.cmu.edu or
0036 % send email to murphy@cmu.edu
0037 
0038 
0039 options = model.nuclearShapeModel.shape_space_options;
0040 
0041 try
0042     verbose = param.verbose;
0043     if isa( verbose, 'logical' )
0044         verbose = false;
0045     end
0046 catch
0047     verbose = false;
0048 end
0049 
0050 %icaoberg 9/28/2012
0051 try
0052     debug = param.debug;
0053     if isa( debug, 'logical' )
0054         debug = false;
0055     end
0056 catch
0057     debug = false;
0058 end
0059 
0060  %set up the framework walk output folder
0061 if ~isfield(param,'framefolder')
0062     framefolder  = ['.' filesep 'frames' filesep];
0063 else
0064     framefolder = param.framefolder;
0065 end
0066 %if the folder doesn't exist yet, make it
0067 if ~isdir(framefolder)
0068     mkdir(framefolder);
0069 end
0070 
0071 
0072 try
0073     %D. Sullivan 6/25/13 added support for random walk
0074     if ~isfield(param,'randomwalk')||param.randomwalk==0
0075         %icaoberg 10/1/2012
0076         %number_shapes = size(y2, 1);
0077         number_shapes = size(model.nuclearShapeModel.positions, 1);
0078         
0079         % Generate a random point inside the convex hull of training shapes in the
0080         % shape space (probably not uniform sampling):
0081         random_weights = rand(number_shapes, 1);
0082         random_weights = random_weights ./ sum(random_weights);
0083         
0084         %icaoberg 10/1/2012
0085         %random_point = random_weights' * y2;
0086         random_point = random_weights' * model.nuclearShapeModel.positions;
0087         
0088         % Generate the shape associated with that point:
0089 %         options = struct();
0090         options.use_compression = false;
0091     
0092     %D. Sullivan 6/25/13 added random walk code below
0093     elseif param.randomwalk==1
0094         
0095        
0096         
0097         %max steps
0098         if ~isfield(param,'walksteps')
0099             warning(['No walkstep number specified (param.walksteps).',...
0100                 'synthesizing only 1 image']);
0101             maxIter = 1;
0102         else
0103             maxIter = param.walksteps;
0104         end
0105 %         maxIter = 50;
0106         
0107         %D. Sullivan 7/16/13
0108         %check if the work has already been done
0109         %In the future we should skip to the last synthed shape instead of
0110         %only skipping when totally done.
0111         if isdir([framefolder filesep 'Cellwalk'])
0112             donecells = ml_ls([framefolder filesep 'Cellwalk' filesep '*.tif']);
0113             donenucs = ml_ls([framefolder filesep 'Nucwalk' filesep '*.tif']);
0114             if length(donecells)==param.walksteps
0115 %                 for i = 1:length(donecells)
0116 %                     nucimg{i} = ml_readimage(donenucs{i});
0117 %                     cellimg{i} = ml_readimage(donecells{i});
0118 
0119                     nucimg = ml_readimage(donenucs{1});
0120                     cellimg = ml_readimage(donecells{1});
0121                     for z = 1:size(cellimg,3)
0122                         nucimg(:,:,z) = imfill(nucimg(:,:,z),'holes');
0123                         cellimg(:,:,z) = imfill(cellimg(:,:,z),'holes');
0124                     end
0125 %                 end
0126                 return
0127             end
0128         end
0129         
0130         %save some intermediate results (may want to delete these in the
0131         %final version)
0132         savepath = [framefolder filesep 'walk'];
0133 %         savepath = '/Users/devinsullivan/MurphyLab/CellShape/walk';
0134         
0135         %set random walk paraemters
0136         %For right now we will hard code all these options while testing
0137         %type of walk
0138         if isfield(param,'walk_type')
0139             walk_type = param.walk_type;
0140         else 
0141 %         walk_type = 'brownian';
0142             walk_type = 'willmore';
0143         end
0144         
0145         if strcmpi(walk_type,'density')
0146             energies = model.cellShapeModel.density';
0147         elseif strcmpi(walk_type,'willmore')
0148             energies = model.cellShapeModel.Willmore_energy;
0149         else
0150             energies = ones(size(model.cellShapeModel.positions,1),1);
0151         end
0152         
0153         
0154         %D. Sullivan 7/22/13 - changed Dc and dt to be params of the
0155         %cellshape motion model
0156         if isfield(model.cellShapeModel,'motion')
0157             if isfield(model.cellShapeModel.motion,'Dc')
0158                 Dc = model.cellShapeModel.motion.Dc;
0159             else
0160                 warning(['No diffusion constant found in motion model,'...
0161                     'defaulting to 10e-3']);
0162                 Dc = 10e-3;
0163             end
0164             if isfield(model.cellShapeModel.motion,'dt')
0165                 dt = model.cellShapeModel.motion.dt;
0166             else
0167                 warning(['No time step found in motion model,'...
0168                     'defaulting to 1']);
0169                 dt = 1;
0170             end
0171         else 
0172             warning(['Shape space walk requested with no motion model.'...
0173                 'defaulting diffusion = 10e-3, time step = 1']);
0174             % diffusion constant
0175             Dc = 10e-3;
0176             
0177             %time step
0178             dt = 1;
0179         end
0180 
0181         
0182         %Set up the shape_space
0183         y2 = model.nuclearShapeModel.positions;
0184         convex_hull = model.nuclearShapeModel.convex_hull;
0185         tes = model.nuclearShapeModel.tesselation;
0186         
0187         shape_space = { y2, convex_hull, tes};
0188         
0189         %make the function call
0190         [result, success] = generate_walk_path(...
0191             shape_space, walk_type, Dc, dt, maxIter, savepath,energies);
0192         
0193         random_point = result.walk_path;
0194         
0195     end
0196 
0197     %icaoberg 11/1/2012
0198     try
0199        default_options.method = param.synthesis.diffeomorphic.method;
0200     catch
0201        %icaoberg
0202        %options.method = 'convnfft';
0203        default_options.method = 'convnfft_fast';
0204     end
0205 
0206     %yy 11/28/2012
0207     try
0208        default_options.convergence_absolute_error = param.synthesis.diffeomorphic.convergence_absolute_error;
0209     catch
0210 
0211     end
0212 
0213     %icaoberg 11/9/2012
0214     %maximum number of iterations
0215     try
0216        default_options.maximum_iterations = param.synthesis.diffeomorphic.maximum_iterations;
0217     catch
0218        default_options.maximum_iterations = 100;
0219     end
0220 
0221     %changes per step
0222     try
0223        default_options.convergence_tolerance = param.synthesis.diffeomorphic.convergence_tolerance;
0224     catch
0225        default_options.convergence_tolerance  = 0;
0226        %options.convergence_tolerance  = 8.279212e-02;
0227     end
0228 
0229     %changes per step
0230     try
0231        default_options.convergence_registration_error_scale = param.synthesis.diffeomorphic.convergence_registration_error_scale;
0232     catch
0233        default_options.convergence_registration_error_scale  = 5e-3;
0234     end
0235 
0236     %tebuck 11/5/2012
0237     %forward euler method as rk integrator
0238     %options.integrator_c = 0;
0239     %options.integrator_b = 1;
0240     %options.integrator_a = 0;
0241     %options.integrator_b2 = [];
0242     %options.integrator_first_same_as_last = false;
0243     
0244     %tebuck 11/9/2012
0245     default_options.alpha = 1;          % coefficient of Lapacian operator
0246     default_options.gamma = 0.04;       % Coefficient of the Identity
0247     default_options.keep_proportion = 0; 
0248     default_options.single_sided = false; 
0249     default_options.max_time = 50000; 
0250     default_options.step_size = 0.01; 
0251     default_options.use_fft = true; 
0252     default_options.use_gaussian_kernel = false; 
0253     default_options.absolute_distance_tolerance = 0; 
0254     default_options.absolute_deformation_tolerance = 0; 
0255     default_options.periodic_space = false; 
0256     default_options.use_compression = false; 
0257     default_options.drop_kernel = true; 
0258     default_options.maximum_registration_error_failures = 0;
0259 
0260     %tebuck 11/25/2012
0261     %note that the following code is intended for use with values of 1 or 2
0262     %downsampling scale
0263     try
0264       default_options.downsampling_scale = param.synthesis.diffeomorphic.downsampling_scale;
0265     catch
0266       default_options.downsampling_scale = 1;
0267     end
0268 
0269     default_options.maximum_deformation_per_step = [4/default_options.downsampling_scale, 4/default_options.downsampling_scale, 0.5];
0270     
0271     %tebuck updated radius to match new model
0272     default_options.window_radius = 192/default_options.downsampling_scale; 
0273     default_options.filter_radius = 16/default_options.downsampling_scale;
0274     
0275     options = process_options_structure(default_options, options);
0276     
0277     
0278     %icaoberg 10/1/2012
0279     %D. Sullivan 6/25/13 - added support for randomwalks
0280     frame = cell(1,size(random_point,1));
0281     nucimg = cell(1,size(random_point,1));
0282     cellimg = cell(1,size(random_point,1));
0283 
0284     energy = zeros(1,size(random_point,1));
0285     if size(random_point,1)>1 && isfield(param,'templateJob')
0286         mkdir walkscripts
0287     end
0288     %D. Sullivan 7/22/13 save walk options so individual jobs can load them
0289     param = ml_initparam(param,struct('tempdir',[pwd filesep 'temp' filesep]));
0290     if size(random_point,1)>1
0291         disp(['Random walk detected,',...
0292             'generating and saving frames of the movie to disk. ',...
0293             'The method will only return the final frame to save memory']);
0294         
0295     end
0296     for i = 1:size(random_point,1)
0297         framefile = [param.tempdir filesep 'genframe' num2str(i)];
0298         [startimage, ~, ~, framefile] = chunk_start(framefile);
0299         
0300         if startimage
0301             try
0302                 currframe = i;
0303                 curr_point = random_point(i,:);
0304 %                 %save an file for the current frame to be read by the script
0305 %                 save(framefile,...
0306 %                     'model','options','curr_point','framefolder','currframe');
0307                 makeFrame(model,options,curr_point,framefolder,currframe);
0308                 chunk_finish(param.tempdir,['genframe' num2str(i)])
0309             catch err
0310                 disp(['Skipping image ' num2str(i) ' due to error']);
0311                 getReport( err, 'extended')
0312             end
0313         else
0314             disp(['Result file found, Skipping image ' num2str(i) ' ']);
0315         end
0316       
0317         chunk_finish(framefile)
0318     end
0319     
0320     %D. Sullivan 6/25/13 temporary conversion back to traditional CO
0321     %output (not cellarray) so that no code is broken by this addition
0322     %     if length(cellimg)==1
0323     %         cellimg = cellimg{1};
0324     %         nucimg = nucimg{1};
0325     %     end
0326     
0327     %D. Sullivan 7/22/13
0328     %initialize counter, this is to hold the job open until the first frame is
0329     %done
0330     if ~isempty(dir([param.tempdir filesep '*.tmp']))
0331         warning(['CellOrganizer has detected .tmp files still exist.',... 
0332             'Either a process is still working or one or more frame',...
0333             'failed to generate. Returning empty cell and nuclear image.']);
0334         cellimg = [];
0335         nucimg = [];
0336     else    
0337         cellimg = ml_readimage([framefolder filesep 'Cellwalk' filesep 'frame1.tif']);
0338         nucimg = ml_readimage([framefolder filesep 'Nucwalk' filesep 'frame1.tif']);
0339     end
0340     
0341 catch err
0342     getReport( err, 'extended')
0343     nucimg = [];
0344     cellimg = [];
0345 end
0346 end