models2report

 MODELS2REPORT Generate a report comparing various SLML models.

 COMMENT: When the method is not deployed the shape of the input arguments is
 slml2report( models, param )

 List Of Input Arguments  Descriptions
 -----------------------  ------------
 models                   A cell array of filenames
 param                    A structure holding the function options

 The shape of param is described

 List Of Parameters        Descriptions
 ------------------        ------------
 targetDirectory           (optional) Directory where the images are going to be saved. Default is current directory.
 prefix                    (optional) Filename prefix for the synthesized images. Default is 'demo'
 numberOfSynthesizedImages (optional) Number of synthesized images. Default is 1.
 compression               (optional) Compression of tiff, i.e. 'none', 'lzw' and 'packbits'
 debug                     (optional) Keeps temporary results and catches errors with full reports
 display                   (optional) Display flag for figures
 verbose                   (optional) Print the intermediate steps to screen. Default is true.
 microscope                (optional) Microscope model from which we select a point spread function. Default is 'none'
 synthesis                 (optional) Synthesis parameter that allows to
 synthesize                           'nucleus', 'framework' or 'all'. Default is 'all'
 sampling.method           (optional) Can be 'disc' or 'sampled'. Default is trimmed
 sampling.method.density   (optional) An integer. Default is empty.
 protein.cytonuclearflag   (optional) Can 'cyto', 'nucleus' or 'all'. Default is all.
 resolution.cell           (optional) The resolution of the cell and nucleus that are being passed in
 resolution.objects        (optional) The resolution of the object model being synthesized
 image_size                (optional) The image size. Default is [1024 1024] for both 2D and 3D in x and y

 output.tifimages           (optional) boolean flag specifying whether to write out tif images
 output.indexedimage        (optional) boolean flag specifying whether to write out indexed image
 output.blenderfile         (optional) boolean flag specifying whether to write out (.obj) files for use in blender
 output.blender.downsample  (optional) downsampling fraction for the creation of object files (1 means no downsampling, 1/5 means 1/5 the size)

 Example
 -------
 instances = { 'model01.xml', 'model02.mat', 'model03.mat' };
 param.targetDirectory = pwd;
 param.prefix = 'demo'
 param.numberOfSynthesizedImages = 100;
 param.compression = 'lzw';
 param.microscope = 'svi';
 param.verbose = true;

 >> slml2img( instances, param );

