diff --git a/ExptWrapper_Display.m b/ExptWrapper_Display.m new file mode 100755 index 0000000..05d8194 --- /dev/null +++ b/ExptWrapper_Display.m @@ -0,0 +1,88 @@ +%exptDir = '~/code/rtAttenPenn/'; +%cd(exptDir) +subjectNum = 500; +projectName = 'rtAttenPenn'; +Screen('Preference', 'SkipSyncTests', 1); +subjectRun = 1; +% **** types of stimuli to train/show to subjects ******* +NEUTRAL = 1; +SAD = 2; +HAPPY = 3; +% ******************************************************* +typeNum = SAD; + +subjectName = [datestr(now,5) datestr(now,7) datestr(now,11) num2str(subjectRun) '_' projectName]; +matchNum = 0; +useButtonBox=1; +realtimeData = 0; +debug=1; +KbName('UnifyKeyNames') +addpath(genpath('/opt/psychtoolbox/')) + +% well the git is set to ignore the data folder so it would just make sense +% to save here +% +%% +runNum=2; +fMRI = 0; +% today's scanning number +realtimeData = 0 +debug = 0 +useButtonBox = 0 +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) +runNum=1; +fMRI = 6; +fMRI = 0; + + +%% +runNum=6; +fMRI = 16; +useButtonBox = 1 +realtimeData = 1 +debug = 1 +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) + + +%% +runNum = 3; +fMRI = 14; +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) + + +%% +runNum = 4; +fMRI = 16; +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) + + +%% +runNum = 5; +fMRI = 18; +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) + +%% +runNum = 6; +fMRI = 20; +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) + + +%% + +runNum = 7; +fMRI = 20; +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) + + +%% + +runNum = 8; +fMRI = 22; +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) + + +%% + +runNum = 9; +fMRI = 24; +[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) diff --git a/ExptWrapper_Display_Subject1.m b/ExptWrapper_Display_Subject1.m deleted file mode 100755 index d704536..0000000 --- a/ExptWrapper_Display_Subject1.m +++ /dev/null @@ -1,85 +0,0 @@ -%exptDir = '~/code/rtAttenPenn/'; -%cd(exptDir) -subjectNum = 3; -projectName = 'rtAttenPenn'; -Screen('Preference', 'SkipSyncTests', 1); -subjectRun = 1; -% **** types of stimuli to train/show to subjects ******* -NEUTRAL = 1; -SAD = 2; -HAPPY = 3; -% ******************************************************* -typeNum = SAD; - -subjectName = [datestr(now,5) datestr(now,7) datestr(now,11) num2str(subjectRun) '_' projectName]; -matchNum = 0; -useButtonBox=1; -realtimeData = 1; -debug=1; -KbName('UnifyKeyNames') -addpath(genpath('/opt/psychtoolbox/')) -dataDirHeader = '/Volumes/norman/amennen/'; -%% -runNum=2; -fMRI = 0; -% today's scanning number -realtimeData = 0 -debug = 0 -useButtonBox = 0 -[blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) -runNum=1; -fMRI = 6; -fMRI = 0; - - -%% -runNum=6; -fMRI = 16; -useButtonBox = 1 -realtimeData = 1 -debug = 1 -[blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) - - -%% -runNum = 3; -fMRI = 14; -[blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) - - -%% -runNum = 4; -fMRI = 16; -[blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) - - -%% -runNum = 5; -fMRI = 18; -[blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) - -%% -runNum = 6; -fMRI = 20; -[blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) - - -%% - -runNum = 7; -fMRI = 20; -[blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) - - -%% - -runNum = 8; -fMRI = 22; -[blockData] = RealTimePunisherDisplay(dataDirHeader,ubjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) - - -%% - -runNum = 9; -fMRI = 24; -[blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,realtimeData,debug) diff --git a/ExptWrapper_FileProcess_Subject1.m b/ExptWrapper_FileProcess.m similarity index 98% rename from ExptWrapper_FileProcess_Subject1.m rename to ExptWrapper_FileProcess.m index 1c8c2e0..ef25d73 100755 --- a/ExptWrapper_FileProcess_Subject1.m +++ b/ExptWrapper_FileProcess.m @@ -9,7 +9,7 @@ %% %conf = loadjson('conf/example.json'); -subjectNum = 200; +subjectNum = 500; projectName = 'rtAttenPenn'; %imgDirHeader = conf.imgDirHeader; subjectRun = 1; @@ -29,7 +29,7 @@ KbName('UnifyKeyNames') addpath(genpath('/opt/psychtoolbox/')) -imgDirHeader = '/Data1/subjects/' +imgDirHeader = '/Data1/subjects/'; %% Generate expt sequence if matchNum ==0 diff --git a/ProcessMask.m b/ProcessMask.m index 1053597..a5df672 100644 --- a/ProcessMask.m +++ b/ProcessMask.m @@ -30,15 +30,15 @@ end setenv('FSLOUTPUTTYPE','NIFTI_GZ'); - +base_path = pwd; if matchNum == 0 - save_dir = ['./data/' num2str(subjectNum)]; + save_dir = fullfile(pwd,['data/' num2str(subjectNum)]); else - save_dir = ['./data/' num2str(matchNum) '_match']; + save_dir = fullfile(pwd,['data/' num2str(matchNum) '_match']); end subjectName = [datestr(now,5) datestr(now,7) datestr(now,11) num2str(runNum) '_' projectName]; dicom_dir = ['/mnt/rtexport/RTexport_Current/' datestr(subjDate,10) datestr(subjDate,5) datestr(subjDate,7) '.' subjectName '.' subjectName '/']; -process_dir = [save_dir 'reg' '/']; +process_dir = fullfile(save_dir,[ 'reg' '/']); roi_name = 'wholebrain_mask'; roi_dir = pwd; % change this path name to wherever you put it on the penn computer! code_dir = pwd; @@ -54,7 +54,7 @@ highresfiles_genstr = sprintf('%s001_0000%s_0*',dicom_dir,num2str(highresScan,'%2.2i')); %general string for ALL mprage files** unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,highresfiles_genstr,highresFN)); unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,highresFN,highresFN_RE)); -unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --2niftihdr -s %s.bxh %s',bxhpath,highresFN_RE,highresFN_RE)) +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,highresFN_RE,highresFN_RE)) unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -R -m',fslpath,highresFN_RE,highresFN_RE)) %unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -r 90 -R',fslpath,highresFN_RE,highresFN_RE)) @@ -173,9 +173,9 @@ mask=mask_brain; %save anatomical mask if matchNum == 0 - save(['./data/' num2str(subjectNum) '/mask_' num2str(subjectNum)],'mask'); + save(fullfile(save_dir, [ '/mask_' num2str(subjectNum)]),'mask'); else - save(['./data/' num2str(matchNum) '_match/mask_' num2str(subjectNum)],'mask'); + save(fullfile(save_dir, [ '/mask_' num2str(subjectNum)]),'mask'); end fprintf('Done with mask creation\n'); diff --git a/ProcessMultiday.m b/ProcessMultiday.m new file mode 100644 index 0000000..b519d21 --- /dev/null +++ b/ProcessMultiday.m @@ -0,0 +1,264 @@ +% ProcessMask +% Written by ACM April 2017 + +% Collect: +% - t1-weighted mp-rage +% - example functional scan +% Complete: +% 1. Register t1 to standard space +% 3. Register epi example to t1, including field map corrections +% 4. Calculate inverse transformation matrices +% 5. Apply to anatomical mask +subjectNum = 5; % multiday test is subject 5, intel demo is subject 3 +matchNum = 0; +runNum = 1; +projectName = 'rtAttenPenn'; +highresScan = 5; +functionalScan=6; +subjDate = '9-22-17'; +%% +startProcess = GetSecs; +% this is for Penn +%biac_dir = '/home/amennen/code/BIAC_Matlab_R2014a/'; +biac_dir = '/Data1/packages/BIAC_Matlab_R2014a/'; + +bxhpath='/opt/BXH/1.11.1/bin/'; +fslpath='/opt/fsl/5.0.9/bin/'; +%add necessary package +if ~exist('readmr','file') + addpath(genpath(biac_dir)); + addpath([biac_dir '/mr/']); + addpath([biac_dir '/general/']) +end + +setenv('FSLOUTPUTTYPE','NIFTI_GZ'); +project_folder = '/Data1/code/rtAttenPenn'; +if matchNum == 0 + save_dir = [project_folder '/data/' num2str(subjectNum)]; + %save(['./data/' num2str(subjectNum) '/mask_' num2str(subjectNum)],'mask'); +else + save_dir = [project_folder '/data/' num2str(matchNum) '_match']; +end + +%subjectName = [datestr(now,5) datestr(now,7) datestr(now,11) num2str(runNum) '_' projectName]; +subjectName = [datestr(subjDate,5) datestr(subjDate,7) datestr(subjDate,11) num2str(runNum) '_' projectName]; + +%dicom_dir = ['/mnt/rtexport/RTexport_Current/' datestr(subjDate,10) datestr(subjDate,5) datestr(subjDate,7) '.' subjectName '.' subjectName '/']; +dicom_dir = ['/Data1/subjects/' datestr(subjDate,10) datestr(subjDate,5) datestr(subjDate,7) '.' subjectName '.' subjectName '/']; + +process_dir = [save_dir '/' 'reg' '/']; +roi_name = 'wholebrain_mask'; +roi_dir = pwd; % change this path name to wherever you put it on the penn computer! +code_dir = pwd; +addpath(genpath(code_dir)); + +if ~exist(process_dir) + mkdir(process_dir) +end +cd(process_dir) +%% Process t1-weighted MPRAGE and check brain extraction! +highresFN = 'highres'; +highresFN_RE = 'highres_re'; +highresfiles_genstr = sprintf('%s001_0000%s_0*',dicom_dir,num2str(highresScan,'%2.2i')); %general string for ALL mprage files** +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,highresfiles_genstr,highresFN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,highresFN,highresFN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,highresFN_RE,highresFN_RE)) +unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -R -m',fslpath,highresFN_RE,highresFN_RE)) +%unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -r 90 -R',fslpath,highresFN_RE,highresFN_RE)) + +% for dcm2niix the command would be 'dcm2niix dicomdir -f test -o dicomdir -s y dicomdir/001_000007_000008.dcm' +fprintf('%sfslview %s.nii.gz\n',fslpath,highresFN_RE) +fprintf('%sfslview %s_brain.nii.gz', fslpath,highresFN_RE) +%% Register standard to nifti +StartReg = GetSecs; +unix(sprintf('%sflirt -in %s_brain.nii.gz -ref $FSLDIR/data/standard/MNI152_T1_2mm_brain.nii.gz -out highres2standard -omat highres2standard.mat -cost corratio -dof 12 -searchrx -30 30 -searchry -30 30 -searchrz -30 30 -interp trilinear',fslpath,highresFN_RE)); +unix(sprintf('%sfnirt --iout=highres2standard_head --in=%s.nii.gz --aff=highres2standard.mat --cout=highres2standard_warp --iout=highres2standard --jout=highres2highres_jac --config=T1_2_MNI152_2mm --ref=$FSLDIR/data/standard/MNI152_T1_2mm.nii.gz --refmask=$FSLDIR/data/standard/MNI152_T1_2mm_brain_mask_dil --warpres=10,10,10', fslpath,highresFN_RE)); +unix(sprintf('%sapplywarp -i %s_brain.nii.gz -r $FSLDIR/data/standard/MNI152_T1_2mm_brain.nii.gz -o highres2standard -w highres2standard_warp',fslpath,highresFN_RE)); +%compute inverse transform (standard to highres) +unix(sprintf('%sconvert_xfm -inverse -omat standard2highres.mat highres2standard.mat', fslpath)); +unix(sprintf('%sinvwarp -w highres2standard_warp -o standard2highres_warp -r %s_brain.nii.gz',fslpath,highresFN_RE)); +t.standard2highres = GetSecs - StartReg; +fprintf('Done with standard2highres registration. Time = %6.2f \n',t.standard2highres); + +%% Process example epi file +startFunctional = GetSecs; + +fileN = 6; % we can choose 10 later2 +functionalFN = 'exfunc'; +functionalFN_RE = 'exfunc_re'; +exfunc_str = sprintf('%s001_0000%s_0000%s.dcm',dicom_dir,num2str(functionalScan,'%2.2i'),num2str(fileN,'%2.2i')); %general string for ALL mprage files** +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,exfunc_str,functionalFN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,functionalFN,functionalFN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,functionalFN_RE,functionalFN_RE)) + +% now register to highres! +t1 = GetSecs; +exfunc2highres_mat='example_func2highres'; +highres2exfunc_mat='highres2example_func'; +unix(sprintf('%sepi_reg --epi=%s --t1=%s --t1brain=%s_brain --out=%s',fslpath,functionalFN_RE,highresFN_RE,highresFN_RE,exfunc2highres_mat)) +timefunc2highres = GetSecs-t1; +unix(sprintf('%sconvert_xfm -inverse -omat %s.mat %s.mat',fslpath,highres2exfunc_mat,exfunc2highres_mat)); + +% now register mask to all data +unix(sprintf('%sapplywarp -i %s/%s.nii.gz -r %s.nii.gz -o %s_exfunc.nii.gz -w standard2highres_warp.nii.gz --postmat=%s.mat',fslpath,roi_dir,roi_name,functionalFN_RE,roi_name,highres2exfunc_mat)); +% check after here that the applied warp is binary and in the right +% orientation so we could just apply to nifti files afterwards +if exist(sprintf('%s_exfunc.nii.gz',roi_name),'file') + unix(sprintf('gunzip %s_exfunc.nii.gz',roi_name)); +end +unix(sprintf('%sbxhabsorb %s_exfunc.nii %s_exfunc.bxh',bxhpath,roi_name,roi_name)); + +% brain extract functional scan to make sure we stay inside the brain of +% the subject! +unix(sprintf('%sbet %s.nii.gz %s_brain -R -m',fslpath,functionalFN_RE,functionalFN_RE)); % check that this is okay! +%CHECK OKAY +fprintf('%sfslview %s_brain_mask.nii.gz', fslpath,functionalFN_RE) + +t.standard2func = GetSecs - startFunctional; +fprintf('Done with standard2func registration, time = %6.2f', t.standard2func); +%% if okay then unzip and make bxh wrapper +if exist(sprintf('%s_brain.nii.gz',functionalFN_RE),'file') + unix(sprintf('gunzip %s_brain.nii.gz',functionalFN_RE)); +end +unix(sprintf('%sbxhabsorb %s_brain.nii %s_brain.bxh',bxhpath,functionalFN_RE,functionalFN_RE)); +%% create mask file for real-time dicom files +% this is where you would make the mask bigger +vol = ReadFile(exfunc_str,64,0); +mask = BrainMask(vol,0,0); +imagesc(mask(:,:,10)); %checks that the dicom files are properly aligned +save(['mask_wholeBrain' '.mat'], 'mask'); + + +startMask = GetSecs; +%load registered anatomical ROI +maskStruct = readmr([roi_name '_exfunc.bxh'],'BXH',{[],[],[]}); +brainExtFunc = readmr([functionalFN_RE '_brain.bxh'], 'BXH',{[],[],[]}); + +%rotate anatomical ROI to be in the same space as the mask - check that this works for different scans/ROIs +anatMaskRot = zeros(size(mask)); +brainExtRot = zeros(size(mask)); +for i = 1:size(maskStruct.data,3) + anatMaskRot(:,:,i) = rot90(maskStruct.data(:,:,i)); %rotates entire slice by 90 degrees + brainExtRot(:,:,i) = rot90(brainExtFunc.data(:,:,i)); +end + +%overwrite whole-brain mask +mask = logical(anatMaskRot); %make it so it's just 1's and 0's +brainExt = logical(brainExtRot); +allinMask = find(anatMaskRot); +allinBrainExt = find(brainExt); +mask_indices = allinMask(find(ismember(allinMask,allinBrainExt))); %these are the good mask indices that are only brain +[gX gY gZ] = ind2sub(size(mask),mask_indices); +mask_brain = zeros(size(mask,1),size(mask,2),size(mask,3)); +for j=1:length(mask_indices) + mask_brain(gX(j),gY(j),gZ(j)) = 1; +end + +checkMask = 1; +if checkMask + plot3Dbrain(mask,[], 'mask') + plot3Dbrain(mask_brain, [], 'mask_brain') +end +% use the mask that has been checked there's nothing outside the functional +% brain +mask=mask_brain; +%save anatomical mask +if matchNum == 0 + save([code_dir '/data/' num2str(subjectNum) '/mask_' num2str(subjectNum)],'mask'); +else + save([code_dir '/data/' num2str(matchNum) '_match/mask_' num2str(subjectNum)],'mask'); +end + +fprintf('Done with mask creation\n'); +t.mask = GetSecs - startMask; +t.total = GetSecs - startProcess; +save(fullfile(process_dir, 'timing'), 't'); +fprintf('Standard2highres time = %7.2f \nStandard2func time = %7.2f \nMask time = %7.2f \n Total time = %7.2f\n', t.standard2highres, t.standard2func,t.mask,t.total); +% if cd into the directory, cd out of it back to the general exp folder +cd(code_dir) + + +%% NOW WE REGISTERED ONE LETS USE +% register flash day 1 to examplefunc day 2 + +flash1Scan = 8; +flash1FN = 'flash'; +flash1FN_RE = 'flash1_re'; +flash1files_genstr = sprintf('%s001_0000%s_0*',dicom_dir,num2str(flash1Scan,'%2.2i')); %general string for ALL mprage files** +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,flash1files_genstr,flash1FN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,flash1FN,flash1FN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,flash1FN_RE,flash1FN_RE)) +unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -R -m',fslpath,flash1FN_RE,flash1FN_RE)) +%unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -r 90 -R',fslpath,highresFN_RE,highresFN_RE)) + + + + +%% TOMORROW- - REDO DAY 2 BECAUSE THE SUBJECT PATH WAS NOT UPDATED! +subjectNum = 5; % multiday test is subject 5, intel demo is subject 3 +matchNum = 0; +runNum = 3; +projectName = 'rtAttenPenn'; +highresScan = 5; +functionalScan=6; +subjDate = '9-22-17'; +subjectName = [datestr(subjDate,5) datestr(subjDate,7) datestr(subjDate,11) num2str(runNum) '_' projectName]; +dicom_dir2 = ['/Data1/subjects/' datestr(subjDate,10) datestr(subjDate,5) datestr(subjDate,7) '.' subjectName '.' subjectName '/']; + +flash2Scan = 8; +flash2FN = 'flash2'; +flash2FN_RE = 'flash2_re'; +flash2files_genstr = sprintf('%s001_0000%s_0*',dicom_dir2,num2str(flash2Scan,'%2.2i')); %general string for ALL mprage files** +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,flash2files_genstr,flash2FN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,flash2FN,flash2FN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,flash2FN_RE,flash2FN_RE)) +unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -R -m',fslpath,flash2FN_RE,flash2FN_RE)) +%% now register flash 1 to flash 2 +unix(sprintf('%sflirt -dof 6 -in %s_brain.nii.gz -ref %s_brain.nii.gz -out flash12flash2 -omat flash12flash2.mat', fslpath, flash1FN_RE,flash2FN_RE)) + +% now apply mask to flash 2 +unix(sprintf('%sflirt -in %s_exfunc.nii -ref %s_brain.nii.gz -out %s_day2 -init flash12flash2.mat -applyxfm', fslpath, roi_name,flash2FN_RE,roi_name)) +%% Mark's suggestion - 2 flirt commands + + +subjectNum = 5; % multiday test is subject 5, intel demo is subject 3 +matchNum = 0; +projectName = 'rtAttenPenn'; +highresScan = 5; +functionalScan=6; +flash_hrScan = 7; +flash_lrScan = 8; + +subjDate = '9-22-17'; +runNum = 1; +subjectName1 = [datestr(subjDate,5) datestr(subjDate,7) datestr(subjDate,11) num2str(runNum) '_' projectName]; +dicom_dir1 = ['/Data1/subjects/' datestr(subjDate,10) datestr(subjDate,5) datestr(subjDate,7) '.' subjectName '.' subjectName '/']; + +runNum = 3; +subjectName2 = [datestr(subjDate,5) datestr(subjDate,7) datestr(subjDate,11) num2str(runNum) '_' projectName]; +dicom_dir2 = ['/Data1/subjects/' datestr(subjDate,10) datestr(subjDate,5) datestr(subjDate,7) '.' subjectName '.' subjectName '/']; + + +% repeat over with flash scans being scan numbers 7 AND OR 8!!! + +% new test: (repeat for each type of flash that is recorded) + +% 1. register exfunc1 --> flash2 +unix(sprintf('%sflirt -dof 6 -in %s_brain.nii.gz -ref %s_brain.nii.gz -out func12flash -omat func12flash.mat', fslpath, functionalFN_RE,flash2FN_RE)) + +% 2. register flash 2 --> exfunc 1 +unix(sprintf('%sflirt -dof 6 -in %s_brain.nii.gz -ref %s_brain.nii.gz -out flash2func2 -omat flash2func2.mat', fslpath, flash2FN_RE,functional2FN_RE)) + +% 3. apply old mask - apply mask to func2flash1 +unix(sprintf('%sflirt -in %MADK1NAME.nii.gz -ref %s_brain.nii.gz -applyxfm -init func2flash.mat -interp nearestneighbor -out mask1-2-flash', fslpath, flash2FN_RE)) + +% 4: apply odl mask - apply mask to flash2func2 +unix(sprintf('%sflirt -in mask-1-2-flash -ref %s_brain.nii.gz -applyxfm -init flash2func2.mat -interp nearestneighbor -out mask1-2-flsah-2-func2', fslpath, functional2FN_RE)) + +% now you have a mask in func 2 space! + + + + +'%flirt -in mask1 -ref t1_flash_brain -applyxfm -init func2flash1.mat -interp nearestneighbor -out mask1-2-flash' +'%flirt -in mask1-2-flash -ref exfunc2 -applyxfm -init flash2func2.mat -interp nearestneighbor -out mask1-2-flash2-func2' \ No newline at end of file diff --git a/ProcessMultiday2.m b/ProcessMultiday2.m new file mode 100644 index 0000000..a7d684a --- /dev/null +++ b/ProcessMultiday2.m @@ -0,0 +1,212 @@ +% declare things that don't change across days + +subjectNum = 5; % multiday test is subject 5, intel demo is subject 3 +matchNum = 0; +projectName = 'rtAttenPenn'; +highresScan = 5; +functionalScan=6; +flash_hrScan = 7; +flash_lrScan = 8; + +biac_dir = '/Data1/packages/BIAC_Matlab_R2014a/'; +bxhpath='/opt/BXH/1.11.1/bin/'; +fslpath='/opt/fsl/5.0.9/bin/'; +%add necessary package +if ~exist('readmr','file') + addpath(genpath(biac_dir)); + addpath([biac_dir '/mr/']); + addpath([biac_dir '/general/']) +end +setenv('FSLOUTPUTTYPE','NIFTI_GZ'); +project_folder = '/Data1/code/rtAttenPenn'; +if matchNum == 0 + save_dir = [project_folder '/data/' num2str(subjectNum)]; +else + save_dir = [project_folder '/data/' num2str(matchNum) '_match']; +end +process_dir = [save_dir '/' 'reg' '/']; +roi_name = 'wholebrain_mask'; +roi_dir = pwd; % change this path name to wherever you put it on the penn computer! +code_dir = pwd; +addpath(genpath(code_dir)); + +if ~exist(process_dir) + mkdir(process_dir) +end +cd(process_dir) +%% +subjDate1 = '9-22-17'; +runNum = 1; +subjectName1 = [datestr(subjDate1,5) datestr(subjDate1,7) datestr(subjDate1,11) num2str(runNum) '_' projectName]; +dicom_dir1 = ['/Data1/subjects/' datestr(subjDate1,10) datestr(subjDate1,5) datestr(subjDate1,7) '.' subjectName1 '.' subjectName1 '/']; + +%% now do steps registration for first day scans to make the mask + +%% Process t1-weighted MPRAGE and check brain extraction! +highresFN = 'highres'; +highresFN_RE = 'highres_re'; +highresfiles_genstr = sprintf('%s001_0000%s_0*',dicom_dir,num2str(highresScan,'%2.2i')); %general string for ALL mprage files** +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,highresfiles_genstr,highresFN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,highresFN,highresFN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,highresFN_RE,highresFN_RE)) +unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -R -m',fslpath,highresFN_RE,highresFN_RE)) +%unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -r 90 -R',fslpath,highresFN_RE,highresFN_RE)) + +% for dcm2niix the command would be 'dcm2niix dicomdir -f test -o dicomdir -s y dicomdir/001_000007_000008.dcm' +fprintf('%sfslview %s.nii.gz\n',fslpath,highresFN_RE) +fprintf('%sfslview %s_brain.nii.gz', fslpath,highresFN_RE) +% Register standard to nifti +StartReg = GetSecs; +unix(sprintf('%sflirt -in %s_brain.nii.gz -ref $FSLDIR/data/standard/MNI152_T1_2mm_brain.nii.gz -out highres2standard -omat highres2standard.mat -cost corratio -dof 12 -searchrx -30 30 -searchry -30 30 -searchrz -30 30 -interp trilinear',fslpath,highresFN_RE)); +unix(sprintf('%sfnirt --iout=highres2standard_head --in=%s.nii.gz --aff=highres2standard.mat --cout=highres2standard_warp --iout=highres2standard --jout=highres2highres_jac --config=T1_2_MNI152_2mm --ref=$FSLDIR/data/standard/MNI152_T1_2mm.nii.gz --refmask=$FSLDIR/data/standard/MNI152_T1_2mm_brain_mask_dil --warpres=10,10,10', fslpath,highresFN_RE)); +unix(sprintf('%sapplywarp -i %s_brain.nii.gz -r $FSLDIR/data/standard/MNI152_T1_2mm_brain.nii.gz -o highres2standard -w highres2standard_warp',fslpath,highresFN_RE)); +%compute inverse transform (standard to highres) +unix(sprintf('%sconvert_xfm -inverse -omat standard2highres.mat highres2standard.mat', fslpath)); +unix(sprintf('%sinvwarp -w highres2standard_warp -o standard2highres_warp -r %s_brain.nii.gz',fslpath,highresFN_RE)); +t.standard2highres = GetSecs - StartReg; +fprintf('Done with standard2highres registration. Time = %6.2f \n',t.standard2highres); + +%% Process example epi file + +startFunctional = GetSecs; + +fileN = 6; % we can choose 10 later2 +functionalFN = 'exfunc'; +functionalFN_RE = 'exfunc_re'; +exfunc_str = sprintf('%s001_0000%s_0000%s.dcm',dicom_dir,num2str(functionalScan,'%2.2i'),num2str(fileN,'%2.2i')); %general string for ALL mprage files** +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,exfunc_str,functionalFN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,functionalFN,functionalFN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,functionalFN_RE,functionalFN_RE)) + +% now register to highres! +t1 = GetSecs; +exfunc2highres_mat='example_func2highres'; +highres2exfunc_mat='highres2example_func'; +unix(sprintf('%sepi_reg --epi=%s --t1=%s --t1brain=%s_brain --out=%s',fslpath,functionalFN_RE,highresFN_RE,highresFN_RE,exfunc2highres_mat)) +timefunc2highres = GetSecs-t1; +unix(sprintf('%sconvert_xfm -inverse -omat %s.mat %s.mat',fslpath,highres2exfunc_mat,exfunc2highres_mat)); + +% now register mask to all data +unix(sprintf('%sapplywarp -i %s/%s.nii.gz -r %s.nii.gz -o %s_exfunc.nii.gz -w standard2highres_warp.nii.gz --postmat=%s.mat',fslpath,roi_dir,roi_name,functionalFN_RE,roi_name,highres2exfunc_mat)); +% check after here that the applied warp is binary and in the right +% orientation so we could just apply to nifti files afterwards +if exist(sprintf('%s_exfunc.nii.gz',roi_name),'file') + unix(sprintf('gunzip %s_exfunc.nii.gz',roi_name)); +end +unix(sprintf('%sbxhabsorb %s_exfunc.nii %s_exfunc.bxh',bxhpath,roi_name,roi_name)); + +% brain extract functional scan to make sure we stay inside the brain of +% the subject! +unix(sprintf('%sbet %s.nii.gz %s_brain -R -m',fslpath,functionalFN_RE,functionalFN_RE)); % check that this is okay! +%CHECK OKAY +fprintf('%sfslview %s_brain_mask.nii.gz', fslpath,functionalFN_RE) + +t.standard2func = GetSecs - startFunctional; +fprintf('Done with standard2func registration, time = %6.2f', t.standard2func); + +%% create mask file for real-time dicom files +% this is where you would make the mask bigger +vol = ReadFile(exfunc_str,64,0); +mask = BrainMask(vol,0,0); +imagesc(mask(:,:,10)); %checks that the dicom files are properly aligned +save(['mask_wholeBrain' '.mat'], 'mask'); + + +startMask = GetSecs; +%load registered anatomical ROI +maskStruct = readmr([roi_name '_exfunc.bxh'],'BXH',{[],[],[]}); +brainExtFunc = readmr([functionalFN_RE '_brain.bxh'], 'BXH',{[],[],[]}); + +%rotate anatomical ROI to be in the same space as the mask - check that this works for different scans/ROIs +anatMaskRot = zeros(size(mask)); +brainExtRot = zeros(size(mask)); +for i = 1:size(maskStruct.data,3) + anatMaskRot(:,:,i) = rot90(maskStruct.data(:,:,i)); %rotates entire slice by 90 degrees + brainExtRot(:,:,i) = rot90(brainExtFunc.data(:,:,i)); +end + +%overwrite whole-brain mask +mask = logical(anatMaskRot); %make it so it's just 1's and 0's +brainExt = logical(brainExtRot); +allinMask = find(anatMaskRot); +allinBrainExt = find(brainExt); +mask_indices = allinMask(find(ismember(allinMask,allinBrainExt))); %these are the good mask indices that are only brain +[gX gY gZ] = ind2sub(size(mask),mask_indices); +mask_brain = zeros(size(mask,1),size(mask,2),size(mask,3)); +for j=1:length(mask_indices) + mask_brain(gX(j),gY(j),gZ(j)) = 1; +end + +checkMask = 1; +if checkMask + plot3Dbrain(mask,[], 'mask') + plot3Dbrain(mask_brain, [], 'mask_brain') +end +% use the mask that has been checked there's nothing outside the functional +% brain +mask=mask_brain; +%save anatomical mask +if matchNum == 0 + save([code_dir '/data/' num2str(subjectNum) '/mask_' num2str(subjectNum)],'mask'); +else + save([code_dir '/data/' num2str(matchNum) '_match/mask_' num2str(subjectNum)],'mask'); +end + +fprintf('Done with mask creation\n'); +t.mask = GetSecs - startMask; +t.total = GetSecs - startProcess; +save(fullfile(process_dir, 'timing'), 't'); +fprintf('Standard2highres time = %7.2f \nStandard2func time = %7.2f \nMask time = %7.2f \n Total time = %7.2f\n', t.standard2highres, t.standard2func,t.mask,t.total); +% if cd into the directory, cd out of it back to the general exp folder +cd(code_dir) + + +%% now register for second day mask +subjDate2 = '9-22-17'; +runNum = 3; +subjectName2 = [datestr(subjDate2,5) datestr(subjDate2,7) datestr(subjDate2,11) num2str(runNum) '_' projectName]; +dicom_dir2 = ['/Data1/subjects/' datestr(subjDate2,10) datestr(subjDate2,5) datestr(subjDate2,7) '.' subjectName2 '.' subjectName2 '/']; + +% get both flashes ready + +flashhrFN = 'flashhr'; +flashhrFN_RE = 'flashhr_re'; +flashhrfiles_genstr = sprintf('%s001_0000%s_0*',dicom_dir,num2str(flash_hrScan,'%2.2i')); +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,flashhrfiles_genstr,flashhrFN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,flashhrFN,flashhrFN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,flashhrFN_RE,flashhrFN_RE)) +unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -R -m',fslpath,flashhrFN_RE,flashhrFN_RE)) + +flashlrFN = 'flashlr'; +flashlrFN_RE = 'flashlr_re'; +flashlrfiles_genstr = sprintf('%s001_0000%s_0*',dicom_dir,num2str(flash_lrScan,'%2.2i')); +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,flashlrfiles_genstr,flashlrFN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,flashlrFN,flashlrFN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,flashlrFN_RE,flashlrFN_RE)) +unix(sprintf('%sbet %s.nii.gz %s_brain.nii.gz -R -m',fslpath,flashlrFN_RE,flashlrFN_RE)) + +% get epi 2 ready + +fileN = 6; % we can choose 10 later2 +functional2FN = 'exfunc2'; +functional2FN_RE = 'exfunc2_re'; +exfunc_str = sprintf('%s001_0000%s_0000%s.dcm',dicom_dir,num2str(functionalScan,'%2.2i'),num2str(fileN,'%2.2i')); %general string for ALL mprage files** +unix(sprintf('%sdicom2bxh %s %s.bxh',bxhpath,exfunc_str,functional2FN)); +unix(sprintf('%sbxhreorient --orientation=LAS %s.bxh %s.bxh',bxhpath,functional2FN,functional2FN_RE)); +unix(sprintf('%sbxh2analyze --overwrite --analyzetypes --niigz --niftihdr -s %s.bxh %s',bxhpath,functional2FN_RE,functional2FN_RE)) + + +% repeat over with flash scans being scan numbers 7 AND OR 8!!! + +% new test: (repeat for each type of flash that is recorded) + +% 1. register exfunc1 --> flash2 +unix(sprintf('%sflirt -dof 6 -in %s_brain.nii.gz -ref %s_brain.nii.gz -out func12flash -omat func12flash.mat', fslpath, functionalFN_RE,flashhrFN_RE)) +% 2. register flash 2 --> exfunc 2 +unix(sprintf('%sflirt -dof 6 -in %s_brain.nii.gz -ref %s_brain.nii.gz -out flash2func2 -omat flash2func2.mat', fslpath, flashhrFN_RE,functional2FN_RE)) +% 3. apply old mask - apply mask to func2flash1 +unix(sprintf('%sflirt -in wholebrain_exfunc.nii.gz -ref %s_brain.nii.gz -applyxfm -init func2flash.mat -interp nearestneighbor -out mask1-2-flash', fslpath, flash2FN_RE)) +% 4: apply odl mask - apply mask to flash2func2 +unix(sprintf('%sflirt -in mask1-2-flash -ref %s_brain.nii.gz -applyxfm -init flash2func2.mat -interp nearestneighbor -out mask1-2-flash-2-func2', fslpath, functional2FN_RE)) + +% now you have a mask in func 2 space! \ No newline at end of file diff --git a/RealTimePunisherDisplay.m b/RealTimePunisherDisplay.m index 8bca097..cf2c744 100755 --- a/RealTimePunisherDisplay.m +++ b/RealTimePunisherDisplay.m @@ -1,4 +1,4 @@ -function [blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,rtData,debug) +function [blockData] = RealTimePunisherDisplay(subjectNum,subjectName,matchNum,runNum,useButtonBox,fMRI,rtData,debug) % function [blockData] = RealTimePunisherDisplay(dataDirHeader,subjectNum,subjectName,runNum,useButtonBox,fMRI,rtData,debug) % % Face/house attention experiment with real-time classifier feedback @@ -106,6 +106,7 @@ fixationSize = 4;% pixels progWidth = 400; % image loading progress bar progHeight = 20; +minimumDisplay = 0.25; % how to average the faces and scenes attImgPropPhase1 = .5; @@ -127,6 +128,7 @@ % skyra: use current design button box (keys 1,2,3,4) KbName('UnifyKeyNames'); LEFT = KbName('1!'); +subj_keycode = LEFT; DEVICENAME = 'Current Designs, Inc. 932'; if useButtonBox && (~debug) [index devName] = GetKeyboardIndices; @@ -138,7 +140,6 @@ else DEVICE = -1; end -KbName('UnifyKeyNames') TRIGGER = '5%'; TRIGGER_keycode = getKeys(TRIGGER); % counterbalancing response mapping based on subject assignment @@ -232,11 +233,11 @@ %% Load or Initialize Real-Time Data & Staircasing Parameters - +dataDirHeader = pwd; if matchNum == 0 - dataHeader = [dataDirHeader 'data/' num2str(subjectNum)]; + dataHeader = fullfile(dataDirHeader,[ 'data/' num2str(subjectNum)]); else - dataHeader = [dataDirHeader 'data/' num2str(matchNum) '_match']; + dataHeader = fullfile(dataDirHeader,['data/' num2str(matchNum) '_match']); % shouldn't this be matchNum?? end runHeader = [dataHeader '/run' num2str(runNum)]; @@ -365,39 +366,60 @@ runInstruct{1} = faceInstruct; runInstruct{2} = sceneInstruct; end -%runInstruct{3} = 'Please press to begin task.'; +runInstruct{3} = '-- Please press to begin once you understand these instructions. --'; for instruct=1:length(runInstruct) tempBounds = Screen('TextBounds',mainWindow,runInstruct{instruct}); + if instruct==3 + textSpacing = textSpacing*3; + end Screen('drawtext',mainWindow,runInstruct{instruct},centerX-tempBounds(3)/2,centerY-tempBounds(4)/5+textSpacing*(instruct-1),textColor); clear tempBounds; end Screen('Flip',mainWindow); - -% wait for experimenter to advance with 'q' key -FlushEvents('keyDown'); -pause; -Screen(mainWindow,'FillRect',backColor); -Screen('Flip',mainWindow); - - +[~,~] = WaitTRPulse(subj_keycode,DEVICE); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% now here we're adding to say waiting for scanner, hold tight! +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +waitMessage = 'Waiting for scanner start, hold tight!'; +tempBounds = Screen('TextBounds', mainWindow, waitMessage); +Screen('drawtext',mainWindow,waitMessage,centerX-tempBounds(3)/2,centerY-tempBounds(4)/5+textSpacing*(instruct-1),textColor); +Screen('Flip', mainWindow); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% now here we're going to say to stay still once the triggers start coming +% in +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Start Experiment +STILLREMINDER = ['The scan is now starting.\n\nMoving your head even a little blurs the image, so '... + 'please try to keep your head totally still until the scanning noise stops.\n\n Do it for science!']; +STILLDURATION = 6; % wait for initial trigger Priority(MaxPriority(screenNum)); -Screen(mainWindow,'FillRect',backColor); -Screen(mainWindow,'FillOval',fixColor,fixDotRect); + if (rtData && ~debug ) % if strcmp(computer,'MACI') % taking out because we're running on a linux! %runStart = WaitTRPulsePTB3_skyra(1); - runStart = WaitTRPulse(TRIGGER_keycode,DEVICE); - % else + timing.trig.wait = WaitTRPulse(TRIGGER_keycode,DEVICE); + runStart = timing.trig.wait; + tempBounds = Screen('TextBounds', mainWindow, STILLREMINDER); + Screen('drawtext',mainWindow,runInstruct{instruct},centerX-tempBounds(3)/2,centerY-tempBounds(4)/5+textSpacing*(instruct-1),textColor); + startTime = Screen('Flip',mainWindow); + elapsedTime = 0; + while (elapsedTime < STILLDURATION) + pause(0.005) + elapsedTime = GetSecs()-startTime; + end + % else % WaitSecs(.5); % runStart = KbWait; % end else runStart = GetSecs; + timing = -1; end +Screen(mainWindow,'FillRect',backColor); +Screen(mainWindow,'FillOval',fixColor,fixDotRect); Screen('Flip',mainWindow); Priority(0); @@ -443,7 +465,7 @@ clear tempBounds; timespec = blockData(iBlock).plannedinstructonset - slack; if rtData - [~,blockData(iBlock).instructpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedinstructonset); + [timing.trig.instructpulses(iBlock),blockData(iBlock).instructpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedinstructonset); end blockData(iBlock).actualinstructonset = Screen('Flip',mainWindow,timespec); %#ok % turn on fprintf('Flip time error = %.4f\n', blockData(iBlock).actualinstructonset-blockData(iBlock).plannedinstructonset) @@ -485,7 +507,7 @@ %wait for pulse if (rtData) && mod(iTrial,nTrialsPerTR)==1 %[~,blockData(iBlock).pulses(iTrial)] = WaitTRPulsePTB3_skyra(1,blockData(iBlock).plannedtrialonsets(iTrial)+allowance); %#ok - [~,blockData(iBlock).pulses(iTrial)] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedtrialonsets(iTrial)); + [timing.trig.pulses(iBlock,iTrial),blockData(iBlock).pulses(iTrial)] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedtrialonsets(iTrial)); timespec = blockData(iBlock).plannedtrialonsets(iTrial) - slack; blockData(iBlock).actualtrialonsets(iTrial) = Screen('Flip',mainWindow,timespec); %#ok % turn on else @@ -541,7 +563,7 @@ Screen(mainWindow,'FillOval',fixColor,fixDotRect); % the IBI goes here! if rtData - [~,blockData(iBlock).IBIpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedIBI); + [timing.trig.IBIpulses(iBlock),blockData(iBlock).IBIpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedIBI); end timespec = blockData(iBlock).plannedIBI -slack; blockData(iBlock).actualIBI = Screen('Flip',mainWindow,timespec); @@ -555,7 +577,7 @@ % show instructions-wait for first trigger if rtData -[~,timing.pulsereston] = WaitTRPulse(TRIGGER_keycode,DEVICE,timing.plannedreston); +[timing.trig.pulserest,timing.pulsereston] = WaitTRPulse(TRIGGER_keycode,DEVICE,timing.plannedreston); end tempBounds = Screen('TextBounds',mainWindow,restInstruct); Screen('drawtext',mainWindow,restInstruct,centerX-tempBounds(3)/2,centerY-tempBounds(4)/5,textColor); @@ -600,7 +622,7 @@ clear tempBounds; timespec = blockData(iBlock).plannedinstructonset - slack; if rtData - [~,blockData(iBlock).instructpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedinstructonset); + [timing.trig.instructpulses(iBlock),blockData(iBlock).instructpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedinstructonset); end blockData(iBlock).actualinstructonset = Screen('Flip',mainWindow,timespec); %#ok % turn on fprintf('Flip time error = %.4f\n', blockData(iBlock).actualinstructonset-blockData(iBlock).plannedinstructonset) @@ -656,7 +678,7 @@ %wait for pulse if (rtData) && mod(iTrial,nTrialsPerTR) % this will be true for every other then %[~,blockData(iBlock).pulses(iTrial)] = WaitTRPulsePTB3_skyra(1,blockData(iBlock).plannedtrialonsets(iTrial)+allowance); %#ok - [~,blockData(iBlock).pulses(iTrial)] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedtrialonsets(iTrial)); + [timing.trig.pulses(iBlock,iTrial),blockData(iBlock).pulses(iTrial)] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedtrialonsets(iTrial)); timespec = blockData(iBlock).plannedtrialonsets(iTrial) - slack; blockData(iBlock).actualtrialonsets(iTrial) = Screen('Flip',mainWindow,timespec); %#ok % turn on else @@ -791,7 +813,7 @@ Screen('FillRect',mainWindow,backColor); Screen(mainWindow,'FillOval',fixColor,fixDotRect); if rtData - [~,blockData(iBlock).IBIpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedIBI); + [timing.trig.IBIpulses(iBlock),blockData(iBlock).IBIpulses] = WaitTRPulse(TRIGGER_keycode,DEVICE,blockData(iBlock).plannedIBI); end timespec = blockData(iBlock).plannedIBI -slack; blockData(iBlock).actualIBI = Screen('Flip',mainWindow,timespec); @@ -803,7 +825,7 @@ %% save % question: do you want to save it to this computer's data or where you % save the data folder??? -save([dataHeader '/blockdata_' num2str(runNum) '_' datestr(now,30)],'blockData','runStart'); +save([dataHeader '/blockdata_' num2str(runNum) '_' datestr(now,30)],'blockData','runStart', 'timing'); % clean up and go home sca; diff --git a/displayText.m b/displayText.m new file mode 100644 index 0000000..6622374 --- /dev/null +++ b/displayText.m @@ -0,0 +1,16 @@ +%% display some text +function startTime = displayText(mainWindow,text,duration,horiz,COLOR,WRAPCHARS,timespec) + + DrawFormattedText(mainWindow,text,'center',horiz,COLOR,WRAPCHARS); + if ~exist('timespec', 'var') + timespec = 0; + end + + startTime = Screen('Flip',mainWindow, timespec); + elapsedTime = 0; + while (elapsedTime < duration) + pause(0.005) + elapsedTime = GetSecs()-startTime; + end + +return diff --git a/imgaussfilt3.m b/imgaussfilt3.m deleted file mode 100755 index 71c88cd..0000000 --- a/imgaussfilt3.m +++ /dev/null @@ -1,413 +0,0 @@ -function B = imgaussfilt3(varargin) -%IMGAUSSFILT3 3-D Gaussian filtering of 3-D images -% -% B = imgaussfilt3(A) filters 3-D image A with a 3-D Gaussian smoothing -% kernel with standard deviation of 0.5. -% -% B = imgaussfilt3(A,SIGMA) filters 3-D image A with a 3-D Gaussian -% smoothing kernel with standard deviation specified by SIGMA. SIGMA can -% be a scalar or a 3-element vector with positive values. If sigma is a -% scalar, a cube Gaussian kernel is used. -% -% B = imgaussfilt3(___,Name,Value,...) filters 3-D image A with a 3-D -% Gaussian smoothing kernel with Name-Value pairs used to control aspects -% of the filtering. -% -% Parameters include: -% -% 'FilterSize' - Scalar or 3-element vector, of positive, odd -% integers that specifies the size of the Gaussian -% filter. If a scalar Q is specified, then a square -% Gaussian filter of size [Q Q Q] is used. -% Default value is 2*ceil(2*SIGMA)+1. -% -% 'Padding' - String or numeric scalar that specifies padding to -% be used on image before filtering. -% -% If a scalar (X) is specified, input image values -% outside the bounds of the image are implicitly -% assumed to have the value X. -% -% If a string is specified, it can be 'replicate', -% 'circular' or 'symmetric'. These options are -% analogous to the padding options provided by -% imfilter. -% -% 'replicate' -% Input image values outside the bounds of the image -% are assumed equal to the nearest image border -% value. -% -% 'circular' -% Input image values outside the bounds of the image -% are computed by implicitly assuming the input image -% is periodic. -% -% 'symmetric' -% Input image values outside the bounds of the image -% are computed by mirror-reflecting the array across -% the array border. -% -% Default value is 'replicate'. -% -% 'FilterDomain' - String that specifies domain in which filtering is -% performed. It can be 'spatial', 'frequency' or -% 'auto'. For 'spatial', convolution is performed in -% the spatial domain, for 'frequency', convolution is -% performed in the frequency domain and for 'auto', -% convolution may be performed in spatial or -% frequency domain based on internal heuristics. -% Default value is 'auto'. -% -% -% Class Support -% ------------- -% The input image A must be a real, non-sparse matrix of 3 dimension of -% the following classes: uint8, int8, uint16, int16, uint32, int32, -% single or double. -% -% Notes -% ----- -% 1. If the image A contains Infs or NaNs, the behavior of imgaussfilt3 -% for frequency domain filtering is undefined. This can happen when -% the 'FilterDomain' parameter is set to 'frequency' or when it is set -% to 'auto' and frequency domain filtering is used internally. To -% restrict the propagation of Infs and NaNs in the output in a manner -% similar to imfilter, consider setting the 'FilterDomain' parameter -% to 'spatial'. -% -% 2. When 'FilterDomain' parameter is set to 'auto', an internal -% heuristic is used to determine whether spatial or frequency domain -% filtering is faster. This heuristic is machine dependent and may -% vary for different configurations. For optimal performance, consider -% comparing 'spatial' and 'frequency' to determine the best filtering -% domain for your image and kernel size. -% -% 3. When the 'Padding' parameter is not specified, 'replicate' padding -% is used as the default, which is different from the default used in -% imfilter. -% -% -% Example 1 -% --------- -% This example smooths an MRI volume with a 3-D Gaussian filter. -% -% vol = load('mri'); -% figure, montage(vol.D), title('Original image volume') -% -% siz = vol.siz; -% vol = squeeze(vol.D); -% sigma = 2; -% -% volSmooth = imgaussfilt3(vol, sigma); -% -% figure, montage(reshape(volSmooth,siz(1),siz(2),1,siz(3))) -% title('Gaussian filtered image volume') -% -% See also imgaussfilt, imfilter. - -% Copyright 2014-2015, The MathWorks, Inc. - -narginchk(1, Inf); - -[A, options] = parseInputs(varargin{:}); - -sigma = options.Sigma; -hsize = options.FilterSize; -padding = options.Padding; -domain = options.FilterDomain; - -[domain, ippFlag] = chooseFilterImplementation(A, hsize, domain); - -switch domain - case 'spatial' - B = spatialGaussianFilter(A, sigma, hsize, padding, ippFlag); - - case 'frequency' - B = frequencyGaussianFilter(A, sigma, hsize, padding); - - otherwise - assert(false, 'Internal Error: Unknown filter domain'); -end - -end - -%-------------------------------------------------------------------------- -% Spatial Domain Filtering -%-------------------------------------------------------------------------- -function A = spatialGaussianFilter(A, sigma, hsize, padding, ippFlag) - -sizeA = size(A); -dtype = class(A); -outSize = [sizeA ones(1,3-ndims(A))]; - -[A, padSize] = padImage(A, hsize, padding); - -[hCol,hRow,hSlc] = createSeparableGaussianKernel(sigma, hsize); - -% cast to double to preserve precision over separable filter calls. -if ~isa(A, 'double') - A = double(A); -end - -% Row filtering -outSize1D = [size(A) ones(1,3-ndims(A))]; -outSize1D(2) = outSize(2); - -startRow = zeros(size(padSize)); -startRow(2) = padSize(2); - -A = imfiltermexCaller(A, outSize1D, hRow, startRow, ippFlag); - -% Column filtering -outSize1D(1) = outSize(1); - -startCol = zeros(size(padSize)); -startCol(1) = padSize(1); - -A = imfiltermexCaller(A, outSize1D, hCol, startCol, ippFlag); - -% Slice filtering -outSize1D(3) = outSize(3); - -startSlc = zeros(size(padSize)); -startSlc(3) = padSize(3); - -% IPP only supports 2-D filtering, disable for filtering along 3rd -% dimension. -ippFlag = false; - -A = imfiltermexCaller(A, outSize1D, hSlc, startSlc, ippFlag); - -if ~isa(A,dtype) - A = cast(A,dtype); -end - -end - -function A = imfiltermexCaller(A, outSize, h, start, ippFlag) - -conn = []; -nz = []; - -convMode = ippFlag; % Set the convolution mode to true when IPP is enabled - % to prevent unnecessary (since Gaussian kernels are - % symmetric) rotation of kernel. Conversly, when IPP is - % disabled, set convolution mode to false so home-brew - % does not rotate the kernel. -sameSize = true; - -if ~ippFlag - conn = h~=0; - nz = h(conn); -end - -A = imfilter_mex(A, outSize, h, nz, conn, start, sameSize, convMode, ippFlag); - -end - -function TF = useIPPL(A, outSize) - -prefFlag = useIPPLibrary(); - -if ~isImageIPPFilterType(class(A)) - TF = false; - return; -end - -tooBig = isImageTooBigForIPPFilter(A, outSize); - -TF = prefFlag && ~tooBig; - -end - -function [hcol,hrow,hslc] = createSeparableGaussianKernel(sigma, hsize) - -isIsotropic = all(sigma==sigma(1)) && all(hsize==hsize(1)); - -hcol = images.internal.createGaussianKernel(sigma(1), hsize(1)); - -if isIsotropic - hrow = hcol; - hslc = hcol; -else - hrow = images.internal.createGaussianKernel(sigma(2), hsize(2)); - hslc = images.internal.createGaussianKernel(sigma(3), hsize(3)); -end - -hrow = reshape(hrow, 1, hsize(2)); -hslc = reshape(hslc, 1, 1, hsize(3)); - -end - -%-------------------------------------------------------------------------- -% Frequency Domain Filtering -%-------------------------------------------------------------------------- -function A = frequencyGaussianFilter(A, sigma, hsize, padding) - -sizeA = [size(A) ones(1,3-ndims(A))]; -dtype = class(A); -outSize = sizeA; - -A = padImage(A, hsize, padding); - -h = images.internal.createGaussianKernel(sigma, hsize); - -% cast to double to preserve precision unless single -if ~isfloat(A) - A = double(A); -end - -fftSize = size(A); -A = ifftn( fftn(A, fftSize) .* fftn(h, fftSize), 'symmetric' ); - -if ~strcmp(dtype,class(A)) - A = cast(A, dtype); -end - -start = 1 + size(A)-outSize; -stop = start + outSize - 1; - -A = A(start(1):stop(1),start(2):stop(2),start(3):stop(3)); - -end - -%-------------------------------------------------------------------------- -% Common Functions -%-------------------------------------------------------------------------- -function [domain, ippFlag] = chooseFilterImplementation(A, hsize, domain) - -ippFlag = useIPPL(A, size(A)); - -if strcmp(domain, 'auto') - domain = chooseFilterDomain3(A, hsize, ippFlag); -end - -end - -function [A, padSize] = padImage(A, hsize, padding) - -padSize = computePadSize(size(A), hsize); - -if ischar(padding) - method = padding; - padVal = []; -else - method = 'constant'; - padVal = padding; -end - -A = padarray_algo(A, padSize, method, padVal, 'both'); - -end - -function padSize = computePadSize(sizeA, sizeH) - -rankA = numel(sizeA); -rankH = numel(sizeH); - -sizeH = [sizeH ones(1,rankA-rankH)]; - -padSize = floor(sizeH/2); - -end - -%-------------------------------------------------------------------------- -% Input Parsing -%-------------------------------------------------------------------------- -function [A, options] = parseInputs(varargin) - -A = varargin{1}; - -supportedClasses = {'uint8','uint16','uint32','int8','int16','int32','single','double'}; -supportedImageAttributes = {'real','nonsparse', '3d'}; -validateattributes(A, supportedClasses, supportedImageAttributes, mfilename, 'A'); - -% Default options -options = struct(... - 'Sigma', [.5 .5 .5],... - 'FilterSize', [3 3 3],... - 'Padding', 'replicate',... - 'FilterDomain', 'auto'); - -beginningOfNameVal = find(cellfun(@isstr,varargin),1); - -if isempty(beginningOfNameVal) && length(varargin)==1 - %imgaussfilt3(A) - return; -elseif beginningOfNameVal==2 - %imgaussfilt3(A,'Name',Value) -elseif (isempty(beginningOfNameVal) && length(varargin)==2) || (~isempty(beginningOfNameVal) && beginningOfNameVal==3) - %imgaussfilt3(A,sigma,'Name',Value,...) - %imgaussfilt3(A,sigma) - options.Sigma = validateSigma(varargin{2}); - options.FilterSize = computeFilterSizeFromSigma(options.Sigma); -else - error(message('images:imgaussfilt:tooManyOptionalArgs')); -end - -numPVArgs = length(varargin) - beginningOfNameVal + 1; -if mod(numPVArgs,2)~=0 - error(message('images:imgaussfilt:invalidNameValue')); -end - -ParamNames = {'FilterSize', 'Padding', 'FilterDomain'}; -ValidateFcn = {@images.internal.validateThreeDFilterSize, @validatePadding, @validateFilterDomain}; - -for p = beginningOfNameVal : 2 : length(varargin)-1 - - Name = varargin{p}; - Value = varargin{p+1}; - - idx = strncmpi(Name, ParamNames, numel(Name)); - - if ~any(idx) - error(message('images:imgaussfilt:unknownParamName', Name)); - elseif numel(find(idx))>1 - error(message('images:imgaussfilt:ambiguousParamName', Name)); - end - - validate = ValidateFcn{idx}; - options.(ParamNames{idx}) = validate(Value); - -end - -end - -function sigma = validateSigma(sigma) - -validateattributes(sigma, {'numeric'}, {'real','positive','finite','nonempty'}, mfilename, 'Sigma'); - -if isscalar(sigma) - sigma = [sigma sigma sigma]; -end - -if numel(sigma)~=3 - error(message('images:imgaussfilt:invalidLength3', 'Sigma')); -end - -sigma = double(sigma); - -end - -function filterSize = computeFilterSizeFromSigma(sigma) - -filterSize = 2*ceil(2*sigma) + 1; - -end - -function padding = validatePadding(padding) - -if ~ischar(padding) - validateattributes(padding, {'numeric','logical'}, {'real','scalar','nonsparse'}, mfilename, 'Padding'); -else - padding = validatestring(padding, {'replicate','circular','symmetric'}, mfilename, 'Padding'); -end - -end - -function filterDomain = validateFilterDomain(filterDomain) - -filterDomain = validatestring(filterDomain, {'spatial','frequency','auto'}, mfilename, 'FilterDomain'); - -end