diff --git a/ExampleInvocation.sh b/ExampleInvocation.sh index f988d3f..d407e70 100755 --- a/ExampleInvocation.sh +++ b/ExampleInvocation.sh @@ -2,6 +2,6 @@ ./RetinotopyAnalysis.sh \ --subject=123456 \ - --movie-files="movie1.mov@movie2.mov@movie3.mov@movie4.mov@movie5.mov" \ + --movie-files="movie1.mat@movie2.mat@movie3.mat@movie4.mat@movie5.mat" \ --image-files="hello.nii.gz@there.nii.gz@dolly.nii.gz" \ --behavior-files="hello_behavior.xml@there_behavior.xml@EMPTY" \ No newline at end of file diff --git a/MATLAB/CHPCsubmit.m b/MATLAB/CHPCsubmit.m new file mode 100644 index 0000000..ec1aaf1 --- /dev/null +++ b/MATLAB/CHPCsubmit.m @@ -0,0 +1,116 @@ +function jobids = CHPCsubmit(userid,jobname,scriptname,numvxs,beginix,endix,flags,varstoexclude) + +% function jobids = CHPCsubmit(userid,jobname,scriptname,numvxs,beginix,endix,flags,varstoexclude) +% +% is a string. this should be the user who is submitting the jobs. +% is a string with the name of the job +% is the name of the script (e.g. 'run_fitnonlinearmodel.sh') +% is the number of voxels to process in each job +% is starting voxel index +% is ending voxel index +% (optional) is a string with some qsub flags. +% Default: '-l nodes=1:ppn=1,walltime=4:00:00,vmem=2gb -q pe1_iD_4hr'. +% (optional) is a cell vector of variables to exclude +% when saving to the temporary .mat file +% +% submit jobs and return a vector of job ids. +% +% some special things we do: +% (1) we make available the caller's workspace to the jobs +% (2) we submit in chunks, being careful not to overload the queue + +% internal constants +maxsize = 10; % warn if workspace is bigger than this number of MB +maxjobs = 1500; % max jobs to have running at one time +jobchunk = 1000; % how many jobs to submit in each iteration +pausetime = 10; % seconds to pause when unexpected exit code +overloadtime = 60; % seconds to wait when overloaded before trying again + +% inputs +if ~exist('flags','var') || isempty(flags) + flags = '-l nodes=1:ppn=1,walltime=4:00:00,vmem=2gb -q pe1_iD_4hr'; % -q dque +end +if ~exist('varstoexclude','var') || isempty(varstoexclude) + varstoexclude = {}; +end + +% setup +mkdirquiet('~/sgeoutput'); + +%%%%% ensure uniqueness + +if exist(sprintf('~/mcc/job_%s.mat',jobname),'file') + error(sprintf('jobname %s already exists!',jobname)); +end + +%%%%% save the caller's workspace to the special file + +% check size of workspace +a = evalin('caller','whos'); +ok = ~ismember(cat(2,{a.name}),varstoexclude); +a = a(ok); +workspacesize = sum(cat(1,a.bytes)); + +% give warning if the workspace is big +if workspacesize/1000/1000 > maxsize + warning(sprintf('We are saving to disk a workspace that is larger than %.1f MB!',maxsize)); +end + +% save caller's workspace to the special .mat file +fprintf('saving workspace to disk...'); +evalin('caller',sprintf('saveexcept(''~/mcc/job_%s.mat'',%s);',jobname,cell2str(varstoexclude))); +fprintf('done.\n'); + +% do in chunks +jobids = []; +totalchunks = ceil((endix-beginix+1)/jobchunk); +for zz=1:totalchunks + fprintf('working on chunk %d of %d.\n',zz,totalchunks); + + %%%%% wait until not overloaded + + fprintf('checking to see if CHPC is overloaded.\n'); + while 1 + [s,r] = unix(sprintf('qstat -t -u %s | wc -l',userid)); + fprintf('exit status was %d.\n',s); + while s~=0 + fprintf('STATUS WASNT 0, WEIRD; TRY AGAIN IN A FEW SECONDS\n'); + pause(pausetime); + [s,r] = unix(sprintf('qstat -t -u %s | wc -l',userid)); + fprintf('exit status was %d.\n',s); + end + if str2double(r) > maxjobs + fprintf('overloaded, so waiting (%s).\n',datestr(now)); + pause(overloadtime); + else + break; + end + end + fprintf('CHPC is not overloaded right now, so let''s proceed.\n'); + + %%%%% issue the qsub call on chpc + + beginix0 = (zz-1)*jobchunk + beginix; + endix0 = min(beginix0 + jobchunk - 1,endix); + cmd = sprintf('qsub -t %d-%d -v MYSCRIPT="%s",MYARG1="~/mcc/job_%s.mat",MYARG2="%d" -N job_%s -o ~/sgeoutput/job_%s.o -e ~/sgeoutput/job_%s.e %s $HCPRETINODIR/mcc/matlabsge.sh', ... + beginix0,endix0,scriptname,jobname,numvxs,jobname,jobname,jobname,flags); + fprintf('this is the qsub command:\n\n%s\n\n',cmd); + fprintf('issuing the qsub command on chpc...'); + [s,r] = unix(cmd); + fprintf('exit status was %d.\n',s); + while s~=0 + fprintf('STATUS WASNT 0, WEIRD; TRY AGAIN IN A FEW SECONDS\n'); + pause(pausetime); + [s,r] = unix(cmd); + end + jobids = [jobids str2double(regexp(r,'\d+','match'))]; + fprintf('done.\n'); + + %%%%% wait to take effect + + if zz ~= totalchunks + fprintf('waiting before trying next chunk.\n'); + pause(overloadtime); + end + +end diff --git a/MATLAB/CHPCwait.m b/MATLAB/CHPCwait.m new file mode 100644 index 0000000..c1cd8f1 --- /dev/null +++ b/MATLAB/CHPCwait.m @@ -0,0 +1,62 @@ +function CHPCwait(jobnames,jobids,userid) + +% function CHPCwait(jobnames,jobids,userid) +% +% is a job name or a cell vector of job names +% is a vector of job IDs +% is a string referring to the user that submitted the jobs +% +% periodically check on job status and report to the command window. +% we exit when all of the jobs are completed. if any job errors are +% detected, we issue an error. + +% internal constants +pausetime = 60; % seconds to wait before checking again +weirdtime = 10; % seconds to wait if unexpected exit code + +% input +if ~iscell(jobnames) + jobnames = {jobnames}; +end + +% wait until jobs are done +fprintf('waiting for jobs to finish.\n'); +while 1 + pause(pausetime); + + % check if jobs are done + fprintf('checking if jobs are done (%s): ',datestr(now)); + temp = catcell(2,cellfun(@(x) [num2str(x) '|'],num2cell(jobids),'UniformOutput',0)); + str = sprintf('qstat -t -u %s | grep -E ''(%s)'' | grep -v '' C '' | wc -l',userid,temp(1:end-1)); + [s,r] = unix(str); + while (s~=0 && s~=1) + fprintf('STATUS WASNT 0 OR 1, WEIRD; TRY AGAIN IN A FEW SECONDS\n'); + pause(weirdtime); + [s,r] = unix(str); + end + + % report the number of jobs found + fprintf('%d jobs. ',str2double(r)); + if str2double(r) ~= 0 + fprintf('qstat is not empty.\n'); + else + fprintf('qstat is empty! checking for errors...'); + for jj=1:length(jobnames) + + strB = sprintf('cat ~/sgeoutput/job_%s.e*',jobnames{jj}); + [s,r] = unix(strB); + while s~=0 + fprintf('STATUS WASNT 0, WEIRD; TRY AGAIN IN A FEW SECONDS\n'); + pause(weirdtime); + [s,r] = unix(strB); + end + + if ~isempty(r) + error('Errors were found in the .e files for job_%s',jobnames{jj}); + end + end + fprintf('done!\n'); + break; + end +end +fprintf('ok, jobs are done!.\n'); diff --git a/MATLAB/RetinotopyAnalysis.m b/MATLAB/RetinotopyAnalysis.m new file mode 100644 index 0000000..1441cb2 --- /dev/null +++ b/MATLAB/RetinotopyAnalysis.m @@ -0,0 +1,236 @@ +function RetinotopyAnalysis(inputfile) + +% function RetinotopyAnalysis(inputfile) +% +% is the location of a .txt file that can be evaluated in order +% to define the following variables: +% is a string +% is a positive integer +% is a string +% is a string +% is a cell vector of .mat files. there should be 5 of them. +% is a cell vector of minimally pre-processed CIFTI files +% is a cell vector of behavioral XML files +% (optional) is +% 0 means regular full processing (involving CHPC machinery) +% 1 means super-fast processing (does not involve CHPC machinery at all) +% 2 means fast processing involving CHPC machinery, but analyzing only a few voxels +% Default: 0. +% +% Analyze the data and write outputs to: +% outpath/results/_XXX.dtseries.nii +% where XXX is the name of a quantity (e.g. angle). +% +% Regarding the image and behavior files, it is assumed that the number of both +% types of files is the same, that there is a one-to-one correspondence between +% the files, and that missing files are indicated by the string 'EMPTY'. +% +% If the image files have fewer than the expected number of time points, we assume +% that these missing time points occur at the end of the run. To handle these cases, +% we just truncate the corresponding stimulus sequence before performing the actual +% analysis. +% +% history: +% - 2015/03/11 - handle the case of data truncation; flexible format for loading the +% stimulus files; rename 'ret_summary' to 'behav_summary'. + +% wrap everything in a try-catch to ensure that MATLAB will exit +try + + %%%%% GET INPUTS + + % evaluate the text file + inputload = loadtext(inputfile); + for p=1:length(inputload) + eval(inputload{p}); + end + + % deal with inputs + if ~exist('debugmode','var') || isempty(debugmode) + debugmode = 0; + end + + %%%%% SETUP + + % internal constants + viewingdistance = 39; % inches + widthofstim = 15; % inches + wbcmd = 'wb_command'; % workbench command + tr = 1; % seconds + stimres = 200; % number of pixels along each dimension of the stimulus + maxiter = 100; % maximum number of iterations + dummycifti = dummy_file; % dummy CIFTI file + + % calc + stimdeg = atan(widthofstim/2/viewingdistance)/pi*180*2; % stimulus diameter in degrees + pxtodeg = stimdeg/stimres; % this converts pixels to degrees + + % debug speed ups + switch debugmode + case 0 + vxs = []; + seedmode = []; + numperjob = 50; + case 1 + vxs = []; + seedmode = [-2]; % super fast seed mode + numperjob = []; + case 2 + vxs = [25138 24995 24818 23601 23862 24994 23160 24987]; % eight voxels only + seedmode = [0]; % only one seed + numperjob = 2; % two voxels in each job + end + + % prepare writecifti function + gifti0 = ciftiopen(dummycifti,wbcmd); + % data is a column vector; loc is a location to write to + writecifti = @(data,loc) ciftisave(setfield(gifti0,'cdata',data),loc,wbcmd); + + %%%%% LOAD STIMULI + + % load stimuli + stimulus = {}; % each element is 200 x 200 x 300, single format + for p=1:length(movie_files) + [~,~,ext] = fileparts(movie_files{p}); + switch lower(ext) + case '.mat' + a1 = load(movie_files{p}); + stimulus{p} = a1.stim; + clear a1; + case '.hdf5' + stimulus{p} = h5read(movie_files{p},'/stim'); + case '.mov' + mobj = VideoReader(movie_files{p}); + stim0 = read(mobj); + stimulus{p} = permute(single(stim0(:,:,1,:)),[1 2 4 3]); + clear mobj stim0; + otherwise + error('Unknown movie file type: %s\n',movie_files{p}); + end + end + + % sanity check + assert(length(stimulus)==5); + if debugmode + stimulus + end + + %%%%% LOAD BEHAVIORAL FILES + + % load behavior files + behaviors = {}; + for p=1:length(behavior_files) + if isequal(behavior_files{p},'EMPTY') + behaviors{p} = []; + else + a1 = xml2struct(behavior_files{p}); % relevant fields: expttype (1-5), ttlStamps + behaviors{p} = a1.behav_summary; + end + end + clear a1; + + % sanity check + if debugmode + behaviors + end + + %%%%% LOAD DATA + + % load fMRI data + data = {}; % each element is ~90000 x 300, single format + for p=1:length(image_files) + if isequal(image_files{p},'EMPTY') + data{p} = []; + else + data{p} = single(getfield(ciftiopen(image_files{p},wbcmd),'cdata')); + end + end + + % sanity check + if debugmode + data + end + + %%%%% DEAL WITH MISSING DATA + + % if any data is EMPTY, kill it all + ok = cellfun(@(x) ~isempty(x),data); + data = data(ok); + behaviors = behaviors(ok); + + %%%%% SANITY CHECK THE BEHAVIORAL DATA AND CONSTRUCT FINAL STIMULUS + + % go through the behavioral files and check the timing. + % as we go, we construct the final stimulus specification. + % we crash if any sanity check fails (see below). + finalstimulus = {}; + for p=1:length(behaviors) + + % if we lack a behavioral file, we crash + if isempty(behaviors{p}) + error('Found an empty behavioral file (subject=%d, filenumber=%d).\n',subject,p); + else + + % pull out the correct stimulus type + fprintf('stimulus type detected: %d\n',behaviors{p}.expttype); + finalstimulus{p} = stimulus{behaviors{p}.expttype}; + + % strict check of the TTL stamps: + % (1) we must have TTL stamps + % (2) the difference between first and last must be in (298.5,299.5) + % (3) the first TTL must have occurred within 100 ms before the stimulus started + ok = ~isempty(behaviors{p}.ttlStamps) && ... + abs(diff(behaviors{p}.ttlStamps([1 end])) - 299) < .5 && ... + behaviors{p}.ttlStamps(1) > -.1 && ... + behaviors{p}.ttlStamps(1) < 0; + if ~ok + error('TTL pulses in behavioral file did not pass sanity checks (subject=%d, filenumber=%d).\n',subject,p); + end + + end + + end + + %%%%% DEAL WITH TRUNCATED DATA + + % we might be missing some data points at the end of a given run. + % in such cases, just truncate the stimulus sequence to match the data. + for p=1:length(data) + if size(data{p},2) < size(finalstimulus{p},3) + finalstimulus{p} = finalstimulus{p}(:,:,1:size(data{p},2)); + end + end + + %%%%% FINALLY, ANALYZE THE DATA + + % analyze the data + results = analyzePRF_chpc(userid,subject,finalstimulus,data,tr, ... + struct('vxs',vxs,'seedmode',seedmode, ... + 'numperjob',numperjob,'maxiter',maxiter,'display','off')); + + %%%%% MASSAGE THE OUTPUTS AND WRITE TO DISK + + % make the output directory + mkdirquiet([outpath '/results']); + + % write out results + writecifti(results.ang, sprintf('%s/results/%d_angle.dtseries.nii',outpath,subject)); + writecifti(results.ecc*pxtodeg, sprintf('%s/results/%d_eccentricity.dtseries.nii',outpath,subject)); + writecifti(results.rfsize*pxtodeg, sprintf('%s/results/%d_size.dtseries.nii',outpath,subject)); + writecifti(results.expt, sprintf('%s/results/%d_exponent.dtseries.nii',outpath,subject)); + writecifti(results.gain./results.meanvol*100, sprintf('%s/results/%d_gain.dtseries.nii',outpath,subject)); + writecifti(results.R2, sprintf('%s/results/%d_varianceexplained.dtseries.nii',outpath,subject)); + writecifti(results.meanvol, sprintf('%s/results/%d_mean.dtseries.nii',outpath,subject)); + +catch me + + fprintf('Error: %s\n',me.message); + if length(me.stack) > 0 + me.stack(1) + end + +end + +% make sure we get out +fprintf('RetinotopyAnalysis.m complete.\n'); +quit; diff --git a/MATLAB/analyzePRF_chpc.m b/MATLAB/analyzePRF_chpc.m new file mode 100644 index 0000000..0592999 --- /dev/null +++ b/MATLAB/analyzePRF_chpc.m @@ -0,0 +1,431 @@ +function results = analyzePRF_chpc(userid,subject,stimulus,data,tr,options) + +% function results = analyzePRF_chpc(userid,subject,stimulus,data,tr,options) +% +% the documentation of this function is identical to that for analyzePRF.m, +% except that we take two additional inputs: +% +% is a string with the user ID (this is used for submitting jobs) +% is a positive integer +% +% history: +% - 2015/02/07 - inherit from analyzePRF.m + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% REPORT + +fprintf('*** analyzePRF_chpc: started at %s. ***\n',datestr(now)); +stime = clock; % start time + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INTERNAL CONSTANTS + +% define +scratchdir = sprintf('/scratch/%s/analyzePRF',userid); +mkdirquiet(scratchdir); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SETUP AND PREPARATION + +% massage cell inputs +if ~iscell(stimulus) + stimulus = {stimulus}; +end +if ~iscell(data) + data = {data}; +end + +% calc +is3d = size(data{1},4) > 1; +if is3d + dimdata = 3; + dimtime = 4; + xyzsize = sizefull(data{1},3); +else + dimdata = 1; + dimtime = 2; + xyzsize = size(data{1},1); +end +numvxs = prod(xyzsize); + +% calc +res = sizefull(stimulus{1},2); +resmx = max(res); +numruns = length(data); + +% deal with inputs +if ~exist('options','var') || isempty(options) + options = struct(); +end +if ~isfield(options,'vxs') || isempty(options.vxs) + options.vxs = 1:numvxs; +end +if ~isfield(options,'wantglmdenoise') || isempty(options.wantglmdenoise) + options.wantglmdenoise = 0; +end +if ~isfield(options,'hrf') || isempty(options.hrf) + options.hrf = []; +end +if ~isfield(options,'maxpolydeg') || isempty(options.maxpolydeg) + options.maxpolydeg = []; +end +if ~isfield(options,'seedmode') || isempty(options.seedmode) + options.seedmode = [0 1 2]; +end +if ~isfield(options,'xvalmode') || isempty(options.xvalmode) + options.xvalmode = 0; +end +if ~isfield(options,'numperjob') || isempty(options.numperjob) + options.numperjob = []; +end +if ~isfield(options,'maxiter') || isempty(options.maxiter) + options.maxiter = 500; +end +if ~isfield(options,'display') || isempty(options.display) + options.display = 'iter'; +end +if ~isfield(options,'typicalgain') || isempty(options.typicalgain) + options.typicalgain = 10; +end + +% massage +wantquick = isequal(options.seedmode,-2); +options.seedmode = union(options.seedmode(:),[]); + +% massage more +if wantquick + opt.xvalmode = 0; + opt.vxs = 1:numvxs; + opt.numperjob = []; +end + +% calc +usecluster = ~isempty(options.numperjob); + +% prepare stimuli +for p=1:length(stimulus) + stimulus{p} = squish(stimulus{p},2)'; % frames x pixels + stimulus{p} = [stimulus{p} p*ones(size(stimulus{p},1),1)]; % add a dummy column to indicate run breaks + stimulus{p} = single(stimulus{p}); % make single to save memory +end + +% deal with data badness (set bad voxels to be always all 0) +bad = cellfun(@(x) any(~isfinite(x),dimtime) | all(x==0,dimtime),data,'UniformOutput',0); % if non-finite or all 0 +bad = any(cat(dimtime,bad{:}),dimtime); % badness in ANY run +for p=1:numruns + data{p}(repmat(bad,[ones(1,dimdata) size(data{p},dimtime)])) = 0; +end + +% calc mean volume +meanvol = mean(catcell(dimtime,data),dimtime); + +% what HRF should we use? +if isempty(options.hrf) + options.hrf = getcanonicalhrf(tr,tr)'; +end +numinhrf = length(options.hrf); + +% what polynomials should we use? +if isempty(options.maxpolydeg) + options.maxpolydeg = cellfun(@(x) round(size(x,dimtime)*tr/60/2),data); +end +if isscalar(options.maxpolydeg) + options.maxpolydeg = repmat(options.maxpolydeg,[1 numruns]); +end +fprintf('using the following maximum polynomial degrees: %s\n',mat2str(options.maxpolydeg)); + +% initialize cluster stuff +if usecluster + filestodelete = {}; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIGURE OUT NOISE REGRESSORS + +if isequal(options.wantglmdenoise,1) + noisereg = analyzePRFcomputeGLMdenoiseregressors(stimulus,data,tr); +elseif isequal(options.wantglmdenoise,0) + noisereg = []; +else + noisereg = options.wantglmdenoise; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE MODEL + +% pre-compute some cache +[d,xx,yy] = makegaussian2d(resmx,2,2,2,2); + +% define the model (parameters are R C S G N) +modelfun = @(pp,dd) conv2run(posrect(pp(4)) * (dd*[vflatten(placematrix(zeros(res),makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0) / (2*pi*abs(pp(3))^2))); 0]) .^ posrect(pp(5)),options.hrf,dd(:,prod(res)+1)); +model = {{[] [1-res(1)+1 1-res(2)+1 0 0 NaN; + 2*res(1)-1 2*res(2)-1 Inf Inf Inf] modelfun} ... + {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 0; + 2*res(1)-1 2*res(2)-1 Inf Inf Inf] @(ss)modelfun}}; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE SEEDS + +% init +seeds = []; + +% generic large seed +if ismember(0,options.seedmode) + seeds = [seeds; + (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5) options.typicalgain 0.5]; +end + +% generic small seed +if ismember(1,options.seedmode) + seeds = [seeds; + (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5)/10 options.typicalgain 0.5]; +end + +% super-grid seed +if any(ismember([2 -2],options.seedmode)) + [supergridseeds,rvalues] = analyzePRFcomputesupergridseeds(res,stimulus,data,modelfun, ... + options.maxpolydeg,dimdata,dimtime, ... + options.typicalgain,noisereg); +end + +% make a function that individualizes the seeds +if exist('supergridseeds','var') + seedfun = @(vx) [[seeds]; + [subscript(squish(supergridseeds,dimdata),{vx ':'})]]; +else + seedfun = @(vx) [seeds]; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PERFORM OPTIMIZATION + +% if this is true, we can bypass all of the optimization stuff! +if wantquick + +else + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE RESAMPLING STUFF + + % define wantresampleruns and resampling + switch options.xvalmode + case 0 + wantresampleruns = []; + resampling = 0; + case 1 + wantresampleruns = 1; + half1 = copymatrix(zeros(1,length(data)),1:round(length(data)/2),1); + half2 = ~half1; + resampling = [(1)*half1 + (-1)*half2; + (-1)*half1 + (1)*half2]; + case 2 + wantresampleruns = 0; + resampling = []; + for p=1:length(data) + half1 = copymatrix(zeros(1,size(data{p},2)),1:round(size(data{p},2)/2),1); + half2 = ~half1; + resampling = cat(2,resampling,[(1)*half1 + (-1)*half2; + (-1)*half1 + (1)*half2]); + end + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE STIMULUS AND DATA + + %%%%% CLUSTER CASE + + if usecluster + + % save stimulus + while 1 + localfile0 = [scratchdir '/' sprintf('stim%d_%s.mat',subject,randomword(5))]; + if ~exist(localfile0,'file') + break; + end + end + save(localfile0,'stimulus'); + filestodelete{end+1} = localfile0; + clear stimulus; + + % define stimulus + stimulus = @() loadmulti(localfile0,'stimulus'); + + % save data + while 1 + % directory name that will contain 001.bin, etc. + localfile0 = [scratchdir '/' sprintf('data%d_%s',subject,randomword(5))]; + if ~exist(localfile0,'dir') + break; + end + end + assert(mkdir(localfile0)); + for p=1:numruns + savebinary([localfile0 sprintf('/%03d.bin',p)],'single',squish(data{p},dimdata)'); % notice squish + end + filestodelete{end+1} = localfile0; + clear data; + + % define data + binfiles = cellfun(@(x) [localfile0 sprintf('/%03d.bin',x)],num2cell(1:numruns),'UniformOutput',0); + data = @(vxs) cellfun(@(x) double(loadbinary(x,'single',[0 numvxs],-vxs)),binfiles,'UniformOutput',0); + + % prepare the output directory + while 1 + localfile0 = [scratchdir '/' sprintf('prfresults%d_%s',subject,randomword(5))]; + if ~exist(localfile0,'dir') + break; + end + end + filestodelete{end+1} = localfile0; + filestodelete{end+1} = [localfile0 '.mat']; % after consolidation + outputdir = localfile0; + + %%%%% NON-CLUSTER CASE + + else + + stimulus = {stimulus}; + data = @(vxs) cellfun(@(x) subscript(squish(x,dimdata),{vxs ':'})',data,'UniformOutput',0); + outputdir = []; + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE OPTIONS + + % last-minute prep + if iscell(noisereg) + noiseregINPUT = {noisereg}; + else + noiseregINPUT = noisereg; + end + + % construct the options struct + opt = struct( ... + 'outputdir',outputdir, ... + 'stimulus',stimulus, ... + 'data',data, ... + 'vxs',options.vxs, ... + 'model',{model}, ... + 'seed',seedfun, ... + 'optimoptions',{{'Display' options.display 'Algorithm' 'levenberg-marquardt' 'MaxIter' options.maxiter}}, ... + 'wantresampleruns',wantresampleruns, ... + 'resampling',resampling, ... + 'metric',@calccod, ... + 'maxpolydeg',options.maxpolydeg, ... + 'wantremovepoly',1, ... + 'extraregressors',noiseregINPUT, ... + 'wantremoveextra',0, ... + 'dontsave',{{'modelfit' 'opt' 'vxsfull' 'modelpred' 'testdata'}}); % 'resnorms' 'numiters' + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIT MODEL + + %%%%% CLUSTER CASE + + if usecluster + + % define job name + jobname = sprintf('%d_%s',subject,makedirid(opt.outputdir,1)); + + % submit jobs + jobids = CHPCsubmit(userid,jobname,'run_fitnonlinearmodel.sh',options.numperjob, ... + 1,ceil(length(options.vxs)/options.numperjob),[], ... + {'data' 'stimulus' 'bad' 'd' 'xx' 'yy' 'modelfun' 'model'}); + + % record additional files to delete + filestodelete{end+1} = sprintf('~/sgeoutput/job_%s.*',jobname); % .o and .e files + filestodelete{end+1} = sprintf('~/mcc/job_%s.mat',jobname); + + % wait for jobs to finish + CHPCwait(jobname,jobids,userid); + + % consolidate and load the results + fitnonlinearmodel_consolidate(outputdir); + a1 = load([outputdir '.mat']); + + %%%%% NON-CLUSTER CASE + + else + + a1 = fitnonlinearmodel(opt); + + end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE OUTPUT + +% depending on which analysis we did (quick or full optimization), +% we have to get the outputs in a common format +if wantquick + paramsA = permute(squish(supergridseeds,dimdata),[3 2 1]); % fits x parameters x voxels + rA = squish(rvalues,dimdata)'; % fits x voxels +else + paramsA = a1.params; % fits x parameters x voxels + rA = a1.trainperformance; % fits x voxels +end + +% calc +numfits = size(paramsA,1); + +% init +clear results; +results.ang = NaN*zeros(numvxs,numfits); +results.ecc = NaN*zeros(numvxs,numfits); +results.expt = NaN*zeros(numvxs,numfits); +results.rfsize = NaN*zeros(numvxs,numfits); +results.R2 = NaN*zeros(numvxs,numfits); +results.gain = NaN*zeros(numvxs,numfits); +results.resnorms = cell(numvxs,1); +results.numiters = cell(numvxs,1); + +% massage model parameters for output and put in 'results' struct +results.ang(options.vxs,:) = permute(mod(atan2((1+res(1))/2 - paramsA(:,1,:), ... + paramsA(:,2,:) - (1+res(2))/2),2*pi)/pi*180,[3 1 2]); +results.ecc(options.vxs,:) = permute(sqrt(((1+res(1))/2 - paramsA(:,1,:)).^2 + ... + (paramsA(:,2,:) - (1+res(2))/2).^2),[3 1 2]); +results.expt(options.vxs,:) = permute(posrect(paramsA(:,5,:)),[3 1 2]); +results.rfsize(options.vxs,:) = permute(abs(paramsA(:,3,:)) ./ sqrt(posrect(paramsA(:,5,:))),[3 1 2]); +results.R2(options.vxs,:) = permute(rA,[2 1]); +results.gain(options.vxs,:) = permute(posrect(paramsA(:,4,:)),[3 1 2]); +if ~wantquick + results.resnorms(options.vxs) = a1.resnorms; + results.numiters(options.vxs) = a1.numiters; +end + +% reshape +results.ang = reshape(results.ang, [xyzsize numfits]); +results.ecc = reshape(results.ecc, [xyzsize numfits]); +results.expt = reshape(results.expt, [xyzsize numfits]); +results.rfsize = reshape(results.rfsize, [xyzsize numfits]); +results.R2 = reshape(results.R2, [xyzsize numfits]); +results.gain = reshape(results.gain, [xyzsize numfits]); +results.resnorms = reshape(results.resnorms, [xyzsize 1]); +results.numiters = reshape(results.numiters, [xyzsize 1]); + +% add some more stuff +results.meanvol = meanvol; +results.noisereg = noisereg; +results.params = paramsA; +results.options = options; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CLEAN UP + +% no clean up necessary in the quick case +if ~wantquick + + %%%%% CLUSTER CASE + + if usecluster + + % delete local files and directories [should make this a function!] + for p=1:length(filestodelete) + if exist(filestodelete{p},'dir') % first dir, then file + rmdir(filestodelete{p},'s'); + elseif exist(filestodelete{p},'file') + delete(filestodelete{p}); + end + end + + %%%%% NON-CLUSTER CASE + + else + + end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% REPORT + +fprintf('*** analyzePRF_chpc: ended at %s (%.1f minutes). ***\n', ... + datestr(now),etime(clock,stime)/60); diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/Contents.m b/MATLAB/utilities/MatlabCIFTI/@gifti/Contents.m new file mode 100755 index 0000000..eaac456 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/Contents.m @@ -0,0 +1,39 @@ +% GIfTI Class for MATLAB +% +% Geometry format under the Neuroimaging Informatics Technology Initiative +% (NIfTI): +% http://www.nitrc.org/projects/gifti/ +% http://nifti.nimh.nih.gov/ +% +% This MATLAB class is part of SPM: +% http://www.fil.ion.ucl.ac.uk/spm/ +% +% It relies on external MATLAB code: +% Base64, by Peter J. Acklam: +% http://home.online.no/~pjacklam/ +% dzip, by Michael Kleder: +% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8899 +% XMLTree, by Guillaume Flandin: +% http://www.artefact.tk/software/matlab/xml/ +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: Contents.m 2076 2008-09-10 12:34:08Z guillaume $ + +% GIfTI file format for MATLAB (The Mathworks, Inc.). +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging +% +% This program 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 2 +% of the License, or any later version. +% +% This program 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 this program; if not, write to the Free Software +% Foundation Inc, 59 Temple Pl. - Suite 330, Boston, MA 02111-1307, USA. \ No newline at end of file diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/display.m b/MATLAB/utilities/MatlabCIFTI/@gifti/display.m new file mode 100755 index 0000000..8ab8c9d --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/display.m @@ -0,0 +1,25 @@ +function display(this) +% Display method for GIfTI objects +% FORMAT display(this) +% this - GIfTI object +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: display.m 4182 2011-02-01 12:29:09Z guillaume $ + +display_name = inputname(1); +if isempty(display_name) + display_name = 'ans'; +end + +if length(this) == 1 && ~isempty(this.data) + eval([display_name ' = struct(this);']); + eval(['display(' display_name ');']); +else + disp(' ') + disp([display_name ' =']); + disp(' '); + eval([display_name ' = this;']); + eval(['disp(' display_name ');']); +end \ No newline at end of file diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/export.m b/MATLAB/utilities/MatlabCIFTI/@gifti/export.m new file mode 100755 index 0000000..0e662b7 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/export.m @@ -0,0 +1,53 @@ +function s = export(this,target) +% Export a GIfTI object into specific MATLAB struct +% FORMAT s = export(this,target) +% this - GIfTI object +% target - string describing target output [default: Matlab] +% s - a structure containing public fields of the object +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: export.m 3578 2009-11-18 16:19:45Z guillaume $ + +if numel(this) > 1, warning('Only handle scalar objects yet.'); end + +if nargin <= 1, target = 'Matlab'; end + +switch lower(target) + case 'matlab' + s = struct(this); + + case 'patch' + if isfield(this,'vertices') + s.vertices = double(subsref(this, substruct('.', 'vertices'))); + end + if isfield(this,'faces') + s.faces = subsref(this, substruct('.', 'faces')); + end + if isfield(this,'cdata') + s.facevertexcdata = double(subsref(this, substruct('.', 'cdata'))); + end + try, s; catch, s = struct([]); end + + case {'fieldtrip', 'ft'} + s = struct('tri',[], 'pnt',[]); + if isfield(this,'vertices') + s.pnt = double(subsref(this, substruct('.', 'vertices'))); + end + if isfield(this,'faces') + s.tri = double(subsref(this, substruct('.', 'faces'))); + end + + case {'spm'} + s = struct('face',[], 'vert',[]); + if isfield(this,'vertices') + s.vert = double(subsref(this, substruct('.', 'vertices'))); + end + if isfield(this,'faces') + s.face = uint16(subsref(this, substruct('.', 'faces'))); + end + + otherwise + error('Unknown target ''%s''.', target); +end diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/fieldnames.m b/MATLAB/utilities/MatlabCIFTI/@gifti/fieldnames.m new file mode 100755 index 0000000..0813696 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/fieldnames.m @@ -0,0 +1,16 @@ +function names = fieldnames(this) +% Fieldnames method for GIfTI objects +% FORMAT names = fieldnames(this) +% this - GIfTI object +% names - field names +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: fieldnames.m 2076 2008-09-10 12:34:08Z guillaume $ + +if numel(this) > 1, warning('Only handle scalar objects yet.'); end + +pfn = {'vertices' 'faces' 'normals' 'cdata' 'mat'}; + +names = pfn(isintent(this,pfn)); \ No newline at end of file diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/gifti.m b/MATLAB/utilities/MatlabCIFTI/@gifti/gifti.m new file mode 100755 index 0000000..3e8dca6 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/gifti.m @@ -0,0 +1,108 @@ +function this = gifti(varargin) +% GIfTI Geometry file format class +% Geometry format under the Neuroimaging Informatics Technology Initiative +% (NIfTI): +% http://www.nitrc.org/projects/gifti/ +% http://nifti.nimh.nih.gov/ +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: gifti.m 3999 2010-07-19 10:54:18Z guillaume $ + +switch nargin + + case 0 + this = giftistruct; + this = class(this,'gifti'); + + case 1 + if isa(varargin{1},'gifti') + this = varargin{1}; + + elseif isstruct(varargin{1}) + f = {'faces', 'face', 'tri' 'vertices', 'vert', 'pnt', 'cdata'}; + ff = {'faces', 'faces', 'faces', 'vertices', 'vertices', 'vertices', 'cdata'}; + [c, ia] = intersect(f,fieldnames(varargin{1})); + if ~isempty(c) + this = gifti; + for i=1:length(c) + this = subsasgn(this,... + struct('type','.','subs',ff{ia(i)}),... + varargin{1}.(c{i})); + end + elseif isempty(setxor(fieldnames(varargin{1}),... + {'metadata','label','data'})) + this = class(varargin{1},'gifti'); + else + error('[GIFTI] Invalid structure.'); + end + + elseif ishandle(varargin{1}) + this = struct('vertices',get(varargin{1},'Vertices'), ... + 'faces', get(varargin{1},'Faces')); + if ~isempty(get(varargin{1},'FaceVertexCData')); + this.cdata = get(varargin{1},'FaceVertexCData'); + end + this = gifti(this); + + elseif isnumeric(varargin{1}) + this = gifti; + this = subsasgn(this,... + struct('type','.','subs','cdata'),... + varargin{1}); + + elseif ischar(varargin{1}) + if size(varargin{1},1)>1 + this = gifti(cellstr(varargin{1})); + return; + end + [p,n,e] = fileparts(varargin{1}); + if strcmpi(e,'.mat') + try + this = gifti(load(varargin{1})); + catch + error('[GIFTI] Loading of file %s failed.', varargin{1}); + end + else + this = read_gifti_file_standalone(varargin{1},giftistruct); + this = class(this,'gifti'); + end + + elseif iscellstr(varargin{1}) + fnames = varargin{1}; + this(numel(fnames)) = giftistruct; + this = class(this,'gifti'); + for i=1:numel(fnames) + this(i) = gifti(fnames{i}); + end + + else + error('[GIFTI] Invalid object construction.'); + end + + otherwise + error('[GIFTI] Invalid object construction.'); +end + +%========================================================================== +function s = giftistruct +s = struct(... + 'metadata', ... + struct(... + 'name', {}, ... + 'value', {} ... + ), ... + 'label', ... + struct(... + 'name', {}, ... + 'index', {} ... + ), ... + 'data', ... + struct(... + 'attributes', {}, ... + 'metadata', struct('name',{}, 'value',{}), ... + 'space', {}, ... + 'data', {} ... + ) ... + ); diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/isfield.m b/MATLAB/utilities/MatlabCIFTI/@gifti/isfield.m new file mode 100755 index 0000000..f257a72 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/isfield.m @@ -0,0 +1,13 @@ +function tf = isfield(this,field) +% Isfield method for GIfTI objects +% FORMAT tf = isfield(this,field) +% this - GIfTI object +% field - string of cell array +% tf - logical array +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: isfield.m 2076 2008-09-10 12:34:08Z guillaume $ + +tf = ismember(field, fieldnames(this)); \ No newline at end of file diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/plot.m b/MATLAB/utilities/MatlabCIFTI/@gifti/plot.m new file mode 100755 index 0000000..1f96725 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/plot.m @@ -0,0 +1,61 @@ +function varargout = plot(varargin) +% plot method for GIfTI objects +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: plot.m 3556 2009-11-11 18:20:34Z guillaume $ + +% if ishandle(varargin{1}) +% h = figure(varargin{1}); +% else +% h = figure; +% %axis equal; +% %axis off; +% %camlight; +% %camlight(-80,-10); +% %lighting phong; +% end +% cameramenu; + + +cdata = []; +if nargin == 1 + this = varargin{1}; + h = gcf; +else + if ishandle(varargin{1}) + h = figure(get(varargin{1},'parent')); + this = varargin{2}; + else + this = varargin{1}; + h = gcf; + cdata = subsref(varargin{2},struct('type','.','subs','cdata')); + end + if nargin > 2 + indc = varargin{3}; + else + indc = 1; + end +end + +hp = patch(struct(... + 'vertices', subsref(this,struct('type','.','subs','vertices')),... + 'faces', subsref(this,struct('type','.','subs','faces'))),... + 'FaceColor', 'b',... + 'EdgeColor', 'none'); + +if ~isempty(cdata) + set(hp,'FaceVertexCData',cdata(:,indc), 'FaceColor','interp') +end + +axis equal; +axis off; +camlight; +camlight(-80,-10); +lighting phong; +cameramenu; + +if nargout + varargout{1} = hp; +end \ No newline at end of file diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/TEST.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/TEST.m new file mode 100755 index 0000000..d89b14c --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/TEST.m @@ -0,0 +1,211 @@ +function this = read_gifti_file(filename, this) +% Low level reader of GIfTI 1.0 files +% FORMAT this = read_gifti_file(filename, this) +% filename - XML GIfTI filename +% this - structure with fields 'metaData', 'label' and 'data'. +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: read_gifti_file.m 8 2008-05-12 09:53:02Z guillaume $ + +% Import XML-based GIfTI file +%-------------------------------------------------------------------------- +try + fid = fopen(filename,'rt'); + xmlstr = fread(fid,'*char')'; + fclose(fid); + t = xml_parser(xmlstr); +catch + error('[GIFTI] Loading of XML file %s failed.', filename); +end + +% Root element of a GIFTI file +%-------------------------------------------------------------------------- +if ~strcmp(xml_get(t,xml_root(t),'name'),'GIFTI') + error('[GIFTI] %s is not a GIFTI 1.0 file.', filename); +end +attr = cell2mat(xml_attributes(t,'get',xml_root(t))); +attr = cell2struct({attr.val},strrep({attr.key},':','___'),2); +if ~all(ismember({'Version','NumberOfDataArrays'},fieldnames(attr))) + error('[GIFTI] Missing mandatory attributes for GIFTI root element.'); +end +if str2double(attr.Version) ~= 1 + warning('[GIFTI] Unknown specification version of GIFTI file (%s).',attr.Version); +end +nbData = str2double(attr.NumberOfDataArrays); + +% Read children elements +%-------------------------------------------------------------------------- +uid = xml_children(t,xml_root(t)); +for i=1:length(uid) + switch xml_get(t,uid(i),'name') + case 'MetaData' + this.metadata = gifti_MetaData(t,uid(i)); + case 'LabelTable' + this.label = gifti_LabelTable(t,uid(i)); + case 'DataArray' + this.data{end+1} = gifti_DataArray(t,uid(i),filename); + otherwise + warning('[GIFTI] Unknown element "%s": ignored.',xml_get(t,uid(i),'name')); + end +end + +if nbData ~= length(this.data) + warning('[GIFTI] Mismatch between expected and effective number of datasets.'); +end + +%========================================================================== +function s = gifti_MetaData(t,uid) +s = struct('name',{}, 'value',{}); +c = xml_children(t,uid); +for i=1:length(c) + for j=xml_children(t,c(i)) + s(i).(lower(xml_get(t,j,'name'))) = xml_get(t,xml_children(t,j),'value'); + end +end + +%========================================================================== +function s = gifti_LabelTable(t,uid) +s = struct('name',{}, 'index',[]); +c = xml_children(t,uid); +for i=1:length(c) + a = xml_attributes(t,'get',c(i)); + s(1).index(i) = str2double(a.val); + s(1).name{i} = xml_get(t,xml_children(t,c(i)),'value'); +end + +%========================================================================== +function s = gifti_DataArray(t,uid,filename) +s = struct(... + 'attributes', {}, ... + 'data', {}, ... + 'metadata', struct([]), ... + 'space', {} ... + ); + +attr = cell2mat(xml_attributes(t,'get',uid)); +s(1).attributes = cell2struct({attr.val},{attr.key},2); +s(1).attributes.Dim = []; +for i=1:str2double(s(1).attributes.Dimensionality) + f = sprintf('Dim%d',i-1); + s(1).attributes.Dim(i) = str2double(s(1).attributes.(f)); + s(1).attributes = rmfield(s(1).attributes,f); +end +s(1).attributes = rmfield(s(1).attributes,'Dimensionality'); +try + s(1).attributes.ExternalFileName = fullfile(fileparts(filename),... + s(1).attributes.ExternalFileName); +catch +end + +c = xml_children(t,uid); +for i=1:length(c) + switch xml_get(t,c(i),'name') + case 'MetaData' + s(1).metadata = gifti_MetaData(t,c(i)); + case 'CoordinateSystemTransformMatrix' + s(1).space(end+1) = gifti_Space(t,c(i)); + case 'Data' + s(1).data = gifti_Data(t,c(i),s(1).attributes); + otherwise + error('[GIFTI] Unknown DataArray element "%s".',xml_get(t,c(i),'name')); + end +end + +%========================================================================== +function s = gifti_Space(t,uid) +s = struct('DataSpace','', 'TransformedSpace','', 'MatrixData',[]); +for i=xml_children(t,uid) + s.(xml_get(t,i,'name')) = xml_get(t,xml_children(t,i),'value'); +end +s.MatrixData = reshape(str2num(s.MatrixData),4,4)'; + +%========================================================================== +function d = gifti_Data(t,uid,s) +tp = getdict; +try + tp = tp.(s.DataType); +catch + error('[GIFTI] Unknown DataType.'); +end + +[unused,unused,mach] = fopen(1); +sb = @deal; %inline('x'); +try + if (strcmp(s.Endian,'LittleEndian') && ~isempty(strmatch('ieee-be',mach))) ... + || (strcmp(s.Endian,'BigEndian') && ~isempty(strmatch('ieee-le',mach))) + sb = @swapbyte; + end +catch + % Byte Order can be absent if encoding is ASCII, assume native otherwise +end + +switch s.Encoding + case 'ASCII' + d = sscanf(xml_get(t,xml_children(t,uid),'value'),tp.format); + + case 'Base64Binary' + d = typecast(sb(base64decode(xml_get(t,xml_children(t,uid),'value'))), tp.cast); + + case 'GZipBase64Binary' + d = typecast(dunzip(sb(base64decode(xml_get(t,xml_children(t,uid),'value')))), tp.cast); + + case 'ExternalFileBinary' + fid = fopen(s.ExternalFileName,'r'); + if fid == -1 + error('[GIFTI] Unable to read binary file %s.',s.ExternalFileName); + end + fseek(fid,str2double(s.ExternalFileOffset),0); + d = sb(fread(fid,prod(s.Dim),['*' tp.class])); + fclose(fid); + + otherwise + error('[GIFTI] Unknown data encoding: %s.',s.Encoding); +end + +if length(s.Dim) == 1, s.Dim(end+1) = 1; end +switch s.ArrayIndexingOrder + case 'RowMajorOrder' + d = permute(reshape(d,fliplr(s.Dim)),length(s.Dim):-1:1); + case 'ColumnMajorOrder' + d = reshape(d,s.Dim); + otherwise + error('[GIFTI] Unknown array indexing order.'); +end + +%========================================================================== +%========================================================================== +function uid = xml_root(tree) +uid = 1; +for i=1:length(tree) + if strcmp(xml_get(tree,i,'type'),'element') + uid = i; + break + end +end + +%-------------------------------------------------------------------------- +function child = xml_children(tree,uid) +if strcmp(tree{uid}.type,'element') + child = tree{uid}.contents; +else + child = []; +end + +%-------------------------------------------------------------------------- +function value = xml_get(tree,uid,parameter) +try + value = tree{uid}.(parameter); +catch + error(sprintf('[XML] Parameter %s not found.',parameter)); +end + +%-------------------------------------------------------------------------- +function varargout = xml_attributes(tree,method,uid) +if ~strcmpi(method,'get'), error('[XML] Unknown attributes method.'); end +if isempty(tree{uid}.attributes) + varargout{1} = {}; +else + varargout{1} = tree{uid}.attributes; +end diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/base64decode.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/base64decode.m new file mode 100755 index 0000000..ca45732 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/base64decode.m @@ -0,0 +1,81 @@ +function y = base64decode(x) +%BASE64DECODE Perform base64 decoding on a string. +% +% BASE64DECODE(STR) decodes the given base64 string STR. +% +% Any character not part of the 65-character base64 subset set is silently +% ignored. +% +% This function is used to decode strings from the Base64 encoding specified +% in RFC 2045 - MIME (Multipurpose Internet Mail Extensions). The Base64 +% encoding is designed to represent arbitrary sequences of octets in a form +% that need not be humanly readable. A 65-character subset ([A-Za-z0-9+/=]) +% of US-ASCII is used, enabling 6 bits to be represented per printable +% character. +% +% See also BASE64ENCODE. + +% Author: Peter J. Acklam +% Time-stamp: 2004-09-20 08:20:50 +0200 +% E-mail: pjacklam@online.no +% URL: http://home.online.no/~pjacklam + +% Modified by Guillaume Flandin, May 2008 + + +% Perform the following mapping +%-------------------------------------------------------------------------- +% A-Z -> 0 - 25 a-z -> 26 - 51 0-9 -> 52 - 61 +% + -> 62 / -> 63 = -> 64 +% anything else -> NaN + +base64chars = NaN(1,256); +base64chars('ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=') = 0:64; +x = base64chars(x); + +% Remove/ignore any characters not in the base64 characters list or '=' +%-------------------------------------------------------------------------- + +x = x(~isnan(x)); + +% Replace any incoming padding ('=' -> 64) with a zero pad +%-------------------------------------------------------------------------- + +if x(end-1) == 64, p = 2; x(end-1:end) = 0; +elseif x(end) == 64, p = 1; x(end) = 0; +else p = 0; +end + +% Allocate decoded data array +%-------------------------------------------------------------------------- + +n = length(x) / 4; % number of groups +x = reshape(uint8(x), 4, n); % input data +y = zeros(3, n, 'uint8'); % decoded data + +% Rearrange every 4 bytes into 3 bytes +%-------------------------------------------------------------------------- +% 00aaaaaa 00bbbbbb 00cccccc 00dddddd +% +% to form +% +% aaaaaabb bbbbcccc ccdddddd + +y(1,:) = bitshift(x(1,:), 2); % 6 highest bits of y(1,:) +y(1,:) = bitor(y(1,:), bitshift(x(2,:), -4)); % 2 lowest bits of y(1,:) + +y(2,:) = bitshift(x(2,:), 4); % 4 highest bits of y(2,:) +y(2,:) = bitor(y(2,:), bitshift(x(3,:), -2)); % 4 lowest bits of y(2,:) + +y(3,:) = bitshift(x(3,:), 6); % 2 highest bits of y(3,:) +y(3,:) = bitor(y(3,:), x(4,:)); % 6 lowest bits of y(3,:) + +% Remove any zero pad that was added to make this a multiple of 24 bits +%-------------------------------------------------------------------------- + +if p, y(end-p+1:end) = []; end + +% Reshape to a row vector +%-------------------------------------------------------------------------- + +y = reshape(y, 1, []); diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/base64encode.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/base64encode.m new file mode 100755 index 0000000..1256731 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/base64encode.m @@ -0,0 +1,157 @@ +function y = base64encode(x, eol) +%BASE64ENCODE Perform base64 encoding on a string. +% +% BASE64ENCODE(STR, EOL) encode the given string STR. EOL is the line ending +% sequence to use; it is optional and defaults to '\n' (ASCII decimal 10). +% The returned encoded string is broken into lines of no more than 76 +% characters each, and each line will end with EOL unless it is empty. Let +% EOL be empty if you do not want the encoded string broken into lines. +% +% STR and EOL don't have to be strings (i.e., char arrays). The only +% requirement is that they are vectors containing values in the range 0-255. +% +% This function may be used to encode strings into the Base64 encoding +% specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions). The +% Base64 encoding is designed to represent arbitrary sequences of octets in a +% form that need not be humanly readable. A 65-character subset +% ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per +% printable character. +% +% Examples +% -------- +% +% If you want to encode a large file, you should encode it in chunks that are +% a multiple of 57 bytes. This ensures that the base64 lines line up and +% that you do not end up with padding in the middle. 57 bytes of data fills +% one complete base64 line (76 == 57*4/3): +% +% If ifid and ofid are two file identifiers opened for reading and writing, +% respectively, then you can base64 encode the data with +% +% while ~feof(ifid) +% fwrite(ofid, base64encode(fread(ifid, 60*57))); +% end +% +% or, if you have enough memory, +% +% fwrite(ofid, base64encode(fread(ifid))); +% +% See also BASE64DECODE. + +% Author: Peter J. Acklam +% Time-stamp: 2004-02-03 21:36:56 +0100 +% E-mail: pjacklam@online.no +% URL: http://home.online.no/~pjacklam + + + % make sure we have the EOL value + if nargin < 2 + eol = ''; %sprintf('\n'); + else + if sum(size(eol) > 1) > 1 + error('EOL must be a vector.'); + end + if any(eol(:) > 255) + error('EOL can not contain values larger than 255.'); + end + end + + if sum(size(x) > 1) > 1 + error('STR must be a vector.'); + end + + x = uint8(x); + eol = uint8(eol); + + ndbytes = length(x); % number of decoded bytes + nchunks = ceil(ndbytes / 3); % number of chunks/groups + nebytes = 4 * nchunks; % number of encoded bytes + + % add padding if necessary, to make the length of x a multiple of 3 + if rem(ndbytes, 3) + x(end+1 : 3*nchunks) = 0; + end + + x = reshape(x, [3, nchunks]); % reshape the data + y = repmat(uint8(0), 4, nchunks); % for the encoded data + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Split up every 3 bytes into 4 pieces + % + % aaaaaabb bbbbcccc ccdddddd + % + % to form + % + % 00aaaaaa 00bbbbbb 00cccccc 00dddddd + % + y(1,:) = bitshift(x(1,:), -2); % 6 highest bits of x(1,:) + + y(2,:) = bitshift(bitand(x(1,:), 3), 4); % 2 lowest bits of x(1,:) + y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4)); % 4 highest bits of x(2,:) + + y(3,:) = bitshift(bitand(x(2,:), 15), 2); % 4 lowest bits of x(2,:) + y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6)); % 2 highest bits of x(3,:) + + y(4,:) = bitand(x(3,:), 63); % 6 lowest bits of x(3,:) + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Now perform the following mapping + % + % 0 - 25 -> A-Z + % 26 - 51 -> a-z + % 52 - 61 -> 0-9 + % 62 -> + + % 63 -> / + % + % We could use a mapping vector like + % + % ['A':'Z', 'a':'z', '0':'9', '+/'] + % + % but that would require an index vector of class double. + % + z = repmat(uint8(0), size(y)); + i = y <= 25; z(i) = 'A' + double(y(i)); + i = 26 <= y & y <= 51; z(i) = 'a' - 26 + double(y(i)); + i = 52 <= y & y <= 61; z(i) = '0' - 52 + double(y(i)); + i = y == 62; z(i) = '+'; + i = y == 63; z(i) = '/'; + y = z; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Add padding if necessary. + % + npbytes = 3 * nchunks - ndbytes; % number of padding bytes + if npbytes + y(end-npbytes+1 : end) = '='; % '=' is used for padding + end + + if isempty(eol) + + % reshape to a row vector + y = reshape(y, [1, nebytes]); + + else + + nlines = ceil(nebytes / 76); % number of lines + neolbytes = length(eol); % number of bytes in eol string + + % pad data so it becomes a multiple of 76 elements + y(nebytes + 1 : 76 * nlines) = 0; + y = reshape(y, 76, nlines); + + % insert eol strings + eol = eol(:); + y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines)); + + % remove padding, but keep the last eol string + m = nebytes + neolbytes * (nlines - 1); + n = (76+neolbytes)*nlines - neolbytes; + y(m+1 : n) = ''; + + % extract and reshape to row vector + y = reshape(y, 1, m+neolbytes); + + end + + % output is a character array + y = char(y); diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/dunzip.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/dunzip.m new file mode 100755 index 0000000..c6ea022 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/dunzip.m @@ -0,0 +1,19 @@ +function M = dunzip(Z) +% DUNZIP - decompress gzipped stream of bytes +% FORMAT M = dzip(Z) +% Z - compressed variable to decompress (uint8 vector) +% M - decompressed output +% +% See also DZIP + +% Carefully tested, but no warranty; use at your own risk. +% Michael Kleder, Nov 2005 +% Modified by Guillaume Flandin, May 2008 + +import com.mathworks.mlwidgets.io.InterruptibleStreamCopier +a = java.io.ByteArrayInputStream(Z); +b = java.util.zip.InflaterInputStream(a); +isc = InterruptibleStreamCopier.getInterruptibleStreamCopier; +c = java.io.ByteArrayOutputStream; +isc.copyStream(b,c); +M = c.toByteArray; diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/dzip.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/dzip.m new file mode 100755 index 0000000..d3905d9 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/dzip.m @@ -0,0 +1,20 @@ +function Z = dzip(M) +% DZIP - losslessly compress data using gzip +% FORMAT Z = dzip(M) +% M - byte stream to compress +% Z - compressed output +% +% See also DUNZIP + +% This function uses the public domain ZLIB Deflater algorithm. +% Carefully tested, but no warranty; use at your own risk. +% Michael Kleder, Nov 2005 +% Modified by Guillaume Flandin, May 2008 + +M = typecast(M(:),'uint8'); +f = java.io.ByteArrayOutputStream(); +g = java.util.zip.DeflaterOutputStream(f); +g.write(M); +g.close; +Z = typecast(f.toByteArray,'uint8'); +f.close; diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/getdict.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/getdict.m new file mode 100755 index 0000000..85986f3 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/getdict.m @@ -0,0 +1,26 @@ +function d = getdict +% Dictionary of GIfTI/NIfTI stuff +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: getdict.m 4505 2011-09-30 11:45:58Z guillaume $ + +persistent dict; +if ~isempty(dict) + d = dict; + return; +end + +table = {... + 'NIFTI_TYPE_UINT8', 'uint8', '%d', @uint8, 'uint8' + 'NIFTI_TYPE_INT32', 'int32', '%d', @int32, 'int32' + 'NIFTI_TYPE_FLOAT32', 'float32', '%f', @single, 'single' + 'NIFTI_TYPE_FLOAT64', 'float64', '%f', @double, 'double'}; + +for i=1:size(table,1) + dict.(table{i,1}) = cell2struct({table{i,2:end}},... + {'class', 'format', 'conv', 'cast'}, 2); +end + +d = dict; \ No newline at end of file diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/isintent.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/isintent.m new file mode 100755 index 0000000..7468e2a --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/isintent.m @@ -0,0 +1,104 @@ +function [a, b] = isintent(this,intent) +% Correspondance between fieldnames and NIfTI intent codes +% FORMAT ind = isintent(this,intent) +% this - GIfTI object +% intent - fieldnames +% a - indices of found intent(s) +% b - indices of dataarrays of found intent(s) +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: isintent.m 4702 2012-03-27 16:40:11Z guillaume $ + +a = []; +b = []; +if ischar(intent), intent = cellstr(intent); end +for i=1:length(this(1).data) + switch this(1).data{i}.attributes.Intent(14:end) + case 'POINTSET' + [tf, loc] = ismember('vertices',intent); + if tf + a(end+1) = loc; + b(end+1) = i; + end + [tf, loc] = ismember('mat',intent); + if tf + a(end+1) = loc; + b(end+1) = i; + end + case 'TRIANGLE' + [tf, loc] = ismember('faces',intent); + if tf + a(end+1) = loc; + b(end+1) = i; + end + case 'VECTOR' + [tf, loc] = ismember('normals',intent); + if tf + a(end+1) = loc; + b(end+1) = i; + end + case cdata + [tf, loc] = ismember('cdata',intent); + if tf + a(end+1) = loc; + b(end+1) = i; + end + otherwise + fprintf('Intent %s is ignored.\n',this.data{i}.attributes.Intent); + end +end +%[d,i] = unique(a); +%if length(d) < length(a) +% warning('Several fields match intent type. Using first.'); +% a = a(i); +% b = b(i); +%end + +function c = cdata + +c = { +'NONE' +'CORREL' +'TTEST' +'FTEST' +'ZSCORE' +'CHISQ' +'BETA' +'BINOM' +'GAMMA' +'POISSON' +'NORMAL' +'FTEST_NONC' +'CHISQ_NONC' +'LOGISTIC' +'LAPLACE' +'UNIFORM' +'TTEST_NONC' +'WEIBULL' +'CHI' +'INVGAUSS' +'EXTVAL' +'PVAL' +'LOGPVAL' +'LOG10PVAL' +'ESTIMATE' +'LABEL' +'NEURONAMES' +'GENMATRIX' +'SYMMATRIX' +'DISPVECT' +'QUATERNION' +'DIMLESS' +'TIME_SERIES' +'RGB_VECTOR' +'RGBA_VECTOR' +'NODE_INDEX' +'SHAPE' +'CONNECTIVITY_DENSE' +'CONNECTIVITY_DENSE_TIME' +'CONNECTIVITY_PARCELLATED' +'CONNECTIVITY_PARCELLATED_TIME' +'CONNECTIVITY_CONNECTIVITY_TRAJECTORY' +}; diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/read_gifti_file.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/read_gifti_file.m new file mode 100755 index 0000000..75c8321 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/read_gifti_file.m @@ -0,0 +1,196 @@ +function this = read_gifti_file(filename, this) +% Low level reader of GIfTI 1.0 files +% FORMAT this = read_gifti_file(filename, this) +% filename - XML GIfTI filename +% this - structure with fields 'metaData', 'label' and 'data'. +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: read_gifti_file.m 4612 2012-01-08 11:54:26Z guillaume $ + +% Import XML-based GIfTI file +%-------------------------------------------------------------------------- +try + t = xmltree(filename); +catch + error('[GIFTI] Loading of XML file %s failed.', filename); +end + +% Root element of a GIFTI file +%-------------------------------------------------------------------------- +if ~strcmp(get(t,root(t),'name'),'GIFTI') + error('[GIFTI] %s is not a GIFTI 1.0 file.', filename); +end +attr = cell2mat(attributes(t,'get',root(t))); +attr = cell2struct({attr.val},strrep({attr.key},':','___'),2); +if ~all(ismember({'Version','NumberOfDataArrays'},fieldnames(attr))) + error('[GIFTI] Missing mandatory attributes for GIFTI root element.'); +end +if str2double(attr.Version) ~= 1 + warning('[GIFTI] Unknown specification version of GIFTI file (%s).',attr.Version); +end +nbData = str2double(attr.NumberOfDataArrays); + +% Read children elements +%-------------------------------------------------------------------------- +uid = children(t,root(t)); +for i=1:length(uid) + switch get(t,uid(i),'name') + case 'MetaData' + this.metadata = gifti_MetaData(t,uid(i)); + case 'LabelTable' + this.label = gifti_LabelTable(t,uid(i)); + case 'DataArray' + this.data{end+1} = gifti_DataArray(t,uid(i),filename); + otherwise + warning('[GIFTI] Unknown element "%s": ignored.',get(t,uid(i),'name')); + end +end + +if nbData ~= length(this.data) + warning('[GIFTI] Mismatch between expected and effective number of datasets.'); +end + +%========================================================================== +function s = gifti_MetaData(t,uid) +s = struct('name',{}, 'value',{}); +c = children(t,uid); +for i=1:length(c) + for j=children(t,c(i)) + s(i).(lower(get(t,j,'name'))) = get(t,children(t,j),'value'); + end +end + +%========================================================================== +function s = gifti_LabelTable(t,uid) +s = struct('name',{}, 'key',[], 'rgba',[]); +c = children(t,uid); +for i=1:length(c) + a = attributes(t,'get',c(i)); + s(1).rgba(i,1:4) = NaN; + for j=1:numel(a) + switch lower(a{j}.key) + case {'key','index'} + s(1).key(i) = str2double(a{j}.val); + case 'red' + s(1).rgba(i,1) = str2double(a{j}.val); + case 'green' + s(1).rgba(i,2) = str2double(a{j}.val); + case 'blue' + s(1).rgba(i,3) = str2double(a{j}.val); + case 'alpha' + s(1).rgba(i,4) = str2double(a{j}.val); + otherwise + end + end + s(1).name{i} = get(t,children(t,c(i)),'value'); +end + +%========================================================================== +function s = gifti_DataArray(t,uid,filename) +s = struct(... + 'attributes', {}, ... + 'data', {}, ... + 'metadata', struct([]), ... + 'space', {} ... + ); + +attr = cell2mat(attributes(t,'get',uid)); +s(1).attributes = cell2struct({attr.val},{attr.key},2); +s(1).attributes.Dim = []; +for i=1:str2double(s(1).attributes.Dimensionality) + f = sprintf('Dim%d',i-1); + s(1).attributes.Dim(i) = str2double(s(1).attributes.(f)); + s(1).attributes = rmfield(s(1).attributes,f); +end +s(1).attributes = rmfield(s(1).attributes,'Dimensionality'); +if isfield(s(1).attributes,'ExternalFileName') && ... + ~isempty(s(1).attributes.ExternalFileName) + s(1).attributes.ExternalFileName = fullfile(fileparts(filename),... + s(1).attributes.ExternalFileName); +end + +c = children(t,uid); +for i=1:length(c) + switch get(t,c(i),'name') + case 'MetaData' + s(1).metadata = gifti_MetaData(t,c(i)); + case 'CoordinateSystemTransformMatrix' + s(1).space(end+1) = gifti_Space(t,c(i)); + case 'Data' + s(1).data = gifti_Data(t,c(i),s(1).attributes); + otherwise + error('[GIFTI] Unknown DataArray element "%s".',get(t,c(i),'name')); + end +end + +%========================================================================== +function s = gifti_Space(t,uid) +s = struct('DataSpace','', 'TransformedSpace','', 'MatrixData',[]); +for i=children(t,uid) + s.(get(t,i,'name')) = get(t,children(t,i),'value'); +end +s.MatrixData = reshape(str2num(s.MatrixData),4,4)'; + +%========================================================================== +function d = gifti_Data(t,uid,s) +tp = getdict; +try + tp = tp.(s.DataType); +catch + error('[GIFTI] Unknown DataType.'); +end + +[unused,unused,mach] = fopen(1); +sb = @deal; %inline('x'); +try + if (strcmp(s.Endian,'LittleEndian') && ~isempty(strmatch('ieee-be',mach))) ... + || (strcmp(s.Endian,'BigEndian') && ~isempty(strmatch('ieee-le',mach))) + sb = @swapbyte; + end +catch + % Byte Order can be absent if encoding is ASCII, assume native otherwise +end + +switch s.Encoding + case 'ASCII' + d = sscanf(get(t,children(t,uid),'value'),tp.format); + + case 'Base64Binary' + d = typecast(sb(base64decode(get(t,children(t,uid),'value'))), tp.cast); + + case 'GZipBase64Binary' + d = typecast(dunzip(sb(base64decode(get(t,children(t,uid),'value')))), tp.cast); + + case 'ExternalFileBinary' + [p,f,e] = fileparts(s.ExternalFileName); + if isempty(p) + s.ExternalFileName = fullfile(pwd,[f e]); + end + if true + fid = fopen(s.ExternalFileName,'r'); + if fid == -1 + error('[GIFTI] Unable to read binary file %s.',s.ExternalFileName); + end + fseek(fid,str2double(s.ExternalFileOffset),0); + d = sb(fread(fid,prod(s.Dim),['*' tp.class])); + fclose(fid); + else + d = file_array(s.ExternalFileName, s.Dim, tp.class, ... + str2double(s.ExternalFileOffset),1,0,'ro'); + end + + otherwise + error('[GIFTI] Unknown data encoding: %s.',s.Encoding); +end + +if length(s.Dim) == 1, s.Dim(end+1) = 1; end +switch s.ArrayIndexingOrder + case 'RowMajorOrder' + d = permute(reshape(d,fliplr(s.Dim)),length(s.Dim):-1:1); + case 'ColumnMajorOrder' + d = reshape(d,s.Dim); + otherwise + error('[GIFTI] Unknown array indexing order.'); +end diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/read_gifti_file_standalone.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/read_gifti_file_standalone.m new file mode 100755 index 0000000..0b25a2b --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/read_gifti_file_standalone.m @@ -0,0 +1,235 @@ +function this = read_gifti_file(filename, this) +% Low level reader of GIfTI 1.0 files +% FORMAT this = read_gifti_file(filename, this) +% filename - XML GIfTI filename +% this - structure with fields 'metaData', 'label' and 'data'. +%__________________________________________________________________________ +% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging + +% Guillaume Flandin +% $Id: read_gifti_file.m 4612 2012-01-08 11:54:26Z guillaume $ + +% Import XML-based GIfTI file +%-------------------------------------------------------------------------- +try + fid = fopen(filename,'rt'); + xmlstr = fread(fid,'*char')'; + fclose(fid); + t = xml_parser(xmlstr); +catch + error('[GIFTI] Loading of XML file %s failed.', filename); +end + +% Root element of a GIFTI file +%-------------------------------------------------------------------------- +if ~strcmp(xml_get(t,xml_root(t),'name'),'GIFTI') + error('[GIFTI] %s is not a GIFTI 1.0 file.', filename); +end +attr = cell2mat(xml_attributes(t,'get',xml_root(t))); +attr = cell2struct({attr.val},strrep({attr.key},':','___'),2); +if ~all(ismember({'Version','NumberOfDataArrays'},fieldnames(attr))) + error('[GIFTI] Missing mandatory attributes for GIFTI root element.'); +end +if str2double(attr.Version) ~= 1 + warning('[GIFTI] Unknown specification version of GIFTI file (%s).',attr.Version); +end +nbData = str2double(attr.NumberOfDataArrays); + +% Read children elements +%-------------------------------------------------------------------------- +uid = xml_children(t,xml_root(t)); +for i=1:length(uid) + switch xml_get(t,uid(i),'name') + case 'MetaData' + this.metadata = gifti_MetaData(t,uid(i)); + case 'LabelTable' + this.label = gifti_LabelTable(t,uid(i)); + case 'DataArray' + this.data{end+1} = gifti_DataArray(t,uid(i),filename); + otherwise + warning('[GIFTI] Unknown element "%s": ignored.',xml_get(t,uid(i),'name')); + end +end + +if nbData ~= length(this.data) + warning('[GIFTI] Mismatch between expected and effective number of datasets.'); +end + +%========================================================================== +function s = gifti_MetaData(t,uid) +s = struct('name',{}, 'value',{}); +c = xml_children(t,uid); +for i=1:length(c) + for j=xml_children(t,c(i)) + s(i).(lower(xml_get(t,j,'name'))) = xml_get(t,xml_children(t,j),'value'); + end +end + +%========================================================================== +function s = gifti_LabelTable(t,uid) +s = struct('name',{}, 'key',[], 'rgba',[]); +c = xml_children(t,uid); +for i=1:length(c) + a = xml_attributes(t,'get',c(i)); + s(1).rgba(i,1:4) = NaN; + for j=1:numel(a) + switch lower(a{j}.key) + case {'key','index'} + s(1).key(i) = str2double(a{j}.val); + case 'red' + s(1).rgba(i,1) = str2double(a{j}.val); + case 'green' + s(1).rgba(i,2) = str2double(a{j}.val); + case 'blue' + s(1).rgba(i,3) = str2double(a{j}.val); + case 'alpha' + s(1).rgba(i,4) = str2double(a{j}.val); + otherwise + end + end + s(1).name{i} = xml_get(t,xml_children(t,c(i)),'value'); +end + +%========================================================================== +function s = gifti_DataArray(t,uid,filename) +s = struct(... + 'attributes', {}, ... + 'data', {}, ... + 'metadata', struct([]), ... + 'space', {} ... + ); + +attr = cell2mat(xml_attributes(t,'get',uid)); +s(1).attributes = cell2struct({attr.val},{attr.key},2); +s(1).attributes.Dim = []; +for i=1:str2double(s(1).attributes.Dimensionality) + f = sprintf('Dim%d',i-1); + s(1).attributes.Dim(i) = str2double(s(1).attributes.(f)); + s(1).attributes = rmfield(s(1).attributes,f); +end +s(1).attributes = rmfield(s(1).attributes,'Dimensionality'); +if isfield(s(1).attributes,'ExternalFileName') && ... + ~isempty(s(1).attributes.ExternalFileName) + s(1).attributes.ExternalFileName = fullfile(fileparts(filename),... + s(1).attributes.ExternalFileName); +end + +c = xml_children(t,uid); +for i=1:length(c) + switch xml_get(t,c(i),'name') + case 'MetaData' + s(1).metadata = gifti_MetaData(t,c(i)); + case 'CoordinateSystemTransformMatrix' + s(1).space(end+1) = gifti_Space(t,c(i)); + case 'Data' + s(1).data = gifti_Data(t,c(i),s(1).attributes); + otherwise + error('[GIFTI] Unknown DataArray element "%s".',xml_get(t,c(i),'name')); + end +end + +%========================================================================== +function s = gifti_Space(t,uid) +s = struct('DataSpace','', 'TransformedSpace','', 'MatrixData',[]); +for i=xml_children(t,uid) + s.(xml_get(t,i,'name')) = xml_get(t,xml_children(t,i),'value'); +end +s.MatrixData = reshape(str2num(s.MatrixData),4,4)'; + +%========================================================================== +function d = gifti_Data(t,uid,s) +tp = getdict; +try + tp = tp.(s.DataType); +catch + error('[GIFTI] Unknown DataType.'); +end + +[unused,unused,mach] = fopen(1); +sb = @deal; %inline('x'); +try + if (strcmp(s.Endian,'LittleEndian') && ~isempty(strmatch('ieee-be',mach))) ... + || (strcmp(s.Endian,'BigEndian') && ~isempty(strmatch('ieee-le',mach))) + sb = @swapbyte; + end +catch + % Byte Order can be absent if encoding is ASCII, assume native otherwise +end + +switch s.Encoding + case 'ASCII' + d = sscanf(xml_get(t,xml_children(t,uid),'value'),tp.format); + + case 'Base64Binary' + d = typecast(sb(base64decode(xml_get(t,xml_children(t,uid),'value'))), tp.cast); + + case 'GZipBase64Binary' + d = typecast(dunzip(sb(base64decode(xml_get(t,xml_children(t,uid),'value')))), tp.cast); + + case 'ExternalFileBinary' + [p,f,e] = fileparts(s.ExternalFileName); + if isempty(p) + s.ExternalFileName = fullfile(pwd,[f e]); + end + if true + fid = fopen(s.ExternalFileName,'r'); + if fid == -1 + error('[GIFTI] Unable to read binary file %s.',s.ExternalFileName); + end + fseek(fid,str2double(s.ExternalFileOffset),0); + d = sb(fread(fid,prod(s.Dim),['*' tp.class])); + fclose(fid); + else + d = file_array(s.ExternalFileName, s.Dim, tp.class, ... + str2double(s.ExternalFileOffset),1,0,'ro'); + end + + otherwise + error('[GIFTI] Unknown data encoding: %s.',s.Encoding); +end + +if length(s.Dim) == 1, s.Dim(end+1) = 1; end +switch s.ArrayIndexingOrder + case 'RowMajorOrder' + d = permute(reshape(d,fliplr(s.Dim)),length(s.Dim):-1:1); + case 'ColumnMajorOrder' + d = reshape(d,s.Dim); + otherwise + error('[GIFTI] Unknown array indexing order.'); +end + +%========================================================================== +%========================================================================== +function uid = xml_root(tree) +uid = 1; +for i=1:length(tree) + if strcmp(xml_get(tree,i,'type'),'element') + uid = i; + break + end +end + +%-------------------------------------------------------------------------- +function child = xml_children(tree,uid) +if strcmp(tree{uid}.type,'element') + child = tree{uid}.contents; +else + child = []; +end + +%-------------------------------------------------------------------------- +function value = xml_get(tree,uid,parameter) +try + value = tree{uid}.(parameter); +catch + error(sprintf('[XML] Parameter %s not found.',parameter)); +end + +%-------------------------------------------------------------------------- +function varargout = xml_attributes(tree,method,uid) +if ~strcmpi(method,'get'), error('[XML] Unknown attributes method.'); end +if isempty(tree{uid}.attributes) + varargout{1} = {}; +else + varargout{1} = tree{uid}.attributes; +end diff --git a/MATLAB/utilities/MatlabCIFTI/@gifti/private/xml_parser.m b/MATLAB/utilities/MatlabCIFTI/@gifti/private/xml_parser.m new file mode 100755 index 0000000..d60f307 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@gifti/private/xml_parser.m @@ -0,0 +1,424 @@ +function tree = xml_parser(xmlstr) +% XML (eXtensible Markup Language) Processor +% FORMAT tree = xml_parser(xmlstr) +% +% xmlstr - XML string to parse +% tree - tree structure corresponding to the XML file +%_______________________________________________________________________ +% +% xml_parser.m is an XML 1.0 (http://www.w3.org/TR/REC-xml) parser +% written in Matlab. It aims to be fully conforming. It is currently not +% a validating XML processor. +% +% A description of the tree structure provided in output is detailed in +% the header of this m-file. +%_______________________________________________________________________ +% Copyright (C) 2002-2008 http://www.artefact.tk/ + +% Guillaume Flandin +% $Id: xml_parser.m 1460 2008-04-21 17:43:18Z guillaume $ + +% XML Processor for MATLAB (The Mathworks, Inc.). +% Copyright (C) 2002-2008 Guillaume Flandin +% +% This program 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 2 +% of the License, or any later version. +% +% This program 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 this program; if not, write to the Free Software +% Foundation Inc, 59 Temple Pl. - Suite 330, Boston, MA 02111-1307, USA. +%----------------------------------------------------------------------- + +% Suggestions for improvement and fixes are always welcome, although no +% guarantee is made whether and when they will be implemented. +% Send requests to +% Check also the latest developments on the following webpage: +% +%----------------------------------------------------------------------- + +% The implementation of this XML parser is much inspired from a +% Javascript parser available at +%----------------------------------------------------------------------- + +% Structure of the output tree: +% There are 5 types of nodes in an XML file: element, chardata, cdata, +% pi and comment. +% Each of them contains an UID (Unique Identifier): an integer between +% 1 and the number of nodes of the XML file. +% +% element (a tag [contents] +% |_ type: 'element' +% |_ name: string +% |_ attributes: cell array of struct 'key' and 'value' or [] +% |_ contents: double array of uid's or [] if empty +% |_ parent: uid of the parent ([] if root) +% |_ uid: double +% +% chardata (a character array) +% |_ type: 'chardata' +% |_ value: string +% |_ parent: uid of the parent +% |_ uid: double +% +% cdata (a litteral string ) +% |_ type: 'cdata' +% |_ value: string +% |_ parent: uid of the parent +% |_ uid: double +% +% pi (a processing instruction ) +% |_ type: 'pi' +% |_ target: string (may be empty) +% |_ value: string +% |_ parent: uid of the parent +% |_ uid: double +% +% comment (a comment ) +% |_ type: 'comment' +% |_ value: string +% |_ parent: uid of the parent +% |_ uid: double +% +%----------------------------------------------------------------------- + +% TODO/BUG/FEATURES: +% - [compile] only a warning if TagStart is empty ? +% - [attribution] should look for " and ' rather than only " +% - [main] with normalize as a preprocessing, CDATA are modified +% - [prolog] look for a DOCTYPE in the whole string even if it occurs +% only in a far CDATA tag, bug even if the doctype is inside a comment +% - [tag_element] erode should replace normalize here +% - remove globals? uppercase globals rather persistent (clear mfile)? +% - xml_findstr is indeed xml_strfind according to Mathworks vocabulary +% - problem with entities: do we need to convert them here? (é) +%----------------------------------------------------------------------- + +%- XML string to parse and number of tags read +global xmlstring Xparse_count xtree; + +%- Check input arguments +error(nargchk(1,1,nargin)); +if isempty(xmlstr) + error('[XML] Not enough parameters.') +elseif ~ischar(xmlstr) || sum(size(xmlstr)>1)>1 + error('[XML] Input must be a string.') +end + +%- Initialize number of tags (<=> uid) +Xparse_count = 0; + +%- Remove prolog and white space characters from the XML string +xmlstring = normalize(prolog(xmlstr)); + +%- Initialize the XML tree +xtree = {}; +tree = fragment; +tree.str = 1; +tree.parent = 0; + +%- Parse the XML string +tree = compile(tree); + +%- Return the XML tree +tree = xtree; + +%- Remove global variables from the workspace +clear global xmlstring Xparse_count xtree; + +%======================================================================= +% SUBFUNCTIONS + +%----------------------------------------------------------------------- +function frag = compile(frag) + global xmlstring xtree Xparse_count; + + while 1, + if length(xmlstring)<=frag.str || ... + (frag.str == length(xmlstring)-1 && strcmp(xmlstring(frag.str:end),' ')) + return + end + TagStart = xml_findstr(xmlstring,'<',frag.str,1); + if isempty(TagStart) + %- Character data + error('[XML] Unknown data at the end of the XML file.'); + xtree{Xparse_count} = chardata; + xtree{Xparse_count}.value = erode(entity(xmlstring(frag.str:end))); + xtree{Xparse_count}.parent = frag.parent; + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + frag.str = ''; + elseif TagStart > frag.str + if strcmp(xmlstring(frag.str:TagStart-1),' ') + %- A single white space before a tag (ignore) + frag.str = TagStart; + else + %- Character data + xtree{Xparse_count} = chardata; + xtree{Xparse_count}.value = erode(entity(xmlstring(frag.str:TagStart-1))); + xtree{Xparse_count}.parent = frag.parent; + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + frag.str = TagStart; + end + else + if strcmp(xmlstring(frag.str+1),'?') + %- Processing instruction + frag = tag_pi(frag); + else + if length(xmlstring)-frag.str>4 && strcmp(xmlstring(frag.str+1:frag.str+3),'!--') + %- Comment + frag = tag_comment(frag); + else + if length(xmlstring)-frag.str>9 && strcmp(xmlstring(frag.str+1:frag.str+8),'![CDATA[') + %- Litteral data + frag = tag_cdata(frag); + else + %- A tag element (empty (<.../>) or not) + if ~isempty(frag.end) + endmk = ['/' frag.end '>']; + else + endmk = '/>'; + end + if strcmp(xmlstring(frag.str+1:frag.str+length(frag.end)+2),endmk) || ... + strcmp(strip(xmlstring(frag.str+1:frag.str+length(frag.end)+2)),endmk) + frag.str = frag.str + length(frag.end)+3; + return + else + frag = tag_element(frag); + end + end + end + end + end + end + +%----------------------------------------------------------------------- +function frag = tag_element(frag) + global xmlstring xtree Xparse_count; + close = xml_findstr(xmlstring,'>',frag.str,1); + if isempty(close) + error('[XML] Tag < opened but not closed.'); + else + empty = strcmp(xmlstring(close-1:close),'/>'); + if empty + close = close - 1; + end + starttag = normalize(xmlstring(frag.str+1:close-1)); + nextspace = xml_findstr(starttag,' ',1,1); + attribs = ''; + if isempty(nextspace) + name = starttag; + else + name = starttag(1:nextspace-1); + attribs = starttag(nextspace+1:end); + end + xtree{Xparse_count} = element; + xtree{Xparse_count}.name = strip(name); + if frag.parent + xtree{Xparse_count}.parent = frag.parent; + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + end + if ~isempty(attribs) + xtree{Xparse_count}.attributes = attribution(attribs); + end + if ~empty + contents = fragment; + contents.str = close+1; + contents.end = name; + contents.parent = Xparse_count; + contents = compile(contents); + frag.str = contents.str; + else + frag.str = close+2; + end + end + +%----------------------------------------------------------------------- +function frag = tag_pi(frag) + global xmlstring xtree Xparse_count; + close = xml_findstr(xmlstring,'?>',frag.str,1); + if isempty(close) + warning('[XML] Tag close || nextspace == frag.str+2 + xtree{Xparse_count}.value = erode(xmlstring(frag.str+2:close-1)); + else + xtree{Xparse_count}.value = erode(xmlstring(nextspace+1:close-1)); + xtree{Xparse_count}.target = erode(xmlstring(frag.str+2:nextspace)); + end + if frag.parent + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + xtree{Xparse_count}.parent = frag.parent; + end + frag.str = close+2; + end + +%----------------------------------------------------------------------- +function frag = tag_comment(frag) + global xmlstring xtree Xparse_count; + close = xml_findstr(xmlstring,'-->',frag.str,1); + if isempty(close) + warning('[XML] Tag ) +% |_ type: 'comment' +% |_ value: string +% |_ parent: uid of the parent +% |_ uid: double +% +%----------------------------------------------------------------------- + +% TODO/BUG/FEATURES: +% - [compile] only a warning if TagStart is empty ? +% - [attribution] should look for " and ' rather than only " +% - [main] with normalize as a preprocessing, CDATA are modified +% - [prolog] look for a DOCTYPE in the whole string even if it occurs +% only in a far CDATA tag, bug even if the doctype is inside a comment +% - [tag_element] erode should replace normalize here +% - remove globals? uppercase globals rather persistent (clear mfile)? +% - xml_findstr is indeed xml_strfind according to Mathworks vocabulary +% - problem with entities: do we need to convert them here? (é) +%----------------------------------------------------------------------- + +%- XML string to parse and number of tags read +global xmlstring Xparse_count xtree; + +%- Check input arguments +error(nargchk(1,1,nargin)); +if isempty(xmlstr) + error('[XML] Not enough parameters.') +elseif ~ischar(xmlstr) || sum(size(xmlstr)>1)>1 + error('[XML] Input must be a string.') +end + +%- Initialize number of tags (<=> uid) +Xparse_count = 0; + +%- Remove prolog and white space characters from the XML string +xmlstring = normalize(prolog(xmlstr)); + +%- Initialize the XML tree +xtree = {}; +tree = fragment; +tree.str = 1; +tree.parent = 0; + +%- Parse the XML string +tree = compile(tree); + +%- Return the XML tree +tree = xtree; + +%- Remove global variables from the workspace +clear global xmlstring Xparse_count xtree; + +%======================================================================= +% SUBFUNCTIONS + +%----------------------------------------------------------------------- +function frag = compile(frag) + global xmlstring xtree Xparse_count; + + while 1, + if length(xmlstring)<=frag.str || ... + (frag.str == length(xmlstring)-1 && strcmp(xmlstring(frag.str:end),' ')) + return + end + TagStart = xml_findstr(xmlstring,'<',frag.str,1); + if isempty(TagStart) + %- Character data + error('[XML] Unknown data at the end of the XML file.'); + xtree{Xparse_count} = chardata; + xtree{Xparse_count}.value = erode(entity(xmlstring(frag.str:end))); + xtree{Xparse_count}.parent = frag.parent; + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + frag.str = ''; + elseif TagStart > frag.str + if strcmp(xmlstring(frag.str:TagStart-1),' ') + %- A single white space before a tag (ignore) + frag.str = TagStart; + else + %- Character data + xtree{Xparse_count} = chardata; + xtree{Xparse_count}.value = erode(entity(xmlstring(frag.str:TagStart-1))); + xtree{Xparse_count}.parent = frag.parent; + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + frag.str = TagStart; + end + else + if strcmp(xmlstring(frag.str+1),'?') + %- Processing instruction + frag = tag_pi(frag); + else + if length(xmlstring)-frag.str>4 && strcmp(xmlstring(frag.str+1:frag.str+3),'!--') + %- Comment + frag = tag_comment(frag); + else + if length(xmlstring)-frag.str>9 && strcmp(xmlstring(frag.str+1:frag.str+8),'![CDATA[') + %- Litteral data + frag = tag_cdata(frag); + else + %- A tag element (empty (<.../>) or not) + if ~isempty(frag.end) + endmk = ['/' frag.end '>']; + else + endmk = '/>'; + end + if strcmp(xmlstring(frag.str+1:frag.str+length(frag.end)+2),endmk) || ... + strcmp(strip(xmlstring(frag.str+1:frag.str+length(frag.end)+2)),endmk) + frag.str = frag.str + length(frag.end)+3; + return + else + frag = tag_element(frag); + end + end + end + end + end + end + +%----------------------------------------------------------------------- +function frag = tag_element(frag) + global xmlstring xtree Xparse_count; + close = xml_findstr(xmlstring,'>',frag.str,1); + if isempty(close) + error('[XML] Tag < opened but not closed.'); + else + empty = strcmp(xmlstring(close-1:close),'/>'); + if empty + close = close - 1; + end + starttag = normalize(xmlstring(frag.str+1:close-1)); + nextspace = xml_findstr(starttag,' ',1,1); + attribs = ''; + if isempty(nextspace) + name = starttag; + else + name = starttag(1:nextspace-1); + attribs = starttag(nextspace+1:end); + end + xtree{Xparse_count} = element; + xtree{Xparse_count}.name = strip(name); + if frag.parent + xtree{Xparse_count}.parent = frag.parent; + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + end + if ~isempty(attribs) + xtree{Xparse_count}.attributes = attribution(attribs); + end + if ~empty + contents = fragment; + contents.str = close+1; + contents.end = name; + contents.parent = Xparse_count; + contents = compile(contents); + frag.str = contents.str; + else + frag.str = close+2; + end + end + +%----------------------------------------------------------------------- +function frag = tag_pi(frag) + global xmlstring xtree Xparse_count; + close = xml_findstr(xmlstring,'?>',frag.str,1); + if isempty(close) + warning('[XML] Tag close || nextspace == frag.str+2 + xtree{Xparse_count}.value = erode(xmlstring(frag.str+2:close-1)); + else + xtree{Xparse_count}.value = erode(xmlstring(nextspace+1:close-1)); + xtree{Xparse_count}.target = erode(xmlstring(frag.str+2:nextspace)); + end + if frag.parent + xtree{frag.parent}.contents = [xtree{frag.parent}.contents Xparse_count]; + xtree{Xparse_count}.parent = frag.parent; + end + frag.str = close+2; + end + +%----------------------------------------------------------------------- +function frag = tag_comment(frag) + global xmlstring xtree Xparse_count; + close = xml_findstr(xmlstring,'-->',frag.str,1); + if isempty(close) + warning('[XML] Tag blah blah +% Now root is the first element node of the tree. + +% Look for the first element in the XML Tree +for i=1:length(tree) + if strcmp(get(tree,i,'type'),'element') + uid = i; + break + end +end diff --git a/MATLAB/utilities/MatlabCIFTI/@xmltree/save.m b/MATLAB/utilities/MatlabCIFTI/@xmltree/save.m new file mode 100755 index 0000000..4e67d5e --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@xmltree/save.m @@ -0,0 +1,133 @@ +function varargout = save(tree, filename) +% XMLTREE/SAVE Save an XML tree in an XML file +% FORMAT varargout = save(tree,filename) +% +% tree - XMLTree +% filename - XML output filename +% varargout - XML string +%_______________________________________________________________________ +% +% Convert an XML tree into a well-formed XML string and write it into +% a file or return it as a string if no filename is provided. +%_______________________________________________________________________ +% Copyright (C) 2002-2008 http://www.artefact.tk/ + +% Guillaume Flandin +% $Id: save.m 1460 2008-04-21 17:43:18Z guillaume $ + +error(nargchk(1,2,nargin)); + +prolog = '\n'; + +%- Return the XML tree as a string +if nargin == 1 + varargout{1} = [sprintf(prolog) ... + print_subtree(tree,'',root(tree))]; +%- Output specified +else + %- Filename provided + if isstr(filename) + [fid, msg] = fopen(filename,'w'); + if fid==-1, error(msg); end + if isempty(tree.filename), tree.filename = filename; end + %- File identifier provided + elseif isnumeric(filename) && numel(filename) == 1 + fid = filename; + prolog = ''; %- With this option, do not write any prolog + else + error('[XMLTree] Invalid argument.'); + end + fprintf(fid,prolog); + save_subtree(tree,fid,root(tree)); + if ischar(filename), fclose(fid); end + if nargout == 1 + varargout{1} = print_subtree(tree,'',root(tree)); + end +end + +%======================================================================= +function xmlstr = print_subtree(tree,xmlstr,uid,order) + if nargin < 4, order = 0; end + xmlstr = [xmlstr blanks(3*order)]; + switch tree.tree{uid}.type + case 'element' + xmlstr = sprintf('%s<%s',xmlstr,tree.tree{uid}.name); + for i=1:length(tree.tree{uid}.attributes) + xmlstr = sprintf('%s %s="%s"', xmlstr, ... + tree.tree{uid}.attributes{i}.key,... + tree.tree{uid}.attributes{i}.val); + end + if isempty(tree.tree{uid}.contents) + xmlstr = sprintf('%s/>\n',xmlstr); + else + xmlstr = sprintf('%s>\n',xmlstr); + for i=1:length(tree.tree{uid}.contents) + xmlstr = print_subtree(tree,xmlstr, ... + tree.tree{uid}.contents(i),order+1); + end + xmlstr = [xmlstr blanks(3*order)]; + xmlstr = sprintf('%s\n',xmlstr,... + tree.tree{uid}.name); + end + case 'chardata' + xmlstr = sprintf('%s%s\n',xmlstr, ... + entity(tree.tree{uid}.value)); + case 'cdata' + xmlstr = sprintf('%s\n',xmlstr, ... + tree.tree{uid}.value); + case 'pi' + xmlstr = sprintf('%s\n',xmlstr, ... + tree.tree{uid}.target, tree.tree{uid}.value); + case 'comment' + xmlstr = sprintf('%s\n',xmlstr,... + tree.tree{uid}.value); + otherwise + warning(sprintf('Type %s unknown: not saved', ... + tree.tree{uid}.type)); + end + +%======================================================================= +function save_subtree(tree,fid,uid,order) + if nargin < 4, order = 0; end + fprintf(fid,blanks(3*order)); + switch tree.tree{uid}.type + case 'element' + fprintf(fid,'<%s',tree.tree{uid}.name); + for i=1:length(tree.tree{uid}.attributes) + fprintf(fid,' %s="%s"',... + tree.tree{uid}.attributes{i}.key, ... + tree.tree{uid}.attributes{i}.val); + end + if isempty(tree.tree{uid}.contents) + fprintf(fid,'/>\n'); + else + fprintf(fid,'>\n'); + for i=1:length(tree.tree{uid}.contents) + save_subtree(tree,fid,... + tree.tree{uid}.contents(i),order+1) + end + fprintf(fid,blanks(3*order)); + fprintf(fid,'\n',tree.tree{uid}.name); + end + case 'chardata' + fprintf(fid,'%s\n',entity(tree.tree{uid}.value)); + case 'cdata' + fprintf(fid,'\n',tree.tree{uid}.value); + case 'pi' + fprintf(fid,'\n',tree.tree{uid}.target, ... + tree.tree{uid}.value); + case 'comment' + fprintf(fid,'\n',tree.tree{uid}.value); + otherwise + warning(sprintf('[XMLTree] Type %s unknown: not saved', ... + tree.tree{uid}.type)); + end + + +%======================================================================= +function str = entity(str) + str = strrep(str,'&','&'); + str = strrep(str,'<','<'); + str = strrep(str,'>','>'); + str = strrep(str,'"','"'); + str = strrep(str,'''','''); diff --git a/MATLAB/utilities/MatlabCIFTI/@xmltree/set.m b/MATLAB/utilities/MatlabCIFTI/@xmltree/set.m new file mode 100755 index 0000000..ea5b118 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@xmltree/set.m @@ -0,0 +1,26 @@ +function tree = set(tree,uid, parameter, value) +% XMLTREE/SET Method (set object properties) +% FORMAT tree = set(tree,uid,parameter,value) +% +% tree - XMLTree object +% uid - array (or cell) of uid's +% parameter - property name +% value - property value +%_______________________________________________________________________ +% +% Set object properties given its uid and pairs parameter/value +% The tree parameter must be in input AND in output +%_______________________________________________________________________ +% Copyright (C) 2002-2008 http://www.artefact.tk/ + +% Guillaume Flandin +% $Id: set.m 1460 2008-04-21 17:43:18Z guillaume $ + +error(nargchk(4,4,nargin)); + +if iscell(uid), uid = [uid{:}]; else uid = uid(:); end + +for i=1:length(uid) + tree.tree{uid(i)} = builtin('subsasgn', tree.tree{uid(i)}, struct('type','.','subs',parameter), value); + %tree.tree{uid(i)} = setfield(tree.tree{uid(i)},parameter,value); +end diff --git a/MATLAB/utilities/MatlabCIFTI/@xmltree/setfilename.m b/MATLAB/utilities/MatlabCIFTI/@xmltree/setfilename.m new file mode 100755 index 0000000..da7ef69 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@xmltree/setfilename.m @@ -0,0 +1,16 @@ +function tree = setfilename(tree,filename) +% XMLTREE/SETFILENAME Set filename method +% FORMAT tree = setfilename(tree,filename) +% +% tree - XMLTree object +% filename - XML filename +%_______________________________________________________________________ +% +% Set the filename linked to the XML tree as filename. +%_______________________________________________________________________ +% Copyright (C) 2002-2008 http://www.artefact.tk/ + +% Guillaume Flandin +% $Id: setfilename.m 1460 2008-04-21 17:43:18Z guillaume $ + +tree.filename = filename; diff --git a/MATLAB/utilities/MatlabCIFTI/@xmltree/view.m b/MATLAB/utilities/MatlabCIFTI/@xmltree/view.m new file mode 100755 index 0000000..e788ddf --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@xmltree/view.m @@ -0,0 +1,17 @@ +function view(tree) +% XMLTREE/VIEW View Method +% FORMAT view(tree) +% +% tree - XMLTree object +%_______________________________________________________________________ +% +% Display an XML tree in a graphical interface +%_______________________________________________________________________ +% Copyright (C) 2002-2008 http://www.artefact.tk/ + +% Guillaume Flandin +% $Id: view.m 1460 2008-04-21 17:43:18Z guillaume $ + +error(nargchk(1,1,nargin)); + +editor(tree); diff --git a/MATLAB/utilities/MatlabCIFTI/@xmltree/xmltree.m b/MATLAB/utilities/MatlabCIFTI/@xmltree/xmltree.m new file mode 100755 index 0000000..916a23b --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/@xmltree/xmltree.m @@ -0,0 +1,61 @@ +function tree = xmltree(varargin) +% XMLTREE/XMLTREE Constructor of the XMLTree class +% FORMAT tree = xmltree(varargin) +% +% varargin - XML filename or XML string +% tree - XMLTree Object +% +% tree = xmltree; % creates a minimal XML tree: '' +% tree = xmltree('foo.xml'); % creates a tree from XML file 'foo.xml' +% tree = xmltree('content') % creates a tree from string +%_______________________________________________________________________ +% +% This is the constructor of the XMLTree class. +% It creates a tree of an XML 1.0 file (after parsing) that is stored +% using a Document Object Model (DOM) representation. +% See http://www.w3.org/TR/REC-xml for details about XML 1.0. +% See http://www.w3.org/DOM/ for details about DOM platform. +%_______________________________________________________________________ +% Copyright (C) 2002-2008 http://www.artefact.tk/ + +% Guillaume Flandin +% $Id: xmltree.m 1460 2008-04-21 17:43:18Z guillaume $ + +switch(nargin) + case 0 + tree.tree{1} = struct('type','element',... + 'name','tag',... + 'attributes',[],... + 'contents',[],... + 'parent',[],... + 'uid',1); + tree.filename = ''; + tree = class(tree,'xmltree'); + case 1 + if isa(varargin{1},'xmltree') + tree = varargin{1}; + elseif ischar(varargin{1}) + % Input argument is an XML string + if (~exist(varargin{1},'file') && ... + ~isempty(xml_findstr(varargin{1},'<',1,1))) + tree.tree = xml_parser(varargin{1}); + tree.filename = ''; + % Input argument is an XML filename + else + fid = fopen(varargin{1},'rt'); + if (fid == -1) + error(['[XMLTree] Cannot open ' varargin{1}]); + end + xmlstr = fread(fid,'*char')'; + %xmlstr = fscanf(fid,'%c'); + fclose(fid); + tree.tree = xml_parser(xmlstr); + tree.filename = varargin{1}; + end + tree = class(tree,'xmltree'); + else + error('[XMLTree] Bad input argument'); + end + otherwise + error('[XMLTree] Too many input arguments'); +end diff --git a/MATLAB/utilities/MatlabCIFTI/README b/MATLAB/utilities/MatlabCIFTI/README new file mode 100755 index 0000000..ebd1ee1 --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/README @@ -0,0 +1,16 @@ +This package contains and requires: + +1) A caret7 caret_command binary for 64bit Linux +2) An updated/bugfixed version of the GIFTI matlab library (based on version 1.1 from http://www.artefact.tk/software/matlab/gifti/, but now properly handles GIFTI label tables/files) +3) ciftiopen.m and ciftisave.m, wrappers that convert CIFTI files into GIFTI external binary files and then read the into matlab, or save a GIFTI external binary file and convert it to a CIFTI file + +Usage: +gifti_object = ciftiopen('path/to/file.dtseries.nii','/path/to/caret7/caret_command'); + +Outputs a GIFTI object. Access CIFTI matrix data with gifti_object.cdata. Access the whole struct with gifti_object.private + +ciftisave(gifti_object,'path/to/file.dtseries.nii','/path/to/caret7/caret_command'); + +Saves a GIFTI object. Need to reuse an input GIFTI object that contains the CIFTI header extension in it (i.e. the original gifti_object from ciftiopen). + + diff --git a/MATLAB/utilities/MatlabCIFTI/ciftiopen.m b/MATLAB/utilities/MatlabCIFTI/ciftiopen.m new file mode 100755 index 0000000..4364c7f --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/ciftiopen.m @@ -0,0 +1,23 @@ +function [ cifti ] = ciftiopen(filename,caret7command) +%Open a CIFTI file by converting to GIFTI external binary first and then +%using the GIFTI toolbox + +rng shuffle +RANDOM = randi(1000000,1); +tic +unix([caret7command ' -cifti-convert -to-gifti-ext ' filename ' ' filename num2str(RANDOM) '.gii']); +toc + +%unix(['/media/1TB/matlabsharedcode/cifticlean.sh ' filename num2str(RANDOM) '.gii ' filename num2str(RANDOM) '_.gii']); + +tic +%cifti = gifti([filename num2str(RANDOM) '_.gii']); +cifti = gifti([filename num2str(RANDOM) '.gii']); +toc + +%unix([' rm ' filename num2str(RANDOM) '.gii ' filename num2str(RANDOM) '.gii.data ' filename num2str(RANDOM) '_.gii']); +unix([' rm ' filename num2str(RANDOM) '.gii ' filename num2str(RANDOM) '.gii.data ']); + + +end + diff --git a/MATLAB/utilities/MatlabCIFTI/ciftisave.m b/MATLAB/utilities/MatlabCIFTI/ciftisave.m new file mode 100755 index 0000000..4d8740f --- /dev/null +++ b/MATLAB/utilities/MatlabCIFTI/ciftisave.m @@ -0,0 +1,19 @@ +function [ output_args ] = ciftisave(cifti,filename,caret7command) +%Save a CIFTI file as a GIFTI external binary and then convert it to CIFTI + +tic +save(cifti,[filename '.gii'],'ExternalFileBinary') +toc + +%unix(['/media/1TB/matlabsharedcode/ciftiunclean.sh ' filename '.gii ' filename '_.gii']); + +%unix(['mv ' filename '_.gii ' filename '.gii']); + +tic +unix([caret7command ' -cifti-convert -from-gifti-ext ' filename '.gii ' filename]); +toc + +unix([' rm ' filename '.gii ' filename '.dat ']); + +end + diff --git a/MATLAB/utilities/loadtext.m b/MATLAB/utilities/loadtext.m new file mode 100644 index 0000000..b80d843 --- /dev/null +++ b/MATLAB/utilities/loadtext.m @@ -0,0 +1,33 @@ +function s = loadtext(file); + +% function s = loadtext(file); +% +% is a file location +% +% read and return a cell vector of strings. +% each string corresponds to one line of . +% see also savetext.m. +% +% note that if consists entirely of numeric text, +% it is much faster to use load.m! +% +% example: +% savetext('temp.txt',{'this' 'is' 'a' 'test'}); +% a = loadtext('temp.txt') + +% open file +fid = fopen(file,'r'); +assert(fid ~= -1,' could not be opened for reading'); + +% do it +s = {}; +while 1 + line = fgetl(fid); + if ~ischar(line) + break; + end + s{end+1} = line; +end + +% close file +assert(fclose(fid)==0); diff --git a/MATLAB/utilities/struct2xml.m b/MATLAB/utilities/struct2xml.m new file mode 100755 index 0000000..fe801fe --- /dev/null +++ b/MATLAB/utilities/struct2xml.m @@ -0,0 +1,206 @@ +function varargout = struct2xml( s, varargin ) +%Convert a MATLAB structure into a xml file +% [ ] = struct2xml( s, file ) +% xml = struct2xml( s ) +% +% A structure containing: +% s.XMLname.Attributes.attrib1 = "Some value"; +% s.XMLname.Element.Text = "Some text"; +% s.XMLname.DifferentElement{1}.Attributes.attrib2 = "2"; +% s.XMLname.DifferentElement{1}.Text = "Some more text"; +% s.XMLname.DifferentElement{2}.Attributes.attrib3 = "2"; +% s.XMLname.DifferentElement{2}.Attributes.attrib4 = "1"; +% s.XMLname.DifferentElement{2}.Text = "Even more text"; +% +% Will produce: +% +% Some text +% Some more text +% Even more text +% +% +% Please note that the following strings are substituted +% '_dash_' by '-', '_colon_' by ':' and '_dot_' by '.' +% +% Written by W. Falkena, ASTI, TUDelft, 27-08-2010 +% On-screen output functionality added by P. Orth, 01-12-2010 +% Multiple space to single space conversion adapted for speed by T. Lohuis, 11-04-2011 +% Val2str subfunction bugfix by H. Gsenger, 19-9-2011 +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Modified by K Jamison, University of Minnesota, 01-30-2015 +% - Now writes data type of content as element attribute +% example: s.XMLname.SomeNumber = [6 7 8]; +% produces: 6 7 8 + + if (nargin ~= 2) + if(nargout ~= 1 || nargin ~= 1) + error(['Supported function calls:' sprintf('\n')... + '[ ] = struct2xml( s, file )' sprintf('\n')... + 'xml = struct2xml( s )']); + end + end + + if(nargin == 2) + file = varargin{1}; + + if (isempty(file)) + error('Filename can not be empty'); + end + + if (isempty(strfind(file,'.xml'))) + file = [file '.xml']; + end + end + + if (~isstruct(s)) + error([inputname(1) ' is not a structure']); + end + + if (length(fieldnames(s)) > 1) + error(['Error processing the structure:' sprintf('\n') 'There should be a single field in the main structure.']); + end + xmlname = fieldnames(s); + xmlname = xmlname{1}; + + %substitute special characters + xmlname_sc = xmlname; + xmlname_sc = strrep(xmlname_sc,'_dash_','-'); + xmlname_sc = strrep(xmlname_sc,'_colon_',':'); + xmlname_sc = strrep(xmlname_sc,'_dot_','.'); + + %create xml structure + docNode = com.mathworks.xml.XMLUtils.createDocument(xmlname_sc); + + %process the rootnode + docRootNode = docNode.getDocumentElement; + + %append childs + parseStruct(s.(xmlname),docNode,docRootNode,[inputname(1) '.' xmlname '.']); + + if(nargout == 0) + %save xml file + xmlwrite(file,docNode); + else + varargout{1} = xmlwrite(docNode); + end +end + +% ----- Subfunction parseStruct ----- +function [] = parseStruct(s,docNode,curNode,pName) + + fnames = fieldnames(s); + for i = 1:length(fnames) + curfield = fnames{i}; + + %substitute special characters + curfield_sc = curfield; + curfield_sc = strrep(curfield_sc,'_dash_','-'); + curfield_sc = strrep(curfield_sc,'_colon_',':'); + curfield_sc = strrep(curfield_sc,'_dot_','.'); + + if (strcmp(curfield,'Attributes')) + %Attribute data + if (isstruct(s.(curfield))) + attr_names = fieldnames(s.Attributes); + for a = 1:length(attr_names) + cur_attr = attr_names{a}; + + %substitute special characters + cur_attr_sc = cur_attr; + cur_attr_sc = strrep(cur_attr_sc,'_dash_','-'); + cur_attr_sc = strrep(cur_attr_sc,'_colon_',':'); + cur_attr_sc = strrep(cur_attr_sc,'_dot_','.'); + + [cur_str,succes] = val2str(s.Attributes.(cur_attr)); + if (succes) + curNode.setAttribute(cur_attr_sc,cur_str); + else + disp(['Warning. The text in ' pName curfield '.' cur_attr ' could not be processed.']); + end + end + else + disp(['Warning. The attributes in ' pName curfield ' could not be processed.']); + disp(['The correct syntax is: ' pName curfield '.attribute_name = ''Some text''.']); + end + elseif (strcmp(curfield,'Text')) + %Text data + [txt,succes] = val2str(s.Text); + if (succes) + curNode.appendChild(docNode.createTextNode(txt)); + else + disp(['Warning. The text in ' pName curfield ' could not be processed.']); + end + else + %Sub-element + if (isstruct(s.(curfield))) + %single element + curElement = docNode.createElement(curfield_sc); + curNode.appendChild(curElement); + parseStruct(s.(curfield),docNode,curElement,[pName curfield '.']) + elseif (iscell(s.(curfield))) + %multiple elements + for c = 1:length(s.(curfield)) + curElement = docNode.createElement(curfield_sc); + curNode.appendChild(curElement); + if (isstruct(s.(curfield){c})) + parseStruct(s.(curfield){c},docNode,curElement,[pName curfield '{' num2str(c) '}.']) + else + disp(['Warning. The cell ' pName curfield '{' num2str(c) '} could not be processed, since it contains no structure.']); + end + end + else + %eventhough the fieldname is not text, the field could + %contain text. Create a new element and use this text + curElement = docNode.createElement(curfield_sc); + curNode.appendChild(curElement); + [txt,succes] = val2str(s.(curfield)); + if (succes) + curElement.appendChild(docNode.createTextNode(txt)); + valtype=class(s.(curfield)); + curElement.setAttribute('datatype',valtype); + else + disp(['Warning. The text in ' pName curfield ' could not be processed.']); + end + end + end + end +end + +%----- Subfunction val2str ----- +function [str,succes] = val2str(val) + + succes = true; + str = []; + + if (isempty(val)) + return; %bugfix from H. Gsenger + elseif (ischar(val)) + %do nothing + elseif (isnumeric(val)) + val = num2str(val); + else + succes = false; + end + + if (ischar(val)) + %add line breaks to all lines except the last (for multiline strings) + lines = size(val,1); + val = [val char(sprintf('\n')*[ones(lines-1,1);0])]; + + %transpose is required since indexing (i.e., val(nonspace) or val(:)) produces a 1-D vector. + %This should be row based (line based) and not column based. + valt = val'; + + remove_multiple_white_spaces = true; + if (remove_multiple_white_spaces) + %remove multiple white spaces using isspace, suggestion of T. Lohuis + whitespace = isspace(val); + nonspace = (whitespace + [zeros(lines,1) whitespace(:,1:end-1)])~=2; + nonspace(:,end) = [ones(lines-1,1);0]; %make sure line breaks stay intact + str = valt(nonspace'); + else + str = valt(:); + end + end +end diff --git a/MATLAB/utilities/xml2struct.m b/MATLAB/utilities/xml2struct.m new file mode 100644 index 0000000..f13c1a7 --- /dev/null +++ b/MATLAB/utilities/xml2struct.m @@ -0,0 +1,250 @@ +function [ s ] = xml2struct( file, parsevals ) +%Convert xml file into a MATLAB structure +% [ s ] = xml2struct( file ) % parse numerical fields +% [ s ] = xml2struct( file, false ) % dont parse numerical fields +% +% A file containing: +% +% Some text +% Some more text +% Even more text +% +% +% Will produce: +% s.XMLname.Attributes.attrib1 = "Some value"; +% s.XMLname.Element.Text = "Some text"; +% s.XMLname.DifferentElement{1}.Attributes.attrib2 = "2"; +% s.XMLname.DifferentElement{1}.Text = "Some more text"; +% s.XMLname.DifferentElement{2}.Attributes.attrib3 = "2"; +% s.XMLname.DifferentElement{2}.Attributes.attrib4 = "1"; +% s.XMLname.DifferentElement{2}.Text = "Even more text"; +% +% Please note that the following characters are substituted +% '-' by '_dash_', ':' by '_colon_' and '.' by '_dot_' +% +% Written by W. Falkena, ASTI, TUDelft, 21-08-2010 +% Attribute parsing speed increased by 40% by A. Wanner, 14-6-2011 +% Added CDATA support by I. Smirnov, 20-3-2012 +% +% Modified by X. Mo, University of Wisconsin, 12-5-2012 +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Modified by K Jamison, University of Minnesota, 01-30-2015 +% - Added text2val() and parsevals functionality +% +% If second argument parsevals=true (default), convert elements with "datatype" +% attribute to a specified numeric datatype. +% +% ex: +% struct2xml(S1,file); +% S2=xml2struct(file,true); %output S2 is identical to input S1 +% +% A file containing: +% +% Some text +% 6 7 8 9 +% +% +% Will produce: +% s = xml2struct(file, true); +% s.XMLName.Element = 'Some text'; +% s.XMLName.DifferentElement = [6 7 8 9]; +% + if (nargin < 1) + clc; + help xml2struct + return + end + + if (nargin < 2 || isempty(parsevals)) + parsevals = true; + end + + if isa(file, 'org.apache.xerces.dom.DeferredDocumentImpl') || isa(file, 'org.apache.xerces.dom.DeferredElementImpl') + % input is a java xml object + xDoc = file; + else + %check for existance + if (exist(file,'file') == 0) + %Perhaps the xml extension was omitted from the file name. Add the + %extension and try again. + if (isempty(strfind(file,'.xml'))) + file = [file '.xml']; + end + + if (exist(file,'file') == 0) + error(['The file ' file ' could not be found']); + end + end + %read the xml file + xDoc = xmlread(file); + end + + %parse xDoc into a MATLAB structure + s = parseChildNodes(xDoc); + if(parsevals) + s = text2val(s); + end + +end + +% ----- Subfunction parseChildNodes ----- +function [children,ptext,textflag] = parseChildNodes(theNode) + % Recurse over node children. + children = struct; + ptext = struct; textflag = 'Text'; + if hasChildNodes(theNode) + childNodes = getChildNodes(theNode); + numChildNodes = getLength(childNodes); + + for count = 1:numChildNodes + theChild = item(childNodes,count-1); + [text,name,attr,childs,textflag] = getNodeData(theChild); + + if (~strcmp(name,'#text') && ~strcmp(name,'#comment') && ~strcmp(name,'#cdata_dash_section')) + %XML allows the same elements to be defined multiple times, + %put each in a different cell + if (isfield(children,name)) + if (~iscell(children.(name))) + %put existsing element into cell format + children.(name) = {children.(name)}; + end + index = length(children.(name))+1; + %add new element + children.(name){index} = childs; + if(~isempty(fieldnames(text))) + children.(name){index} = text; + end + if(~isempty(attr)) + children.(name){index}.('Attributes') = attr; + end + else + %add previously unknown (new) element to the structure + children.(name) = childs; + if(~isempty(text) && ~isempty(fieldnames(text))) + children.(name) = text; + end + if(~isempty(attr)) + children.(name).('Attributes') = attr; + end + end + else + ptextflag = 'Text'; + if (strcmp(name, '#cdata_dash_section')) + ptextflag = 'CDATA'; + elseif (strcmp(name, '#comment')) + ptextflag = 'Comment'; + end + + %this is the text in an element (i.e., the parentNode) + if (~isempty(regexprep(text.(textflag),'[\s]*',''))) + if (~isfield(ptext,ptextflag) || isempty(ptext.(ptextflag))) + ptext.(ptextflag) = text.(textflag); + else + %what to do when element data is as follows: + %Text More text + + %put the text in different cells: + % if (~iscell(ptext)) ptext = {ptext}; end + % ptext{length(ptext)+1} = text; + + %just append the text + ptext.(ptextflag) = [ptext.(ptextflag) text.(textflag)]; + end + end + end + + end + end +end + +% ----- Subfunction getNodeData ----- +function [text,name,attr,childs,textflag] = getNodeData(theNode) + % Create structure of node info. + + %make sure name is allowed as structure name + name = toCharArray(getNodeName(theNode))'; + name = strrep(name, '-', '_dash_'); + name = strrep(name, ':', '_colon_'); + name = strrep(name, '.', '_dot_'); + + attr = parseAttributes(theNode); + if (isempty(fieldnames(attr))) + attr = []; + end + + %parse child nodes + [childs,text,textflag] = parseChildNodes(theNode); + + if (isempty(fieldnames(childs)) && isempty(fieldnames(text))) + %get the data of any childless nodes + % faster than if any(strcmp(methods(theNode), 'getData')) + % no need to try-catch (?) + % faster than text = char(getData(theNode)); + text.(textflag) = toCharArray(getTextContent(theNode))'; + end + +end + +% ----- Subfunction parseAttributes ----- +function attributes = parseAttributes(theNode) + % Create attributes structure. + + attributes = struct; + if hasAttributes(theNode) + theAttributes = getAttributes(theNode); + numAttributes = getLength(theAttributes); + + for count = 1:numAttributes + %attrib = item(theAttributes,count-1); + %attr_name = regexprep(char(getName(attrib)),'[-:.]','_'); + %attributes.(attr_name) = char(getValue(attrib)); + + %Suggestion of Adrian Wanner + str = toCharArray(toString(item(theAttributes,count-1)))'; + k = strfind(str,'='); + attr_name = str(1:(k(1)-1)); + attr_name = strrep(attr_name, '-', '_dash_'); + attr_name = strrep(attr_name, ':', '_colon_'); + attr_name = strrep(attr_name, '.', '_dot_'); + attributes.(attr_name) = str((k(1)+2):(end-1)); + end + end +end + +% ----- Subfunction text2val ----- +function s_new = text2val(s) + % If element has "datatype" attribute, convert element.Text into + % numeric field with that data type. + s_new = []; + fields = fieldnames(s); + + for f = 1:numel(fields) + name = fields{f}; + val=s.(fields{f}); + valtype=[]; + if(isstruct(val)) + if(isfield(val,'Attributes') && isfield(val.Attributes,'datatype')) + valtype=val.Attributes.datatype; + end + + if(isfield(val,'Text')) + val = val.Text; + else + val = text2val(val); + end + end + + numtypes='single|double|int8|int16|int32|int64|uint8|uint16|uint32|uint64|logical'; + if(isempty(val)) + val = []; + elseif(~isempty(valtype) ... + && ~isempty(regexp(valtype,numtypes,'ONCE'))) + num = str2num(val); %#ok + if(~isempty(num)) + val=cast(num,valtype); + end + end + s_new.(name) = val; + end +end diff --git a/README.md b/README.md new file mode 100644 index 0000000..95171fe --- /dev/null +++ b/README.md @@ -0,0 +1,55 @@ +# HCPretinotopyanalysis + +History of major code changes: +- 2015/02/07 - Version 1.0. + +## CONTENTS + +Contents: +- ExampleInvocation.sh - example script for invoking the analysis +- MATLAB - This folder contains MATLAB source code + - analyzePRF_chpc.m - a version of analyzePRF.m specifically adapted for HCP use on CHPC + - CHPCsubmit.m - utility for submitting jobs to CHPC + - CHPCwait.m - utility for waiting for jobs to complete on CHPC + - RetinotopyAnalysis.m - entry point for the MATLAB processing + - utilities - various MATLAB utility functions, including MatlabCIFTI +- mcc - This folder contains various files related to the compiled MATLAB code and the CHPC setup + - compilescript.m - a script that compiles the necessary MATLAB functions + - fitnonlinearmodel and + run_fitnonlinearmodel.sh - compiled version of fitnonlinearmodel.m + - RetinotopyAnalysis and + run_RetinotopyAnalysis.sh - compiled version of RetinotopyAnalysis.m + - matlabsge.sh - a shell script used by CHPC jobs in order to invoke MATLAB +- README.md - this file +- RetinotopyAnalysis.sh - top-level analysis script + +There are several external dependencies: +1. Connectome Workbench. This is required by MatlabCIFTI. It is required that + the command-line utility "wb_command" is accessible from the shell. +2. analyzePRF (http://github.com/kendrickkay/analyzePRF/). This is the MATLAB toolbox + that actually analyzes the data. We use tagged version 1.2. If you wish to + compile the MATLAB functions yourself, you will need to download this toolbox. + +Note that the files in the MATLAB directory are not actually used during the +analysis. This is because the processing relies on compiled MATLAB code (as contained +in the mcc directory). + +Installation and setup issues: +1. The CHPC jobs write .o and .e files to ~/sgeoutput +2. It is assumed that the mcc directory is installed in the home directory (~/) +3. It is assumed that a dummy CIFTI file is available at ~/dummy.dtseries.nii +4. It is assumed that scratch space is available at /scratch// + +The code workflow is as follows: +1. RetinotopyAnalysis.sh + [this is the top-level script that is called by the user. + this script generates a .txt file that is really a MATLAB + script that defines some input variables.] +2. -> run_RetinotopyAnalysis.sh + [this is a compiled version of RetinotopyAnalysis.m.] + [this function accepts the location of the .txt file created in the previous step, + evaluates the file, prepares inputs, calls analyzePRF_chpc.m with these inputs, + takes the outputs, massages these outputs, and writes them to CIFTI files.] +3. -> analyzePRF_chpc.m + [this does the real work, and is an alternative version of analyzePRF.m.] + [note that this function encodes a lot of specific details about the CHPC setup.] diff --git a/RetinotopyAnalysis.sh b/RetinotopyAnalysis.sh index 30c3ef0..3583ac9 100755 --- a/RetinotopyAnalysis.sh +++ b/RetinotopyAnalysis.sh @@ -64,7 +64,7 @@ show_version() { usage() { local scriptName=$(basename ${0}) echo "" - echo " Perform Retinotophy Analysis ... TO BE WRITTEN" + echo " Perform Retinotopy Analysis ... TO BE WRITTEN" echo "" echo " Usage: ${scriptName} " echo "" @@ -80,7 +80,8 @@ usage() { echo " : id of subject for the data being processed" echo "" echo " --movie-files=" - echo " : @ symbol separated list of movie files (QuickTime files) used as stimuli for the retinotopy task" + echo " : @ symbol separated list of movie files used as stimuli for the retinotopy task" + echo " Expects Quicktime .mov files, Matlab .mat files, or .hdf5 echo "" echo " --image-files=" echo " : @ symbol separated list of minimally preprocessed functional \(fMRI\) image files" @@ -124,7 +125,7 @@ usage() { # Global output variables # ${userid} - input - user login id # ${subject} - input - subject id -# ${movie_files} - input - @ symbol separated list of movie files (QuickTime files) used as stimuli for the +# ${movie_files} - input - @ symbol separated list of movie files (mov,mat,or hdf5) used as stimuli for the # retinotopy task # ${image_files} - input - @ symbol separated list of minimally preprocessed functional (fMRI) image files # ${behavior_files} - input - @ symbol separated list of behavior files from which time offsets can be obtained @@ -238,7 +239,7 @@ main() { # Global Variables Set # ${userid} - input - user login id # ${subject} - input - subject id - # ${movie_files} - input - @ symbol separated list of movie files (QuickTime files) used as stimuli for the + # ${movie_files} - input - @ symbol separated list of movie files (mov,mat,or hdf5) used as stimuli for the # retinotopy task # ${image_files} - input - @ symbol separated list of minimally preprocessed functional (fMRI) image files # ${behavior_files} - input - @ symbol separated list of behavior files from which time offsets can be obtained @@ -341,7 +342,7 @@ EOF # Run a compiled Matlab script, passing it the variables file we just created export MCR_CACHE_ROOT=/tmp export MATLAB_HOME="/export/matlab/R2012b" - ./run_RetinotopyAnalysis.sh ${MATLAB_HOME}/MCR ${subject}_matlab_variables.txt + ~/mcc/run_RetinotopyAnalysis.sh ${MATLAB_HOME}/MCR ${subject}_matlab_variables.txt } # diff --git a/mcc/RetinotopyAnalysis b/mcc/RetinotopyAnalysis new file mode 100755 index 0000000..d8a4bab Binary files /dev/null and b/mcc/RetinotopyAnalysis differ diff --git a/mcc/compilescript.m b/mcc/compilescript.m new file mode 100755 index 0000000..c79b8dd --- /dev/null +++ b/mcc/compilescript.m @@ -0,0 +1,27 @@ +% make sure the startup.m file does not exist or has just a dummy line + +% ensure a vanilla path +restoredefaultpath; + +% add RetinotopySkeleton to path +addpath(genpath('path/to/RetinotopySkeleton')); + +% add analyzePRF version 1.2 to path +addpath(genpath('path/to/analyzePRF-1.2')); + +% compile fitnonlinearmodel.m into an executable +mcc -a makegaussian2d.m ... % extra functions in anonymous functions + -a loadbinary.m ... + -a subscript.m ... + -a loadmulti.m ... + -a conv2run.m ... + -a posrect.m ... + -a vflatten.m ... + -a placematrix.m ... + -a calccod.m ... + -a analyzePRF_chpc.m ... + -v -m -R -nojvm -R -nodisplay -R singleCompThread ... + fitnonlinearmodel.m; + +% compile RetinotopyAnalysis.m into an executable +mcc -v -m -R nojvm -R -nodisplay -R singleCompThread RetinotopyAnalysis.m; diff --git a/mcc/fitnonlinearmodel b/mcc/fitnonlinearmodel new file mode 100755 index 0000000..3e73fdf Binary files /dev/null and b/mcc/fitnonlinearmodel differ diff --git a/mcc/matlabsge.sh b/mcc/matlabsge.sh new file mode 100755 index 0000000..e68c0e3 --- /dev/null +++ b/mcc/matlabsge.sh @@ -0,0 +1,7 @@ +#!/bin/sh + +MCR_CACHE_ROOT=/tmp +export MCR_CACHE_ROOT + +echo Now issuing MATLAB command. +~/mcc/$MYSCRIPT /export/matlab/R2012b/MCR $MYARG1 $MYARG2 $PBS_ARRAYID diff --git a/mcc/run_RetinotopyAnalysis.sh b/mcc/run_RetinotopyAnalysis.sh new file mode 100755 index 0000000..cbc8ade --- /dev/null +++ b/mcc/run_RetinotopyAnalysis.sh @@ -0,0 +1,39 @@ +#!/bin/sh +# script for execution of deployed applications +# +# Sets up the MCR environment for the current $ARCH and executes +# the specified command. +# +exe_name=$0 +exe_dir=`dirname "$0"` +echo "------------------------------------------" +if [ "x$1" = "x" ]; then + echo Usage: + echo $0 \ args +else + echo Setting up environment variables + MCRROOT="$1" + echo --- + LD_LIBRARY_PATH=.:${MCRROOT}/runtime/glnxa64 ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/bin/glnxa64 ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/sys/os/glnxa64; + MCRJRE=${MCRROOT}/sys/java/jre/glnxa64/jre/lib/amd64 ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/native_threads ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/server ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/client ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE} ; + XAPPLRESDIR=${MCRROOT}/X11/app-defaults ; + export LD_LIBRARY_PATH; + export XAPPLRESDIR; + echo LD_LIBRARY_PATH is ${LD_LIBRARY_PATH}; + shift 1 + args= + while [ $# -gt 0 ]; do + token=$1 + args="${args} ${token}" + shift + done + "${exe_dir}"/RetinotopyAnalysis $args +fi +exit + diff --git a/mcc/run_fitnonlinearmodel.sh b/mcc/run_fitnonlinearmodel.sh new file mode 100755 index 0000000..3d4b04e --- /dev/null +++ b/mcc/run_fitnonlinearmodel.sh @@ -0,0 +1,39 @@ +#!/bin/sh +# script for execution of deployed applications +# +# Sets up the MCR environment for the current $ARCH and executes +# the specified command. +# +exe_name=$0 +exe_dir=`dirname "$0"` +echo "------------------------------------------" +if [ "x$1" = "x" ]; then + echo Usage: + echo $0 \ args +else + echo Setting up environment variables + MCRROOT="$1" + echo --- + LD_LIBRARY_PATH=.:${MCRROOT}/runtime/glnxa64 ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/bin/glnxa64 ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/sys/os/glnxa64; + MCRJRE=${MCRROOT}/sys/java/jre/glnxa64/jre/lib/amd64 ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/native_threads ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/server ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/client ; + LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE} ; + XAPPLRESDIR=${MCRROOT}/X11/app-defaults ; + export LD_LIBRARY_PATH; + export XAPPLRESDIR; + echo LD_LIBRARY_PATH is ${LD_LIBRARY_PATH}; + shift 1 + args= + while [ $# -gt 0 ]; do + token=$1 + args="${args} ${token}" + shift + done + "${exe_dir}"/fitnonlinearmodel $args +fi +exit + diff --git a/run_RetinotopyAnalysis.sh b/run_RetinotopyAnalysis.sh deleted file mode 100755 index fc59f46..0000000 --- a/run_RetinotopyAnalysis.sh +++ /dev/null @@ -1,9 +0,0 @@ -#!/bin/bash - -SCRIPT_NAME=$0 -MCR_DIR=$1 -VAR_FILE=$2 - -echo "SCRIPT_NAME: ${SCRIPT_NAME}" -echo "MCR_DIR: ${MCR_DIR}" -echo "VAR_FILE: ${VAR_FILE}" \ No newline at end of file