%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This program is a translation of the Cochlea.r script written in R
% language and distributed by the Kral Lab here: 
% http://www.neuroprostheses.com/AK/CochleaModel.html

% The paper it referes to, and should be cited if used is:
% Pietsch, M., Aguirre Dávila, L., Erfurt, P., Avci, E., Lenarz, T., 
% & Kral, A. (2017). Spiral Form of the Human Cochlea Results from Spatial
% Constraints. Scientific Reports, 7(1), 7500. 
% https://doi.org/10.1038/s41598-017-07795-4
% For the PDF directly:
% http://www.neuroprostheses.com/AK/CochleaModel_files/Sci_Rep_2017.pdf

% For any questions, bugs and queries please contact the author 
% Lukas Aguirre Davila via e-mail: aguirredavila.lukas@mh-hannover.de

% The translation from R to Matlab was done by Andreas Prokopiou. 
% Contact: prokopiou.a@me.com
%
% Date: 9th Nov 2017
%
% Few modifications by A. Kral (making 3D plot proportional along all 
% 3 axes, adding the option to switch it off, adding the grid, adding the 
% turnturn)
% 
% Date: 10th Nov 2017
%
% Further modifications suggested by D. Schurzig, MedEl Comp., Hannover,
% who used a quicker 3D plotting routine. So we replaced the old routine
% with the faster one.
%
% Date: 30th Nov 2017
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cochlea(varargin)

% cochlea.m generates functions to predict the form of a cochlea from base
% measures from a CT scan (see paper for more details)
% Input: 
%   - no input (i.e. cochlea()) = generates a cochlea with the mean
%   parameter values of the study.
%   - cochlea('mean') = same as no input call.
%   - cochlea('random') = generates a random cochlea.
%   - cochlea(A), = generates a cochlea with the parameter set A, 
%     where A is a 1x4 vector representing the base measures of a cochlea (see paper).
%     Example: cochlea([5.97 3.95 3.26 2.85])

% Check the input arguments

movie=0;
if nargin == 0
    A = generate_cochlea_parameters('mean');
    disp('No input chooses the mean cochlea values')
elseif nargin == 1
    if ischar(varargin{1})
        switch varargin{1}
            case 'mean'
                A = generate_cochlea_parameters('mean');
                disp('Mean cochlea values')
            case 'random'
                % generate random base measures
                A = generate_cochlea_parameters('random');
                disp('Random cochlea values')
            otherwise
                error('The input options are ''mean'' and ''random'' only.')
        end
    else
        A = varargin{1};
        if length(A) ~= 4
            error('Input array must be 4 elements.')
        end        
    end
    
else
    A = varargin{1};
    if length(A) ~= 4
       error('Input array must be 4 elements.')
    end
    if nargin==2
        error('Too many inputs.')
    end
end

% The A measurment is the distance from modiolus at the base for 0, 90, 180
% and 270 degrees. Check out paper for details.

% the following are the coefficients for parameter estimation of the model
% functions. Check our paper for details.
coefs.height=[-2.025562262381,	0.308684774745994, 		-0.0245970612199774,	0.00093743547596123,	-6.49142944919421e-06 ;...
    -0.483933701821246,	0.150263587271025, 		-0.0249700764778651,	0.00184112837255666,	-4.68696577096306e-05;...
    0.0936906993853392,	0.0485929242340616,		-0.0111388827166082,	0.000852542615496957,	-2.39666478178502e-05;...
    0.437215616236419,	-0.609467318753365,		0.147019129524084,		-0.012345472153677,		0.000341806092537278;...
    0.19884386428529,		0.0706538846984097,		-0.0490269768729332,	0.00553612473205947,	-0.000178036145575138];

coefs.modiolus=[-0.0972007477234853,	0.0652777719428317,		0.00579010295961996,	-0.000410053336041606 ;...
    0.745407586062517,	-0.297795950930632,		0.0353079060718429,		-0.00125774660897636;...
    0.349940654425921,	-0.0488173880945223,	0.00235713329581525,	-3.38189250648794e-05;...
    0.0567728210841444,	0.115930983938798,		-0.0170603835499829,	0.000605305012135639;...
    -0.0613995454924607,	0.170502985259421,		-0.0267552502171719,	0.00111572946871731];

coefs.turns = [963.166310413576,	6.37772934525638,	-26.9585473096045,	-36.3953582023656,	66.9416454453684];


% Generate the cartesian functions (x , y ,z)
c_length = turn_number_estimation(A,coefs.turns)*2*pi;

syms phi v

r_modiolus = subs(radius_modiolus_estimation(A,coefs.modiolus),'x',phi);
r_scala_vent_dor = (c_length-phi)/c_length*0.5 + 0.6;
p = subs(height_estimation(A,coefs.height),'x',phi);

x = (r_modiolus+r_scala_vent_dor*(cos(v)-1))*cos(phi);
y = (r_modiolus+r_scala_vent_dor*(cos(v)-1))*sin(phi);
z = r_scala_vent_dor*sin(v)+p;



