clear all close all grating = displayPattern(360/22); velocity = displayCosine(12.5); [leftdisk,leftparam] = displayDiskStimuli('left'); [rightdisk,rightparam] = displayDiskStimuli('right'); singledisk = displayDiskStimuli('single'); [diskhandle,diskdata,fitdata] = gainAnalysisDisk('Resolution',50); variant = 3; % 1 = regular, 2 = rotated, 3 = combined/corrected %%%%% !!!! %%%%% hd.collage.fg = figure(); hd.collage.fg.Name = 'Sphere manuscript (Fig 6); OKR gain dependence on stimulus location'; hd.collage.fg.Color = [1 1 1]; hd.collage.fg.Position = [0 40 1070 960]; % Populate the combined figure from previously generated partial figures... hd.collage.ax(1,1) = copyobj(grating.pattern.ax(1), hd.collage.fg); hd.collage.ax(1,2) = copyobj(velocity.cosine.ax(1), hd.collage.fg); hd.collage.ax(1,3) = copyobj(leftdisk.default.ax(1), hd.collage.fg); hd.collage.ax(1,4) = copyobj(rightdisk.default.ax(1), hd.collage.fg); hd.collage.ax(1,5) = copyobj(singledisk.default.ax(1), hd.collage.fg); % hd.collage.ax(2,1) = copyobj(diskdata.gain.ax(7), hd.collage.fg); hd.collage.ax(2,1) = axes; for subpanel = 2:4 hd.collage.ax(2,subpanel) = copyobj(diskhandle.gain(variant).ax(subpanel-1), hd.collage.fg); end for subpanel = 5:7 hd.collage.ax(2,subpanel) = copyobj(hd.collage.ax(2,subpanel-3), hd.collage.fg); end % ...get rid of the partial figures once pillaged... % close(1:7) % ...then adjust the size, position and other properties of those objects. [hd.collage.ax(1,1:5).Units, hd.collage.ax(2,1:7).Units] = deal('pixels'); hd.collage.ax(1,1).Position = [100 850 180 70]; hd.collage.ax(1,2).Position = [130 730 140 55]; hd.collage.ax(1,3).Position = [350 730 200 200]; hd.collage.ax(1,4).Position = hd.collage.ax(1,3).Position + [250 0 0 0]; hd.collage.ax(1,5).Position = hd.collage.ax(1,3).Position + [520 0 0 0]; hd.collage.ax(2,1).Position = [ 50 240 50 160]; hd.collage.ax(2,2).Position = [ 90 280 350 350]; for subpanel = 3:4 hd.collage.ax(2,subpanel).Position = hd.collage.ax(2,subpanel-1).Position + [320 0 0 0]; end for subpanel = 5:7 hd.collage.ax(2,subpanel).Position = hd.collage.ax(2,subpanel-3).Position - [0 280 0 0]; end darkgreyscale = (.8:.01:.99)' * [1 1 1]; lightgreyscale = (.9:.01:.99)' * [1 1 1]; [hd.collage.ax(1,3:4).Colormap] = deal(lightgreyscale); [hd.collage.ax(1,5).Colormap] = deal(darkgreyscale); % % Shared properties of all axis objects anywhere: [hd.collage.ax(1,1:5).FontSize, ... hd.collage.ax(2,1:7).FontSize] = deal(15); [hd.collage.ax(1,1:5).LabelFontSizeMultiplier, ... hd.collage.ax(1,1:5).TitleFontSizeMultiplier, ... hd.collage.ax(2,1:7).LabelFontSizeMultiplier, ... hd.collage.ax(2,1:7).TitleFontSizeMultiplier] = deal(1); % Create content of first subpanel of panel b (colorbar): hd.collage.ax(2,1).CLim = hd.collage.ax(2,2).CLim; hd.collage.ob(2,1,1) = colorbar; hd.collage.ob(2,1,1).Label.String = 'OKR gain'; hd.collage.ob(2,1,1).Label.Units = 'normalized'; hd.collage.ob(2,1,1).AxisLocation = 'in'; hd.collage.ax(2,1).Visible = 'off'; hd.collage.ob(2,1,1).Units = 'pixels'; hd.collage.ob(2,1,1).Position(3) = 9; hd.collage.ob(2,1,1).FontSize = hd.collage.ax(2,1).FontSize; hd.collage.ob(2,1,1).Label.Position = [-.8 1 1] .* hd.collage.ob(2,1,1).Label.Position; % Add and format titles to panel a: for panel = 1 hd.collage.ax(panel,1).Title.String = 'stimulus pattern'; hd.collage.ax(panel,3).Title.String = '----------------------- stimulus crop -----------------------'; hd.collage.ax(panel,3).Title.Interpreter = 'Latex'; % hd.collage.ax(panel,4).Title.String = {'stimulus crop','(right hemisphere)'}; hd.collage.ax(panel,2).XLabel.String = 'time (sec)'; hd.collage.ax(panel,2).YLabel.String = {'velocity','(deg/sec)'}; hd.collage.ax(panel,3).XLabel.String = 'left hemisphere'; hd.collage.ax(panel,4).XLabel.String = 'right hemisphere'; [hd.collage.ax(panel,3).XLabel.Visible, ... hd.collage.ax(panel,4).XLabel.Visible] = deal('on'); hd.collage.ax(panel,5).Title.String = 'stimulus'; for subpanel = 1:5 hd.collage.ax(panel,subpanel).Title.FontWeight = 'normal'; [hd.collage.ax(panel,subpanel).Title.Units, ... hd.collage.ax(panel,subpanel).XLabel.Units] = deal('pixels'); end % hd.collage.ax(panel,1).Title.Position(2) = hd.collage.ax(panel,1).Position(4) - 14; hd.collage.ax(panel,1).XLabel.Position(2) = hd.collage.ax(panel,1).XLabel.Position(2) + 25; % hd.collage.ax(panel,2).XLabel.Position(2) = hd.collage.ax(panel,2).XLabel.Position(2) - 65; hd.collage.ax(panel,2).XLabel.Position(2) = hd.collage.ax(panel,2).XLabel.Position(2) + 25; hd.collage.ax(panel,3).Title.Position(1) = hd.collage.ax(1,3).Title.Position(1) + 126; for subpanel = 3:4 hd.collage.ax(panel,subpanel).Title.Position(2) = hd.collage.ax(panel,subpanel).Position(4) - 2; hd.collage.ax(panel,subpanel).XLabel.Position(2) = hd.collage.ax(panel,subpanel).XLabel.Position(2) + 23; end % hd.collage.ax(panel,5).Title.Position(2) = hd.collage.ax(panel,4).Position(4) + 2; end % axes(hd.collage.ax(1,3)) % leftline = text(hd.collage.ax(1,3).Title.Position(1) - 50, ... % hd.collage.ax(1,3).Title.Position(2) + 20, ... % '-----------------------','Interpreter','Latex','HorizontalAlignment','right','Units','pixels') % Force visibility of axis ticks on bottom left (sine) plot of panel a: % hd.collage.ax(1,2).XLim = [-13 13]; hd.collage.ax(1,2).YLim = [-13 13]; % Format titles on top row of panels b to d: for subpanel = 2:4 hd.collage.ax(2,subpanel).Title.Units = 'pixels'; hd.collage.ax(2,subpanel).Title.Position(2) = hd.collage.ax(2,subpanel).Title.Position(2) - 25; end % Correct for odd treatment of indices (or their absence): hd.collage.ax(2,3).Title.Position(2) = hd.collage.ax(2,3).Title.Position(2) + 9; % Remove titles from bottom row of panels b to d: for subpanel = 5:7 hd.collage.ax(2,subpanel).Title.Visible = 'off'; end % % Make fish on panels b to d more visible by adjusting sphere alpha: % [hd.collage.ax( % Adjust camera perspective on the lower row of panels b to d: [hd.collage.ax(2,5:7).CameraPosition] = deal([19.5 0 0]); % Display panel lettering: hd.collage.ax(4,1) = axes; hd.collage.ax(4,1).Color = 'none'; hd.collage.ax(4,1).Position = [0 0 1 1]; hd.collage.ax(4,1).Units = 'pixels'; hd.collage.ax(4,1).Visible = 'off'; hold on hd.collage.lt(1) = text(hd.collage.ax(1,1).Position(1) - 85 , ... hd.collage.ax(1,1).Position(2) + hd.collage.ax(1,1).Position(4) + 20, ... 'a', 'Units', 'pixels'); hd.collage.lt(2) = text(hd.collage.ax(2,2).Position(1) - 50, ... hd.collage.ax(2,2).Position(2) + hd.collage.ax(2,2).Position(4) - 10, ... 'b', 'Units', 'pixels'); hd.collage.lt(3) = text(hd.collage.ax(2,3).Position(1) + 35, ... hd.collage.ax(2,3).Position(2) + hd.collage.ax(2,3).Position(4) - 10, ... 'c', 'Units', 'pixels'); hd.collage.lt(4) = text(hd.collage.ax(2,4).Position(1) + 35, ... hd.collage.ax(2,4).Position(2) + hd.collage.ax(2,4).Position(4) - 10, ... 'd', 'Units', 'pixels'); hold off [hd.collage.lt(:).FontSize] = deal(15); [hd.collage.lt(:).FontWeight] = deal('bold'); % plotAsymmetry(diskdata(variant)) % mercatorPanelOldStimuli20190527 mercatorPanelOldStimuli20190616 delete(hd.collage.ax(1,1:2)) % remove panels about basic stimulus pattern close(1:8) for subpanel = 3:4 hd.collage.ax(1,subpanel).Position = hd.collage.ax(1,subpanel).Position - [200 0 0 0]; end hd.collage.ax(1,5).Position = hd.collage.ax(1,5).Position - [100 0 0 0]; % [hd.collage.ax(1:2,1:2).Units,hd.collage.ax(3:5,:).Units] = deal('pixels'); % [hd.collage.ax(1:2,1:2).FontSize,hd.collage.ax(3:5,:).FontSize] = deal(15); % [hd.collage.ax(1:2,1:2).LabelFontSizeMultiplier, ... % hd.collage.ax(1:2,1:2).TitleFontSizeMultiplier, ... % hd.collage.ax(3:5,:).LabelFontSizeMultiplier, ... % hd.collage.ax(3:5,:).TitleFontSizeMultiplier] = deal(1); function plotAsymmetry(diskdata) col.barface = [.2 .8 .4]; col.neutralbarface = .6*[1 1 1]; col.baredge = 'none'; col.errorbar = .3*[1 1 1]; hd.asymindex.fg = figure(); hd.asymindex.fg.Color = [1 1 1]; hd.asymindex.fg.Position = [0 100 1900 390]; hd.asymindex.ax(1,1) = axes(); % hd.asymindex.ax(1,2) = axes(); [hd.asymindex.ax(:,:).Units] = deal('pixels'); hd.asymindex.ax(1,1).Position = [90 50 1780 300]; y5a = diskdata.asymmetry.A1; y5b = diskdata.asymmetry.A2; % e5a = mscdisk.botheyes.stdgain; % e5b = upside.botheyes.stdgain; e5a = NaN(size(y5a)); e5b = NaN(size(y5b)); pvalue.asymmetry = signrank(y5a,y5b); x5a = (1:numel(y5a))' - .2; x5b = (1:numel(y5b))' + .2; axes(hd.asymindex.ax(1,1)) hold on hd.asymindex.barregular = bar(x5a,y5a); hd.asymindex.barcontrol = bar(x5b,y5b); hd.asymindex.errorregular = errorbar(x5a,y5a,e5a,'LineStyle','none'); hd.asymindex.errorcontrol = errorbar(x5b,y5b,e5b,'LineStyle','none'); hd.asymindex.explanation = text(24,.351,{'biological asymmetry shown in green', ... '(error bars show s.e.m.)'},'HorizontalAlignment','center'); hd.asymindex.pvalue = text(36,.351,{['p = ',num2str(pvalue.asymmetry)], ... '(Wilcoxon signed-rank test)'},'HorizontalAlignment','center'); hold off [hd.asymindex.ax(1,:).XTick] = deal([]); hd.asymindex.barregular.FaceColor = col.neutralbarface; hd.asymindex.barcontrol.FaceColor = col.barface; [hd.asymindex.barregular.BarWidth, ... hd.asymindex.barcontrol.BarWidth] = deal(.35); [hd.asymindex.errorregular.Color, ... hd.asymindex.errorcontrol.Color] = deal(col.errorbar); [hd.asymindex.errorregular.LineWidth, ... hd.asymindex.errorcontrol.LineWidth] = deal(1); hd.asymindex.ax(1,1).XLim = [7 46]; % hd.asymindex.ax(1,1).YLim = hd.asymindex.ax(1,2).YLim; hd.asymindex.ax(1,1).YLabel.String = 'mean OKR gain (pooled across eyes)'; % hd.asymindex.ax(1,2).YLabel.String = 'mean OKR gain (right eye)'; for panel = 1:1 hd.asymindex.ax(1,panel).YLabel.Units = 'pixels'; hd.asymindex.ax(1,panel).YLabel.Position(1) = hd.asymindex.ax(1,panel).YLabel.Position(1) - 20; end hd.asymindex.ax(1,1).XLabel.String = 'stimulus types'; % hd.asymindex.ax(1,2).XLabel.String = 'stimulus types'; [hd.asymindex.explanation.FontSize, ... hd.asymindex.pvalue.FontSize] = deal(13); for panel = 1:1 hd.asymindex.ax(1,panel).XLabel.FontSize = 13; hd.asymindex.ax(1,panel).FontSize = 13; hd.asymindex.ax(1,panel).LabelFontSizeMultiplier = 1; hd.asymindex.ax(1,panel).YGrid = 'on'; hd.asymindex.ax(1,panel).TickLength = [0 0]; end end function hd = displayPattern(barwidth) spatialfreq = 180/barwidth; grid.colour = .3*[1 1 1]; grid.width = 360; grid.height = 180; grid.x0 = -180; grid.y0 = -90; hd.pattern.fg = figure(); hd.pattern.ax = axes(); hold on dy = grid.height; numbar = ceil(grid.width*spatialfreq); dx = grid.width / spatialfreq / 2; for barindex = 1:numbar x0 = grid.x0 + (barindex-1) * dx * 2; xmax = grid.x0 + grid.width; y0 = grid.y0; darkbar(barindex) = patch([x0*[1 1],min(x0+dx,xmax)*[1 1]], ... [y0,y0+dy,y0+dy,y0], grid.colour); end [darkbar(:).EdgeColor] = deal('none'); outline = plot(grid.x0*[1 1 1 1 1] + grid.width*[0 0 1 1 0], ... grid.y0*[1 1 1 1 1] + grid.height*[0 1 1 0 0], ... 'Color', grid.colour); hold off hd.pattern.ax.PlotBoxAspectRatio = [grid.width, grid.height, grid.height]; % % [hd.pattern.ax.XAxis.Visible, hd.pattern.ax.YAxis.Visible] = deal('off'); hd.pattern.ax.XLim = [-180 180]; hd.pattern.ax.YLim = [-90 90]; hd.pattern.ax.XTick = [-180 180]; hd.pattern.ax.YTick = [-90 0 90]; hd.pattern.ax.XLabel.String = 'azimuth'; hd.pattern.ax.YLabel.String = 'elevation'; end function hd = displayCosine(maxvelocity) time = 0:.2:25; velocity = sin(time*2*pi/10) * maxvelocity; hd.cosine.fg = figure(); hd.cosine.ax = axes(); hold on hd.cosine.ob(1) = plot(time([1 end]),[0 0]); hd.cosine.ob(2) = plot(time,velocity); hold off hd.cosine.ob(1).LineStyle = ':'; hd.cosine.ob(1).Color = .0*[1 1 1]; hd.cosine.ob(2).Color = .2*[1 1 1]; hd.cosine.ax.Box = 'off'; hd.cosine.ax.XGrid = 'on'; % hd.cosine.ax.XMinorGrid = 'on'; % hd.cosine.ax.YGrid = 'on'; % hd.cosine.ax.XAxisLocation = 'origin'; hd.cosine.ax.XTick = [0 25]; % hd.cosine.ax.XTick = 25; % hd.cosine.ax.XTickLabel = hd.cosine.ax.XTick; hd.cosine.ax.YTick = [-1 1] * maxvelocity; end function [handle,data,datafit] = gainAnalysisDisk(varargin) %#ok % GAINANALYSIS and all subfunctions were written by Florian Alexander % Dehmelt in 2018. This is version 20180911, a.k.a. gainAnalysis20180911.m. % Apologies for the lack of polish (and comments) at this point. % % All input arguments are optional: % % Resolution - scalar, indicating your choice of resolution for the % rendering of the spherical surface as part % of the % figures. Actual data fitting is completely unaffected. % CamAzimuth - scalar, horizontal angle at which sphere is shown % CamElevation - scalar, vertical angle at which sphere is shown % Normalisation - options: nothing, or 'mean gain per hemisphere', % or 'max gain per hemisphere', % or 'ratio of mean gains per hemisphere' % % Output arguments are: % % handle - structure containing handles to all graphics objects % (figures, axes, plots, decorations) % data - structure containing pre-processed data, with those fields: % .left - numerical array of means of OKR gain for the left eye % .right - numerical array of means of OKR gain for the right eye % .both - numerical array of mean of OKR gain means across both eyes % datafit - structure containing the fits to data, with following fields: % .centre - 2x2 numerical array of azimuths (1st row) and elevation % angles (2nd row) of the two individial von Mises-Fisher % functions fit to data (1st column = 1st function, etc.). % These are NOT the peaks of the actual fit (= sum of both % functions)! I still have to implement that feature. % .value - 36x1 numerical array containing the value of the actual fit % (= sum of both functions) at the locations sampled by data % .param - 1x9 numerical array containing all best-fit parameters % (offset, amplitude1, azimuth1, elevation1, concentration1, % amplitude2, azimuth2, elevation2, concentration2) % % The recommended default is calling "hd = gainAnalysis20180911(200);". % % To specify camera position before analysis, call the function as follows: % hd = gainAnalysis20180911(200,'CamAzimuth',-30,'CamElevation',15); % % To edit camera position post-hoc, adjust and execute the following lines: % az = 0; el = 0; % substitute your desired camera coordinates here % [hd.gain(end).ax(1:6).CameraPosition] = ... % deal([cosd(az)*cosd(el),sind(az)*cosd(el),sind(el)]/19.0526); % figure(21) % MANAGE VARIABLE INPUT p = inputParser; valid.resolution = @(x) isnumeric(x) && numel(x)==1 && x>0; valid.camazimuth = @(x) isnumeric(x) && numel(x)==1; valid.camelevation = @(x) isnumeric(x) && numel(x)==1; valid.normalisation = @(x) ischar(x); % valid.rotation = @(x) ischar(x); default.resolution = 200; default.camazimuth = 0; default.camelevation = 45; default.normalisation = ''; % default.rotation = 'none'; addOptional(p, 'Resolution', default.resolution, valid.resolution); addOptional(p, 'CamAzimuth', default.camazimuth, valid.camazimuth); addOptional(p, 'CamElevation', default.camelevation, valid.camelevation); addOptional(p, 'Normalisation', default.normalisation, valid.normalisation); % addOptional(p, 'Rotation', default.rotation, valid.rotation); parse(p,varargin{:}); resolution = p.Results.Resolution; cam.azimuth = p.Results.CamAzimuth; cam.elevation = p.Results.CamElevation; normalisation = p.Results.Normalisation; % rotation = p.Results.Rotation; % GENERAL SETTINGS warning('off','MATLAB:interp1:UsePCHIP') % The line above deactivates 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) % Identify which stimulus entries correspond to the disk-mask stimuli. xlsrange = 2:37; % xlsrange = 9:46; % % xlsrange = 9:44; % !!! IMPORTANT !!! These are the Excel row numbers, not stimulus types. % READ DATA % Query user to manually select the main data regularset. Other files have % standard names and are expected to be in the same folder. [mfilefolder,~,~] = fileparts(mfilename('fullpath')); % % [regularset.result,regularset.folder] = uigetfile([mfilefolder,'\*.mat']); % regularset.file = 'result_20180110_125537.mat'; % regularset.folder = ['..\..\data\Fig3 default disk data\Data D Series (neue Experimente)\']; % rotatedset.file = 'result_20181001_114504.mat'; % rotatedset.folder = ['..\..\data\Fig3 default disk data\Control Data Rotation (arena flip)\']; % upside down embedding, no rotation of arena % (both identical to bypass correction!) regularset.file = 'result_20180913_151415.mat'; regularset.folder = ['..\..\data\Fig4 rotated arena data\Control Data Upside Down (Kontrolle ventral)\']; rotatedset.file = regularset.file; rotatedset.folder = regularset.folder; regularset.structure = load([regularset.folder,regularset.file],'result'); rotatedset.structure = load([rotatedset.folder,rotatedset.file],'result'); regularset.result = regularset.structure.result; rotatedset.result = rotatedset.structure.result; % Discard unwanted trials, based on .xlsx documenting manual assessment. regularset.result = discardTrial(regularset.result,regularset.folder); rotatedset.result = discardTrial(rotatedset.result,rotatedset.folder); % Pool gain data in various ways (this is a significant subfunction!). regularset.gain = poolGainData(regularset.result); rotatedset.gain = poolGainData(rotatedset.result); % READ STIMULUS INFORMATION regularset.stimcentre = readStimulus(regularset.folder,xlsrange); rotatedset.stimcentre = readStimulus(rotatedset.folder,xlsrange); % reordering = [1;2;3;5;4;7;6]; % hemisphere stimuli, U vs. D and L vs. R flipped; now add disks: % for regularindex = 1:size(regularset.stimcentre,1) % % Find rotated-arena stimulus appearing at same position relative to fish: % rotatedindex = find(sum(regularset.stimcentre(regularindex,:)==rotatedset.stimcentre,2)==2); % reordering = [reordering; rotatedindex]; %#ok % end % % Correction above is unnecessary if stimulus coordinates were already in % fish-centred coordinates. In this case, do not reorder data: % reordering = 1:45; reordering = 1:43; % CHOOSE SUBSET OF DATA TO ANALYSE, AND PROCEED TO ANALYSE IT variantlist = {'regular','rotated','corrected'}; for variant = 1:numel(variantlist) switch variantlist{variant} case 'regular' data(variant).left = regularset.gain.mean.FR{1}; data(variant).right = regularset.gain.mean.FR{2}; case 'rotated' data(variant).left = rotatedset.gain.mean.FR{1}(reordering); data(variant).right = rotatedset.gain.mean.FR{1}(reordering); case 'corrected' % data(variant).left = (regularset.gain.mean.FR{1} + rotatedset.gain.mean.FR{1}) / 2; % data(variant).right = (regularset.gain.mean.FR{2} + rotatedset.gain.mean.FR{1}) / 2; data(variant).left = (regularset.gain.mean.FR{1} + rotatedset.gain.mean.FR{1}(reordering)) / 2; data(variant).right = (regularset.gain.mean.FR{2} + rotatedset.gain.mean.FR{1}(reordering)) / 2; for selection = xlsrange-1 gLreg = regularset.gain.mean.FR{1}(selection); gRreg = regularset.gain.mean.FR{2}(selection); gLrot = rotatedset.gain.mean.FR{1}(reordering(selection)); gRrot = rotatedset.gain.mean.FR{2}(reordering(selection)); data(variant).asymmetry.A1(selection) = ((gLreg-gRreg)/(gLreg+gRreg) - (gLrot-gRrot)/(gLrot+gRrot)) / 2; data(variant).asymmetry.A2(selection) = ((gLreg-gRreg)/(gLreg+gRreg) + (gLrot-gRrot)/(gLrot+gRrot)) / 2; end end data(variant).both = [data(variant).left(xlsrange(1:end/2)-1);data(variant).right(xlsrange(end/2+1:end)-1)]; % data(variant).both = mean([data(variant).left,data(variant).right],2); % important change! no more averaging % Create figure(s). hd.gain(variant).fg = figure(); hd.gain(variant).fg.Name = 'OKR gain'; hd.gain(variant).fg.Color = [1 1 1]; hd.gain(variant).fg.Position = [50 50 1000 650]; hd.gain(variant).fg.Colormap = flipud(cbrewer('div','RdBu',100)); % hd.gain(variant).fg.Colormap = cbrewer('div','PRGn',100); choicelist = {'left','both','right'}; for choice = 1:numel(choicelist) switch choicelist{choice} case 'right' response = data(variant).right(xlsrange-1); % -1 is offset by Excel title row case 'left' response = data(variant).left(xlsrange-1); case 'both' response = data(variant).both; otherwise error(['Input argument must be one of the following strings: ' ... '''left'', ''right'', or ''both''.']) end [x,y,z] = sphere(resolution); [datafit(choice).matrix,datafit(choice).param] = ... vonMisesFisherFit(regularset.stimcentre,response,x,y,z); % set 3D fish parameters scale = .1; yaw = 180; pitch = 0; roll = 0; % set decorative arc parameters; compute their cartesian coordinates equator.base = -1:.01:1; equator.x = [equator.base, fliplr(equator.base)]; equator.y = [sqrt(1-equator.base.^2), -sqrt(1-equator.base.^2)]; equator.z = 0*ones(1,numel(equator.x)); meridian.x = equator.x; meridian.y = equator.z; meridian.z = equator.y; aura.x = equator.z; aura.y = equator.y; aura.z = equator.x; [sample.x,sample.y,sample.z] = ... geo2car(regularset.stimcentre(:,1), regularset.stimcentre(:,2), ... ones(size(regularset.stimcentre,1),1),'CCW'); % plot all of the above onto the correct panel of the existing figure panel = choice; figure(hd.gain(variant).fg.Number) hd.gain(variant).ax(panel) = axes; hd.gain(variant).ax(panel).Position = [.13 .45 .35 .45] + (panel-1)*[.25 0 0 0]; hold on hd.gain(variant).ob(panel,1) = surf(x,y,z,datafit(panel).matrix); shading interp colour.eye = 0*[1 1 1]; colour.body = 1*[1 1 1]; hd.gain(variant).ob(panel,2:4) = plot3dFish(scale,yaw,pitch,roll,colour); hd.gain(variant).ob(panel,5) = plot3(equator.x,equator.y,equator.z,'k','LineWidth',2); hd.gain(variant).ob(panel,6) = plot3(meridian.x,meridian.y,meridian.z,'--k'); hd.gain(variant).ob(panel,7) = plot3(aura.x,aura.y,aura.z,':k'); hd.gain(variant).ob(panel,8) = scatter3(1.03*sample.x,1.03*sample.y,1.03*sample.z, ... 180*ones(size(sample.x)),response,'Marker','o'); hold off end % Apply full formatting to the figure(s). Subfunction included below. hd = fishbowlFigure(hd,cam,choicelist); % % Compute yoking/lateralisation index per individual stimulus location. % convention = 'CCW'; % switch convention % case 'CCW' % stim.isleft = stim.azimuth > 0; % stim.isright = stim.azimuth < 0; % case 'CW' % stim.isleft = stim.azimuth < 0; % stim.isright = stim.azimuth > 0; % otherwise % error('Sign convention must either be ''CCW'' or ''CW''.') % end end % Collect output arguments, unless already done (e.g., "datafit"). handle = hd; end function stimcentre = readStimulus(datafolder,xlsrange) stimulus = 'Stimulus.xlsx'; [~,tableconvention,~] = xlsread([datafolder,stimulus],'C1:C1'); azimuth = xlsread([datafolder,stimulus],['C',num2str(xlsrange(1)),':C',num2str(xlsrange(end))]); elevation = xlsread([datafolder,stimulus],['D',num2str(xlsrange(1)),':D',num2str(xlsrange(end))]); % Stimulus tables may use CW convention; default code uses CCW. switch tableconvention{:} case 'azimuth (CCW)' case 'azimuth (CW)' azimuth = -azimuth; otherwise error(['Field C1 of Stimulus.xlsx should contain either ' ... 'string ''azimuth (CCW)'' or string ''azimuth (CW)''.']) end stimcentre = [azimuth,elevation]; end function gain = poolGainData(result) %POOLGAINDATA converts a raw result structure into plottable data by %pooling OKR gain data in various ways, and storing those results in a new %structure called "gain". % Pool data in various ways (e.g., .FR = pooled across all [F]ish and % all [R]epetitions, but with a separate array for each [E]ye). for stimtype = 1:size(result,3) gainpool(stimtype).FER = [result(:,:,stimtype,:).gain]; for eye = 1:size(result,2) gainpool(stimtype).FR{eye} = [result(:,eye,stimtype,:).gain]; end for fish = 1:size(result,1) gainpool(stimtype).ER{fish} = [result(fish,:,stimtype,:).gain]; end for repetition = 1:size(result,4) gainpool(stimtype).FE{repetition} = [result(:,:,stimtype,repetition).gain]; end for fish = 1:size(result,1) for eye = 1:size(result,2) gainpool(stimtype).R{fish,eye} = [result(fish,eye,stimtype,:).gain]; end end end % Make sure that in following section, if a stimulus type is not present, % the array contain a NaN at that position, instead of a zero. nantemplate = NaN(numel(gainpool),1); gain.mean .FER = nantemplate; gain.median.FER = nantemplate; gain.std .FER = nantemplate; gain.sem .FER = nantemplate; for eye = 1:size(result,2) gain.mean .FR{eye} = nantemplate; gain.median.FR{eye} = nantemplate; gain.std .FR{eye} = nantemplate; gain.sem .FR{eye} = nantemplate; end for fish = 1:size(result,1) gain.mean .ER{fish} = nantemplate; gain.median.ER{fish} = nantemplate; gain.std .ER{fish} = nantemplate; gain.sem .ER{fish} = nantemplate; end for fish = 1:size(result,1) for eye = 1:size(result,2) gain.mean .R{fish,eye} = nantemplate; gain.median.R{fish,eye} = nantemplate; gain.std .R{fish,eye} = nantemplate; gain.sem .R{fish,eye} = nantemplate; end end % mean, median, stdev, sem across all pooled data, one per stimulus type for stimtype = 1:numel(gainpool) selection = gainpool(stimtype).FER; gain.mean .FER(stimtype) = mean(selection); gain.median.FER(stimtype) = median(selection); gain.std .FER(stimtype) = std(selection); gain.sem .FER(stimtype) = std(selection)/sqrt(numel(selection)); for eye = 1:size(result,2) selection = gainpool(stimtype).FR{eye}; gain.mean .FR{eye}(stimtype) = mean(selection); gain.median.FR{eye}(stimtype) = median(selection); gain.std .FR{eye}(stimtype) = std(selection); gain.sem .FR{eye}(stimtype) = std(selection)/sqrt(numel(selection)); end for fish = 1:size(result,1) selection = gainpool(stimtype).ER{fish}; gain.mean .ER{fish}(stimtype) = mean(selection); gain.median.ER{fish}(stimtype) = median(selection); gain.std .ER{fish}(stimtype) = std(selection); gain.sem .ER{fish}(stimtype) = std(selection)/sqrt(numel(selection)); end for fish = 1:size(result,1) for eye = 1:size(result,2) selection = gainpool(stimtype).R{fish,eye}; gain.mean .R{fish,eye}(stimtype) = mean(selection); gain.median.R{fish,eye}(stimtype) = median(selection); gain.std .R{fish,eye}(stimtype) = std(selection); gain.sem .R{fish,eye}(stimtype) = std(selection)/sqrt(numel(selection)); end end end end function [fitresult,fitparam] = vonMisesFisherFit(stimcentre,response,x,y,z) predictor = stimcentre; % Initialise fit parameters. offset = 0; amplitude(1) = 1; amplitude(2) = 1; meanazimuth(1) = 70; meanelevation(1) = 20; % azimuth, elevation meanazimuth(2) = -70; meanelevation(2) = 20; % azimuth, elevation concentration(1) = 2.5; concentration(2) = 2.5; initialparam = [offset, ... amplitude(1), ... meanazimuth(1), ... meanelevation(1), ... concentration(1), ... amplitude(2), ... meanazimuth(2), ... meanelevation(2), ... concentration(2)]; % Fit. option.MaxIter = 1e4; datafit.param = nlinfit(predictor,response,@vonmisesfisher9param,initialparam,option); datafit.value = vonmisesfisher9param(datafit.param,predictor); % Identify centres. azimuth(1) = datafit.param(3); elevation(1) = datafit.param(4); azimuth(2) = datafit.param(7); elevation(2) = datafit.param(8); for k = 1:2 % azimuth(k) = datafit.param(3); % This was incorrect - same value read twice. % elevation(k) = datafit.param(4); uniqueazimuth(k) = mod(azimuth(k)+180,360) - 180; uniqueelevation(k) = mod(elevation(k)+90,180) - 90; datafit.centre(:,k) = [uniqueazimuth(k), uniqueelevation(k)]; end % compute geographic coordinates of sphere, and its colour data values [azimuth,elevation,radius] = car2geo(x,y,z,'CCW'); predictor = [azimuth(:),elevation(:),radius(:)]; fitresult = reshape(vonmisesfisher9param(datafit.param,predictor),size(radius)); fitparam = datafit.param; end function hd = fishbowlFigure(hd,cam,choicelist) %FISHBOWLFIGURE configures an existing figure of fishbowl plots, adds decorations % % HD = fishbowlFigure sets the parameters of an existing fishbowl figure, % including axis limits, labels, colours, axes object sizes, etc., and % adds additional decorations such as colorbars. cmin = min([min(min(hd.gain(end).ob(1,1).CData)), ... min(min(hd.gain(end).ob(2,1).CData)), ... min(min(hd.gain(end).ob(3,1).CData)), ... min(min(hd.gain(end).ob(1,8).CData)), ... min(min(hd.gain(end).ob(2,8).CData)), ... min(min(hd.gain(end).ob(3,8).CData))]); cmax = max([max(max(hd.gain(end).ob(1,1).CData)), ... max(max(hd.gain(end).ob(2,1).CData)), ... max(max(hd.gain(end).ob(3,1).CData)), ... max(max(hd.gain(end).ob(1,8).CData)), ... max(max(hd.gain(end).ob(2,8).CData)), ... max(max(hd.gain(end).ob(3,8).CData))]); clim = max(abs([cmin cmax]))*[-1 1]; for panel = 1:3 switch choicelist{panel} case 'left' string.title = 'left eye gain g_L'; case 'right' string.title = 'right eye gain g_R'; case 'both' string.title = 'merge stimulated'; % string.title = 'mean = (g_L+ g_R)/2'; end hd.gain(end).ax(panel).Title.String = string.title; % title(string.title,'Interpreter','LaTeX') hd.gain(end).ob(panel,8).MarkerFaceColor = hd.gain(end).ob(panel,8).MarkerEdgeColor; hd.gain(end).ob(panel,8).MarkerEdgeColor = [0 0 0]; hd.gain(end).ax(panel).CLim = [0 cmax]; % hd.gain(end).ax(panel).CLim = [0 max(max(datafit(choice).matrix))]; % hd.gain(end).ob(panel,1).FaceAlpha = .85; hd.gain(end).ob(panel,1).FaceAlpha = .75; hd.gain(end).ob(panel,1).EdgeColor = 'none'; % hd.gain(end).ob(panel,8).Shading = 'flat'; shading interp axis square xlabel('x') ylabel('y') zlabel('z') hd.gain(end).ax(panel).XLim = [-1.1 1.1]; hd.gain(end).ax(panel).YLim = [-1.1 1.1]; hd.gain(end).ax(panel).ZLim = [-1.1 1.1]; % hd.gain(end).ax(panel).CameraPosition = [12.2934, -9.0967,8.1315]; cam.radius = 19.0526; % MATLAB automatically reverts to this values while rotating, anyway. Might as well choose it as default. [cam.x,cam.y,cam.z] = geo2car(cam.azimuth,cam.elevation,cam.radius,'CCW'); hd.gain(end).ax(panel).CameraPosition = [cam.x cam.y cam.z]; % FD20190507: this hides 3rd fish!!! hd.gain(end).ax(panel).CameraTarget = [0 0 0]; hd.gain(end).ax(panel).CameraUpVector = [0 0 1]; % hd.gain(end).ax(panel).CameraViewAngle = 10.134; hd.gain(end).ax(panel).CameraViewAngle = 9; end hd.gain(end).ax(4) = axes; hd.gain(end).ax(4).Position = hd.gain(end).ax(1).Position - [.13 -.1 .2 .2]; hd.gain(end).ax(4).CLim = hd.gain(end).ax(1).CLim; hd.gain(end).ob(4,1) = colorbar; hd.gain(end).ob(4,1).Label.String = 'OKR gain'; hd.gain(end).ob(4,1).Label.Units = 'normalized'; hd.gain(end).ob(4,1).Label.Position = [-.8 1 1] .* hd.gain(end).ob(4,1).Label.Position; [hd.gain(end).ob(4,1).Label.FontSize, hd.gain(end).ob(4,1).FontSize, ... hd.gain(end).ax.FontSize] = deal(14); [hd.gain(end).ax(1).Title.FontSize, ... hd.gain(end).ax(2).Title.FontSize, ... hd.gain(end).ax(3).Title.FontSize] = deal(14); [hd.gain(end).ax(1).Title.FontWeight, ... hd.gain(end).ax(2).Title.FontWeight, ... hd.gain(end).ax(3).Title.FontWeight] = deal('normal'); [hd.gain(end).ax.Visible] = deal('off'); % [hd.gain(end).ax.Title.Visible] = deal('on'); set(([hd.gain(end).ax.Title]),'Visible','on') set(([hd.gain(end).ax.Title]),'Position',[1 1 .7] .* hd.gain(end).ax(1).Title.Position) end function [retained,original] = discardTrial(original,folderstring) % DISCARDTRIAL removes data obtained from trials consider "poorly fit" etc. % % This function takes a "result" structure as its only input argument, and % requires an XLSX file containing the indices of trials to be discarded. % It then removes those bad trials, and returns a reduced structure as its % primary output argument, "retained". For easy access, it also returns its % original input. % % This function can safely be executed multiple times, as "discarded" % trials are in fact just set to empty, so their array indices are not % reassigned to other trials. % % Written by Florian Alexander Dehmelt, Tuebingen University, in 2018. % Identify and load an Excel file containing the indices of all unwanted % trials to be discarded. This file must contain four columns of numbers % (from left to right: fish no., eye no., stimulus type, and repetition. % It may contain as many text comments as desired; these will be ignored. folder.badtrial = folderstring; % folder.badtrial = 'C:\Users\fdehmelt\Dropbox\ThinkPad Dropbox\MATLAB\201808 SphereAnalysis\Sample Data\D series\'; regularset.badtrial = 'DiscardTrial.xlsx'; badtrial = xlsread([folder.badtrial,regularset.badtrial]); % The "retained" structure is the same as the "original" one... retained = original; % ...except that unwanted trials are overwritten by empty values. % They thus appear the same way as non-existent trials already did. for trial = 1:size(badtrial,1) % Find the position of the bad trial within the structure. fish = badtrial(trial,1); eye = badtrial(trial,2); stimtype = badtrial(trial,3) + 1; % The +1 converts index >=0 to >=1. repetition = badtrial(trial,4); % Overwrite all field entries with "[]". fieldname = fieldnames(retained); for field = 1:numel(fieldname) retained(fish,eye,stimtype,repetition).(fieldname{field}) = []; end end end function [altered,original] = zeroTrial(original,folderstring) % ZEROTRIAL sets gains to zero based on visual inspection of trials. % % This function takes a "result" structure as its only input argument, and % requires an XLSX file containing the indices of trials to be set to zero. % It then adjusts those gains, and returns an altered structure as its % primary output argument, "retained". For easy access, it also returns its % original input. % % This function can safely be executed multiple times, as "altered" % trials are in fact just set to zero, so their array indices are not % reassigned to other trials. % % Written by Florian Alexander Dehmelt, Tuebingen University, in 2018. % Identify and load an Excel file containing the indices of all unwanted % trials to be discarded. This file must contain four columns of numbers % (from left to right: fish no., eye no., stimulus type, and repetition. % It may contain as many text comments as desired; these will be ignored. folder.badtrial = folderstring; % folder.badtrial = 'C:\Users\fdehmelt\Dropbox\ThinkPad Dropbox\MATLAB\201808 SphereAnalysis\Sample Data\D series\'; regularset.badtrial = 'ZeroTrial.xlsx'; badtrial = xlsread([folder.badtrial,regularset.badtrial]); % The "altered" structure is the same as the "original" one... altered = original; % ...except that gains of visually identified trials are set to zero. % They thus appear the same way as non-existent trials already did. for trial = 1:size(badtrial,1) % Find the position of the bad trial within the structure. fish = badtrial(trial,1); eye = badtrial(trial,2); stimtype = badtrial(trial,3) + 1; % The +1 converts index >=0 to >=1. repetition = badtrial(trial,4); % Overwrite the gain field entries with "0", retain all other values. altered(fish,eye,stimtype,repetition).gain = 0; end end function combinedvalue = vonmisesfisher9param(param,predictor) %VONMISESFISHER9PARAM Evaluate two overlapping von Mises-Fisher distributions. % % VALUE = VONMISESFISHER9PARAM computes the value of two(!) overlapping, % 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 = '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 = 2; value = NaN(size(predictor,1),numberofdistributions); for k = 1:numberofdistributions offset = param(1); amplitude = param(4*k-2); meanazimuth = param(4*k-1); meanelevation = param(4*k); concentration = param(4*k+1); [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; % 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 function [x,y,z] = geo2car(azimuth,elevation,radius,convention) switch convention case 'CCW' x = radius .* cosd(elevation).*cosd(azimuth); y = radius .* cosd(elevation).*sind(azimuth); z = radius .* sind(elevation); case 'CW' x = radius .* cosd(elevation).*cosd(azimuth); y = radius .* cosd(elevation).*-sind(azimuth); z = radius .* sind(elevation); otherwise error('Convention not supported.') end end function [azimuth,elevation,radius] = car2geo(x,y,z,convention) azimuth = atand(y./x) + 180*((x<0).*(y>=0) - (x<0).*(y<0)); switch convention case 'CCW' case 'CW' azimuth = -azimuth; otherwise error('Convention not supported.') end radius = sqrt(x.^2+y.^2+z.^2); elevation = asind(z./radius); end function [lefteyehandle,righteyehandle,bodyhandle] = ... plot3dFish(scale,yaw,pitch,roll,colour) % This corresponds to version 20170505a, a.k.a. plot3dFish20170505a.m. set(gcf,'Color',[1 1 1]) resolution = 10; % colour.eye = .1*[1 1 1]; % colour.body = .7*[1 1 1]; eye.x = scale * (cos(-pi:pi/resolution:3*pi-pi/resolution)'); eye.y = scale * (0.7*sin(0:pi/resolution:4*pi-pi/resolution)'); eye.z = scale * ([.7*ones(2*resolution,1); ... 1*ones(2*resolution,1)]); N = 2*resolution; lefteye.vertices = [eye.x, eye.y, eye.z]; righteye.vertices = [eye.x, -eye.y, eye.z]; lefteye.faces = [[1:2*resolution 1]; ... [1+2*resolution:4*resolution 1+2*resolution]]; for k = 1:N lefteye.faces = [lefteye.faces; ... [mod([k-1 k],N)+1, ... mod([k k-1],N)+1+N, ... mod(k-1,N)+1, ... NaN(1,N-4)]]; end lefteye.facevertexcdata = repmat(colour.eye,[size(lefteye.faces,1) 1]); righteye.faces = lefteye.faces; righteye.facevertexcdata = lefteye.facevertexcdata; lefteyehandle = patch(lefteye); righteyehandle = patch(righteye); set([lefteyehandle,righteyehandle],'FaceColor','flat','EdgeColor','none') mainbody.x = scale * (-.1 + [1.0*cos(-pi:pi/resolution:pi-pi/resolution), ... 0.4*cos(-pi:pi/resolution:pi-pi/resolution)]'); mainbody.y = scale * ([0.7*sin(0:pi/resolution:2*pi-pi/resolution), ... 0.2*sin(0:pi/resolution:2*pi-pi/resolution)]'); mainbody.z = scale * ([1.1*ones(2*resolution,1); ... -7*ones(2*resolution,1)]); N = 2*resolution; body.vertices = [mainbody.x, mainbody.y, mainbody.z]; body.faces = [[1:2*resolution 1]; ... [1+2*resolution:4*resolution 1+2*resolution]]; for k = 1:N body.faces = [body.faces; ... [mod([k-1 k],N)+1, ... mod([k k-1],N)+1+N, ... mod(k-1,N)+1, ... NaN(1,N-4)]]; end body.facevertexcdata = repmat(colour.body,[size(body.faces,1) 1]); bodyhandle = patch(body); set(bodyhandle,'FaceColor','flat','EdgeColor','none') rotate(lefteyehandle,[1 0 0],80) rotate(lefteyehandle,[0 0 1],-10) rotate(righteyehandle,[1 0 0],-80) rotate(righteyehandle,[0 0 1],10) rotate(bodyhandle,[0 1 0],-90) rotate(lefteyehandle,[0 0 1],yaw) rotate(righteyehandle,[0 0 1],yaw) rotate(bodyhandle,[0 0 1],yaw) rotate(lefteyehandle,[0 1 0],pitch) rotate(righteyehandle,[0 1 0],pitch) rotate(bodyhandle,[0 1 0],pitch) rotate(lefteyehandle,[1 0 0],roll) rotate(righteyehandle,[1 0 0],roll) rotate(bodyhandle,[1 0 0],roll) end function [hd,stim] = displayDiskStimuli(choice) warning('off','MATLAB:interp1:UsePCHIP') % The line above deactivates 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) % close all % clf % regularset.folder = '..\..\data\Fig3 default disk data\Data D Series (neue Experimente)\'; % xlsrange = 9:46; regularset.folder = ['..\..\data\Fig4 rotated arena data\Control Data Upside Down (Kontrolle ventral)\']; xlsrange = 2:37; % regularset.folder = '..\..\data\Fig5 size data\Data A Series (A-Series)\'; % xlsrange = 2:43; % xlsrange = 2:8; % regularset.folder = '..\..\data\Fig5 frequency data\Data V Series (V-Serie)\'; % xlsrange = 2:43; regularset.stimulus = 'Stimulus.xlsx'; [~,stim.tableconvention,~] = xlsread([regularset.folder,regularset.stimulus],'C1:C1'); stim.azimuth = xlsread([regularset.folder,regularset.stimulus],['C',num2str(xlsrange(1)),':C',num2str(xlsrange(end))]); stim.elevation = xlsread([regularset.folder,regularset.stimulus],['D',num2str(xlsrange(1)),':D',num2str(xlsrange(end))]); % stim.radius = xlsread([regularset.folder,regularset.stimulus],['E',num2str(xlsrange(1)),':E',num2str(xlsrange(end))]); % stim.azimuth = stim.azimuth(1:2:end); % stim.elevation = stim.elevation(1:2:end); stim.radius = 64*ones(size(stim.elevation)); % % [surface.]backgroundSphere % % compute background sphere with shading % [surface.x,surface.y,surface.z] = sphere(100); % reference.geo = [30;20;1]; % desired azimuth, elevation, radius of maximum colour value % reference.azim = reference.geo(1); % reference.elev = reference.geo(2); % reference.rads = ones(size(reference.geo(1))); % [reference.x,reference.y,reference.z] = geo2car(reference.azim,reference.elev,reference.rads,'CW'); % scalarproduct = [surface.x(:),surface.y(:),surface.z(:)] * [reference.x;reference.y;reference.z]; % surface.c = reshape(scalarproduct,size(surface.x)); % surface.c = 10.^(surface.c+1); % % surface.c = 2.^(surface.c+1); % compute disk boundaries diskcentre.geo = [stim.azimuth, stim.elevation]; diskcentre.angle = stim.radius; % diskcentre.geo = [ 45 20; ... % 30 5; ... % 60 -45]; % diskcentre.angle = 40 * [ones(size(diskcentre.geo,1),1)]; diskcentre.azim = diskcentre.geo(:,1); diskcentre.elev = diskcentre.geo(:,2); diskcentre.rads = ones(size(diskcentre.geo(:,1))); resolution = .2; numdisk = size(diskcentre.geo,1); for disk = 1:numdisk [boundary.azim{disk}, ... boundary.elev{disk}] = diskBoundary(diskcentre.azim(disk), ... diskcentre.elev(disk), diskcentre.angle(disk), resolution); [boundary.x{disk}, ... boundary.y{disk}, ... boundary.z{disk}] = geo2car(boundary.azim{disk}, boundary.elev{disk}, ... ones(size(boundary.azim{disk})),'CW'); end [diskcentre.x,diskcentre.y,diskcentre.z] = ... geo2car(diskcentre.azim,diskcentre.elev,diskcentre.rads,'CW'); [textcentre.x,textcentre.y,textcentre.z] = ... geo2car(diskcentre.azim,.91*diskcentre.elev,1.10*diskcentre.rads,'CW'); % create figure hd.default.fg = figure(); hd.default.fg.Name = 'This is a sphere.'; hd.default.fg.Position = [50 50 500 800]; hd.default.fg.Color = [1 1 1]; hd.default.ax(1) = axes; hd.default.ax(1).Position = [0 0 1 1]; hold on % plot background sphere with shading hd.default.globe = displaySphere; % hd.default.globe = surf(surface.x,surface.y,surface.z,surface.c); % shading flat % axis square % greyscale = (.8:.01:.99)' * [1 1 1]; % colormap(greyscale) fish.scale = .05; fish.yaw = 180; fish.pitch = 0; fish.roll = 0; colour.eye = 0*[1 1 1]; colour.body = .5*[1 1 1]; plot3dFish(fish.scale,fish.yaw,fish.pitch,fish.roll,colour); hd.default.globe.FaceAlpha = .5; cropdisk = 8; switch choice case 'left' disklist = 1 : numdisk/2; case 'right' disklist = numdisk/2+1 : numdisk; case 'single' disklist = cropdisk; end % plot disk boundaries % diskcolour = cbrewer('seq','Greens',numdisk*2); cbrewergreen = cbrewer('seq','Greens',numdisk); diskcolour = repmat(cbrewergreen(round(numdisk/1.4),:),[numdisk,1]); % diskcolour = flipud(cbrewer('seq','Greys',numdisk*2)); % diskcolour = diskcolour(randperm(numdisk),:); for disk = disklist if sum(disk==cropdisk) plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',2,'Color',0*[1 1 1]) disktext(disk) = text(textcentre.x(disk),textcentre.y(disk),textcentre.z(disk), ... ['',num2str(disk)], ... 'FontSize', 13, 'HorizontalAlignment','center','VerticalAlignment','middle'); else plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',1,'Color',.6*[1 1 1]) if sum(disk==disklist) disktext(disk) = text(textcentre.x(disk),textcentre.y(disk),textcentre.z(disk), ... ['',num2str(disk)], ... 'FontSize', 13, 'HorizontalAlignment','center','VerticalAlignment','middle'); end % plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',2,'Color',1*[1 1 1]) % plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',1,'Color',diskcolour(disk,:)) % plot3(boundary.x{disk},boundary.y{disk},boundary.z{disk},'-','LineWidth',1,'Color',.2*[1 1 1]) end end hold off switch choice case 'single' disktext(cropdisk).Visible = 'off'; barwidth = 360/22/2; % gratingazimuth = -180:barwidth:(180-barwidth); gratingazimuth = -180:barwidth:-80; for bar = 1:2:numel(gratingazimuth)-1 % if bar < numel(gratingazimuth) azimlimit = gratingazimuth([bar bar+1]); % else % azimlimit = gratingazimuth([bar 1]); % end displaySphereBar(azimlimit(1),azimlimit(2),boundary.azim{cropdisk},boundary.elev{cropdisk}) end case {'left','right'} disktext(disklist(ceil(end/2))).Position = ... disktext(disklist(ceil(end/2))).Position - [0 0 .16]; end % hd.default.ax(1).CameraPosition = [2.3 17 2.2]; switch choice case {'left','single'} hd.default.ax(1).CameraPosition = [0 20 0]; case 'right' hd.default.ax(1).CameraPosition = [0 -20 0]; end hd.default.ax(1).CameraViewAngle = 6; [hd.default.ax(1).XAxis.Visible, ... hd.default.ax(1).YAxis.Visible, ... hd.default.ax(1).ZAxis.Visible] = deal('off'); end function globe = displaySphere [surface.x,surface.y,surface.z] = sphere(300); reference.geo = [30;20;1]; % desired azimuth, elevation, radius of maximum colour value reference.azim = reference.geo(1); reference.elev = reference.geo(2); reference.rads = ones(size(reference.geo(1))); [reference.x,reference.y,reference.z] = geo2car(reference.azim,reference.elev,reference.rads,'CW'); scalarproduct = [surface.x(:),surface.y(:),surface.z(:)] * [reference.x;reference.y;reference.z]; surface.c = reshape(scalarproduct,size(surface.x)); surface.c = 10.^(surface.c); globe = surf(surface.x,surface.y,surface.z,surface.c); shading interp axis square greyscale = (.8:.01:.99)' * [1 1 1]; % greyscale = (.9:.01:.99)' * [1 1 1]; colormap(greyscale) end function [azimuthOut, elevationOut] = diskBoundary(azimuthIn, elevationIn, angleIn, 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; % This equation is used to compute border.azimuth below. centre.azimuth = azimuthIn; centre.elevation = elevationIn; border.angle = angleIn / 2; % convert solid angle ("diameter") to "radius" border.resolution = resolutionIn; intended = -90 : border.resolution : 90; % limit = cosd([centre.elevation-border.angle centre.elevation+border.angle]); % intended = acosd(min(limit) : border.resolution : max(limit)); % "auxiliary" elevations are explicitly added to guarantee coverage near top/bottom of circle: auxiliary = centre.elevation + border.angle*[-1+1e-9 1-1e-9]; border.elevation = unique([intended,auxiliary]); 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)-168; azimuthOut = border.azimuth'; elevationOut = border.elevation'; end function hd = displaySphereBar(azimlimit1,azimlimit2,boundaryazim,boundaryelev) resolution = .01; if azimlimit2 < azimlimit1 error('sdasdsa') end crop.index = logical((azimlimit1 <= boundaryazim) .* (boundaryazim <= azimlimit2)); if sum(crop.index) > 0 crop.azim = boundaryazim(crop.index)'; crop.elev = boundaryelev(crop.index)'; crop.azim1 = crop.azim(diff(crop.azim) > 0); crop.azim2 = crop.azim(diff(crop.azim) < 0); crop.elev1 = crop.elev(diff(crop.azim) > 0); crop.elev2 = crop.elev(diff(crop.azim) < 0); sequence = @(a,b) a:(sign(b-a)*resolution):b; edge.elev1 = sequence(crop.elev1(1), crop.elev2(end)); edge.elev2 = sequence(crop.elev2(1), crop.elev1(end)); % edge.elev1 = sequence(crop.elev1(1), crop.elev2(end)); % edge.elev2 = sequence(crop.elev2(1), crop.elev1(end)); edge.azim1 = azimlimit1*ones(size(edge.elev1)); edge.azim2 = azimlimit2*ones(size(edge.elev2)); azim = [edge.azim1,fliplr(crop.azim2),edge.azim2,fliplr(crop.azim1)]; elev = [edge.elev1,fliplr(crop.elev2),edge.elev2,fliplr(crop.elev1)]; rads = 1.009*ones(size(elev)); [bar.x,bar.y,bar.z] = geo2car(azim,elev,rads,'CW'); hold on, p = patch(bar.x,bar.y,bar.z,'k'); hold off % p.FaceColor = [.2 .8 .4]; p.FaceColor = .2*[1 1 1]; p.EdgeColor = 'none'; % hold on, plot3(bar.x,bar.y,bar.z,'r'); hold off end end