|
|
@@ -0,0 +1,771 @@
|
|
|
+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
|
|
|
+
|