0001 function models2report( varargin )
0002 % MODELS2REPORT Generate a report comparing various SLML models.
0003 %
0004 % COMMENT: When the method is not deployed the shape of the input arguments is
0005 % slml2report( models, param )
0006 %
0007 % List Of Input Arguments  Descriptions
0008 % -----------------------  ------------
0009 % models                   A cell array of filenames
0010 % param                    A structure holding the function options
0011 %
0012 % The shape of param is described
0013 %
0014 % List Of Parameters        Descriptions
0015 % ------------------        ------------
0016 % targetDirectory           (optional) Directory where the images are going to be saved. Default is current directory.
0017 % prefix                    (optional) Filename prefix for the synthesized images. Default is 'demo'
0018 % numberOfSynthesizedImages (optional) Number of synthesized images. Default is 1.
0019 % compression               (optional) Compression of tiff, i.e. 'none', 'lzw' and 'packbits'
0020 % debug                     (optional) Keeps temporary results and catches errors with full reports
0021 % display                   (optional) Display flag for figures
0022 % verbose                   (optional) Print the intermediate steps to screen. Default is true.
0023 % microscope                (optional) Microscope model from which we select a point spread function. Default is 'none'
0024 % synthesis                 (optional) Synthesis parameter that allows to
0025 % synthesize                           'nucleus', 'framework' or 'all'. Default is 'all'
0026 % sampling.method           (optional) Can be 'disc' or 'sampled'. Default is trimmed
0027 % sampling.method.density   (optional) An integer. Default is empty.
0028 % protein.cytonuclearflag   (optional) Can 'cyto', 'nucleus' or 'all'. Default is all.
0029 % resolution.cell           (optional) The resolution of the cell and nucleus that are being passed in
0030 % resolution.objects        (optional) The resolution of the object model being synthesized
0031 % image_size                (optional) The image size. Default is [1024 1024] for both 2D and 3D in x and y
0032 %
0033 % output.tifimages           (optional) boolean flag specifying whether to write out tif images
0034 % output.indexedimage        (optional) boolean flag specifying whether to write out indexed image
0035 % output.blenderfile         (optional) boolean flag specifying whether to write out (.obj) files for use in blender
0036 % output.blender.downsample  (optional) downsampling fraction for the creation of object files (1 means no downsampling, 1/5 means 1/5 the size)
0037 %
0038 % Example
0039 % -------
0040 % instances = { 'model01.xml', 'model02.mat', 'model03.mat' };
0041 % param.targetDirectory = pwd;
0042 % param.prefix = 'demo'
0043 % param.numberOfSynthesizedImages = 100;
0044 % param.compression = 'lzw';
0045 % param.microscope = 'svi';
0046 % param.verbose = true;
0047 %
0048 % >> slml2img( instances, param );
0049 
0050 % Author: Robert F. Murphy
0051 % Created: June 29, 2013
0052 
0053 % Copyright (C) 2013 Murphy Lab
0054 % Carnegie Mellon University
0055 %
0056 % This program is free software; you can redistribute it and/or modify
0057 % it under the terms of the GNU General Public License as published
0058 % by the Free Software Foundation; either version 2 of the License,
0059 % or (at your option) any later version.
0060 %
0061 % This program is distributed in the hope that it will be useful, but
0062 % WITHOUT ANY WARRANTY; without even the implied warranty of
0063 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0064 % General Public License for more details.
0065 %
0066 % You should have received a copy of the GNU General Public License
0067 % along with this program; if not, write to the Free Software
0068 % Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
0069 % 02110-1301, USA.
0070 %
0071 % For additional information visit http://murphylab.web.cmu.edu or
0072 % send email to murphy@cmu.edu
0073 
0074 % if ~isdeployed
0075 %     if nargin ~= 2
0076 %         error( ['CellOrganizer: Wrong number of input arguments. If slml2img ' ...
0077 %             ' is not deployed, then the number of input arguments must be 2.'] );
0078 %     end
0079 %
0080 % end
0081 
0082 %INPUT ARGUMENTS
0083 param = varargin{2};
0084 
0085 param = ml_initparam(param,struct('verbose',true,'includenuclear',true,'includecell',true,'includeprot',true));
0086 
0087 models = varargin{1};
0088 
0089 n=length(models);
0090 
0091 try
0092     indvarlabels = varargin{3};
0093     nolabels=false;
0094     legendentry = regexp(indvarlabels, ';', 'split');
0095 catch
0096     nolabels=true;
0097 end
0098 
0099 try
0100     indvar = varargin{4};
0101     if nolabels indvarlabels = 'Indepedent variable'; end
0102     for i = 1:length(indvar)
0103         legendentry{i} = num2str(indvar(i));
0104     end
0105 catch
0106     indvar = [1:n];
0107     if nolabels
0108         indvarlabels = 'Model number';
0109         for i = 1:length(indvar)
0110             legendentry{i} = num2str(indvar(i));
0111         end
0112     end
0113 end
0114 
0115 colors = jet(length(indvar));
0116 
0117 is3D = strcmpi(models{1}.dimensionality, '3D');
0118 
0119 if is3D
0120     protfield = 'proteinShape';
0121 else
0122     protfield = 'proteinModel';
0123 end
0124 
0125 [ nuclearparams, cellparams, proteinparams ] = getModelReportParams( models{1}, param );
0126 
0127 %
0128 % compare protein model parameters
0129 %
0130 
0131 for j = 1:n
0132     if is3D
0133     compartments{j} = eval(['models{j}.' protfield '.cytonuclearflag']);
0134     else
0135         compartments{j} = 'all';
0136     end
0137 end
0138 
0139 
0140 y = []; 
0141 if param.includeprot
0142     
0143     for j=1:n
0144         
0145         
0146         
0147         for i=1:length(proteinparams)
0148             y(j,i)=eval(['models{j}.' protfield '.' proteinparams{i}]');
0149         end
0150         
0151         if is3D
0152             [positionBeta, fractdist] = beta2posMap(models{j}.proteinShape.position.beta);
0153         else
0154             [positionBeta, fractdist] = beta2posMap(models{j}.proteinModel.positionModel.beta);
0155         end
0156         
0157         objpdf = squeeze(sum(sum(positionBeta,1),3));
0158 
0159         if strcmpi(compartments{j}, 'cyto')
0160             domain = fractdist >= 0;
0161             objpdf(~domain) = [];
0162             objpdf = [zeros(1,1+sum(~domain)), objpdf ] ;
0163             fractdists{j} = [fractdist(~domain) 0 fractdist(domain)];
0164         
0165         elseif strcmpi(compartments{j}, 'nuc')
0166             domain = fractdist <= 0;
0167             objpdf(~domain) = [];
0168             objpdf = [objpdf zeros(1,1+sum(~domain))];
0169             fractdists{j} = [fractdist(domain) 0 fractdist(~domain)];
0170             
0171         else
0172             fractdists{j} = fractdist;
0173         end
0174         
0175         %renormalize
0176         objpdf = objpdf./sum(objpdf);
0177         
0178         objpdfs{j} = objpdf;
0179     end
0180 
0181     
0182 end
0183 %
0184 % compare nuclear parameters
0185 %
0186 y2 = []; 
0187 y3 = []; cellparams = [];
0188 if param.includenuclear | param.includecell
0189     if param.includenuclear
0190        
0191         for i=1:length(nuclearparams)
0192             for j=1:n
0193                 y2(j,i)=eval(['models{j}.nuclearShapeModel.' nuclearparams{i}]');
0194             end
0195         end
0196 
0197     end
0198 
0199     tparam = [];
0200     tparam.generatemeanshape = true;
0201     tparam.display = false;
0202     tparam.debug = false;
0203 
0204     havealready = length(nuclearparams);
0205     for j=1:n
0206         model = models{j};
0207         tparam.resolution.cell = model.cellShapeModel.resolution;
0208         tparam.resolution.objects = eval(['model.' protfield '.resolution']);
0209         % generate average nuclear instance
0210         
0211         if is3D
0212             instance = tp_genspsurf(model.nuclearShapeModel,tparam);
0213             [nucimg,nucsurf,outres] = tp_gennucshape(instance,tparam);
0214 
0215             nucleus.nucimgsize = size(nucimg);
0216             nucleus.nucsurf = nucsurf; clear nucsurf;
0217             nucleus.nucimg = nucimg; clear nucimg;
0218             % generate cell shape instance from it
0219             [cellimg,nucimg] = ml_gencellshape3d( ...
0220                         model.cellShapeModel, nucleus,tparam );
0221 
0222             if param.includenuclear
0223                 stats=rm_get3Dstats(nucimg);
0224                 statnames = fieldnames(stats);
0225                 for k=1:length(statnames)
0226                     y2(j,k+havealready) = eval(['stats.' statnames{k}]);
0227                     nuclearparams{k+havealready} = statnames{k};
0228                 end
0229 
0230                 clear nucimg;
0231             end
0232 
0233             if param.includecell
0234                 stats=rm_get3Dstats(cellimg);
0235                 statnames = fieldnames(stats);
0236                 for k=1:length(statnames)
0237                     y3(j,k) = eval(['stats.' statnames{k}]);
0238                 end
0239                 cellparams=statnames';
0240 
0241                 clear cellimg;
0242             end
0243         
0244         end
0245     end
0246 end
0247 
0248 paramnames = [proteinparams nuclearparams cellparams];
0249 yall = [y y2 y3];
0250 
0251 %
0252 % plot all of the parameters for all of the models individually, and plot
0253 % the ones that change the most together
0254 %
0255 
0256 %save testing
0257 
0258 %end
0259 
0260 figure
0261 npanel = length(paramnames);
0262 ncolumns = 3; nrows = 5; ntot = ncolumns*nrows;
0263 for iplot=1:ceil(npanel/ntot);
0264     nleft = npanel-(iplot-1)*ntot;
0265     for i=(iplot-1)*ntot+1:(iplot-1)*ntot+min(nleft,ntot)
0266         subplot((ceil(min(nleft,ntot)/ncolumns)),ncolumns,i-(iplot-1)*ntot);
0267         plot(indvar,yall(:,i),'rd-');
0268         xlabel(indvarlabels);
0269         ylabel(paramnames{i});
0270     end
0271     figure
0272 end
0273 
0274 normy = yall./repmat(mean(yall),n,1);
0275 cv=std(normy);
0276 [sorted,order]=sort(cv,'descend');
0277 % find the ones whos CV is greater than 50%
0278 idx = find(sorted>0.5);
0279 mintodisplay = 3;
0280 if length(idx)> mintodisplay
0281     order = order(find(sorted>0.5));
0282 %    order = order(end:-1:1);
0283 else
0284     order = order(1:mintodisplay);
0285 %    order = order(end:-1:end-mintodisplay+1);
0286 end
0287 
0288 colorstr = {'rd-','g+-','b^-','co-','mv-','yx-'};
0289 for i=1:min(length(colorstr),length(order))
0290     plot(indvar,normy(:,order(i)),colorstr{i});
0291     legends{i}=paramnames{order(i)};
0292     hold on;
0293 end
0294 legend(legends);
0295 xlabel(indvarlabels);
0296 ylabel('Coefficient of variation');
0297 title('Parameters ordered by extent of variation');
0298 hold off;
0299 
0300 if param.includeprot
0301     figure, 
0302     for i = 1:length(models)
0303         %marginalize out the anglular dimensions.
0304 
0305         plot(fractdists{i}, objpdfs{i} , 'color', colors(i,:), 'LineWidth',2)
0306         hold on;
0307     end
0308 
0309     xlabel('Fractional distance between nuclear and plasma membranes')
0310     ylabel('Relative probability density')
0311     legend(legendentry);
0312     set(gca,...
0313     'YTickLabel','')
0314     
0315     
0316     
0317     if is3D
0318         
0319 
0320         objmu = y(:,strcmpi(proteinparams, 'frequency.mu'));
0321         objsigma = y(:,strcmpi(proteinparams, 'frequency.sigma'));
0322         xrange = [0 max(exp(objmu+objsigma))];
0323 
0324         figure, 
0325         for i = 1:length(models)
0326             semilogx(xrange(1):xrange(2), lognpdf(xrange(1):xrange(2), objmu(i), objsigma(i)), 'color', colors(i,:), 'LineWidth',2)
0327             hold on
0328         end
0329         xlabel('Number of objects')
0330         ylabel('Relative probability density')
0331         legend(legendentry);
0332         set(gca,...
0333         'YTickLabel','')
0334     else
0335         objalpha = y(:,strcmpi(proteinparams, 'objectModel.numStatModel.alpha'));
0336         objbeta = y(:,strcmpi(proteinparams, 'objectModel.numStatModel.beta'));
0337 
0338         xrange = 0:1000;
0339         
0340         for i = 1:length(models)
0341             gamcdfs(i,:) = gamcdf(xrange, objalpha(i), objbeta(i));
0342             gampdfs(i,:) = gampdf(xrange, objalpha(i), objbeta(i));
0343         end
0344         
0345         inds = any(gamcdfs < 0.95,1);
0346         
0347         figure, 
0348         
0349         for i = 1:length(models)
0350             plot(xrange(inds), gampdfs(i,inds)','color', colors(i,:))
0351             legendentry{i} = num2str(indvar(i));
0352             hold on
0353         end
0354         
0355         xlabel('Number of objects')
0356         ylabel('Relative probability density')        
0357         legend(legendentry);
0358         set(gca,...
0359         'YTickLabel','')
0360     end
0361 end
0362 
0363 end