rest_powerspectrum.m

While resting-state fMRI is drawing more and more attention, there has not been a software for its d

function [varargout]=rest_powerspectrum(AOperation, varargin)
%Visualize the Power Spectrum and the Time course of user specified voxel whose coordinates could be set with SliceViewer By Xiao-Wei Song
%------------------------------------------------------------------------------------------------------------------------------
%	Copyright(c) 2007~2010
%	State Key Laboratory of Cognitive Neuroscience and Learning in Beijing Normal University
%	Written by Xiao-Wei Song 
%	http://resting-fmri.sourceforge.net
%------------------------------------------------------------------------------------------------------------------------------
% Draw the time course fluctions at the specific position in the 3D+time Dataset
% 	<a href="[e-mail]">Mail to Author</a>: Xiaowei Song
%	Version=1.0;
%	Release=20070903;

if nargin<1, help rest_powerspectrum; return; end

persistent ReHo_WaveGraph_Cfg; % run-time persistent config
%ADataDir, AVoxelPosition, ABandLow, ABandHigh
if ~mislocked(mfilename),mlock; end

switch upper(AOperation),
case 'SHOWFLUCTUATION',	%ShowFluctuation
	if nargin~=5, error('Usage: rest_powerspectrum(''ShowFluctuation'', ADataDir, AVoxelPosition, ASamplePeriod, ABandRange);'); end
	ADataDir 		=varargin{1};
    AVoxelPosition 	=varargin{2};
	ASamplePeriod  	=varargin{3};
	ABandRange  	=varargin{4};
	theFig =ExistDisplayFigure(ReHo_WaveGraph_Cfg, ADataDir);
	isExistFig =rest_misc( 'ForceCheckExistFigure' , theFig);	%Force check whether the figure exist
	if ~isExistFig
		%the specific image didn't exist, so I create one
		ReHo_WaveGraph_Cfg.Config(1+GetDisplayCount(ReHo_WaveGraph_Cfg)) =InitControls(ADataDir, AVoxelPosition, ASamplePeriod, ABandRange);
		%To Force display 
		theFig =ReHo_WaveGraph_Cfg.Config(GetDisplayCount(ReHo_WaveGraph_Cfg)).hFig;
	else
		%Update the VoxelPosition, SamplePeriod, BandRange
		theCardinal =ExistDisplay(ReHo_WaveGraph_Cfg, ADataDir);
		if theCardinal>0,			
			ReHo_WaveGraph_Cfg.Config(theCardinal).VoxelPosition =AVoxelPosition;
			%Add the voxel to the voxel array
			ReHo_WaveGraph_Cfg.Config(theCardinal).VoxelArray =[AVoxelPosition; ReHo_WaveGraph_Cfg.Config(theCardinal).VoxelArray];
			
			ReHo_WaveGraph_Cfg.Config(theCardinal).SamplePeriod	 =ASamplePeriod;
			ReHo_WaveGraph_Cfg.Config(theCardinal).BandRange     =ABandRange;			
			% if get(ReHo_WaveGraph_Cfg.Config(theCardinal).hDetrendBeforeFFT, 'Value'),
				% ReHo_WaveGraph_Cfg.Config(theCardinal).DetrendBeforeFFT ='Yes';
			% else
				% ReHo_WaveGraph_Cfg.Config(theCardinal).DetrendBeforeFFT ='No';
			% end
			PlotFluctuation(ReHo_WaveGraph_Cfg.Config(theCardinal));
		end
	end
	figure(theFig);
	varargout{1} =theFig;	
	
case 'ONFIGUREDELETE',				%OnFigureDelete
	if nargin~=2, error('Usage: rest_powerspectrum(''OnFigureDelete'', ADataDir);'); end
	
	ADataDir =varargin{1};
	ReHo_WaveGraph_Cfg =DeleteFigure(ReHo_WaveGraph_Cfg, ADataDir);	

