Hi

I've attached an updated version of the lmgs detrending procedure described
in P.M. Macey, K.E. Macey, R. Kumar, and R.M. Harper. A method for removal
of global effects from fMRI time series. NeuroImage 22 (2004) 360-366. This
follows a request for a version that works with SPM5. It's otherwise
unchanged from the earlier one posted; type "help cspm_lmgs" for details.

Best wishes,
Paul Macey
Dept. Neurobiology, UCLA



--'8Lmkg-vOV/CzF5wf-z830'aKChstm4KwEh'u8,hA0LV8STADRB,14NfouoMhWc?KCVeAy
Content-Transfer-Encoding: 7bit
Content-Type: TEXT/PLAIN;
 name="cspm_lmgs.m"

function cspm_lmgs(P,overwrite,prefix,display)
% CSPM_LMGS - Detrend fMRI image series using LMGS method.
% CSPM_LMGS(P, OVERWRITE, PREFIX, DISPLAY )
% P is array of images returned by spm, e.g. (SPM2):
% P = spm_get(Inf, '.img',{'Please select images for detrending'});
% You will be prompted for the files the argument is omitted.
% OVERWRITE (optional) - defaults to true; if false and detrended
% images exist, detrending is skipped.
% PREFIX (optional) - prefix for detrended files; defaults to "d".
% DISPLAY (optional) - if true, global trends of raw and detrended
% data are plotted; defaults to false.
%
% The procedure applies the detrending to the time-series of files P, and writes
% detrended files prepended with PREFIX (defaults to "d"). To run, type
% >>cspm_lmgs
% To see the result of the detrending, set the DISPLAY flag to 1, e.g.:
% >>P = spm_get(Inf, '.img',{'Please select images for detrending'});
% >>cspm_lmgs(P, 1, 'd', 1)
%
% LMGS (Linear Model of Global Signal) detrending is described in:
% P.M. Macey, K.E. Macey, R. Kumar, and R.M. Harper. A method for removal of
% global effects from fMRI time series. NeuroImage 22 (2004) 360-366.
%
% Works with SPM5 and SPM2.
%

