happy
Description:
happy is a new(er) addition to the rapidtide suite (it was added in 2019, so it’s hardly new at this point). It’s complementary to rapidtide - it’s focussed on fast, cardiac signals in fMRI, rather than the slow, LFO signals we are usually looking at. It’s sort of a Frankenprogram - it has three distinct jobs, which are related, but are very distinct.
The first thing happy does is try to extract a cardiac waveform from the fMRI data. This is something I’ve been thinking about for a long time. It occurred to me that while the TR for most scans is long compared to the required sample time for recording a cardiac waveform, the scanner is actually recording data at a much faster rate than that - each slice, or in the case of multiband data, stack of slices, is acquired somewhere during each TR, so the _effective_ samplerate is TR/(number of acquisitions within a TR).
The second task is to take this raw estimate of the cardiac waveform, and clean it up using a deep learning filter. The original signal is useful, but pretty gross, but I figured you should be able to exploit the pseudoperiodic nature of the signal to greatly improve it. This was also a testbed to work on using neural nets to process time domain signals. It seemed like a worthwhile project, so it got grafted in.
The final task (which was actually the initial task, and the reason I wrote happy to begin with) is to implement Henning Voss’ totally cool hypersampling with analytic phase projection [1] (guess where the name “happy” comes from). This is fairly straightforward, as Voss describes his method very clearly. But I have lots of data with no simultaneously recorded cardiac signals, and I was too lazy to go find datasets with pleth data to play with, so that’s why I did the cardiac waveform extraction part. In retrospect, that’s part is pretty cool in it’s own right, if I do say so myself.
The paper describing the theory of operation, development, and testing of this program can be found here [2].
Inputs:
Happy needs a 4D BOLD fMRI data file (space by time) as input. This can be Nifti1 or Nifti2. If you have a simultaneously recorded cardiac waveform, it will happily (heh heh) use it, otherwise it will try to construct (and refine) one. NOTE: the 4D input dataset needs to be completely unpreprocessed - gradient distortion correction and motion correction can destroy the relationship between slice number and actual acquisition time, and slice time correction does not behave as expected for aliased signals (which the cardiac component in fMRI most certainly is), and in any case we need the slice time offsets to construct our waveform.
Outputs:
Outputs are space or space by time Nifti or text files, depending on what the input data file was, and some text files containing textual information, histograms, or numbers. File formats and naming follow BIDS conventions for derivative data for fMRI input data. Output spatial dimensions and file type match the input dimensions and file type (Nifti1 in, Nifti1 out). Depending on the file type of map, there can be no time dimension, a time dimension that matches the input file, or something else, such as a time lag dimension for a correlation map.
BIDS Outputs:
Name |
Extension(s) |
Content |
When present |
|---|---|---|---|
XXX_commandline |
.txt |
The command line used to run happy |
Always |
XXX_formattedcommandline |
.txt |
The command line used to run happy, attractively formatted |
Always |
XXX_desc-rawapp_info |
.nii.gz |
The analytic phase projection map of the cardiac waveform |
Always |
XXX_desc-app_info |
.nii.gz |
The analytic phase projection map of the cardiac waveform, voxelwise minimum subtracted |
Always |
XXX_desc-normapp_info |
.nii.gz |
The analytic phase projection map of the cardiac waveform, voxelwise minimum subtracted and normalized |
Always |
XXX_desc-apppeaks_hist |
.tsv.gz, .json |
Not sure |
Always |
XXX_desc-apppeaks_hist_centerofmass |
.txt |
Not sure |
Always |
XXX_desc-apppeaks_hist_peak |
.txt |
Not sure |
Always |
XXX_desc-slicerescardfromfmri_timeseries |
.tsv.gz, .json |
Cardiac timeseries at the time resolution of slice acquisition ((1/TR * number of slices / multiband factor |
Always |
XXX_desc-stdrescardfromfmri_timeseries |
.tsv.gz, .json |
Cardiac timeseries at standard time resolution (25.O Hz) |
Always |
XXX_desc-cardpulsefromfmri_timeseries |
.tsv.gz, .json |
The average (over time from minimum) of the cardiac waveform over all voxels |
Always |
XXX_desc-cardiaccyclefromfmri_timeseries |
.tsv.gz, .json |
The average (over a single cardiac cycle) of the cardiac waveform over all voxels |
Always |
XXX_desc-cine_info |
.nii.gz |
Average image of the fMRI data over a single cardiac cycle |
Always |
XXX_desc-cycleaverage_timeseries |
.tsv.gz, .json |
Not sure |
Always |
XXX_desc-maxphase_map |
.nii.gz |
Map of the average phase where the maximum amplitude occurs for each voxel |
Always |
XXX_desc-minphase_map |
.nii.gz |
Map of the average phase where the minimum amplitude occurs for each voxel |
Always |
XXX_desc-processvoxels_mask |
.nii.gz |
Map of all voxels used for analytic phase projection |
Always |
XXX_desc-vessels_map |
.nii.gz |
Amplitude of variance over a cardiac cycle (large values are assumed to be vessels) |
Always |
XXX_desc-vessels_mask |
.nii.gz |
Locations of voxels with variance over a cardiac cycle that exceeds a threshold (assumed to be vessels) |
Always |
XXX_desc-arteries_map |
.nii.gz |
High variance vessels with early maximum values within the cardiac cycle |
Always |
XXX_desc-veins_map |
.nii.gz |
High variance vessels with late maximum values within the cardiac cycle |
Always |
XXX_info |
.json |
Run parameters and derived values found during the run (quality metrics, derived thresholds, etc.) |
Always |
XXX_memusage |
.csv |
Memory statistics at multiple checkpoints over the course of the run |
Always |
XXX_runtimings |
.txt |
Detailed timing information |
Always |
Usage:
Hypersampling by Analytic Phase Projection - Yay!.
usage: happy [-h] [--cardcalconly] [--skipdlfilter] [--slicetimesareinseconds]
[--model MODELNAME] [--mklthreads NTHREADS] [--nprocs NPROCS]
[--numskip SKIP] [--motskip SKIP] [--motionfile MOTFILE]
[--motionhp HPFREQ] [--motionlp LPFREQ] [--nomotorthogonalize]
[--nomotderiv] [--motfiltorder N] [--discardmotionfiltered]
[--gmsfilt] [--estweights WEIGHTSNAME] [--minhr MINHR]
[--maxhr MAXHR] [--minhrfilt MINHR] [--maxhrfilt MAXHR]
[--hilbertcomponents NCOMPS] [--envcutoff CUTOFF]
[--notchwidth WIDTH] [--invertphysiosign] [--teoffset TE]
[--cardiacfile FILE[:COL]]
[--cardiacfreq FREQ | --cardiactstep TSTEP]
[--cardiacstart START] [--forcehr BPM]
[--respirationfile FILE[:COL]]
[--respirationfreq FREQ | --respirationtstep TSTEP]
[--respirationstart START] [--forcerr BreathsPM]
[--spatialregression | --temporalregression] [--stdfreq FREQ]
[--outputbins BINS] [--gridbins BINS]
[--gridkernel {old,gauss,kaiser}] [--projmask MASKNAME]
[--projectwithraw] [--fliparteries] [--arteriesonly]
[--pulsatilitysigma SIGMAINMM] [--pulsatilitythreshold THRESHPCT]
[--version] [--detailedversion] [--processmask MASKNAME]
[--aliasedcorrelation] [--aliasedcorrelationwidth WIDTH]
[--upsample] [--estimateflow] [--noprogressbar]
[--wrightiterations WRIGHTITERATIONS] [--useoldvesselmethod]
[--infotag tagkey tagvalue] [--debug] [--nodetrend DETRENDORDER]
[--normvesselmap] [--noorthog] [--disablenotch] [--nomask]
[--nocensor] [--noappsmooth] [--nophasefilt] [--nocardiacalign]
[--saveinfoastext] [--saveintermediate] [--increaseoutputlevel]
[--decreaseoutputlevel] [--nompdetrend] [--nompphaseproject]
[--noprefillcongrid] [--nocongridcache] [--focaldebug]
fmrifilename slicetimename outputroot
Positional Arguments
- fmrifilename
The input data file (BOLD fmri file or NIRS text file)
- slicetimename
Text file containing the offset time in seconds of each slice relative to the start of the TR, one value per line, OR the BIDS sidecar JSON file.NB: FSL slicetime files give slice times in fractions of a TR, BIDS sidecars give slice times in seconds. Non-json files are assumed to be the FSL style (fractions of a TR) UNLESS the –slicetimesareinseconds flag is used.
- outputroot
The root name for the output files
Processing steps
- --cardcalconly
Stop after all cardiac regressor calculation steps (before phase projection).
Default:
False- --skipdlfilter
Disable deep learning cardiac waveform filter.
Default:
True- --slicetimesareinseconds
If a non-json slicetime file is specified, happy assumes the file is FSL style (slice times are specified in fractions of a TR). Setting this flag overrides this assumption, and interprets the slice time file as being in seconds. This does nothing when the slicetime file is a .json BIDS sidecar.
Default:
False- --model
Use model MODELNAME for dl filter (default is model_ppgattention_pytorch_w128_fulldata - from the revised NeuroImage paper.)
Performance
- --mklthreads
Use NTHREADS MKL threads to accelerate processing (defaults to 1 - more threads up to the number of cores can accelerate processing a lot, but can really kill you on clusters unless you’re very careful. Use at your own risk
Default:
1- --nprocs
Use NPROCS CPUs to accelerate processing (defaults to 1 - more CPUs up to the number of cores can accelerate processing a lot, but you need to remember to ask for this many CPUs on clusters.) Entering a mulitprocessor routine disables mklthreads (otherwise there’s chaos).
Default:
1
Preprocessing
- --numskip
Skip SKIP tr’s at the beginning of the fMRI file (default is 0).
Default:
0- --motskip
Skip SKIP tr’s at the beginning of the motion regressor file (default is 0).
Default:
0- --motionfile
Read 6 columns of motion regressors out of MOTFILE file (.par or BIDS .json) (with timepoints rows) and regress them, their derivatives, and delayed derivatives out of the data prior to analysis.
- --motionhp
Highpass filter motion regressors to HPFREQ Hz prior to regression.
- --motionlp
Lowpass filter motion regressors to LPFREQ Hz prior to regression.
- --nomotorthogonalize
Do not orthogonalize motion regressors prior to regressing them out of the data.
Default:
True- --nomotderiv
Do not use motion derivative regressors.
Default:
True- --motfiltorder
Include powers of each confound regressor up to order N. Default is 1 (no expansion).
Default:
1- --discardmotionfiltered
Do not save data after motion filtering.
Default:
True- --gmsfilt
Remove low frequency variation (experimental).
Default:
False
Cardiac estimation tuning
- --estweights
Generation of cardiac waveform from data will be restricted to voxels in MASKNAME and weighted by the mask intensity. If this is selected, happy will only make a single pass through the data (the initial vessel mask generation pass will be skipped).
- --minhr
Limit lower cardiac frequency search range to MINHR BPM (default is 40).
Default:
40.0- --maxhr
Limit upper cardiac frequency search range to MAXHR BPM (default is 140).
Default:
140.0- --minhrfilt
Highpass filter cardiac waveform estimate to MINHR BPM (default is 40).
Default:
40.0- --maxhrfilt
Lowpass filter cardiac waveform estimate to MAXHR BPM (default is 1000).
Default:
1000.0- --hilbertcomponents
Retain NCOMPS components of the cardiac frequency signal to Hilbert transform (default is 1).
Default:
1- --envcutoff
Lowpass filter cardiac normalization envelope to CUTOFF Hz (default is 0.4 Hz).
Default:
0.4- --notchwidth
Set the width of the notch filter, in percent of the notch frequency (default is 1.5).
Default:
1.5- --invertphysiosign
Invert the waveform extracted from the physiological signal. Use this if there is a contrast agent in the blood.
Default:
False- --teoffset
Specify the echo time in seconds. This is used when combining multiecho data. Default is 0.
External cardiac waveform options
- --cardiacfile
Read the cardiac waveform from file FILE. If COL is an integer, and FILE is a text file, use the COL’th column. If FILE is a BIDS format json file, use column named COL. If no file is specified, estimate the cardiac signal from the fMRI data.
- --cardiacfreq
Cardiac waveform in cardiacfile has sample frequency FREQ (default is 32Hz). NB: –cardiacfreq and –cardiactstep are two ways to specify the same thing.
Default:
-32.0- --cardiactstep
Cardiac waveform in cardiacfile has time step TSTEP (default is 1/32 sec). NB: –cardiacfreq and –cardiactstep are two ways to specify the same thing.
Default:
-32.0- --cardiacstart
The time delay in seconds into the cardiac file, corresponding to the first TR of the fMRI file (default is 0.0)
- --forcehr
Force heart rate fundamental detector to be centered at BPM (overrides peak frequencies found from spectrum). Usefulif there is structured noise that confuses the peak finder.
External respiration waveform options
- --respirationfile
Read the respiration waveform from file FILE. If COL is an integer, and FILE is a text file, use the COL’th column. If FILE is a BIDS format json file, use column named COL.
- --respirationfreq
Respiration waveform in respirationfile has sample frequency FREQ (default is 32Hz). NB: –respirationfreq and –respirationtstep are two ways to specify the same thing.
Default:
-32.0- --respirationtstep
Respiration waveform in respirationfile has time step TSTEP (default is 1/32 sec). NB: –respirationfreq and –respirationtstep are two ways to specify the same thing.
Default:
-32.0- --respirationstart
The time delay in seconds into the respiration file, corresponding to the first TR of the fMRI file (default is 0.0)
- --forcerr
Force respiratory rate fundamental detector to be centered at BreathsPM (overrides peak frequencies found from spectrum). Usefulif there is structured noise that confuses the peak finder.
Output processing - mutually exclusive options
Optional cardiac noise removal. You can do either spatial or temporal regression, but not both.
- --spatialregression
Generate framewise cardiac signal maps and filter them out of the input data.
Default:
False- --temporalregression
Generate voxelwise aliased synthetic cardiac regressors and filter them out of the input data.
Default:
False
Output options
- --stdfreq
Frequency to which the physiological signals are resampled for output. Default is 25.
Default:
25.0
Phase projection tuning
- --outputbins
Number of output phase bins (default is 32).
Default:
32- --gridbins
Width of the gridding kernel in output phase bins (default is 3.0).
Default:
3.0- --gridkernel
Possible choices: old, gauss, kaiser
Convolution gridding kernel. Default is kaiser
Default:
'kaiser'- --projmask
Phase projection will be restricted to voxels in MASKNAME (overrides normal intensity mask.)
- --projectwithraw
Use fMRI derived cardiac waveform as phase source for projection, even if a plethysmogram is supplied.
Default:
False- --fliparteries
Attempt to detect arterial signals and flip over the timecourses after phase projection (since relative arterial blood susceptibility is inverted relative to venous blood).
Default:
False- --arteriesonly
Restrict cardiac waveform estimation to putative arteries only.
Default:
False- --pulsatilitysigma
Split pulsatility map spatial frequencies with a gaussian filter with sigma=SIGMAINMM. Default is 6.0mm.
Default:
6.0- --pulsatilitythreshold
Consider voxels with a high spatial frequency pulsatility above THRESHPCT percent to be vessels. Default is 0.5mm.
Default:
6.0
Version options
- --version
Show simplified version information and exit
- --detailedversion
Show detailed version information and exit
Miscellaneous options.
- --processmask
Instead of dynamically generating a processing mask from the data, use the mask specified in the file MASKNAME. The mask must be a 3D NIFTI file with x, y, z dimensions matching the fMRI data.
- --aliasedcorrelation
Attempt to calculate absolute delay using an aliased correlation (experimental).
Default:
False- --aliasedcorrelationwidth
Width of the aliased correlation calculation (default is 5.0).
Default:
5.0- --upsample
Attempt to temporally upsample the fMRI data (experimental).
Default:
False- --estimateflow
Estimate blood flow using optical flow (experimental).
Default:
False- --noprogressbar
Will disable showing progress bars (helpful if stdout is going to a file).
Default:
True- --wrightiterations
Number of iterations for calculating Wright map. Set to 0 to disable.
Default:
0- --useoldvesselmethod
Will use the old method to find blood vessels.
Default:
False- --infotag
Additional key, value pairs to add to the info json file (useful for tracking analyses).
Debugging options (probably not of interest to users)
- --debug
Turn on debugging information.
Default:
False- --nodetrend
Disable data detrending.
Default:
3- --normvesselmap
Find vessels using normalized phase projection.
Default:
True- --noorthog
Disable orthogonalization of motion confound regressors.
Default:
True- --disablenotch
Disable subharmonic notch filter.
Default:
False- --nomask
Disable data masking for calculating cardiac waveform.
Default:
True- --nocensor
Bad points will not be excluded from analytic phase projection.
Default:
True- --noappsmooth
Disable smoothing app file in the phase direction.
Default:
True- --nophasefilt
Disable the phase trend filter (probably not a good idea).
Default:
True- --nocardiacalign
Disable alignment of pleth signal to fMRI derived cardiac signal.
Default:
True- --saveinfoastext
Save the info file in text format rather than json.
Default:
True- --saveintermediate
Save some data from intermediate passes to help debugging.
Default:
False- --increaseoutputlevel
Increase the number of intermediate output files.
Default:
0- --decreaseoutputlevel
Decrease the number of intermediate output files.
Default:
0- --nompdetrend
Disable multiproc detrending.
Default:
True- --nompphaseproject
Disable multiproc phase projection.
Default:
True- --noprefillcongrid
Don’t prefill the congrid value cache.
Default:
True- --nocongridcache
Disable the congrid value cache completely.
Default:
True- --focaldebug
Enable targeted additional debugging output (used during development).
Default:
False
Example:
Extract the cardiac waveform and generate phase projections
Case 1: When you don’t have a pleth recording
There are substantial improvements to the latest versions of happy. In the old versions, you actually had to run happy twice - the first time to estimate the vessel locations, and the second to actually derive the waveform. Happy now combines these operations interpolation a single run with multiple passes - the first pass locates voxels with high variance, labels them as vessels, then reruns the derivation, restricting the cardiac estimation to these high variance voxels. This gives substantially better results.
Using the example data in the example directory, try the following:
happy \ rapidtide/data/examples/src/sub-HAPPYTEST.nii.gz \ rapidtide/data/examples/src/sub-HAPPYTEST.json \ rapidtide/data/examples/dst/happytest
This will perform a happy analysis on the example dataset. To see the extracted cardiac waveform (original and filtered), you can use showtc (also part of them rapidtide package):
showtc \ rapidtide/data/examples/src/happytest_desc-slicerescardfromfmri_timeseries.json:cardiacfromfmri,cardiacfromfmri_dlfiltered \ --format separate
Case 2: When you DO have a pleth recording
If you do have a pleth recording, then by all mean use it. You simply tell happy where it is (it can be in pretty much any sort of texty kind of file - .txt, .csv, .tsv, .tsv.gz, properly BIDS compliant, etc.), and tell it the sample rate and start time (you don’t have to cut up a continuous physio recording). If the data is BIDS, you don’t need to specify samplerate and start time, since they are in the file already.:
happy \ --cardiacfile FILE[:COL] \ --cardiacfreq FREQ \ --cardiacstart STARTTIME
- Where:
FILE[:COL] is the regressor text filename - if multicolumn, use :COL to specify the column to use either numerically, with 0 being the first column, or by name if the file is a BIDS compliant .tsv file.
FREQ is the sample rate of the text file in Hz. Alternately, you can use –cardiactstep TSTEP to give the sample period - these are two ways to specify the same thing.
STARTTIME is the time delay into the file, in seconds, where the part of the waveform you want to use begins. This is useful if you have a long pleth recording and don’t want to have to chop it up into multiple files.
NB: If you have a high quality plethysmogram, then a lot of the constraints on happy go away. The requirement of a combination of TR, multiband factor and slice number that allows you to successfully extract a cardiac waveform aren’t really relevant to analytic phase projection. That should work regardless of any of those factors, so if you want to make a spiffy movie of pulsatility, and you have a pleth recording, you’re good to go.
Performance tuning
If your data is fairly high quality, then happy should Just Work without any tweaking. However, if you have bad motion, low SNR, longish TR or lowish multiband factors, it may need a little help finding the signal, so you can play with the many command line parameters that limit the heartrate search range, etc. to see a list of options, along with fairly self explanatory descriptions, type:
happy --help
Multiecho data
Does this work on multiecho data? That’s a very good question, and I’m glad you asked me that!
Yes, it does. There are a couple of things to consider though. The first is that as far as I can tell, the majority of the cardiac signal that you see in fMRI data is a non-BOLD signal, which is to say that one of the cool things about multiecho data, the fact that you can separate BOLD and non-BOLD data doesn’t really help you, since if you use tedana or ME-ICA on your data, you’re going to strip out the signal you’re looking for, which would be bad. But in fact, that’s not really an issue, since you want to work on the fully unprocessed data anyway, so happy works on data long before you’ve even gotten to the multiecho processing stage.
Which brings us to the next question - there’s additional timing information from the fact that each echo is recorded at a slightly different time - how do we smartly combine the data from all the echoes? The answer for now is - you don’t, because I’m not that smart. What I can confirm, through my limited testing, is that in a 60 slice, TR=1.33s MB=4 4 echo acquisition (a trendy set of parameters in the circles I travel in) you can extract the cardiac waveform from _any_ of the individual echoes, although the first echo has the highest SNR - it seems that the short echo time enhances the cardiac variation. The noise in the echoes seems largely uncorrelated (other than motion effects), so combining echoes is almost certainly a good idea. That said, averaging doesn’t make things look noticeably better. Maybe use PCA to extract the most important timecourse from the set of derived cardiac regressors? Or maybe actually using the timing information of the individual echoes to enhance the time resolution? I await your PR implementing this with bated breath, dear reader…
BTW if you use the cardiac noise removal steps below, you should apply them to each echos data separately prior to doing the multiecho processing. That works fine with an fMRIprep workflow, since happy processing would be done before you start fMRIprep (replace the individual raw echo data in the source directory).
Cardiac noise removal
Since we know the phase of the plethysmogram signal at any given time over the course of the experiment, and what the waveform looks like across the cardiac cycle in every voxel, this suggests that we should be able to model the cardiac signal at every voxel in the brain at every timepoint, and then remove it from the data. To this end, I’ve implemented a cardiac noise removal strategy in happy (well, two actually).
The first is a standard temporal linear regression of the cardiac signal on the data. We know the phase of the plethysmogram
signal at every time during the fMRI acquisition, and we know what the signal looks like in each voxel at any given
cardiac phase (on average), so we can calculate the expected value of the cardiac signal in each voxel at each
timepoint of the fMRI acquisition. This won’t necessarily look like a cardiac signal, since it’s sampled at way
below the Nyquist rate, so it’s extremely aliased, but
if we have the timing right, that shouldn’t really matter - we only care what the value is at the time the fMRI data
is sampled, which we can calculate, and we can regress the signal out of each voxel
using the --temporalregression option. What this fails to account for is the fact that the amplitude of the cardiac
signal varies over the length of the experiment - or at least the amplitude of the plethysmogram does. The signal in
the brain seems to vary far less then in the fingertip. So maybe not a big deal.
The second approach is to do a spatial regression - for each timepoint, we simply calculate the cardiac signal in
each voxel given our known phase of the plethysmogram and the analytic projection at each voxel at that phase, and regress
that spatial pattern out of the image at that timepoint. This allows the absolute amplitude of the signal to vary
at each timepoint, so it might(?) work better. This allows you to ignore the concept of aliasing, since you’re just
saying “there will be a certain pattern of response across the brain at each timepoint”. To do this, use
--spatialregression.
Which of these is better? I haven’t done enough systematic testing to say for sure, so I don’t know. My gut reaction was that spatial regression was going to be better, since it allows the size of the cardiac response to vary over time even if we don’t know the form of the variation (and also, since you’ll get a timecourse of signal amplitude over time, it gives you a way to measure the variation over time, which probably tells you something interesting about sympathetic nervous system function [3]).
Be that as it may, it appears, in my (extremely limited) testing, that in the current implementation, temporal
regression works substantially better. It certainly removes a lot more variance than the spatial regression, so that’s
a strong argument for it. The R2 of the temporal fits are pretty high, indicating in some voxels that it’s taking out
around 50% of the variance. It’s certainly worth getting to the bottom of this, but for now,
use --temmporalregression.
WHOCARES
If you want to lean hard into cardiac noise removal, you should also look at WHOCARES [4] (https://github.com/gferrazzi/WHOCARES). It’s a more sophisticated approach to cardiac noise removal that uses happy as a starting point, but adds a number of additional steps to improve the cardiac regression. You would use this instead of, rather than in addition to, happy. Happy’s regression has improved somewhat since the original WHOCARES paper; I’m not sure what the relative performance is now.