CNC
/
fMRI_scripts
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

The alpha version. The first working version.

This commit is contained in:
p147685 2020-04-16 16:40:31 +02:00
parent 3b9a1b9444
commit 2399e5095f
2 changed files with 325 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

322
GIFT_Clean/GiftIcaReport.m
View File

@ -1,21 +1,49 @@
classdef GiftIcaReport<handle
%GIFTICAREPORT
%GIFTICAREPORT
% This Class will be a little tool to print to file images and timecourses
% to quickly inspect the result of GIFT for a cleanup of fMRI data
% some guidlines can be found below!
% https://doi.org/10.1016/j.neuroimage.2016.12.036
% to do: auto scale images
% to do: check slices to be plotted actually exist!!
% to do: create the output
properties
ICAfile = {}; % .nii file containing the ICA images
TCfile = {}; % .nii file containing the TimeCourse data
Orientation = 'sag';% orientation of images to create
SliceNumbers = []; % slice numbers to plot
ICAfile = {}; % .nii file containing the ICA images
TCfile = {}; % .nii file containing the TimeCourse data
Orientation = 'sag'; % orientation of images to create
DirOut = '.' %default uses the current directory
FileOut = 'Report.png'; %default report name
Layout = 'north'; %define the location of the TC plot
Colormap=gray; % colormap for images
end
properties (Dependent=true)
PlotSliceNumbers % slice numbers to plot
NPlots % number of slices to plot
ImageLimits %defines the minium and max value for colorbar
end
properties (Access=private)
OrientationOptions={'tra','sag','cor','all'};
SPMPath
FigureObject
ImageData
TimeCourseData
VData
VTimeCourse
PData
PTimeCourse
SliceNumbersInternal
NSlicesInternal = 10;
TCPlot %subplot containing the time course
SLPlot %subplot containing the
CLim = [0,1]; %limits for plotting
%% the matlab reporting gives an licence error.
% ReportObj % the place holder for the Report object
end
properties (Constant,Access=private)
OrientationOptions={'tra','sag','cor','all'}; %option all still needs to be implemented
end
methods %constructor, set and get functions
@ -23,6 +51,18 @@ classdef GiftIcaReport<handle
%GIFTICAREPORT Construct an instance of this class
% Detailed explanation goes here
obj.SPMPath=[]; %check if SPM is in the path
obj.FigureObject=figure; %create a figure
%% the matlab reporting gives an licence error at the CNC
%import mlreportgen.report.*
%obj.ReportObj = Report(fullfile(obj.DirOut,obj.FileOut));
end
function delete(obj)
if ~isempty(obj.FigureObject)
if isa(obj.FigureObject,'matlab.ui.Figure')
close(obj.FigureObject);
end
end
end
function set.SPMPath(obj,~)
@ -30,48 +70,294 @@ classdef GiftIcaReport<handle
if isempty(SPMFound) %SPM not found
id='GiftIcaReport:set_SPMPath:SPM_not_found';
msg='SPM not found';
throw(MException(id,msg));
throw(MException(id,msg));
end
%to do add a check on spm version
obj.SPMPath=fileparts(SPMFound);
end
function set.TCfile(obj,val)
obj.TCfile=obj.setfile(val);
obj.TCfile=obj.setfile(val);
end
function set.ICAfile(obj,val)
obj.ICAfile=obj.setfile(val);
obj.ICAfile=obj.setfile(val);
end
function set.Orientation(obj,val)
switch lower(val)
case obj.OrientationOptions %#ok<MCSUP>
case obj.OrientationOptions
obj.Orientation=val;
otherwise
id = 'GiftIcaReport:set_orientation:invalid_orientation';
msg = sprintf('%s\n%s\n%s\t%s\t%s\n',...
'invalid option for orientation',...
'pick one of the following options',...
obj.OrientationOptions{:}); %#ok<MCSUP>
obj.OrientationOptions{:});
throw(MException(id,msg));
end
end
function set.SliceNumbers(obj,val)
if isnumeric(val) && ~any(val<1)
obj.SliceNumbers=val;
function set.PlotSliceNumbers(obj,val)
if isempty(val)
if ~isempty(obj.ImageData) && ~isempty(obj.NPlots)
%autofill the plotslicenumbers
switch lower(obj.Orientation)
case 'tra'
nmax=size(obj.ImageData,3)
case 'cor'
nmax=size(obj.ImageData,2)
case 'sag'
nmax=size(obj.ImageData,1);
otherwise
return;
end
obj.SliceNumbersInternal=round(linspace(1,nmax,obj.NPlots));
end
else
id = 'GiftIcaReport:set_slicenumbers:invalid_input';
if isnumeric(val) && ~any(val<1)
obj.SliceNumbersInternal=val;
obj.NSlicesInternal=numel(val);
else
id = 'GiftIcaReport:set_slicenumbers:invalid_input';
msg = sprintf('invalid input for slicenumbers');
throw(MException(id,msg));
end
end
end
function val=get.PlotSliceNumbers(obj)
val=obj.SliceNumbersInternal;
end
function set.NPlots(obj,val)
if isnumeric(val) && numel(val)==1 && val>0
tmp=obj.SliceNumbersInternal;
if ~isempty(tmp) && val<=numel(tmp)
obj.PlotSliceNumbers=tmp(1:val);
else
obj.NSlicesInternal=val;
obj.PlotSliceNumbers=[];
end
else
id = 'GiftIcaReport:set_NSlices:invalid_input';
msg = sprintf('invalid input for slicenumbers');
throw(MException(id,msg));
end
end
function val=get.NPlots(obj)
val=obj.NSlicesInternal;
end
function set.ImageLimits(obj,val)
if isnumeric(val) && numel(val)==2
obj.CLim=transpose(sort(val(:))); %guarante it is a row vector
else
id = 'GiftIcaReport:setImageLimits:invalid_input';
msg = sprintf('Invalid input, set limits to [min max]');
throw(MException(id,msg));
end
end
function val=get.ImageLimits(obj)
val=obj.CLim;
end
function set.Layout(obj,val)
switch lower(val)
case {'north','south','east','west',''}
obj.Layout=val;
otherwise
id = 'GiftIcaReport:setLayout:invalid_input';
msg = sprintf('Invalid input for Layout.\n Use north east south, west or empty value');
throw(MException(id,msg));
end
end
end
methods %user callable functions
function LoadData(obj)
% read ICA file
[d,f,~]=fileparts(obj.ICAfile{1});
obj.PData=spm_select('ExtFPList',d,f,1:10000);
obj.VData=spm_vol(obj.PData);
obj.ImageData=spm_read_vols(obj.VData);
% read TimeCourse file
[d,f,~]=fileparts(obj.TCfile{1});
obj.PTimeCourse=spm_select('ExtFPList',d,f,1:10000);
obj.VTimeCourse=spm_vol(obj.PTimeCourse);
obj.TimeCourseData=spm_read_vols(obj.VTimeCourse);
if isempty(obj.PlotSliceNumbers) %auto set the slice numbers
obj.PlotSliceNumbers([]);
end
end
function BuildFigure(obj)
if isempty(obj.FigureObject)
obj.FigureObject=figure;
end
figure(obj.FigureObject);
clf;
[nrow,ncol] = obj.getrowcol;
obj.TCPlot={};
obj.SLPlot=cell(nrow,ncol);
switch lower(obj.Layout)
case {'north','n'}
r=nrow+1;
c=ncol;
PlRow=setdiff(1:r,1);
PlCol=1:c;
obj.TCPlot={subplot(r,1,1)};
case {'east','e'}
r=nrow;
c=ncol+1;
PlRow=1:r;
PlCol=setdiff(1:c,1);
obj.TCPlot={subplot(1,c,c)};
case {'south','s'}
r=nrow+1;
c=ncol;
PlRow=setdiff(1:r,r);
PlCol=1:c;
obj.TCPlot={subplot(r,1,r)};
case {'west','w'}
r=nrow;
c=ncol+1;
PlRow=1:r;
PlCol=setdiff(1:c,1);
obj.TCPlot={subplot(1,c,1)};
otherwise %center option
r=nrow+1;
c=ncol;
PlRow=setdiff(1:r,floor(r/2));
PlCol=1:c;
obj.TCPlot={subplot(r,1,floor(r/2))};
end
[x,y]=ndgrid(PlRow,PlCol);
x=transpose(x);
y=transpose(y);
for k=1:obj.NPlots
p=(x(k)-1)*c+y(k);
obj.SLPlot{k}=subplot(r,c,p);
end
end
function PlotTimeCourse(obj,IC)
if isempty(obj.TimeCourseData)
id = 'GiftIcaReport:PlotTC:NoDataToPlot';
msg = 'No time data to plot';
throw(MException(id,msg));
end
if IC>size(obj.TimeCourseData,2)
id = 'GiftIcaReport:PlotTC:InvalidComponent';
msg = 'invalid component';
throw(MException(id,msg));
end
axes(obj.TCPlot{1});
plot(obj.TimeCourseData(:,IC));
end
function PlotSlices(obj,IC)
if isempty(obj.ImageData)
id = 'GiftIcaReport:PlotTC:NoDataToPlot';
msg = 'No time data to plot';
throw(MException(id,msg));
end
if IC>size(obj.ImageData,4)
id = 'GiftIcaReport:PlotTC:InvalidComponent';
msg = 'invalid component';
throw(MException(id,msg));
end
for k=1:obj.NPlots
axes(obj.SLPlot{k}); %#ok<LAXES>
switch lower(obj.Orientation)
case 'cor'
if obj.PlotSliceNumbers(k)>size(obj.ImageData,2);continue;end
im=squeeze(obj.ImageData(:,obj.PlotSliceNumbers(k),:,IC));
case 'tra'
if obj.PlotSliceNumbers(k)>size(obj.ImageData,3);continue;end
im=squeeze(obj.ImageData(:,:,obj.PlotSliceNumbers(k),IC));
case 'sag'
if obj.PlotSliceNumbers(k)>size(obj.ImageData,1);continue;end
im=squeeze(obj.ImageData(obj.PlotSliceNumbers(k),:,:,IC));
case 'all'
end
imagesc(im,obj.CLim);
camroll(90)
colormap(obj.Colormap);
axis image
axis off
end
end
function GenReport(obj)
% load all info if not already present
if isempty(obj.ImageData)||isempty(obj.TimeCourseData)
obj.LoadData;
end
obj.BuildFigure;
% loop over IC
for IC=1:size(obj.TimeCourseData,2)
obj.PlotTimeCourse(IC);
obj.PlotSlices(IC);
obj.TCPlot{1}.Title.String= sprintf('IC:%.4d',IC);
[~,f,e]=fileparts(obj.FileOut);
switch lower(e)
case '.png'
fout=fullfile(obj.DirOut,sprintf('%s_%.4d.png',f,IC));
type=-'png';
otherwise
id = 'GiftIcaReport:GenReport:invalidFileType';
msg = sprintf('Filetype not implemented');
throw(MException(id,msg));
end
saveas(gcf,fout,type);
end
%% this block of code gives a licence error. weird as I could create the object ?
% if isempty(obj.ReportObj)
% id = 'GiftIcaReport:GenReport:EmptyReportObject';
% msg = 'No report object set in constructor of this object';
% throw(MException(id,msg));
% end
% if ~isa(obj.ReportObj,'mlreportgen.report.Report')
% id = 'GiftIcaReport:GenReport:WrongClassReportObject';
% msg = 'No report object set in constructor of this object';
% throw(MException(id,msg));
% end
% import mlreportgen.report.*
% %open(obj.ReportObj)
% % generate the figure
% surf(peaks) %just for testing
% add(obj.ReportObj,Figure)
% %close(obj.ReportObj)
% disp(obj.ReportObj)
%%
end
end
methods %private methods
function out = setfile(~,val)
if isempty(val)
function [nrow,ncol]=getrowcol(obj)
p=factor(obj.NPlots);
while numel(p)>2
p(1)=p(1)*p(end);
p(end)=[];
p=sort(p);
end
if numel(p)==1 % catch primes
p(2)=1;
p=sort(p);
end
if 9*p(2)>16*p(1) %very elongated plot
ncol=ceil(sqrt(obj.NPlots));
nrow=ceil(obj.NPlots/ncol);
else
nrow=p(1);
ncol=p(2);
end
end
function out = setfile(~,val)
if isempty(val)
out={};
return
end
@ -98,7 +384,7 @@ classdef GiftIcaReport<handle
return
end
end
end
end
end
end