% Verbose output
disp(['Parameters: A1 = ' num2str(A(1))   ', B1 = ' num2str(A(2)) ...
    ', A2 = ' num2str(A(3)) ', B2 = ' num2str(A(4)) ])
disp(['estimated number of turns: ' num2str(c_length/(2*pi))])
disp('estimated length: ')
disp(['   with height: ' num2str(length_calculation_estimation_with_height(c_length,A,coefs)) ' mm'])
disp(['   without height: ' num2str(length_calculation_estimation(c_length,A,coefs)) ' mm'])


% Create plots
resolution = 0.1;

phi = 0:resolution:c_length;
V = 0:resolution:2*pi;


% Plot the Distance from modiolus and Height for the instance
turns = phi/ (2*pi);
A_turns = [0, 90, 180, 270]/360;

figure_handle = figure('color','w','WindowStyle','docked');
subplot(2,1,1)
r_modiolus = double(subs(radius_modiolus_estimation(A,coefs.modiolus),'x',phi));
plot(turns,r_modiolus)
hold on 
plot(A_turns,A,'o')
xlabel('turns')
ylabel('mm')
title('Distance from Modiolus - estimation')
axis([-0.25 ceil(turns(end)) 0 1+ceil(max(r_modiolus))])
drawnow;

annotation(figure_handle,'textbox',...
    [0.332142857142857 0.814285714285715 0.526785714285714 0.0966189857142875],...
    'String',{'modiolus distance polynomial',...
    char(vpa(radius_modiolus_estimation(A,coefs.modiolus),4))},...
    'EdgeColor',[1 1 1]);

subplot(2,1,2)
height = double(subs(height_estimation(A,coefs.height),'x',phi));
plot(turns,height)
xlabel('turns')
ylabel('mm')
title('Height - estimation')
axis([-0.25 ceil(turns(end)) -4 0])
annotation(figure_handle,'textbox',...
        [0.451785714285714 0.142857142857143 0.4375 0.111904761904762],...
        'String',{'height polynomial', ...
        char(vpa(height_estimation(A,coefs.height),3))},...
        'FitBoxToText','off',...
        'EdgeColor',[1 1 1]);
drawnow;

%% Modified version after Daniel Schurzig, MedEl Comp.

inp.meas = A;

% Angular resolution of output spiral in deg
% (e.g. a value of 1 corresponds to points placed in 1 deg steps)
inp.res = 5;

% Show modiolus in plot
% ('show' || 'hide')
inp.modiolus = 'show';

% Show scala representation in plot
% ('show' || 'hide')
inp.scala = 'show';

%% Computation of cochlea spiral