case 'ONFIGURERESIZE', 		%OnFigureResize
	if nargin~=2, error('Usage: rest_powerspectrum(''ResizeFigure'', ADataDir);'); end
	ADataDir 		=varargin{1};
	theFig =ExistDisplayFigure(ReHo_WaveGraph_Cfg, ADataDir);
	isExistFig =rest_misc( 'ForceCheckExistFigure' , theFig);	%Force check whether the figure exist
	if isExistFig
		theCardinal =ExistDisplay(ReHo_WaveGraph_Cfg, ADataDir);
		if theCardinal>0			
			OnFigureResize(ReHo_WaveGraph_Cfg.Config(theCardinal));
		end
	end
	
case 'ONOPTIONDETRENDBEFOREFFT', 		%OnOptionDetrendBeforeFFT	
	if nargin~=2, error('Usage: rest_powerspectrum(''OnOptionDetrendBeforeFFT'', ADataDir);'); end
	ADataDir 		=varargin{1};
	theFig =ExistDisplayFigure(ReHo_WaveGraph_Cfg, ADataDir);
	isExistFig =rest_misc( 'ForceCheckExistFigure' , theFig);	%Force check whether the figure exist
	if isExistFig
		theCardinal =ExistDisplay(ReHo_WaveGraph_Cfg, ADataDir);
		if theCardinal>0,
			if get(ReHo_WaveGraph_Cfg.Config(theCardinal).hDetrendBeforeFFT, 'Value'),
				ReHo_WaveGraph_Cfg.Config(theCardinal).DetrendBeforeFFT ='Yes';
			else
				ReHo_WaveGraph_Cfg.Config(theCardinal).DetrendBeforeFFT ='No';
			end
			PlotFluctuation(ReHo_WaveGraph_Cfg.Config(theCardinal));
		end
	end
	varargout{1} =theFig;
	
case 'QUITALLPOWERSPECTRUM',		%QuitAllPowerSpectrum
	if nargin~=1, error('Usage: rest_powerspectrum(''QuitAllPowerSpectrum'');'); end
	for x=1:GetDisplayCount(ReHo_WaveGraph_Cfg), 
		rest_powerspectrum('OnFigureDelete', ReHo_WaveGraph_Cfg.Config(x).DataDir);
	end
	clear ReHo_WaveGraph_Cfg;
		
	
otherwise
end