26
GIFT_Clean/UnityTestGiftIcaReport.m
View File

@ -8,8 +8,9 @@ P={'X:\My Documents\MATLAB\Projects\fMRI_script\GIFT_Clean\TestData\g_sub01_comp
if not(exist(P{1},'file'));error('Test file not found,unity test not performed');end
P2={'X:\My Documents\MATLAB\Projects\fMRI_script\GIFT_Clean\TestData\g_sub01_timecourses_ica_s1_.nii'};
if not(exist(P2{1},'file'));error('Test file not found,unity test not performed');end
%% basic call
%% basic call and close
Q = GiftIcaReport;
clear('Q')
%% basic call with no SPM error
try
tmp=fileparts(which('spm'));
@ -25,6 +26,7 @@ catch ME
end
end
%% test set ICAfile
Q = GiftIcaReport;
Q.ICAfile=[]; %set to empty
try
Q.ICAfile = 'junk';
@ -41,7 +43,7 @@ Q.ICAfile=P{1};
%Q.ICAfile='ui'; % request a user interface
%% test set TC file
%Q.TCfile='ui';
Q.ICAfile=P2;
Q.TCfile=P2;
%% test set orientation
try
Q.Orientation='junk';
@ -53,13 +55,14 @@ catch ME
disp(ME)
end
end
Q.Orientation='all';
Q.Orientation='tra';
Q.Orientation='cor';
Q.Orientation='sag';
%% test setting slice numbers I
try
Q.SliceNumbers=0:5;
Q.PlotSliceNumbers=0:5;
catch ME
if ~strcmpi(ME.identifier,'GiftIcaReport:set_slicenumbers:invalid_input')
rethrow(ME);
@ -70,7 +73,7 @@ catch ME
end
%% test setting slice numbers II
try
Q.SliceNumbers='junk';
Q.PlotSliceNumbers='junk';
catch ME
if ~strcmpi(ME.identifier,'GiftIcaReport:set_slicenumbers:invalid_input')
rethrow(ME);
@ -79,4 +82,17 @@ catch ME
disp(ME);
end
end
Q.SliceNumbers=1:2:100;
%% test ImageLimits
Q.ImageLimits = [0 1]; %to do check on properly working error
%% test NPlots set
Q.NPlots=25;
disp(Q.PlotSliceNumbers);
%
%Q.PlotSliceNumbers=1:2:100;
Q.Layout='north'; % to do build error handling check
Q.LoadData;
Q.NPlots=9; %if PlotSliceNumbers is empty, auto set the slices to plot
Q.PlotSliceNumbers=[]; %auto set the slices to plot
Q.BuildFigure;
Q.GenReport;