Delete 'code/Fig6 cone code/figureCone20190604.m'

code/Fig6 cone code/figureCone20190604.m

@@ -1,771 +0,0 @@
-function hd = coneAnalysis
-  
-%   datachoice = 'visual field';  % for main figure
-  datachoice = 'retina';  % for supplementary figure
-
-
-  warning('off','MATLAB:interp1:UsePCHIP')
-  % deactivating the following cbrewer warning:
-  % Warning: INTERP1(...,'CUBIC') will change in a future release. Use INTERP1(...,'PCHIP') instead. 
-  % > In interp1>sanitycheckmethod (line 265)
-  %   In interp1>parseinputs (line 401)
-  %   In interp1 (line 80)
-  %   In interpolate_cbrewer (line 31)
-  %   In cbrewer (line 101)
-  %   In heatmap20180823 (line 21) 
-  
-%   global globalparam
-  close all
-  
-  cd(['C:\Users\dehmelt\Dropbox\Documents\MATLAB\', ...
-      '201808 SphereAnalysis\fromTakeshiBadenlab_180615\5dpf_180628_correctedbug\'])
-  
-    
-  % Define which character sequences identify the different types of data file.
-  filestring.all   = 'allCones';
-  filestring.blue  = 'Blue';
-  filestring.green = 'Green';
-  filestring.red   = 'Red';
-  filestring.uv    = 'UV';
-  
-  
-  % Manually set max. cone density per square mm, as in Zimmermann et al.
-  maxdensity = 35000;
-  
-  % Manually select data files, and automatically check their names.
-  switch datachoice
-    case 'visual field'
-      fileselection = {'VisualFieldProjection_Red_8dpf.txt', ...
-                       'VisualFieldProjection_Green_8dpf.txt','VisualFieldProjection_Blue_8dpf.txt', ...
-                       'VisualFieldProjection_UV_8dpf.txt','VisualFieldProjection_allCones_8dpf.txt'};
-    case 'retina'
-      fileselection = {'DensityinRetina_Red_8dpf.txt', ...
-                       'DensityinRetina_Green_8dpf.txt','DensityinRetina_Blue_8dpf.txt', ...
-                       'DensityinRetina_UV_8dpf.txt','DensityinRetina_allCones_8dpf.txt'};
-  end
-  occurrence = struct('all',[],'blue',[],'green',[],'red',[],'uv',[]);
-  for file = 1:numel(fileselection)
-    occurrence.all   = [occurrence.all,   ~isempty(strfind(fileselection{file},filestring.all))];
-    occurrence.blue  = [occurrence.blue,  ~isempty(strfind(fileselection{file},filestring.blue))];
-    occurrence.green = [occurrence.green, ~isempty(strfind(fileselection{file},filestring.green))];
-    occurrence.red   = [occurrence.red,   ~isempty(strfind(fileselection{file},filestring.red))];
-    occurrence.uv    = [occurrence.uv,    ~isempty(strfind(fileselection{file},filestring.uv))];
-  end
-  
-  
-  % Check whether each type of file was selected exactly once.
-  if sum(occurrence.all) == 0
-    disp(['None of the file names contain the characters ''',filestring.all,'''.'])
-  elseif sum(occurrence.all) > 1
-    error(['More than one file name contains the characters ''',filestring.all,'''!'])
-  else
-    datafile.all = fileselection{find(occurrence.all)}; %#ok<FNDSB>
-  end
-  if sum(occurrence.blue) == 0
-    disp(['None of the file names contain the characters ''',filestring.blue,'''.'])
-  elseif sum(occurrence.blue) > 1
-    error(['More than one file name contains the characters ''',filestring.blue,'''!'])
-  else
-    datafile.blue = fileselection{find(occurrence.blue)}; %#ok<FNDSB>
-  end
-  if sum(occurrence.green) == 0
-    disp(['None of the file names contain the characters ''',filestring.green,'''.'])
-  elseif sum(occurrence.green) > 1
-    error(['More than one file name contains the characters ''',filestring.green,'''!'])
-  else
-    datafile.green = fileselection{find(occurrence.green)}; %#ok<FNDSB>
-  end
-  if sum(occurrence.red) == 0
-    disp(['None of the file names contain the characters ''',filestring.red,'''.'])
-  elseif sum(occurrence.red) > 1
-    error(['More than one file name contains the characters ''',filestring.red,'''!'])
-  else
-    datafile.red = fileselection{find(occurrence.red)}; %#ok<FNDSB>
-  end
-  if sum(occurrence.uv) == 0
-    disp(['None of the file names contain the characters ''',filestring.uv,'''.'])
-  elseif sum(occurrence.uv) > 1
-    error(['More than one file name contains the characters ''',filestring.uv,'''!'])
-  else
-    datafile.uv = fileselection{find(occurrence.uv)}; %#ok<FNDSB>
-  end
-
-
-
-
-%   % Fit centres from FitAzimElev20181217.txt:
-%   leftnormmaxazim = 74.1675;
-%   leftnormmaxelev = 5.9175;
-%   leftnormmeanazim = 74.9488;
-%   leftnormmeanelev = 5.7009;
-%   leftnormratioazim = 74.9490;
-%   leftnormratioelev = 5.7007;
-%   rightnormmaxazim = -72.0633;
-%   rightnormmaxelev = 7.3858;
-%   rightnormmeanazim = -71.1931;
-%   rightnormmeanelev = 7.5766;
-%   rightnormratioazim = -71.1933;
-%   rightnormratioelev = 7.5767;
-  leftazim_body = -80.2647;
-  leftelev_body = 6.0557;
-  rightazim_body = +77.0306;
-  rightelev_body = -1.9776;
-
-  % Median eye position, pooled once, from positionAnalysis20181217.m:
-  leftmedianazim =  +6.2450;
-  rightmedianazim = -9.8760;
-  leftmeanazim =  +5.2455;
-  rightmeanazim = -9.8840;
-  leftstdazim = 6.1703;
-  rightstdazim = 6.5227;
-  leftsemazim = 0.2472;
-  rightsemazim = 0.2613;
-
-  % Median vertical eye position, from figureVertical20190225.m:
-  leftmeanelev = 3.4937;  % meanmeanposition.UUH(:,:,1);
-  rightmeanelev = 4.8735; % meanmeanposition.UUH(:,:,2);
-  
-  % Corrected to eye-centred visual field angles:
-  leftazim_eye = abs(leftazim_body - leftmeanazim);
-  leftelev_eye = leftelev_body - leftmeanelev;
-  rightazim_eye = abs(rightazim_body - rightmeanazim);
-  rightelev_eye = rightelev_body - rightmeanelev;
-
-
-
-
-  % Create one figure panel per file, and all plots on it.
-  datapresent = fieldnames(datafile);
-  % Reorder field names to reflect desired arrangement of panels on figure.
-  datapresent = {datapresent{[4 3 2 5 1]}};
-  
-  for typeindex = 1:numel(datapresent)
-    
-    conetype = datapresent{typeindex};
-    textfile = datafile.(conetype);
-
-    switch datachoice
-      case 'visual field'
-        conedata = importdata(textfile) * maxdensity;
-      case 'retina'
-        conedata = importdata(textfile);
-    end
-    switch datachoice
-      case 'visual field'
-        conedata = rot90(conedata,3); % so up is at top, anterior to the left
-      case 'retina'
-        conedata = rot90(conedata,1); % so up is at bottom, anterior to the right  
-    end
-    conedata(conedata==0) = NaN; % consider nonzero data only, i.e. from the true area of the retina
-
-    
-    switch datachoice
-      case 'visual field'
-        % Undo approximate correction by Takeshi, then apply a more precise one:
-        offsetAzimTakeshi = 8;
-        offsetElevTakeshi = 5;
-        offsetAzimAccurate = abs(leftmedianazim);
-        offsetElevAccurate = leftmeanelev;
-        resolutionAzim = size(conedata,2)/360;
-        resolutionElev = size(conedata,1)/180;
-        discreteShiftAzim = round((offsetAzimTakeshi-offsetAzimAccurate)/resolutionAzim);
-        discreteShiftElev = round((offsetElevTakeshi-offsetElevAccurate)/resolutionElev);
-        nandata = NaN(size(conedata));
-        % Execute shift (min/max sanitise indices in case shift is negative):
-        nandata(max(1,1-discreteShiftElev):min(end,end-discreteShiftElev), ...
-                max(1,1-discreteShiftAzim):min(end,end-discreteShiftAzim)) = ...
-          conedata(max(1,1+discreteShiftElev):min(end,end+discreteShiftElev), ...
-                   max(1,1+discreteShiftAzim):min(end,end+discreteShiftAzim));
-        % Overwrite original data:
-        conedata = nandata;
-      case 'retina'
-        % Cannot apply correction, as VF filter used by Takeshi is unknown.
-    end
-
-
-  %   cd(['C:\Users\dehmelt\Dropbox\Documents\MATLAB\', ...
-  %       '201808 SphereAnalysis'])
-  % 
-  %   numelevation = size(conedata,1);
-  %   numazimuth = size(conedata,2);
-  % 
-  %   % The following four values represent the coordinates of the four corners
-  %   % of the 2D data sets provided by Baden lab (and are currently unknown).
-  %   elevation.max = 90;
-  %   elevation.min = -90;
-  %   azimuth.max = 180;
-  %   azimuth.min = -180;
-
-  %   % Build an (NX*NY x 2) array, where NX is the number of azimuth values,
-  %   % NY is the number of elevation values, NX*NY is the total number of data
-  %   % points, and for each data point (each row), the azimuth (column 1) and
-  %   % elevation (column 2) are listed.
-  %   elevation.step = (elevation.max-elevation.min)/(numelevation-1);
-  %   azimuth.step = (azimuth.max-azimuth.min)/(numazimuth-1);
-  %   elevationlist = elevation.max:-elevation.step:elevation.min;
-  %   azimuthlist = azimuth.max:-azimuth.step:azimuth.min;
-  %   predictor = [];
-  %   for k = 1:numazimuth
-  %     predictor = [predictor; ...
-  %                  repmat(azimuthlist(k),[numelevation 1]), elevationlist'];
-  %   end
-
-  %   % One raw data value per data point, might have to be reformatted to match
-  %   % the order in which the list of data azimuths/elevations in "predictor" is
-  %   % being created. Predictor goes top-to-bottom, left-to-right.
-  %   response = conedata(:);
-  % 
-  %   offset = 0;
-  %   amplitude = 1;
-  %   meanazimuth = 90;
-  %   meanelevation = 0;
-  %   concentration = 1;
-  %   
-  %   globalparam.offset = offset;
-  %   globalparam.amplitude = amplitude;
-  %   globalparam.concentration = concentration;
-  % 
-  % 
-  %   figure(1)
-  % %   imagesc(conedata)
-  %   surf(conedata)
-  %   shading flat
-  %   axis square
-  %   colormap(cbrewer('div','PRGn',100))
-  %   hold on
-  %   plot([1,numazimuth],[numelevation/2 numelevation/2],'--w','LineWidth',2)
-  %   plot([numazimuth/2, numazimuth/2],[1,numelevation],'--w','LineWidth',2)
-  %   hold off
-  %   set(gca,'XTick',[1,numazimuth*[1/4,1/2,3/4,1]],'XTickLabel',{180 90 0 -90 -180})
-  %   set(gca,'YTick',[1,numelevation*[1/4,1/2,3/4,1]],'YTickLabel',{90 45 0 -45 -90})
-  %   xlabel('azimuth'),ylabel('elevation'),set(gcf,'Color',[1 1 1])
-  %   dataclim = get(gca,'CLim');
-
-
-    hd.cone.fg = figure(1);
-    hd.cone.fg.Position = [400 160 1000 700];
-%     hd.cone.fg.Position = [500 50 540 380];
-    hd.cone.fg.Color = [1 1 1];
-
-    hd.cone.ax(typeindex,1) = axes;
-    if typeindex < 4
-      hd.cone.ax(typeindex,1).Position = [.11+(typeindex-1)*.29 .61 hd.cone.fg.Position([4 3])/3500];
-    else
-      hd.cone.ax(typeindex,1).Position = [.19+(typeindex-4)*.42 .11 hd.cone.fg.Position([4 3])/3500];
-    end
-    
-    az = 0:180/(size(conedata,1)-1):180;
-    el = -90:180/(size(conedata,2)-1):90;
-  %   conedatanonnan = fliplr(flipud(conedata));
-    conedatanonnan = conedata;
-    conedatanonnan(isnan(conedatanonnan)) = 0;
-    density = conv2(conedatanonnan,fspecial('gaussian',[100 100],3),'same');
-  %   surf(az,el,conv2(conedatanonnan,fspecial('gaussian',[100 100],10),'same'))
-  %   shading flat
-    hold on
-    if strcmp(datachoice,'retina') && typeindex == 5 % FD20190829 added as simple correction
-      density = fliplr(density)
-    end
-    imagesc(az,el,density)
-%     if strcmp(conetype,'all')
-%       contour(az,el,density,'LevelStep',1e-1,'Color',.5*[1 1 1])
-%     else
-%       contour(az,el,density,'LevelStep',3e-2,'Color',.5*[1 1 1])
-%     end
-    contourlevels = (.1:.1:1) * max(max(density));  % increments = 10%-of-max.
-    contour(az,el,density,contourlevels,'Color',.5*[1 1 1]); 
-    hold off
-%     axis square
-%     uistack(hd.cone.ax(typeindex,2),'bottom')
-%     hd.cone.ax(2,2) = copyobj(hd.cone.ax(2,1),hd.cone.fg)
-
-    hd.cone.ax(typeindex,1).XLim = [0,180];
-    hd.cone.ax(typeindex,1).YLim = [-90,90];
-    oldclim = hd.cone.ax(typeindex,1).CLim;
-%     climcutoff = .8;
-    climcutoff = .60;
-%     climcutoff = 0;
-    hd.cone.ax(typeindex,1).CLim = [min(oldclim)+climcutoff*(max(oldclim)-min(oldclim)), ...
-                            max(oldclim)];
-    switch datachoice
-      case 'visual field'
-        hd.cone.ax(typeindex,1).XLabel.String = 'azimuth';
-        hd.cone.ax(typeindex,1).YLabel.String = 'elevation';
-      case 'retina'
-        hd.cone.ax(typeindex,1).XLabel.String = 'x  (a.u.)';
-        hd.cone.ax(typeindex,1).YLabel.String = 'y  (a.u.)';
-    end
-    textfilesplit = strsplit(textfile,'_');
-    switch textfilesplit{2}
-      case 'allCones'
-        densitylabel = 'all cones combined';
-      case 'Green'
-        densitylabel = 'green cones (mm^{-2})';
-      case 'Red'
-        densitylabel = 'red cones (mm^{-2})';
-      case 'Blue'
-        densitylabel = 'blue cones (mm^{-2})';
-      case 'UV'
-        densitylabel = 'ultraviolet cones';
-    end
-%     switch textfilesplit{2}
-%       case 'allCones'
-%         densitylabel = {'retinal density','(any cones)'};
-%       case 'Green'
-%         densitylabel = {'retinal density / mm^{-2}','(green cones)'};
-%       case 'Red'
-%         densitylabel = {'retinal density','(red cones)'};
-%       case 'Blue'
-%         densitylabel = {'retinal density','(blue cones)'};
-%       case 'UV'
-%         densitylabel = {'retinal density','(ultraviolet cones)'};
-%     end
-%     switch textfilesplit{2}
-%       case 'allCones'
-%         densitylabel = 'density (any cones)';
-%       case 'Green'
-%         densitylabel = 'density (green cones)';
-%       case 'Red'
-%         densitylabel = 'density (red cones)';
-%       case 'Blue'
-%         densitylabel = 'density (blue cones)';
-%       case 'UV'
-%         densitylabel = 'density (ultraviolet cones)';
-%     end
-
-    hd.cone.ax(typeindex,1).Title.String = densitylabel;
-    hd.cone.ax(typeindex,1).Box = 'on';
-%     hd.cone.ax(typeindex,1).XTick = 0:15:180;
-%     hd.cone.ax(typeindex,1).YTick = -90:15:90;
-%     hd.cone.ax(typeindex,1).XTickLabel = {0,[],[],45,[],[],90,[],[],135,[],[],180};
-%     hd.cone.ax(typeindex,1).YTickLabel = {-90,[],[],-45,[],[],0,[],[],45,[],[],90};
-%     hd.cone.ax(typeindex,1).XTick = 0:45:180;
-%     hd.cone.ax(typeindex,1).YTick = -90:45:90;
-%     hd.cone.ax(typeindex,1).XTickLabel = {0,45,90,135,180};
-%     hd.cone.ax(typeindex,1).YTickLabel = {-90,-45,0,45,90};
-    hd.cone.ax(typeindex,1).XTick = 0:30:180;
-    hd.cone.ax(typeindex,1).YTick = -90:30:90;
-    switch datachoice
-      case 'visual field'
-        hd.cone.ax(typeindex,1).XTickLabel = {0,[],[],90,[],[],180};
-        hd.cone.ax(typeindex,1).YTickLabel = {-90,[],[],0,[],[],90};
-      case 'retina'
-        hd.cone.ax(typeindex,1).XTickLabel = {0,[],[],[],[],[],1};
-        hd.cone.ax(typeindex,1).YTickLabel = {0,[],[],[],[],[],1};
-    end
-    hd.cone.ax(typeindex,1).FontSize = 14;
-    hd.cone.ax(typeindex,1).TitleFontSizeMultiplier = 1;
-    hd.cone.ax(typeindex,1).TitleFontWeight = 'bold';
-%     hd.cone.ax(typeindex,1).Title.Interpreter = 'none';
-    hd.cone.ax(typeindex,1).Title.Interpreter = 'tex';
-  %   [hd.cone.ax(typeindex,1).XLabel.FontSize, ...
-  %    hd.cone.ax(typeindex,1).YLabel.FontSize, ...
-  %    hd.cone.ax(typeindex,1).Title.FontSize, ...
-  %    ha.cone.ax(typeindex,1).FontSize,
-    hd.cone.ax(typeindex,1).XGrid = 'on';
-    hd.cone.ax(typeindex,1).YGrid = 'on';
-    hd.cone.ax(typeindex,1).Layer = 'top';
-    
-    
-    hold on
-    plot([0,180],[0,0],'--k')
-    plot([90,90],[-90,90],'--k')
-    hold off
-
-   
-    
-% %     leftnormmaxazimVF = abs(leftnormmaxazim - leftmedianazim);
-% %     leftnormmeanazimVF = abs(leftnormmeanazim - leftmedianazim);
-% %     leftnormratioazimVF = abs(leftnormratioazim - leftmedianazim);
-% %     rightnormmaxazimVF = abs(rightnormmaxazim - rightmedianazim);
-% %     rightnormmeanazimVF = abs(rightnormmeanazim - rightmedianazim);
-% %     rightnormratioazimVF = abs(rightnormratioazim - rightmedianazim);
-
-
-
-    hold on
-    hd.cone.ob(1,1) = plot([153,165],[-65,-65],'k');
-    hd.cone.ob(1,2) = plot([159,159],[-59,-71],'k');
-    switch datachoice
-      case 'visual field'
-        hd.cone.ob(1,3) = text(153,-64.5,'A','HorizontalAlignment','right','VerticalAlignment','middle');
-        hd.cone.ob(1,4) = text(166,-64.5,'P','HorizontalAlignment','left','VerticalAlignment','middle');
-        hd.cone.ob(1,5) = text(159.2,-59.5,'U','HorizontalAlignment','center','VerticalAlignment','bottom');
-        hd.cone.ob(1,6) = text(159.2,-69.5,'D','HorizontalAlignment','center','VerticalAlignment','top');
-      case 'retina'
-        hd.cone.ob(1,3) = text(153,-64.5,'P','HorizontalAlignment','right','VerticalAlignment','middle');
-        hd.cone.ob(1,4) = text(166,-64.5,'A','HorizontalAlignment','left','VerticalAlignment','middle');
-        hd.cone.ob(1,5) = text(159.2,-59.5,'D','HorizontalAlignment','center','VerticalAlignment','bottom');
-        hd.cone.ob(1,6) = text(159.2,-69.5,'U','HorizontalAlignment','center','VerticalAlignment','top');
-    end
-    hold off
-    
-%     circle.radius = 40/2;
-%     circle.xcentre = leftnormmaxazimVF; % pick preferred measure here
-%     circle.ycentre = leftnormmaxelev; % pick preferred measure here
-%     circle.x = -circle.radius:.1:circle.radius;
-%     circle.y = sqrt(circle.x(1)^2-circle.x.^2);
-%     hold on
-%     plot([circle.xcentre + circle.x, fliplr(circle.xcentre + circle.x)], ...
-%          [circle.ycentre + circle.y, circle.ycentre - circle.y],'-','LineWidth',3,'Color',0*[1 1 1])    
-%     plot([circle.xcentre + circle.x, fliplr(circle.xcentre + circle.x)], ...
-%          [circle.ycentre + circle.y, circle.ycentre - circle.y],'-','LineWidth',2,'Color',1*[1 1 1])    
-%     hold off
-    
-%     boundary.angle = 64;
-% %     centre.azimuth = leftnormmaxazimVF;
-% %     centre.elevation = leftnormmaxelev;
-%     centre.azimuth = rightazim;
-%     centre.elevation = rightelev;
-%     boundary.resolution = .01;
-%     [boundary.azimuth,boundary.elevation] = ...
-%       diskBoundary(centre.azimuth,centre.elevation,boundary.angle,boundary.resolution);
-%     hold on
-%     plot(boundary.azimuth,boundary.elevation,'-','LineWidth',3,'Color',0*[1 1 1])    
-%     plot(boundary.azimuth,boundary.elevation,'-','LineWidth',2,'Color',1*[1 1 1])    
-%     hold off
-
-    boundary.angle = 40;
-%     centre.azimuth = leftnormmaxazimVF;
-%     centre.elevation = leftnormmaxelev;
-%     centre.azimuth = leftazim_body;
-%     centre.elevation = leftelev_body;
-    centre.azimuth = abs(leftazim_body);
-    centre.elevation = leftelev_body;
-    switch datachoice
-      case 'retina'
-        centre.azimuth = 180-centre.azimuth;
-        centre.elevation = -centre.elevation;
-    end
-    switch datachoice
-      case 'visual field'
-        boundary.resolution = .01;
-        [boundary.azimuth,boundary.elevation] = ...
-          diskBoundary(centre.azimuth,centre.elevation,boundary.angle,boundary.resolution);
-        hold on
-        plot(boundary.azimuth,boundary.elevation,'-','LineWidth',3,'Color',0*[1 1 1])    
-        plot(boundary.azimuth,boundary.elevation,'-','LineWidth',2,'Color',1*[1 1 1])    
-        hold off
-    end
-
-
-    mtp = 'o';
-    msz = 100;
-    mec = [0 0 0];
-%     mfc = [.2 .8 .2];
-%     mfa = .2;
-    mfc = 1*[1 1 1];
-    mfa = .5;
-    hold on
-%     hd.cone.ob(1,7) = scatter(rightazim,rightelev,msz,mtp, ...
-%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
-    hd.cone.ob(1,8) = scatter(centre.azimuth,centre.elevation,msz,mtp, ...
-                               'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',.8);
-%     hd.cone.ob(1,9) = scatter(leftnormratioazimVF,leftnormratioelev,msz,mtp, ...
-%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
-%     hd.cone.ob(1,10) = scatter(rightnormmaxazimVF,rightnormmaxelev,msz,mtp, ...
-%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
-%     hd.cone.ob(1,11) = scatter(rightnormmeanazimVF,rightnormmeanelev,msz,mtp, ...
-%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
-%     hd.cone.ob(1,12) = scatter(rightnormratioazimVF,rightnormratioelev,msz,mtp, ...
-%                                 'MarkerEdgeColor',mec,'MarkerFaceColor',mfc,'MarkerFaceAlpha',mfa);
-    hold off
-    
-    
-    
-    if typeindex == 2 || typeindex == 3
-      hd.cone.ax(typeindex,1).YLabel.String = [];
-      hd.cone.ax(typeindex,1).YTickLabel = [];
-    end
-    switch conetype
-      case 'blue'
-        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Blues',100))
-      case 'green'
-        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Greens',100))
-      case 'red'
-        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Reds',100))
-      case 'uv'
-        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Purples',100))
-      case 'all'
-        colormap(hd.cone.ax(typeindex,1),cbrewer('seq','Greys',100))
-    end
-    
-    hd.cone.ax(typeindex,2) = axes;
-    hd.cone.ax(typeindex,2).Position = hd.cone.ax(typeindex,1).Position + ...
-                                       hd.cone.ax(typeindex,1).Position(3) * .28*[1 0 0 0];
-    hd.cone.cbar(typeindex) = colorbar;
-
-    hd.cone.ax(typeindex,2).Colormap = hd.cone.ax(typeindex,1).Colormap;
-    hd.cone.ax(typeindex,2).CLim     = hd.cone.ax(typeindex,1).CLim;
-    hd.cone.ax(typeindex,2).FontSize = hd.cone.ax(typeindex,1).FontSize;
-    hd.cone.ax(typeindex,2).Color = 'none';
-    [hd.cone.ax(typeindex,2).XAxis.Visible,hd.cone.ax(typeindex,2).YAxis.Visible] = deal('off');
-        
-        
-%     colormap(cbrewer('div','PRGn',100))
-
-
-
-  end
-  
-  hd.cone.ax(1,3) = axes;
-  hd.cone.ax(1,3).Position = [8 0 1 1];
-  hd.cone.ax(1,3).Color = 'none';
-  [hd.cone.ax(1,3).XAxis.Visible,hd.cone.ax(1,3).YAxis.Visible] = deal('off');
-  for typeindex = 1:numel(datapresent)
-    parentpos = hd.cone.ax(typeindex,1).Position;
-    xpos = parentpos(1)+parentpos(3)/2;
-    ypos = parentpos(2)+parentpos(4)*1.128;
-    hd.cone.ob(3,typeindex) = text(xpos, ypos, 'retinal density / mm^{-2}');
-%     hd.cone.ob(3,typeindex) = text(xpos, ypos, 'retinal density (µm^{-2})');
-    hd.cone.ob(3,typeindex).FontWeight = 'normal';
-    hd.cone.ob(3,typeindex).Interpreter = 'tex';
-    hd.cone.ob(3,typeindex).FontSize = hd.cone.ax(typeindex,1).FontSize;
-    hd.cone.ob(3,typeindex).HorizontalAlignment = 'center';
-    hd.cone.ob(3,typeindex).VerticalAlignment = 'bottom';
-  end
-  uistack(hd.cone.ax(1,3),'bottom')
-
-  hd.cone.ax(1,4) = axes;
-  hd.cone.ax(1,4).Position = [0 0 1 1];
-  hd.cone.ax(1,4).Color = 'none';
-  [hd.cone.ax(1,4).XAxis.Visible,hd.cone.ax(1,4).YAxis.Visible] = deal('off');
-  for typeindex = 1:numel(datapresent)
-    parentpos = hd.cone.ax(typeindex,1).Position;
-%     xpos = parentpos(1)+parentpos(3)*-.18;
-    xpos = parentpos(1)+parentpos(3)*-.23;
-    ypos = parentpos(2)+parentpos(4)*1.30;
-%     hd.cone.ob(3,typeindex) = text(xpos, ypos, 'retinal density / mm^{-2}');
-    switch typeindex
-      case 1
-        letterstring = 'a';
-      case 4
-        letterstring = 'b';
-      case 5
-        letterstring = 'c';
-      otherwise
-        letterstring = '';  
-    end
-    hd.cone.ob(4,typeindex) = text(xpos, ypos, letterstring);
-    hd.cone.ob(4,typeindex).FontWeight = 'bold';
-    hd.cone.ob(4,typeindex).Interpreter = 'tex';
-    hd.cone.ob(4,typeindex).FontSize = hd.cone.ax(typeindex,1).FontSize * 18/14;
-    hd.cone.ob(4,typeindex).HorizontalAlignment = 'right';
-    hd.cone.ob(4,typeindex).VerticalAlignment = 'top';
-  end
-  uistack(hd.cone.ax(1,3),'bottom')
-
-  [hd.cone.ax(:,1).TickLength] = deal([0 0]);
-
-  
-%   for k = 1:5, hd.cone.cbar(k,1).Position(1) = hd.cone.ax(k,1).Position(1) + .83*hd.cone.ax(k,1).Position(3); end
-  
-    runtime = clock;
-    timestring = ['_', ...
-                num2str(runtime(1),'%.4u'), ...
-                num2str(runtime(2),'%.2u'), ...
-                num2str(runtime(3),'%.2u'), ...
-                '_', ...
-                num2str(runtime(4),'%.2u'), ...
-                num2str(runtime(5),'%.2u'), ...
-                num2str(runtime(6),'%02.0f')];
-                            
-%   savefig(2,['conefigure/',textfilesplit{2},timestring])
-%   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-depsc')
-% %   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-dtiff')
-%   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-dtiffn')
-%   figure(2),print(['conefigure/',textfilesplit{2},timestring],'-dpng')
-  
-end
-
-
-
-function [azimuthOut, elevationOut] = diskBoundary(azimuthIn, elevationIn, angleIn, resolutionIn)
-
-  centre.azimuth = azimuthIn;
-  centre.elevation = elevationIn;
-  border.angle = angleIn / 2; % convert solid angle ("diameter") to "radius"
-  border.resolution = resolutionIn;
-
-  % The following equation describes orthodromic (or "great circle") distance as an angle:
-  % arccos(sin(centre.elevation)*sin(border.elevation) + ...
-  %        cos(centre.elevation)*cos(border.elevation)*cos(border.azimuth-centre.azimuth)) = angle;
-
-%   border.elevation = centre.elevation-border.angle/2 : .1 : centre.elevation+border.angle/2;
-  intended  = -90 : border.resolution : 90;
-%   intended  = (centre.elevation-border.angle) : border.resolution : (centre.elevation+border.angle);
-  % "auxiliary" elevations are explicitly added to guarantee coverage near top/bottom of circle:
-  auxiliary = centre.elevation + border.angle*[-1 1];
-%   border.elevation = unique([intended,auxiliary]);
-  
-  limit = cosd([centre.elevation-border.angle centre.elevation+border.angle]);
-%   intended = acosd(min(limit) : border.resolution : max(limit));
-
-  %   % XXX Not sure why the following line leads to half of the curve segments not being computed...
-%   intended = [-acosd(0:border.resolution:1),acosd(1-border.resolution:-border.resolution:0)];
-%   auxiliary = centre.elevation + border.angle*[-1 1];
-
-  border.elevation = unique([intended,auxiliary]);
-%   figure(6)
-%   plot(intended)
-
-  border.azimuth = centre.azimuth + ...
-                   acosd((cosd(border.angle) - sind(centre.elevation)*sind(border.elevation)) ./ ...
-                         (cosd(centre.elevation)*cosd(border.elevation)));
-
-
-  % Eliminate complex elements...
-  border.azimuth(imag(border.azimuth)~=0) = NaN; % This statement sometimes fails horribly.
-%   border.azimuth(imag(border.azimuth)>1e-5) = NaN; % This statement is silly, but more stable.
-
-  % ...and compute corresponding values for the second half of the boundary.
-  border.azimuth = [border.azimuth, 2*centre.azimuth-border.azimuth];
-  border.elevation = [border.elevation, border.elevation];
-  
-  border.azimuth = mod(border.azimuth+180,360)-180;
-%   jump = diff(border.azimuth) > 10*border.resolution;
-%   border.azimuth(jump)   = NaN;
-%   border.elevation(jump) = NaN;
-  
-  azimuthOut   = border.azimuth;
-  elevationOut = border.elevation;
-
-end
-
-
-
-% function combinedvalue = vonmisesfisher2param_centre(param,predictor)
-% 
-%   global globalparam
-%   
-%   convention = 'CCW';
-%   dimension = 3;
-%   radius = 1;
-%   if size(predictor,2) == 2
-%     predictor(:,3) = radius*ones(size(predictor(:,1)));
-%   end
-%   [x,y,z] = geo2car(predictor(:,1),predictor(:,2),predictor(:,3),convention);
-%   direction = [x,y,z];
-%   
-%   numberofdistributions = 1; % <---
-%   value = NaN(size(predictor,1),numberofdistributions);
-%   
-%   for k = 1:numberofdistributions
-%     
-%     meanazimuth   = param(2*k-1);
-%     meanelevation = param(2*k);
-%     
-%     offset        = globalparam.offset(k);
-%     amplitude     = globalparam.amplitude(k);
-%     concentration = globalparam.concentration(k);
-%     
-%     [meanx,meany,meanz] = geo2car(meanazimuth,meanelevation,radius,convention);
-%     meandirection = [meanx,meany,meanz];
-% 
-%     bessel = concentration / ...
-%              (2*pi * (exp(concentration) - exp(-concentration)));
-% 
-%     normalisation = concentration^(dimension/2 - 1) / ...
-%                     ((2*pi)^(dimension/2) * bessel);
-%                 
-%     value(:,k) = normalisation * exp(concentration*meandirection*direction')';
-%     value(:,k) = amplitude * value(:,k) + offset;
-%    
-%     badness = fitpenalty(meanazimuth,-180,180) + ...
-%               fitpenalty(meanelevation,-90,90) + ...
-%               fitpenalty(amplitude,0,Inf)      + ...
-%               fitpenalty(concentration,0,Inf);
-%             
-%     value(:,k) = value(:,k) + 1e3*badness;
-%     
-%   end
-%   
-%   combinedvalue = sum(value,2);
-%   
-% end
-
-
-
-% function penalty = fitpenalty(value,lowerlimit,upperlimit)
-%   
-%   factor = 1;
-%   
-%   if value < lowerlimit
-%     
-%     penalty = factor*(lowerlimit-value).^1;
-%     
-%   elseif value > upperlimit
-%     
-%     penalty = factor*(value-upperlimit).^1;
-%     
-%   else
-%     
-%     penalty = 0;
-%     
-%   end
-%   
-% end
-
-
-
-% function combinedvalue = oneVonMisesFisher(params,predictor)
-% %TWOVONMISESFISHER Evaluate two overlapping von Mises-Fisher distributions.
-% %
-% %   VALUE = VONMISESFISHER computes the value of only one(!) 
-% %   spherical (p=3) von Mises-Fisher distributions at the chosen location.
-% %
-% %   It is formatted to match the input/output structure required by the
-% %   built-in non-linear regression function nlinfit.
-% 
-%   convention = 'CW';
-%   dimension = 3;
-%   radius = 1;
-%   if size(predictor,2) == 2
-%     predictor(:,3) = radius*ones(size(predictor(:,1)));
-%   end
-%   [x,y,z] = geo2car(predictor(:,1),predictor(:,2),predictor(:,3),convention);
-%   direction = [x,y,z];
-%   
-%   numberofdistributions = 1; % <-- This is the only difference to "twoVonMisesFisher.m".
-%   value = NaN(size(predictor,1),numberofdistributions);
-%   
-%   for k = 1:numberofdistributions
-%     
-%     offset        = params(5*k-4);
-%     amplitude     = params(5*k-3);
-%     meanazimuth   = params(5*k-2);
-%     meanelevation = params(5*k-1);
-%     
-%         meanelevation = 25; % HACK!
-% 
-%     
-%     [meanx,meany,meanz] = geo2car(meanazimuth,meanelevation,radius,convention);
-%     meandirection = [meanx,meany,meanz];
-% %     meandirection = geo2car([params((-2:-1)+4*k),radius]);
-%     concentration = params(5*k);
-%     
-%     bessel = concentration / ...
-%              (2*pi * (exp(concentration) - exp(-concentration)));
-% 
-%     normalisation = concentration^(dimension/2 - 1) / ...
-%                     ((2*pi)^(dimension/2) * bessel);
-%                 
-%     value(:,k) = normalisation * exp(concentration*meandirection*direction')';
-%     value(:,k) = amplitude * value(:,k) + offset;
-%     
-%     % HACK to constrain values
-%     if abs(meanazimuth) > 180  || ...
-%        abs(meanelevation) > 90 || ...
-%        amplitude < 0           || ...
-%        concentration < 0
-%      
-%       value(:,k) = 1e10*ones(size(value(:,k)));
-%       
-%     end
-% 
-%   end
-%   
-%   combinedvalue = sum(value,2);
-%   
-% end
-