function [nucimg,nucsurf,outres] = tp_gennucshape( instance,param )
% TP_GENNUCSHAPE generates a nuclear shape from the tensor-product spline
% surface structure

% Author: Ivan E. Cao-Berg
%
% Copyright (C) 2012-2013 Murphy Lab
% Lane Center for Computational Biology
% School of Computer Science
% Carnegie Mellon University
%
% July 23, 2012 R.F. Murphy Fix bug that caused approximately only the 
%                   bottom half of nucleus to be returned
% August 1, 2012 R.F.Murphy Don't shrink nuclear size - throws off scaling
%                   of protein objects relative to nucleus and cell
% Feb 22, 2013 D. Sullivan made scaling factor resolution dependent
% March 5, 2013 D. Sullivan f is now 1x3 use to account for resolution
%               adjustment in each dimension
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published
% by the Free Software Foundation; either version 2 of the License,
% or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA.
%
% For additional information visit http://murphylab.web.cmu.edu or
% send email to murphy@cmu.edu


if nargin < 1
    error('1 or 2 arguments are required');
end

if ~exist('param','var')
    param.verbose = false;
    param.debug = false;
    param.display = false;
end

%G. Johnson 7/9/2013
%if no resolution is specified assume everything is the same resolution
if ~isfield(param, 'resolution')
    param.resolution.objects = [1,1,1];
    param.resolution.cell = [1,1,1];
    param.resolution.nucleus = [1,1,1];
end

%G. Johnson 7/9/2013
%assume framework synthesis by default
if ~isfield(param, 'synthesis')
    param.synthesis = 'framework';
end
    
% downsampling factor used in generating the models
% f = 5;
%D. Sullivan 2/22/13 changing this to be resolution dependent
%note, objects should have 
%D. Sullivan 3/5/13 f is now 1x3 use to account for resolution adjustment
%in each dimension

%icaoberg 8/7/2013
%this change was done to allow users to synthesize only nuclear shape model
%files
%G. Johson 9/7/2013 
%Bug fixes
if strcmpi(param.synthesis,'all')
    f = param.resolution.cell./param.resolution.objects;
    outres = param.resolution.objects;
outres = param.resolution.objects;
elseif strcmpi( param.synthesis, 'framework' )
    f = [1,1,1];
    outres = param.resolution.cell;
else
    f = [1,1,1];
    outres = param.resolution.nucleus;
end

instance.coefs = shiftdim(instance.coefs,-1);

%D. Sullivan 3/5/13 f is now 1x3 use f(3) for z resolution adjustment
%icaober 7/8/2013
if strcmpi( param.synthesis, 'framework' ) || ...
        strcmpi( param.synthesis, 'all' )
    H = f(3)*instance.height;
else
    H = instance.height;
end

%icaoberg 8/2/2012
%param = ml_initparam(param,...
%    struct('xsize',1024,'ysize',1024,'zsize',H+1,'samp_rate',360,'debug',true));
factor = 1.25;
% factor = 1;
%D. Sullivan 7/6/13 - this should not be hardcoded at factor*1024. it makes
%no sense. 
param = ml_initparam(param,...
    struct('xsize',factor*1024,'ysize',factor*1024,'zsize',H+1,'samp_rate',360,'debug',true));

xcenter = param.xsize / 2;
ycenter = param.ysize / 2;

delta = 2*pi/param.samp_rate;
Phi = -pi:delta:pi;

% generate one nucleus surface at the specified height
Z = 0:(1/H):1;
[Phi_grid, Z_grid] = meshgrid(Phi,Z);
mesh_data = [Z_grid(:), Phi_grid(:)]';
%D. Sullivan 3/5/13 f is now 1x3 use f(1) for x and y resolution adjustment
nucsurf = f(1) * reshape(fnval(instance,mesh_data),[length(Z),length(Phi)]);
%nucsurf = nucsurf/2;

%icaoberg 7/1/2013
if param.debug && param.display
    plotcylsurf(nucsurf,delta);
end

% generate another one at the mean height for later use in finding the
% nuclear position
MEAN_HEIGHT = 85;
Z = 0:(1/MEAN_HEIGHT):1;
[Phi_grid, Z_grid] = meshgrid(Phi,Z);
mesh_data = [Z_grid(:), Phi_grid(:)]';
%D. Sullivan 3/5/13 f is now 1x3 use f(1) for x and y resolution adjustment
nucsurf2 = f(1) * reshape(fnval(instance,mesh_data),[length(Z),length(Phi)]);
%nucsurf2 = nucsurf2/2;

nucimg = uint8(zeros(param.ysize,param.xsize,param.zsize));
for i = 1:(H+1)
    [x,y] = pol2cart(Phi,nucsurf(i,:));
    x = x + xcenter;
    y = y + ycenter;
    sliceimg = zeros(param.ysize,param.xsize);
    for t = 1:length(x)-1
        rpts = round(linspace(y(t),y(t+1),100));
        cpts = round(linspace(x(t),x(t+1),100));
        index = sub2ind([param.ysize,param.xsize],rpts,cpts);
        sliceimg(index) = 255;
    end
    nucimg(:,:,i) = imfill(sliceimg,'holes');
end

%icaoberg 7/1/2013
if param.debug && param.display
    try
        for i=1:size(nucimg,3)
            imshow(nucimg(:,:,i),[0 255]);
            pause(0.1)
        end
    catch
        disp('Unable to display image');
    end
end

%blankimg = zeros(param.ysize,param.xsize);
%nucimg = cat(3,blankimg,nucimg,blankimg);
