00001 function [err] = m3dReorder (Input, Prefix, Mapfile, Opt)
00002 %
00003 %   [err] = m3dReorder (Input, Prefix, Mapfile, Opt)
00004 %
00005 %Purpose:
00006 %   Reorders a time series data set (3D+time) a la AFNI plugin Reorder   
00007 %   
00008 %   
00009 %Input Parameters:
00010 %
00011 %   Input: Input of 3d+time brick
00012 %   Prefix: prefix of output data set
00013 %   Mapfile: Name of map file
00014 %   Opt is the options structure with the following fields
00015 %       .Dup : [Col]/Ave Collate(default) or Average the multiple
00016 %             instances in the map file. 
00017 %     .Detrend : 0/1/[2] Linear trend removal. 0 for none, 
00018 %                1 for mean only, 2 for linear trend (default)  
00019 %     .Verbose : [0]/1 verbosity of function ...
00020 %
00021 %Output Parameters:
00022 %   err : 0 No Problem
00023 %       : 1 Mucho Problems
00024 %   
00025 %   
00026 %      
00027 %Key Terms:
00028 %   
00029 %More Info :
00030 %   
00031 %   help button in Reorder plugin
00032 %   
00033 %
00034 %     Author : Ziad Saad
00035 %     Date : Tue Sep 18 10:39:36 EDT 2001
00036 %     LBC/NIMH/ National Institutes of Health, Bethesda Maryland
00037 
00038 %Define the function name for easy referencing
00039 FuncName = 'm3dReorder';
00040 
00041 %Debug Flag
00042 DBG = 1;
00043 
00044 %initailize return variables
00045 err = 1;
00046 
00047 %check for options      
00048         if (nargin < 4), Opt.Dup = 'Col'; Opt.Detrend = 2; Opt.Verbose = 0; end
00049         
00050         if (~isfield(Opt,'Dup') | isempty(Opt.Dup)), Opt.Dup = 'Col'; end
00051         if (~isfield(Opt,'Detrend') | isempty(Opt.Detrend)), Opt.Detrend = 2; end
00052         if (~isfield(Opt,'Verbose') | isempty(Opt.Verbose)), Opt.Verbose = 0; end
00053         
00054 %check for valid parameters
00055         if (~eq_str(Opt.Dup,'Col') & ~eq_str(Opt.Dup,'Ave')),
00056                 err = ErrEval(FuncName,'Err_Bad value for Opt.Dup');
00057                 return;
00058         end
00059 
00060 %check the prefix
00061         [status, Prefix] = PrefixStatus(Prefix); 
00062         if (status == 0), err = ErrEval(FuncName,'Err_Bad Prefix'); return; end
00063         
00064 %check and read the input file
00065         if (Opt.Verbose), fprintf(1,'\nLoading & Sorting Mapfile ...'); end
00066         fid = fopen (Mapfile,'ro');
00067         sall = fscanf(fid,'%c');
00068         fclose (fid);
00069         
00070 %load the file into a cell string
00071         %search for Newlines
00072         inl = find (sall == 10);
00073         N_inl = length(inl);
00074         %fill in sall to call, ignoring #
00075         icell = 1;
00076         ipos = 1;
00077         vsall = '';
00078         for (i = 1:1:N_inl),
00079                 iend = inl(i) -1;
00080                 sword = sall(ipos:iend);
00081                 if (sall(ipos) ~= '#' & ~isempty(sword)), %Not a comment 
00082                                 %vertically concatenate the strings
00083                                 vsall = strvcat(vsall, sword);
00084                 end
00085                 ipos = inl(i)+1;
00086         end
00087 
00088 % sort the results
00089         [vsall_sort,i_sort] = sortrows (vsall);
00090 
00091 %find the ignore points
00092         ikeep = find (vsall_sort(:,1) ~= '-');
00093         imap = i_sort(ikeep);
00094         %vsall_sort(1:20,:)
00095         %ikeep(1:5)
00096         %length(imap) 
00097         
00098 %load the data set
00099         %make sure that this is a 3D+time
00100         [err, Info] = BrikInfo(Input);
00101         if (err), 
00102                 err = ErrEval(FuncName,'Err_error in BrikInfo. Check Input Brick Name'); return; 
00103         end
00104         if (Info.SCENE_DATA(2) ~= 2), 
00105                 err = ErrEval(FuncName,sprintf('Err_%s is not an EPI type', Input)); return;
00106         end
00107         %make sure length of time series is OK
00108         if (Info.DATASET_RANK(2) ~= size(vsall,1)),
00109                 err = ErrEval(FuncName,...
00110                 sprintf('Err_Length mismatch between time series (%g) and Mapfile(%g)',...
00111                  Info.DATASET_RANK(2),size(vsall,1))); return;  
00112         end
00113         
00114         
00115         %OK, load data
00116         if (Opt.Verbose), fprintf(1,'\nLoading Brick ...'); end
00117         OptR.Format = 'vector';
00118         [err, V, Info, ErrMessage] = BrikLoad (Input, OptR);
00119         if (err), err = ErrEval(FuncName,'Err_error in BrikLoad'); return; end
00120         
00121         %figure(1);clf; subplot (211);  plot (V(1,:), 'b'); 
00122         
00123         %detrend ?
00124         if (Opt.Detrend == 1),
00125                 if (Opt.Verbose), fprintf(1,'\nDetrending ...'); end
00126                 V = detrend(V','constant'); V = V';
00127         elseif (Opt.Detrend == 2),
00128                 if (Opt.Verbose), fprintf(1,'\nDetrending ...'); end
00129                 V = detrend(V','linear'); V = V';
00130         elseif (Opt.Detrend ~= 0),
00131                 err = ErrEval(FuncName, 'Err_Bad value for Opt.Detrend'); return;
00132         end     
00133         %plot (V(1,:), 'r'); hold on; drawnow
00134         %pause
00135 
00136 %Now the output brick is:
00137         Vo = V(:,imap);
00138         %subplot (212); plot (Vo(1,:)); drawnow
00139         %pause
00140 %update the header
00141         Info_o = Info;
00142         Info_o.RootName = '';
00143         Info_o.IDCODE_STRING = '';
00144         Info_o.IDCODE_DATE = '';
00145         Info_o.TAXIS_OFFSETS = [];
00146         Info_o.BRICK_FLOAT_FACS = [];
00147         Info_o.BRICK_STATS = [];
00148         Info_o.BRICK_TYPES = Info.BRICK_TYPES(imap);
00149         Info_o.DATASET_RANK(2) = length(imap);
00150         Info_o.TAXIS_NUMS(1:2) = [length(imap) 0]; %remove time offset because it's meaningless when you scramble the data
00151 
00152 %Write out results
00153         if (Opt.Verbose), fprintf(1,'\nWriting Brick ...'); end
00154         OptW.Scale = 1;
00155         OptW.Prefix = Prefix;
00156         OptW.AppendHistory = 1;
00157         OptW.NoCheck = 1;
00158         [err, ErrMessage, Info] = WriteBrik (Vo, Info_o, OptW);
00159         if (err), err = ErrEval(FuncName,'Err_Error writing brick to disk'); return; end
00160 
00161 err = 0;
00162 return;
00163