/trunk/examples/toolboxes/fieldtrip/ft_sliceinterp.m
MATLAB | 517 lines | 364 code | 32 blank | 121 comment | 90 complexity | b8765d3acf0d4923a089281b8a231cab MD5 | raw file
Possible License(s): BSD-3-Clause, AGPL-1.0, LGPL-2.0, GPL-3.0, GPL-2.0, LGPL-2.1, LGPL-3.0, BSD-2-Clause
- function [outim]=ft_sliceinterp(cfg, ininterp)
- % FT_SLICEINTERP plots a 2D-montage of source reconstruction and anatomical MRI
- % after these have been interpolated onto the same grid.
- %
- % Use as
- % ft_sliceinterp(cfg, interp)
- % or
- % [rgbimage] = ft_sliceinterp(cfg, interp), rgbimage is the monatage image
- %
- % where interp is the output of sourceinterpolate and cfg is a structure
- % with any of the following fields:
- %
- % cfg.funparameter string with the functional parameter of interest (default = 'source')
- % cfg.maskparameter parameter used as opacity mask (default = 'none')
- % cfg.clipmin value or 'auto' (clipping of source data)
- % cfg.clipmax value or 'auto' (clipping of source data)
- % cfg.clipsym 'yes' or 'no' (default) symmetrical clipping
- % cfg.colormap colormap for source overlay (default is jet(128))
- % cfg.colmin source value mapped to the lowest color (default = 'auto')
- % cfg.colmax source value mapped to the highest color (default = 'auto')
- % cfg.maskclipmin value or 'auto' (clipping of mask data)
- % cfg.maskclipmax value or 'auto' (clipping of mask data)
- % cfg.maskclipsym 'yes' or 'no' (default) symmetrical clipping
- % cfg.maskmap opacitymap for source overlay
- % (default is linspace(0,1,128))
- % cfg.maskcolmin mask value mapped to the lowest opacity, i.e.
- % completely transparent (default ='auto')
- % cfg.maskcolmin mask value mapped to the highest opacity, i.e.
- % non-transparent (default = 'auto')
- % cfg.alpha value between 0 and 1 or 'adaptive' (default)
- % cfg.nslices integer value, default is 20
- % cfg.dim integer value, default is 3 (dimension to slice)
- % cfg.spacemin 'auto' (default) or integer (first slice position)
- % cfg.spacemax 'auto' (default) or integer (last slice position)
- % cfg.resample integer value, default is 1 (for resolution reduction)
- % cfg.rotate number of ccw 90 deg slice rotations (default = 0)
- % cfg.title optional title (default is '')
- % cfg.whitebg 'yes' or 'no' (default = 'yes')
- % cfg.flipdim flip data along the sliced dimension, 'yes' or 'no'
- % (default = 'no')
- % cfg.marker [Nx3] array defining N marker positions to display
- % cfg.markersize radius of markers (default = 5);
- % cfg.markercolor [1x3] marker color in RGB (default = [1 1 1], i.e. white)
- % cfg.interactive 'yes' or 'no' (default), interactive coordinates
- % and source values
- %
- % if cfg.alpha is set to 'adaptive' the opacity of the source overlay
- % linearly follows the source value: maxima are opaque and minima are
- % transparent.
- %
- % if cfg.spacemin and/or cfg.spacemax are set to 'auto' the sliced
- % space is automatically restricted to the evaluated source-space
- %
- % if cfg.colmin and/or cfg.colmax are set to 'auto' the colormap is mapped
- % to source values the following way: if source values are either all
- % positive or all negative the colormap is mapped to from
- % min(source) to max(source). If source values are negative and positive
- % the colormap is symmetrical mapped around 0 from -max(abs(source)) to
- % +max(abs(source)).
- %
- % If cfg.maskparameter specifies a parameter to be used as an opacity mask
- % cfg.alpha is not used. Instead the mask values are maped to an opacitymap
- % that can be specified using cfg.maskmap. The mapping onto that
- % opacitymap is controlled as for the functional data using the
- % corresponding clipping and min/max options.
- %
- % if cfg.whitebg is set to 'yes' the function estimates the head volume and
- % displays a white background outside the head, which can save a lot of black
- % printer toner.
- %
- % if cfg.interactive is set to 'yes' a button will be displayed for
- % interactive data evaluation and coordinate reading. After clicking the
- % button named 'coords' you can click on any position in the slice montage.
- % After clicking these coordinates and their source value are displayed in
- % a text box below the button. The coordinates correspond to indeces in the
- % input data array:
- %
- % f = interp.source(coord_1,coord_2,coord_3)
- %
- % The coordinates are not affected by any transformations used for displaying
- % the data such as cfg.dim, cfg.rotate,cfg.flipdim or cfg.resample.
- % Copyright (C) 2004, Markus Siegel, markus.siegel@fcdonders.kun.nl
- %
- % This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
- % for the documentation and details.
- %
- % FieldTrip is free software: you can redistribute it and/or modify
- % it under the terms of the GNU General Public License as published by
- % the Free Software Foundation, either version 3 of the License, or
- % (at your option) any later version.
- %
- % FieldTrip is distributed in the hope that it will be useful,
- % but WITHOUT ANY WARRANTY; without even the implied warranty of
- % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- % GNU General Public License for more details.
- %
- % You should have received a copy of the GNU General Public License
- % along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
- %
- % $Id: ft_sliceinterp.m 948 2010-04-21 18:02:21Z roboos $
- fieldtripdefs
- % check if the input data is valid for this function
- ininterp = checkdata(ininterp, 'datatype', 'volume', 'feedback', 'yes');
- if ~isfield(cfg, 'clipmin'); cfg.clipmin = 'auto'; end
- if ~isfield(cfg, 'clipmax'); cfg.clipmax = 'auto'; end
- if ~isfield(cfg, 'clipsym'); cfg.clipsym = 'no'; end
- if ~isfield(cfg, 'alpha'); cfg.alpha = 'adaptive'; end
- if ~isfield(cfg, 'nslices'); cfg.nslices = 20; end
- if ~isfield(cfg, 'dim'); cfg.dim = 3; end
- if ~isfield(cfg, 'colormap'); cfg.colormap = jet(128); end
- if ~isfield(cfg, 'spacemin'); cfg.spacemin = 'auto'; end
- if ~isfield(cfg, 'spacemax'); cfg.spacemax = 'auto'; end
- if ~isfield(cfg, 'colmin'); cfg.colmin = 'auto'; end
- if ~isfield(cfg, 'colmax'); cfg.colmax = 'auto'; end
- if ~isfield(cfg, 'resample'); cfg.resample = 1; end
- if ~isfield(cfg, 'rotate'); cfg.rotate = 0; end
- if ~isfield(cfg, 'title'); cfg.title = ''; end
- if ~isfield(cfg, 'whitebg'); cfg.whitebg = 'no'; end
- if ~isfield(cfg, 'flipdim'); cfg.flipdim = 'no'; end
- if ~isfield(cfg, 'marker'); cfg.marker = []; end
- if ~isfield(cfg, 'markersize'); cfg.markersize = 5; end
- if ~isfield(cfg, 'markercolor'); cfg.markercolor = [1,1,1]; end
- if ~isfield(cfg, 'interactive'); cfg.interactive = 'no'; end
- if ~isfield(cfg, 'maskclipmin'); cfg.maskclipmin = 'auto'; end
- if ~isfield(cfg, 'maskclipmax'); cfg.maskclipmax = 'auto'; end
- if ~isfield(cfg, 'maskclipsym'); cfg.maskclipsym = 'no'; end
- if ~isfield(cfg, 'maskmap'); cfg.maskmap = linspace(0,1,128); end
- if ~isfield(cfg, 'maskcolmin'); cfg.maskcolmin = 'auto'; end
- if ~isfield(cfg, 'maskcolmax'); cfg.maskcolmax = 'auto'; end
- if ~isfield(cfg, 'maskparameter');cfg.maskparameter = []; end
- % perform some checks on the configuration for backward compatibility
- if ~isfield(cfg, 'funparameter') && isfield(ininterp, 'source')
- % if present, the default behavior should be to use this field for plotting
- cfg.funparameter = 'source';
- end
- % make the selection of functional and mask data consistent with the data
- cfg.funparameter = parameterselection(cfg.funparameter, ininterp);
- cfg.maskparameter = parameterselection(cfg.maskparameter, ininterp);
- % only a single parameter should be selected
- try, cfg.funparameter = cfg.funparameter{1}; end
- try, cfg.maskparameter = cfg.maskparameter{1}; end
- % check anatomical data
- if isfield(ininterp,'anatomy');
- interp.anatomy = reshape(ininterp.anatomy, ininterp.dim);
- else
- error('no anatomical data supplied');
- end
- % check functional data
- if ~isempty(cfg.funparameter)
- interp.source = double(reshape(getsubfield(ininterp, cfg.funparameter), ininterp.dim));
- else
- error('no functional data supplied');
- end
- % check mask data
- if ~isempty(cfg.maskparameter)
- interp.mask = double(reshape(getsubfield(ininterp,cfg.maskparameter), ininterp.dim));
- maskdat = 1;
- else
- fprintf('no opacity mask data supplied\n');
- interp.mask = [];
- maskdat = 0;
- end
- % only work with the copy of the relevant parameters in "interp"
- clear ininterp;
- % convert anatomy data type and optimize contrast
- if isa(interp.anatomy, 'uint8') || isa(interp.anatomy, 'uint16')
- fprintf('converting anatomy to floating point values...');
- interp.anatomy = double(interp.anatomy);
- fprintf('done\n');
- end
- fprintf('optimizing contrast of anatomical data ...');
- minana = min(interp.anatomy(:));
- maxana = max(interp.anatomy(:));
- interp.anatomy = (interp.anatomy-minana)./(maxana-minana);
- fprintf('done\n');
- % store original data if 'interactive' mode
- if strcmp(cfg.interactive,'yes')
- data.source = interp.source;
- end
- % place markers
- marker = zeros(size(interp.anatomy));
- if ~isempty(cfg.marker)
- fprintf('placing markers ...');
- [x,y,z] = ndgrid([1:size(interp.anatomy,1)],[1:size(interp.anatomy,2)],[1:size(interp.anatomy,3)]);
- for imarker = 1:size(cfg.marker,1)
- marker(find(sqrt((x-cfg.marker(imarker,1)).^2 + (y-cfg.marker(imarker,2)).^2 + (z-cfg.marker(imarker,3)).^2)<=cfg.markersize)) = 1;
- end
- fprintf('done\n');
- end
- % shift dimensions
- fprintf('sorting dimensions...');
- interp.anatomy = shiftdim(interp.anatomy,cfg.dim-1);
- interp.source = shiftdim(interp.source,cfg.dim-1);
- interp.mask = shiftdim(interp.mask,cfg.dim-1);
- marker = shiftdim(marker,cfg.dim-1);
- fprintf('done\n');
- % flip dimensions
- if strcmp(cfg.flipdim,'yes')
- fprintf('flipping dimensions...');
- interp.anatomy = flipdim(interp.anatomy,1);
- interp.source = flipdim(interp.source,1);
- interp.mask = flipdim(interp.mask,1);
- marker = flipdim(marker,1);
- fprintf('done\n');
- end
- % set slice space
- if ischar(cfg.spacemin)
- fprintf('setting first slice position...');
- spacemin = min(find(~isnan(max(max(interp.source,[],3),[],2))));
- fprintf('%d...done\n',spacemin);
- else
- spacemin = cfg.spacemin;
- end
- if ischar(cfg.spacemax)
- fprintf('setting last slice position...');
- spacemax = max(find(~isnan(max(max(interp.source,[],3),[],2))));
- fprintf('%d...done\n',spacemax);
- else
- spacemax = cfg.spacemax;
- end
- % clip funtional data
- if ~ischar(cfg.clipmin)
- fprintf('clipping functional minimum...');
- switch cfg.clipsym
- case 'no'
- interp.source(find(interp.source<cfg.clipmin)) = nan;
- case 'yes'
- interp.source(find(abs(interp.source)<cfg.clipmin)) = nan;
- end
- fprintf('done\n');
- end
- if ~ischar(cfg.clipmax)
- fprintf('clipping functional maximum...');
- switch cfg.clipsym
- case 'no'
- interp.source(find(interp.source>cfg.clipmax)) = nan;
- case 'yes'
- interp.source(find(abs(interp.source)>cfg.clipmax)) = nan;
- end
- fprintf('done\n');
- end
- % clip mask data
- if maskdat
- if ~ischar(cfg.maskclipmin)
- fprintf('clipping mask minimum...');
- switch cfg.maskclipsym
- case 'no'
- interp.mask(find(interp.mask<cfg.maskclipmin)) = nan;
- case 'yes'
- interp.mask(find(abs(interp.mask)<cfg.maskclipmin)) = nan;
- end
- fprintf('done\n');
- end
- if ~ischar(cfg.maskclipmax)
- fprintf('clipping mask maximum...');
- switch cfg.maskclipsym
- case 'no'
- interp.mask(find(interp.mask>cfg.maskclipmax)) = nan;
- case 'yes'
- interp.mask(find(abs(interp.mask)>cfg.maskclipmax)) = nan;
- end
- fprintf('done\n');
- end
- end
- % scale functional data
- fprintf('scaling functional data...');
- fmin = min(interp.source(:));
- fmax = max(interp.source(:));
- if ~ischar(cfg.colmin)
- fcolmin = cfg.colmin;
- else
- if sign(fmin)==sign(fmax)
- fcolmin = fmin;
- else
- fcolmin = -max(abs([fmin,fmax]));
- end
- end
- if ~ischar(cfg.colmax)
- fcolmax = cfg.colmax;
- else
- if sign(fmin)==sign(fmax)
- fcolmax = fmax;
- else
- fcolmax = max(abs([fmin,fmax]));
- end
- end
- interp.source = (interp.source-fcolmin)./(fcolmax-fcolmin);
- if ~ischar(cfg.colmax)
- interp.source(find(interp.source>1)) = 1;
- end
- if ~ischar(cfg.colmin)
- interp.source(find(interp.source<0)) = 0;
- end
- fprintf('done\n');
- % scale mask data
- if maskdat
- fprintf('scaling mask data...');
- fmin = min(interp.mask(:));
- fmax = max(interp.mask(:));
- if ~ischar(cfg.maskcolmin)
- mcolmin = cfg.maskcolmin;
- else
- if sign(fmin)==sign(fmax)
- mcolmin = fmin;
- else
- mcolmin = -max(abs([fmin,fmax]));
- end
- end
- if ~ischar(cfg.maskcolmax)
- mcolmax = cfg.maskcolmax;
- else
- if sign(fmin)==sign(fmax)
- mcolmax = fmax;
- else
- mcolmax = max(abs([fmin,fmax]));
- end
- end
- interp.mask = (interp.mask-mcolmin)./(mcolmax-mcolmin);
- if ~ischar(cfg.maskcolmax)
- interp.mask(find(interp.mask>1)) = 1;
- end
- if ~ischar(cfg.maskcolmin)
- interp.mask(find(interp.mask<0)) = 0;
- end
- fprintf('done\n');
- end
- % merge anatomy, functional data and mask
- fprintf('constructing overlay...');
- if ischar(cfg.colormap)
- % replace string by colormap using standard Matlab function
- cfg.colormap = colormap(cfg.colormap);
- end
- cmap = cfg.colormap;
- cmaplength = size(cmap,1);
- maskmap = cfg.maskmap(:);
- maskmaplength = size(maskmap,1);
- indslice = round(linspace(spacemin,spacemax,cfg.nslices));
- nvox1 = length(1:cfg.resample:size(interp.anatomy,2));
- nvox2 = length(1:cfg.resample:size(interp.anatomy,3));
- if mod(cfg.rotate,2)
- dummy = nvox1;
- nvox1 = nvox2;
- nvox2 = dummy;
- end
- out = zeros(nvox1,nvox2,3,cfg.nslices);
- for islice = 1:cfg.nslices
- dummy1 = squeeze(interp.anatomy(indslice(islice),1:cfg.resample:end,1:cfg.resample:end));
- dummy2 = squeeze(interp.source(indslice(islice),1:cfg.resample:end,1:cfg.resample:end));
- indmarker = find(squeeze(marker(indslice(islice),1:cfg.resample:end,1:cfg.resample:end)));
- indsource = find(~isnan(dummy2));
- if maskdat
- dummymask = squeeze(interp.mask(indslice(islice),1:cfg.resample:end,1:cfg.resample:end));
- indsource = find(~isnan(dummy2) & ~isnan(dummymask));
- end
- for icol = 1:3
- dummy3 = dummy1;
- if not(maskdat)
- if ~ischar(cfg.alpha)
- try
- dummy3(indsource) = ...
- (1-cfg.alpha) * dummy3(indsource) + ...
- cfg.alpha * cmap(round(dummy2(indsource)*(cmaplength-1))+1,icol);
- end
- else
- try
- dummy3(indsource) = ...
- (1-dummy2(indsource)) .* dummy3(indsource) + ...
- dummy2(indsource) .* cmap(round(dummy2(indsource)*(cmaplength-1))+1,icol);
- end
- end
- else
- dummy3(indsource) = ...
- (1-maskmap(round(dummymask(indsource)*(maskmaplength-1))+1)).* ...
- dummy3(indsource) + ...
- maskmap(round(dummymask(indsource)*(maskmaplength-1))+1) .* ...
- cmap(round(dummy2(indsource)*(cmaplength-1))+1,icol);
- end
- dummy3(indmarker) = cfg.markercolor(icol);
- out(:,:,icol,islice) = rot90(dummy3,cfg.rotate);
- end
- if strcmp(cfg.whitebg,'yes')
- bgmask = zeros(nvox1,nvox2);
- bgmask(find(conv2(mean(out(:,:,:,islice),3),ones(round((nvox1+nvox2)/8))/(round((nvox1+nvox2)/8).^2),'same')<0.1)) = 1;
- for icol = 1:3
- out(:,:,icol,islice) = bgmask.*ones(nvox1,nvox2) + (1-bgmask).* out(:,:,icol,islice);
- end
- end
- end
- fprintf('done\n');
- clf;
- fprintf('plotting...');
- axes('position',[0.9 0.3 0.02 0.4]);
- image(permute(cmap,[1 3 2]));
- set(gca,'YAxisLocation','right');
- set(gca,'XTick',[]);
- set(gca,'YDir','normal');
- set(gca,'YTick',linspace(1,cmaplength,5));
- set(gca,'YTickLabel',linspace(fcolmin,fcolmax,5));
- set(gca,'Box','on');
- axes('position',[0.01 0.01 0.88 0.90]);
- [h,nrows,ncols]=slicemon(out);
- xlim=get(gca,'XLim');
- ylim=get(gca,'YLim');
- text(diff(xlim)/2,-diff(ylim)/100,cfg.title,'HorizontalAlignment','center','Interpreter','none');
- drawnow;
- fprintf('done\n');
- if nargout > 0
- outim=get(h,'CData');
- end
- if strcmp(cfg.interactive,'yes')
- data.sin = size(interp.source);
- data.nrows = nrows;
- data.ncols = ncols;
- data.out = out;
- data.indslice = indslice;
- data.cfg = cfg;
- data.hfig = gcf;
- uicontrol('Units','norm', 'Position', [0.9 0.2 0.08 0.05], 'Style','pushbutton', 'String','coords',...
- 'Callback',@getcoords,'FontSize',7);
- data.hcoords = uicontrol('Units','norm', 'Position', [0.9 0.05 0.08 0.13], 'Style','text', 'String','','HorizontalAlign','left','FontSize',7);
- guidata(data.hfig,data);
- end
- % ---------------- subfunctions ----------------
- function getcoords(h,eventdata,handles,varargin)
- data = guidata(gcf);
- [xi,yi] = ginput(1);
- co(2,1) = round(mod(yi,size(data.out,1)));
- co(3,1) = round(mod(xi,size(data.out,2)));
- switch mod(data.cfg.rotate,4)
- case 1,
- t1 = co(2);
- co(2) = co(3);
- co(3) = data.sin(3)-t1;
- case 2,
- co(2) = data.sin(2)-co(2);
- co(3) = data.sin(3)-co(3);
- case 3,
- t1 = co(3);
- co(3) = co(2);
- co(2) = data.sin(2)-t1;
- end
- try
- co(1) = data.indslice(fix(xi/size(data.out,2)) + fix(yi/size(data.out,1))*data.ncols + 1);
- catch
- co(1) = NaN;
- end
- if strcmp(data.cfg.flipdim, 'yes')
- co(1) = data.sin(1) - co(1) + 1;
- end
- co = co(:);
- co(2:3) = round(co(2:3)*data.cfg.resample);
- for ishift = 1:data.cfg.dim-1
- co = [co(3);co(1);co(2)];
- end
- set(data.hcoords,'String',sprintf('1: %d\n2: %d\n3: %d\nf: %0.4f',co(1),co(2),co(3),data.source(co(1),co(2),co(3))));
- function [h,nrows,ncols] = slicemon(a) % display the montage w/o image_toolbox
- siz = [size(a,1) size(a,2) size(a,4)];
- nn = sqrt(prod(siz))/siz(2);
- mm = siz(3)/nn;
- if (ceil(nn)-nn) < (ceil(mm)-mm),
- nn = ceil(nn); mm = ceil(siz(3)/nn);
- else
- mm = ceil(mm); nn = ceil(siz(3)/mm);
- end
- b = a(1,1);
- b(1,1) = 0;
- b = repmat(b, [mm*siz(1), nn*siz(2), size(a,3), 1]);
- rows = 1:siz(1); cols = 1:siz(2);
- for i=0:mm-1,
- for j=0:nn-1,
- k = j+i*nn+1;
- if k<=siz(3),
- b(rows+i*siz(1),cols+j*siz(2),:) = a(:,:,:,k);
- end
- end
- end
- hh = image(b);
- axis image;
- box off;
- set(gca,'XTick',[],'YTick',[],'Visible','off');
- if nargout > 0
- h = hh;
- nrows = mm;
- ncols = nn;
- end