% Script for plotting the Jansonius retinal nerve fiber bundle trajectory map.


% Described in and to be cited as:
%
% Jansonius NM, Nevalainen J, Selig B, Zangwill LM, Sample PA, Budde WM, Jonas JB,
% Lagreze WA, Airaksinen PJ, Vonthein R, Levin LA, Paetzold J, Schiefer U (2009)
% A mathematical description of nerve fiber bundle trajectories and their
% variability in the human retina. Vision Res 49:2157-­2163.
%
% ERRATUM in: Vision Res 50:1501 (2010) {this corrects an errror made by the publisher}
%
% Jansonius NM, Schiefer J, Nevalainen J, Paetzold J, Schiefer U (2012) A 
% mathematical model for describing the retinal nerve fiber bundle trajectories
% in the human eye: average course, variability, and influence of refraction,
% optic disc size and optic disc position. Exp Eye Res 105:70­-78.
%
% Qiu K, Schiefer J, Nevalainen J, Schiefer U, Jansonius NM (2015) Influence of the
% retinal blood vessel topography on the variability of the retinal nerve fiber bundle
% trajectories in the human retina. Invest Ophtahlmol Vis Sci 56:6320-5


% Latest change: July 10, 2017 by NJ


clear all;
close;

r0 = 4;
r = 4:0.1:30;

hemisphere = 1; % 1 = superior; 0 = inferiror

if (hemisphere == 1)

% superior

betas = -1.9; % -1.9 = mean value of ln b
ampls = 3.9;
fi0 = 60:10:180;

[rr,ffi0] = meshgrid(r',fi0');

ffihulp = (rr-r0).^(1.9+1.4*tanh((ffi0-121)/14));
ffi = ffihulp.*exp(betas+ampls*tanh(-(ffi0-121)/14))+ffi0;

xm = rr.*cos(ffi*2*pi/360);
ym = rr.*sin(ffi*2*pi/360);

x = xm+15;
yhulp1 = (x.^2)/225.0;
yhulp2 = (x.^2)/225.0;
yhulp3 = yhulp2.*sign(x);
yhulp4 = yhulp1.+yhulp3;
y = yhulp4.+ym;

figure(1)
plot(x',y'); axis([-20,20,0,20]); axis equal;

else

% inferior

betai = 0.7; % 0.7 = mean value of ln(-b)
ampli = 1.5;
fi0 = -60:-10:-180;

[rr,ffi0] = meshgrid(r',fi0');

ffihulp1 = (rr-r0).^(1.0+0.5*tanh((-1.0*ffi0-90)/25));
ffihulp2 = ffihulp1.*exp(betai+ampli*tanh(-(-1.0*ffi0-90)/25));
ffi = -1.0*ffihulp2+ffi0;

xm = rr.*cos(ffi*2*pi/360);
ym = rr.*sin(ffi*2*pi/360);

x = xm+15;
yhulp1 = (x.^2)/225.0;
yhulp2 = (x.^2)/225.0;
yhulp3 = yhulp2.*sign(x);
yhulp4 = yhulp1.+yhulp3;
y = yhulp4.+ym;

figure(1)
plot(x',y'); axis([-20,20,-20,0]); axis equal;

endif

