00001 %script MatAFNI_Demo
00002 %
00003 %
00004 %
00005 %Purpose:
00006 % Demonstrate the use of the matlab AFNI library
00007 % The script is echoed onto the screen but it is best
00008 % you read through it with a text editor.
00009 %
00010 %
00011 %Input:
00012 % You must have the data sets
00013 % ARzs_CW_avvr+orig, ARzs_CW_avvr.DEL+orig and ARzsspgrax+orig
00014 % in the working directory
00015 %
00016 %Output:
00017 %
00018 %
00019 %
00020 %
00021 %
00022 %Key Terms:
00023 %
00024 %More Info :
00025 %
00026 %
00027 %
00028 %
00029 % Author : Ziad Saad
00030 % Date : Tue Aug 28 16:21:54 PDT 2001
00031 % SSCC/NIMH/ National Institutes of Health, Bethesda Maryland
00032
00033 echo on
00034
00035 clear all
00036
00037 %check for data sets before proceeding
00038 if (exist('ARzs_CW_avvr+orig.HEAD') ~= 2 | exist('ARzs_CW_avvr+orig.BRIK') ~= 2),
00039 fprintf ('\nERROR: Missing ARzs_CW_avvr+orig data.\nCtrl+c to Abort.\n');
00040 while (1),
00041 pause;
00042 end
00043 end
00044
00045 if (exist('ARzs_CW_avvr.DEL+orig.HEAD') ~= 2 | exist('ARzs_CW_avvr.DEL+orig.BRIK') ~= 2),
00046 fprintf ('\nERROR: Missing ARzs_CW_avvr.DEL+orig data.\nCtrl+c to Abort.\n');
00047 while (1),
00048 pause;
00049 end
00050 end
00051 if (exist('ARzsspgrax+orig.HEAD') ~= 2 | exist('ARzsspgrax+orig.BRIK') ~= 2),
00052 fprintf ('\nERROR: Missing ARzsspgrax+orig data.\nCtrl+c to Abort.\n');
00053 while (1),
00054 pause;
00055 end
00056 end
00057
00058 %%%%%%%%%%%%%% Read time series and access time series %%%%%%%%%%%%%%%%%%%%%%%%
00059 fprintf (1,'\tpatience...\n');
00060 %read in a brick and display voxel with AFNI i j k: 29, 33, 3
00061 BrikName = 'ARzs_CW_avvr+orig';
00062 Vox_Ind = [29 33 3];
00063
00064 %Do it using matrix format for V.
00065 Opt.Format = 'matrix';
00066 [err, Vm, Info, ErrMessage] = BrikLoad (BrikName, Opt);
00067 %plot it
00068 figure(1); clf; subplot (211);
00069 plot (squeeze(Vm(Vox_Ind(1)+1,Vox_Ind(2)+1, Vox_Ind(3)+1, :)));
00070
00071
00072 %Do it using vector format for V. Note that this example does not highlight
00073 %the advantage of using the vector format.
00074 Opt.Format = 'vector';
00075 [err, Vv, Info, ErrMessage] = BrikLoad (BrikName, Opt);
00076 %figure out the AFNI 1D index of the voxel
00077 [err, Indx] = AfniXYZ2AfniIndex (Vox_Ind, Info.DATASET_DIMENSIONS(1), Info.DATASET_DIMENSIONS(2));
00078 %plot it. Remember that indexing in matlab is always augmented by 1 relative to AFNI's
00079 figure(1); subplot (212);
00080 plot (Vv(Indx+1,:));
00081 plotsign2(1);
00082
00083 drawnow
00084 msgbox(sprintf('Showing Time Series From %s.\nFrom AFNI''s viewer, use jump to: %g %g %g to see the same time series',...
00085 BrikName, Vox_Ind(1),Vox_Ind(2), Vox_Ind(3)));
00086 fprintf (1,'\nHit Enter to Proceed with Demo (ctrl+c to abort)...\n'); pause
00087
00088 %%%%%%%%%%%%%% Read Anatomy and grab a slice %%%%%%%%%%%%%%%%%%%%%%%%
00089 fprintf (1,'\tpatience...\n');
00090 %read in anatomy brick
00091 BrikNameAnat = 'ARzsspgrax+orig';
00092 [err, Vanat, InfoAnat, ErrMessage] = BrikLoad (BrikNameAnat);
00093
00094 %Grab a few slices in the same plane
00095 OptDisp.plane = 'Ax';
00096 OptDisp.iSlc = [107 79 45]; %Those are AFNI indices, see help GetAfniSlice for more info
00097 [err, slc, slcDisp] = GetAfniSlice (Vanat, InfoAnat, OptDisp);
00098 figure(2); clf
00099 subplot 321; imagesc (slc.M(:,:,1)); title (sprintf('Ax #%g, As In Brick', OptDisp.iSlc(1)));
00100 subplot 322; imagesc (slcDisp.M(:,:,1)); title (sprintf('Ax #%g, As Displayed', OptDisp.iSlc(1)));
00101 subplot 323; imagesc (slc.M(:,:,2));title (sprintf('Ax #%g, As In Brick', OptDisp.iSlc(2)));
00102 subplot 324; imagesc (slcDisp.M(:,:,2)); title (sprintf('Ax #%g, As Displayed', OptDisp.iSlc(2)));
00103 subplot 325; imagesc (slc.M(:,:,3));title (sprintf('Ax #%g, As In Brick', OptDisp.iSlc(3)));
00104 subplot 326; imagesc (slcDisp.M(:,:,3)); title (sprintf('Ax #%g, As Displayed', OptDisp.iSlc(3)));
00105 colormap gray %you need to use some scaling to get the display to work best
00106 plotsign2(2);
00107
00108 drawnow
00109 msgbox(sprintf('Showing Anatomy Brick %s.\nFrom AFNI''s axial viewer, use slider to view slices %g, %g and %g that are displayed in the right column of the figure.',...
00110 BrikNameAnat, OptDisp.iSlc(1),OptDisp.iSlc(2), OptDisp.iSlc(3)));
00111 fprintf (1,'\nHit Enter to Proceed with Demo (ctrl+c to abort)...\n'); pause
00112
00113 %%%%%%%%%%%%%% Read Functional brick and grab a slice triplet %%%%%%%%%%%%%%%%%%%%%%%%
00114 fprintf (1,'\tpatience...\n');
00115 %grab a triplet of slices from .DEL functional Brick
00116 BrikNameFunc = 'ARzs_CW_avvr.DEL+orig';
00117 [err, Vfunc, Infofunc, ErrMessage] = BrikLoad (BrikNameFunc);
00118
00119 DM = AFNI_SliceDispManip (Infofunc);
00120 OptDispFunc.iSlc = [31 32 6];
00121 OptDispFunc.index = 1;
00122 [err, slc] = GetAfniSliceTriplet (Vfunc, Infofunc, DM, OptDispFunc);
00123 figure(3); clf
00124 subplot 321; imagesc (slc(1).M(:,:)); title (sprintf('Ax #%g, As In Brick', OptDispFunc.iSlc(1)));
00125 subplot 322; imagesc (slc(1).Mdisp(:,:)); title (sprintf('Ax #%g, As Displayed', OptDispFunc.iSlc(1)));
00126 subplot 323; imagesc (slc(2).M(:,:));title (sprintf('Sa #%g, As In Brick', OptDispFunc.iSlc(2)));
00127 subplot 324; imagesc (slc(2).Mdisp(:,:)); title (sprintf('Sa #%g, As Displayed', OptDispFunc.iSlc(2)));
00128 subplot 325; imagesc (slc(3).M(:,:));title (sprintf('Co #%g, As In Brick', OptDispFunc.iSlc(3)));
00129 subplot 326; imagesc (slc(3).Mdisp(:,:)); title (sprintf('Co #%g, As Displayed', OptDispFunc.iSlc(3)));
00130 colormap gray %you need to use some scaling to get the display to work best
00131 plotsign2(3);
00132
00133 drawnow
00134 msgbox(sprintf('Showing Slice Triplet from Functional Brick %s.\nFrom AFNI''s axial viewer, use slider to view slices %g, %g and %g that are displayed in the right column of the figure.\nSet Func underlay and set Func sub-brick to #%g',...
00135 BrikNameFunc, OptDispFunc.iSlc(1),OptDispFunc.iSlc(2), OptDispFunc.iSlc(3), OptDispFunc.index));
00136 fprintf (1,'\nHit Enter to Proceed with Demo (ctrl+c to abort)...\n'); pause
00137
00138
00139 %%%%%%%%%%%%%%%%% Show a cool graph %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00140 fprintf (1,'\tpatience...\n');
00141 %Extract the delay and cross correlation coefficient sub-bricks
00142 %see Infofunc.BRICK_LABS for indexing info
00143 Vdel = squeeze(Vfunc(:,:,:,1)); %that's not necessary but it makes annotation and referencing clear
00144 Vxcor = squeeze(Vfunc(:,:,:,3));
00145
00146 %find the cross correlation coefficients < 0.5
00147 inoact = find(Vxcor < 0.5);
00148
00149 %mask those voxels with -1
00150 Vdel(inoact) = -1;
00151
00152 %make sure no voxels having a delay outside [0 .. 40] are masked
00153 iout = find(Vdel < 0 | Vdel > 40);
00154 Vdel(iout) = -1; %Note that for efficiency, the previous two steps can be combined into 1
00155
00156 %also mask the Cross Correlation Subbrick for later use
00157 Vxcor(iout) = 0;
00158
00159 %Extract axial slice 29
00160 OptFunc2.iSlc = 29;
00161 OptFunc2.plane = 'Ax';
00162 [err, slc, slcDisp] = GetAfniSlice (Vdel, Infofunc, OptFunc2);
00163
00164 %interpolate for pretty pictures
00165 SlcShow = interp2(slcDisp.M, 3, 'cubic');
00166
00167 %display a surface mesh of delay values for slice
00168 figure(4); clf
00169 subplot (211);
00170 hs = surf(SlcShow, 'EdgeColor', 'none'); title ('Surface Plot of Delay Values');
00171 subplot (223);
00172 hc = contour (SlcShow); title ('Contour Plot');
00173 subplot (224);
00174 ipos = find(Vdel >= 0);
00175 hist(Vdel(ipos),20); title ('Histogram');
00176 plotsign2(4);
00177 drawnow;
00178 fprintf (1,'\nHit Enter to Proceed with Demo (ctrl+c to abort)...\n'); pause
00179
00180 %%%%%%%%%%%%%%%%% Writing Out Functional Data %%%%%%%%%%%%%%%%%%%%%%%%%%%
00181 %Write out delay brick with only the masked voxels showing
00182
00183 fprintf (1,'\tpatience...\n');
00184 %We need to setup the header info for the new data set to write out
00185 %Copy the header info from the .DEL brick
00186 InfoDelOut = Infofunc;
00187 %modify the necessary fields
00188 %YOU MUST READ THE HELP for the function WriteBrik and BrikInfo
00189 InfoDelOut.RootName = ''; %that'll get set by WriteBrik
00190 InfoDelOut.DATASET_RANK(2) = 2; %two sub-bricks
00191 InfoDelOut.BRICK_TYPES = [1 1]; %store data as shorts
00192 InfoDelOut.BRICK_STATS = []; %automatically set
00193 InfoDelOut.BRICK_FLOAT_FACS = [];%automatically set
00194 InfoDelOut.BRICK_LABS = 'Masked Delay~Masked Cross Correlation';
00195 InfoDelOut.IDCODE_STRING = '';%automatically set
00196 InfoDelOut.BRICK_STATAUX = [1 2 3 160 2 2];
00197
00198 OptDelOut.Scale = 1;
00199 OptDelOut.Prefix = 'Demo1_func';
00200 OptDelOut.verbose = 0;
00201 !rm Demo1_func*
00202 %To write the two sub-bricks, I placed them in an Nx2 matrix [Vdel(:) Vxcor(:)]
00203 [err, ErrMessage, InfoDelOut] = WriteBrik ([Vdel(:) Vxcor(:)], InfoDelOut, OptDelOut);
00204
00205 if (err),
00206 errordlg(sprintf('%s\nAbort now (ctrl+c).', ErrMessage));
00207 while (1),
00208 pause;
00209 end
00210 end
00211
00212 msgbox(sprintf('Data set %s+orig was written to disk, check it out with AFNI.',...
00213 OptDelOut.Prefix));
00214
00215 fprintf (1,'\nHit Enter to Proceed with Demo (ctrl+c to abort)...\n'); pause
00216
00217 %%%%%%%%%%%%%%%%% Writing Out time series %%%%%%%%%%%%%%%%%%%%%%%%%%%
00218 fprintf (1,'\tpatience...\n');
00219 %Write out random time series where functional voxels passed the threshold above
00220 %find voxel indices that are not masked
00221 [i_in] = setdiff([1:1:prod(Info.DATASET_DIMENSIONS(1:3))],iout);
00222 %create that many random time series
00223 Mrand = randn(length(i_in), Info.DATASET_RANK(2));
00224 %create an empty volume
00225 Vnew = zeros(prod(Info.DATASET_DIMENSIONS(1:3)), Info.DATASET_RANK(2));
00226 %replace zeros with random time series
00227 Vnew(i_in,:) = Mrand;
00228
00229 %prepare Header for writing
00230 InfoNewTSOut = Info;
00231 InfoNewTSOut.RootName = '';
00232 InfoNewTSOut.BRICK_STATS = []; %automatically set
00233 InfoNewTSOut.BRICK_FLOAT_FACS = []; %automatically set
00234 InfoNewTSOut.IDCODE_STRING = '';%automatically set
00235
00236 OptTSOut.Scale = 1;
00237 OptTSOut.Prefix = 'Demo1_TS';
00238 OptTSOut.verbose = 1;
00239 !rm Demo1_TS*
00240 %write it
00241 [err, ErrMessage, InfoNewTSOut] = WriteBrik (Vnew, InfoNewTSOut, OptTSOut);
00242 msgbox(sprintf('Data set %s+orig was written to disk, check it out with AFNI.',...
00243 OptTSOut.Prefix));
00244
00245 if (err),
00246 errordlg(sprintf('%s\nAbort now (ctrl+c).', ErrMessage));
00247 while (1),
00248 pause;
00249 end
00250 end
00251
00252 fprintf (1,'\nDemo Done\n');
00253
00254 echo off