function Result =InitControls(ADataDir, AVoxelPosition, ASamplePeriod, ABandRange)	
	%Construct 3D+time dataset	
	[theDataset, theVoxelSize, theImgFileList, Origin]=rest_to4d(ADataDir);
	[nDim1, nDim2, nDim3, nDim4] =size(theDataset);
	rest_waitbar;
	
	%Initialization 20070525
	theFig =figure('Units', 'pixel', 'MenuBar', 'none', 'Toolbar', 'figure', ...
				'NumberTitle', 'off', 'Name', ADataDir, ...
				'DeleteFcn', sprintf('rest_powerspectrum(''OnFigureDelete'', ''%s'');', rest_misc('ReplaceSingleQuota', ADataDir)) , ...
				'ResizeFcn', sprintf('rest_powerspectrum(''OnFigureResize'', ''%s'');', rest_misc('ReplaceSingleQuota', ADataDir)) );
	MarginX =10; MarginY =10;
	OffsetX =3*MarginX; 	OffsetY =MarginY;
		
	%Create Axes and lines and images			
	hVoxelPosition  =uicontrol(theFig, 'Style','radiobutton', 'Units','pixels', ...
							  'String', sprintf('(%d,%d,%d)',AVoxelPosition), ...
							  'BackgroundColor', get(theFig,'Color'), ...
							  'HorizontalAlignment', 'left', ...
							  'Position',[OffsetX, OffsetY, 150,15]);
		
	%Create options to draw the Power Spectrum
	OffsetX =3*MarginX; 	OffsetY =MarginY +25;				
	hDetrendBeforeFFT =uicontrol(theFig, 'Style','checkbox', 'Units','pixels', ...
							  'String', 'Remove Linear Trend Before Power Spectrum', ...
							  'BackgroundColor', get(theFig,'Color'), ...
							  'HorizontalAlignment', 'left', ...
							  'Value', 1, ...
							  'Callback', ...
							  sprintf('rest_powerspectrum(''OnOptionDetrendBeforeFFT'', ''%s'');', rest_misc('ReplaceSingleQuota', ADataDir)), ...
							  'Position',[OffsetX, OffsetY, 250,15]);
		
	%Create Axes
	OffsetX =6*MarginX; 	OffsetY =3*MarginY +15 +3*MarginY;				
	hAxesTimeCourse	=axes('Parent', theFig, ...
						  'Units', 'pixel', 'DrawMode','fast', ...
						  'Position', [OffsetX OffsetY 3*nDim4 150]);
						  
	OffsetX =6*MarginX; 	OffsetY =3*MarginY +15 +6*MarginY +150 +6*MarginY;
	hAxesAmplitude	=axes('Parent', theFig, ...
						  'Units', 'pixel', 'DrawMode','fast', ...
						  'Position', [OffsetX OffsetY 3*nDim4 150], ...
						  'NextPlot', 'replacechildren');						  
	%Save to config
	AConfig.hFig			=theFig;			%handle of the config
	%Save Axes's handles
	AConfig.hAxesTimeCourse 	=hAxesTimeCourse;
	AConfig.hAxesAmplitude 		=hAxesAmplitude;
	
	%Save Voxel position label handle
	AConfig.hVoxelPosition 		=hVoxelPosition;
	% Save Options' handles
	AConfig.hDetrendBeforeFFT	=hDetrendBeforeFFT;
		
	%Save important variables
	AConfig.DataDir 		=ADataDir;		
	AConfig.VoxelPosition 	=AVoxelPosition;		%[x y z] the Voxel position in Volume current showing its amplitude and its time course
	AConfig.VoxelArray		=AVoxelPosition;		%20070718
	
	AConfig.Dataset 		=theDataset;
	AConfig.SamplePeriod 	=ASamplePeriod;		%TR
	AConfig.BandRange 		=ABandRange;		%[BandLow BandHigh] or [HighCutoff LowCutoff]
	AConfig.DetrendBeforeFFT='Yes';				% Remove linear trend before FFT
    Result =AConfig;
	
	
	%%Display Images
	PlotFluctuation(AConfig);
	
	%Resize figure width and height
	FigWidth  =6*MarginX + 3*nDim4 + 6*MarginX;
	FigHeight =3*MarginY +15 +6*MarginY +150 +6*MarginY + 150+ 3*MarginY ;
	thePos =get(theFig, 'Position');	
	theScreenSize =get(0,'ScreenSize');
	if thePos(1)>= theScreenSize(3)
		thePos(1) =theScreenSize(1);
	end	
	if 	(thePos(2) +FigHeight) +100>= theScreenSize(4)
		thePos(2) =theScreenSize(4) -FigHeight -100;
	end
	thePos =[thePos(1), thePos(2), FigWidth,FigHeight];	
	set(theFig, 'Position', thePos);	
	%Force resize the figure	
	OnFigureResize(AConfig);
	
	return;



function Result =DeleteFigure(AGlobalConfig, ADataDir)
	x =ExistDisplay(AGlobalConfig, ADataDir);
	if x>0
		theDisplayCount =GetDisplayCount(AGlobalConfig);
		isExistFig =rest_misc( 'ForceCheckExistFigure' , AGlobalConfig.Config(x).hFig);
		if isExistFig
			delete(AGlobalConfig.Config(x).hFig);
			if theDisplayCount>x
				for y=x:theDisplayCount-1
					AGlobalConfig.Config(x) =AGlobalConfig.Config(x+1);
                end
            end	
            AGlobalConfig.Config(theDisplayCount)=[];
		end	
	end
	Result =AGlobalConfig;
function Result =GetDisplayCount(AGlobalConfig)
%Get the Count of display, this program allow multi-view of brain like MRIcro
	if isempty(AGlobalConfig) || isempty(AGlobalConfig.Config),
		Result =0;		
	else
		Result =length(AGlobalConfig.Config);