|
@@ -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
|
|
|
|
-
|
|
|