% Computation of coefficients for radial and height polynomial functions
spir.rcoeff = coefs.modiolus'*[1;inp.meas'];
spir.hcoeff = coefs.height'*[1;inp.meas'];
spir.nturns = coefs.turns *[1;inp.meas']/360;

% Computation of corresponding point coordinates
spir.phi = (0:inp.res:spir.nturns*360);
spir.lwpts = [polyval(flipud(spir.rcoeff),spir.phi*pi/180).*cosd(spir.phi);...
              polyval(flipud(spir.rcoeff),spir.phi*pi/180).*sind(spir.phi);...
              polyval(flipud(spir.hcoeff)',spir.phi*pi/180)];

% Computation of scala profile
spir.scala.x = zeros(floor(360/inp.res),length(spir.phi));
spir.scala.y = zeros(floor(360/inp.res),length(spir.phi));
spir.scala.z = zeros(floor(360/inp.res),length(spir.phi));
for ii = 1:length(spir.phi)
    % Generate basic shape
    rtmp = (spir.nturns*360-spir.phi(ii))/(spir.nturns*360)*0.5+0.6;
    stmp = [norm(spir.lwpts(1:2,ii))+rtmp*(cosd(0:inp.res:360-inp.res)-1);...
            zeros(1,floor(360/inp.res));...
            spir.lwpts(3,ii)+rtmp*sind(0:inp.res:360-inp.res)];
    % Turn to corect position along spiral, save
    stmp = [cosd(spir.phi(ii)),-sind(spir.phi(ii)),0;...
            sind(spir.phi(ii)), cosd(spir.phi(ii)),0;...
            0,0,1]*stmp;
    spir.scala.x(:,ii) = stmp(1,:)';
    spir.scala.y(:,ii) = stmp(2,:)';
    spir.scala.z(:,ii) = stmp(3,:)';
end

clear rtmp stmp ii

%% 3D Plotting

figure;         %('WindowStyle','docked');
%figure('Position',[100 100 1380 820],'Color','black');      % Make the size for movie

%title('Towards Individualized Cochlear Implants','Color','white');
hold on
if strcmp(inp.modiolus,'show')
    plot3([0,0],[0,0],[spir.lwpts(3,1)-1,spir.lwpts(3,end)+1],...
        'color',[0.6,0.6,0.6],'LineWidth',5)
end
if strcmp(inp.scala,'show')
    mesh(spir.scala.x,spir.scala.y,spir.scala.z,...
        repmat(spir.phi,floor(360/inp.res),1),'FaceColor','none')
    colormap(flipud(winter))
end
plot3(spir.lwpts(1,:),spir.lwpts(2,:),spir.lwpts(3,:),'k','LineWidth',2)
xlabel('x [mm]'),ylabel('y [mm]'),zlabel('z [mm]')
axis equal
axis off

view(15,15)
turnturn(movie);
axis on
grid on

%axis off
%view(0,90);
%saveas(gcf,'figure.bmp');
end


function length = length_calculation_estimation_with_height(v,A,coefs)

syms z
x_deriv = subs(diff(radius_modiolus_estimation(A,coefs.modiolus)),'x',z)*cos(z) ...
    - subs(radius_modiolus_estimation(A,coefs.modiolus),'x',z)*sin(z);

y_deriv = subs(diff(radius_modiolus_estimation(A,coefs.modiolus)),'x',z)*sin(z) ...
    + subs(radius_modiolus_estimation(A,coefs.modiolus),'x',z)*cos(z);

z_deriv = subs(diff(height_estimation(A,coefs.height)),'x',z);

length_deriv = sqrt(x_deriv^2 + y_deriv^2 + z_deriv^2);

length = double(subs(int(length_deriv,0,z),z,v));

end


function length = length_calculation_estimation(v,A,coefs)

syms z
x_deriv = subs(diff(radius_modiolus_estimation(A,coefs.modiolus)),'x',z)*cos(z) ...
    - subs(radius_modiolus_estimation(A,coefs.modiolus),'x',z)*sin(z);

y_deriv = subs(diff(radius_modiolus_estimation(A,coefs.modiolus)),'x',z)*sin(z) ...
    + subs(radius_modiolus_estimation(A,coefs.modiolus),'x',z)*cos(z);

length_deriv = sqrt(x_deriv^2 + y_deriv^2);

length = double(subs(int(length_deriv,0,z),z,v));

end

function radius_modiolus = radius_modiolus_estimation(A,modiolus_coefs)

radius_modiolus = [ [1 A] * modiolus_coefs(:,1) ,...
    [1 A] * modiolus_coefs(:,2) ,...
    [1 A] * modiolus_coefs(:,3) ,...
    [1 A] * modiolus_coefs(:,4) ] ;

radius_modiolus = poly2sym(fliplr(radius_modiolus));

end

function height = height_estimation(A,height_coefs)

height = [ [1 A] * height_coefs(:,1) ,...
    [1 A] * height_coefs(:,2) ,...
    [1 A] * height_coefs(:,3) ,...
    [1 A] * height_coefs(:,4) ,...
    [1 A] * height_coefs(:,5) ] ;

height = poly2sym(fliplr(height));

end

function turns = turn_number_estimation(A,turns_coefs)

turns = ([1 A] * turns_coefs')/360 ;

end


function A = generate_cochlea_parameters(kind_of_cochlea)

switch kind_of_cochlea
    case 'mean'
        % these are the mean base measures in the study.
        mean_cochlea = [5.97,3.95,3.26,2.85];
        A = mean_cochlea;
    case 'random'
        % correlation matrix for random generation of base measures
        a1b1=0.53476;
        a1a2=-0.12441;
        a1b2=-0.07296;
        
        b1a2=0.11668;
        b1b2=-0.43748;
        
        a2b2=0.57846;
        
        Cor = [  1 , a1b1, a1a2, a1b2 ;...
            a1b1,   1 , b1a2, b1b2 ;...
            a1a2, b1a2,  1  , a2b2;...
            a1b2, b1b2, a2b2,  1 ];
        
        mean = 1;
        A = mvnrnd(mean*ones(1,size(Cor,1)),Cor);
        X = normcdf(A);
        
        A(1) = norminv(X(1),5.97,0.36);
        A(2) = norminv(X(2),3.95,0.35);
        A(3) = norminv(X(3),3.26,0.28);
        A(4) = norminv(X(4),2.85,0.33);
    otherwise
        error('Undefined kind of cochlea in generate_cochlea_parameters')
end
end

function turnturn(movie)
%
%   Turns the current figure 360° in the given elevation
%

fh=figure(gcf);
[azi,elev] = view;

l=gca;
set(l,'CameraViewAngle',10);

if movie==1
    cm=colormap;
    % 'compression','Indeo5' does not work in Windows 64bit
    %film=avifile('film.avi','colormap',cm,...
    %   'fps',10,'keyframe',10,'quality',100,'videoname','film');
    v=VideoWriter('film.avi');
    open(v);
end

for j=1:360
    
    view(j,elev);
    
    if (movie==1)
         M=getframe(fh);
         writeVideo(v,M);
    else
        pause(0.02)
    end
end

if movie==1
    close(v);
end
view(azi,elev)
end