% @(#)cspm_lmgs.m 1.1 Paul Macey 2006-04-01


% Check inputs
% Images
if nargin < 1
    if strcmp(spm('ver'),'SPM5')
        P = spm_select([3 Inf],'image','Select images for detrending');
    else
        P = spm_get(Inf, '.img','Select images for detrending');
    end
    if isempty(P), return, end
end
if iscell(P)
    P = char(P);
end
if length(P) < 3
    msgbox('Detrending requires a time-series of images (multiple volumes)',...
        'LMGS',...
        'warn')
    return
end
% Flags, etc.
if nargin < 2 | isempty(overwrite)
    overwrite = 1;
end
if nargin < 3 | isempty(prefix)
    prefix = 'd';
else
    if ~ischar(prefix)
        msgbox('Invalid "Prefix" argument - need character',...
            'LMGS','warn')
        help cspm_lmgs
        return
    end
end
if nargin < 4
    display = 0;
end

% Number of scans
nscan = size(P,1);

% Finish now if overwrite flag is false and files exist
if ~overwrite
    skip = 1;
    % Check for existing files
    for i = 1:nscan

        [pth,nm,ext] = fileparts(P(i,:));

        % Detrended
        newfname = [pth,'\',prefix,nm,ext];

        if ~exist(newfname,'file')
            skip = 0;
            break
        end
    end
    if skip, return, end
end

% Map volumes
V = spm_vol(P);
% Check similar dimensions, etc. (Taken from spm_imcalc_ui.)
if (strcmp(spm('ver'),'SPM5') & any(any(diff(cat(1,V.dim),1,1),1))) | ...
        (~strcmp(spm('ver'),'SPM5') & any(any(diff(cat(1,V.dim),1,1),1)&[1,1,1,0])) | ...
        any(any(any(diff(cat(3,V.mat),1,3),3)))
    msgbox('Images don''t have same dimmensions, orientation and voxel size - detrending cancelled.',...
        'LMGS','error')
    return
end

% Number of voxels
num_vox = prod(V(1).dim(1:3));

% Initialize - create independent variable X for model
% Global values
gt = cspm_globaltrend(P,0,0,0)';
% Remove offset
X = detrend(gt,'constant');
% Add Constant
X = [X ones(nscan,1)];

% Create output images
for i = 1:nscan
    Vnew = V(i);
    [pth,nm,ext] = fileparts(Vnew.fname);
    Vnew.fname = [pth,filesep,prefix,nm,ext];
    if strcmp(spm('ver'),'SPM99')
        Vout(i) = spm_create_image(Vnew);
    else
        Vout(i) = spm_create_vol(Vnew);
    end
end

adjustcount = 0;
num_vox_per_plane = num_vox / V(1).dim(3);
dim1 = V(1).dim(1);
dim2 = V(1).dim(2);
dim3 = V(1).dim(3);
dim12 = V(1).dim(1:2);
% Loop through planes
fprintf('\n%-60s','Detrending')
for z = 1:dim3
    fprintf('%s%-60s',char(sprintf('\b')*ones(1,60)),['Detrending plane ',int2str(z),' / ',int2str(dim3)])

    % Get raw data for plane across all volumes.
    for i = 1:nscan
        yn(:,:,i) = spm_slice_vol(V(i),spm_matrix([0 0 z]),dim12,0);
    end

    % Use 2-D array to speed processing.
    yadj = reshape(yn,num_vox_per_plane,nscan)';

    warning off
    % Perform detrending voxel-by-voxel in this plane
    for v = 1:num_vox_per_plane
        yvox = yadj(:,v);
        if any(yvox)
            b = regression(yvox,X);
            % The following would require the statistics toolbox:
            % [b,bint,r,rint,stats] = regress(yvox,X,0.05);
            % Test for positive effect - b(1);
            % Test for significance as p-value - - stats(3);
            % Test for correlation (R^2) - stats(1)
            % Optional: skip voxels that do not meet thresholds
            % if b(1) <= 0 | stats(3) > 0.8 | stats(1) < 0.3
            % continue
            % end
            yvoxmodel = b(1)*X(:,1);

            % Adjust time-series
            adjustcount = adjustcount + 1;
            yadj(:,v) = yvox - yvoxmodel;

        end
    end
    warning on
    yadj = reshape(yadj',dim1,dim2,nscan);

    % Write adjusted data
    for i = 1:nscan
        Vout(i) = spm_write_plane(Vout(i),yadj(:,:,i),z);
    end

end % End loop through planes
fprintf('%s%-60s',char(sprintf('\b')*ones(1,60)),'Detrending')
fprintf('%-60s\n','Done')

disp(['Adjusted ', int2str(adjustcount),...
    ' non-zero voxels (',int2str(num_vox),' total).'])

% Save file names if displaying
if display
    for i = 1:nscan
        Pin{i} = V(i).fname;
        Pd{i} = Vout(i).fname;
    end
end

% Close volumes, if SPM2
if strcmp(spm('ver'),'SPM2')
    spm_close_vol(V);
    spm_close_vol(Vout);
end

% Show in percent change the pre- and post-detrending global timetrends
if display
    t1 = gt'; % Variable gt was calculated earlier.
    t2 = cspm_globaltrend(Pd,0,0,0)';

    figure('NumberTitle','off','Name','Effect of adjustment')

    b1 = mean(t1(1:end));
    b2 = mean(t2(1:end));

    % Percent change
    t1pc = 100*(t1 - b1)/b1;
    t2pc = 100*(t2 - b2)/b2;

    plot(t1pc,'b.-')
    a = gca;
    title ('Un-adjusted images (%)')
    YLim = get(a,'YLim');
    range = YLim(2) - YLim(1);
    hold on
    plot(t2pc,'r.-')
    title ('Adjusted images (%)')
    set(a,'YLim',YLim)
    legend('Raw (blue)','Adjusted (red)')
end

% =========================================================================
function gt = cspm_globaltrend( P, pc, usespm, display )
% CSPM_GLOBALTREND - Calculate and display global time trend of images.
% T = CSPM_GLOBALTREND( P, PC, SPM, DISPLAY)
% P - cell array of files, e.g.,
% P = spm_get(Inf, '.img',{'Please select images for detrending'});
% PC - flag 1 = percent change relative to mean (default),
% 0 = absolute value
% SPM - flag 1 = use spm_global (all voxels above 1/8 of maximum
% intensity; not very accurate)
% 0 = avearge of all voxels in volume (default)
% DISPLAY - flag 1 = make figure (default)
% 0 = no figure
% Output: T - global time trend

% @(#)cspm_globaltrend.m 1.0 Paul Macey 2005-08-01

% Calculate global trend of series of images
if nargin < 1
    if strcmp(spm('ver'),'SPM5')
        P = spm_select(Inf,'image','Select images for global trend calculation');
    else
        P = spm_get(Inf, '.img','Select images for global trend calculation');
    end
    if isempty(P), return, end
end
if iscell(P)
    P = char(P);
end
if nargin < 2
    pc = 1;
end
if nargin < 3
    usespm = 0;
end
if nargin < 4
    display = 1;
end

V = spm_vol(P);
for i = 1:length(V)
    % Note: spm_global estimates the mean after discounting voxels outside
    % the object using a criteria of greater than > (global mean)/8.
    % However, for large global signal changes, this does not accurately
    % estimate the global signal and LMGS detrending based on spm_global
    % does not remove all gobal components.
    if usespm
        gt(i) = spm_global(V(i));
    else
        Y = spm_read_vols(V(i));
        gt(i) = mean(mean(mean(Y)));
    end
end
if strcmp(spm('ver'),'SPM2')
    spm_close_vol(V);
end

if pc
    bl = mean(gt);
    gt = 100*(gt-bl)/bl;
    pstr = ' (% change)';
else
    pstr = ' (raw)';
end

if display
    figure('Name',['Global timetrend for ',int2str(length(V)),' files',pstr])
    plot(gt)
    ylabel(['Signal Intensity',pstr])
    xlabel(['Scan Number'])
end

if nargout == 0
    clear gt
end

% =========================================================================
function b = regression(y,X)
% Check that matrix (X) and left hand side (y) have compatible dimensions
% [n,p] = size(X);
% [n1,collhs] = size(y);
% if n ~= n1,
% error('The number of rows in Y must equal the number of rows in X.');
% end
%
% if collhs ~= 1,
% error('Y must be a vector, not a matrix');
% end

% Remove missing values, if any
wasnan = (isnan(y) | any(isnan(X),2));
if (any(wasnan))
    y(wasnan) = [];
    X(wasnan,:) = [];
    n = length(y);
end

% Find the least squares solution.
[Q, R]=qr(X,0);
b = R\(Q'*y);



--'8Lmkg-vOV/CzF5wf-z830'aKChstm4KwEh'u8,hA0LV8STADRB,14NfouoMhWc?KCVeAy--