00001 /***********************************************************************
00002   retroicor.h: Implementation of Retrospective Image
00003   Correction for physiological motion effects, using a slightly
00004   modified version of the RETROICOR algorithm described in:
00005 
00006     Glover, G. H., Li, T., & Ress, D. (2000). Image-based method for
00007   retrospective correction of physiological motion effects in fMRI:
00008   RETROICOR. Magnetic Resonance in Medicine, 44, 162-167.
00009 
00010   Fred Tam
00011   Sunnybrook & Women's College Health Sciences Centre
00012   April 15, 2002
00013 
00014   Copyright (C) 2002 Sunnybrook & Women's College Health Sciences Centre,
00015   Toronto, ON, Canada. As part of the AFNI software package, this software
00016   is distributed under the GNU General Public License, Version 2.
00017   See the file README.copyright for details.
00018 ************************************************************************/
00019 
00020 #ifndef _SW__RETROICOR_H_
00021 #define _SW__RETROICOR_H_
00022 
00023 #include "afni.h"
00024 
00025 #define RIC_HISTSIZE 100
00026 #define RIC_HISTFUDGE 0.0000001
00027 
00028 /* Transform cardiac waveform to cardiac phase:
00029 
00030      P_c(t) = 2 * PI * (t - t_1) / (t_2 - t_1)
00031 
00032    where t_1 is the time of the last R-wave peak before t
00033      and t_2 is the time of the next R-wave peak after t
00034 
00035    PRE:
00036    card is the cardiac waveform 1D data;
00037    card has at least 2 timepoints;
00038    card contains float data;
00039    threshold is the value that card data must exceed to be considered part of
00040    an R-wave peak;
00041 
00042    POST:
00043    return value is NULL on error,
00044    else, return value is the phase representation of the entire card dataset
00045    with the same number of timepoints as card;
00046 
00047    It is the caller's responsibility to call mri_free() on the returned struct.
00048 
00049    Note that in this implementation the start and end of the card data are
00050    considered to be R-wave peaks (see above P_c formula for why).
00051 */
00052 MRI_IMAGE * RIC_ToCardiacPhase(const MRI_IMAGE * card, float threshold);
00053 
00054 /* Transform resp waveform to resp phase:
00055 
00056      P_r(t) = PI * sum(H(b), b = 1 .. rnd(R(T)/R_max)) /
00057                    sum(H(b), b = 1 .. 100) * sgn(dR/dt)
00058 
00059    where R(t) is the resp signal amplitude normalized to (0, R_max),
00060          H(b) is the number of occurances in bin b of histogram H
00061               which divides (0, R_max) into 100 bins, each centred at
00062               b * R_max / 100
00063      and dR/dt is computed by convolving with a 39-point kernel (quadratic)
00064 
00065    PRE:
00066    resp is the resp waveform 1D data;
00067    resp has at least 2 timepoints;
00068    resp contains float data;
00069    winsize is the size of the window on the resp data used to estimate slope;
00070    winsize >= 2;
00071 
00072    POST:
00073    return value is NULL on error,
00074    else, return value is the phase representation of the entire resp dataset
00075    with the same number of timepoints as resp;
00076 
00077    It is the caller's responsibility to call mri_free() on the returned struct.
00078 
00079    Implementation note: Convolution is overkill just to get sgn(dR/dt), so
00080    this implementation estimates the slope with the difference between 
00081    winsize-point means on either side of each point. Try setting winsize to
00082    1/2 the sampling rate of resp in Herz. The size of the histogram H is
00083    defined in RIC_HISTSIZE.
00084 */
00085 MRI_IMAGE * RIC_ToRespPhase(const MRI_IMAGE * resp, int winsize);
00086 
00087 /* Calculates the average value for each voxel in dset across all timepoints.
00088 
00089    PRE:
00090    dset is a valid 3D+T dataset with at least 1 timepoint;
00091    ignore is the number of timepoints to ignore at the beginning of dset;
00092    0 <= ignore < number of timepoints in dset;
00093 
00094    POST:
00095    return value is NULL on error,
00096    else, return value is an array of floats with the same dimensions as the
00097    subbricks of dset, containing the voxel value averages;
00098 
00099    Note: The caller is responsible for free()ing the returned block of memory
00100    when done.
00101 
00102    Note2: The complex datatype is not supported, and any such bricks will
00103    result in an error (NULL return value).
00104 */
00105 double * RIC_CalcVoxelMeans(const THD_3dim_dataset * dset, int ignore);
00106 
00107 /* Calculate the a and b per-voxel coefficients for the correction:
00108 
00109    a_x_m = sum((y(t_n) - y_u) * cos(m * P_x(t_n)), n = 1..N) /
00110            sum((cos(m * P_x(t_n))) ^ 2, n = 1..N)
00111 
00112    b_x_m = sum((y(t_n) - y_u) * sin(m * P_x(t_n)), n = 1..N) /
00113            sum((sin(m * P_x(t_n))) ^ 2, n = 1..N)
00114 
00115    where x is from {c, r}, for cardiac, resp, resp.
00116          y_u is the average of y over the time series
00117          m is the order, a value [1, M]
00118 
00119    PRE:
00120    dset is a valid 3D+T dataset with at least 1 timepoint and valid time axis;
00121    phase is the cardiac or resp phase data, spanning the same time as dset;
00122    avg points to an array of voxel averages corresponding to dset (generate
00123     avg with something like this: "avg = RIC_CalcVoxelMeans(dset, ignore);")
00124    a and b are non-NULL pointers to double*;
00125    M is the maximal order of the correction to be made;
00126    M > 0;
00127    ignore is the number of timepoints to ignore at the beginning of dset;
00128    0 <= ignore < number of timepoints in dset;
00129 
00130    POST:
00131    return value is -1 on error,
00132    else, return value is 0 and a and b are pointers to arrays containing the
00133    a and b coefficients described above. The coefficients are stored in the
00134    same order as the voxels in the dset subbricks, offset by m-1 times the
00135    size of the subbricks such that a[i + (m - 1) * DSET_NVALS(dset)] is the
00136    a coefficient for the i'th voxel in dset (zero-based) with order m. Nothing
00137    is done to dset.
00138 
00139    Note: The caller is responsible for free()ing the returned blocks of memory
00140    for a and b when done.
00141 
00142    The complex datatype is not supported, and any such bricks will result in
00143    an error (return value -1).
00144 */
00145 int RIC_CalcCoeffAB(THD_3dim_dataset * dset, const MRI_IMAGE * phase,
00146                     const double * avg, double ** a, double ** b,
00147                     int M, int ignore);
00148 
00149 /* Apply the correction to a dataset _in_place_:
00150 
00151      y_c(t) = y(t) - y_p(t)
00152 
00153      y_p(t) = sum(a_x_m * cos(m * P_x(t)) + b_x_m * sin(m * P_x(t)), m = 1..M)
00154 
00155    where M = 2 is typical, and x is from {c, r} for cardiac, resp, resp.
00156 
00157    PRE:
00158    dset is a valid 3D+T dataset with at least 1 timepoint and valid time axis;
00159    phase is the cardiac or resp phase data, spanning the same time as dset;
00160    a and b are a and b coefficient arrays, as returned by RIC_CalcCoeffAB();
00161    a and b have M * (number of voxels in dset subbrick) entries;
00162    M is the maximal order of the correction to be made;
00163    M > 0;
00164    ignore is the number timepoints to ignore at the beginning of dset;
00165    0 <= ignore < number of timepoints in dset;
00166 
00167    POST:
00168    return value is -1 on error,
00169    else, return value is 0 and dset has been corrected as described above.
00170 
00171    The complex datatype is not supported, and any such bricks will result in
00172    an error (return value -1).
00173 
00174    WARNING WARNING WARNING
00175    =======================
00176    The correction is made to dset "in-place" so if you want to keep the
00177    original dataset, COPY it and pass the copy to this function.
00178 */
00179 int RIC_CorrectDataset(THD_3dim_dataset * dset, const MRI_IMAGE * phase,
00180                        const double * a, const double * b,
00181                        int M, int ignore);
00182 
00183 #endif