Using Subject Space ROIs from Freesurfer#
An example using the AFQ API to find bundles as defined by endpoint ROIs from freesurfer. This example can be modified to work with ROIs in subject space from pipelines other than freesurfer.
import os.path as op
import nibabel as nib
import plotly
import numpy as np
from AFQ.api.group import GroupAFQ
import AFQ.data.fetch as afd
from AFQ.definitions.image import RoiImage
import AFQ.api.bundle_dict as abd
2026-05-26 22:58:36,475 INFO util.py:154 -- Missing packages: ['ipywidgets']. Run `pip install -U ipywidgets`, then restart the notebook server for rich notebook output.
Get some example data#
Retrieves High angular resolution diffusion imaging (HARDI) dataset from Stanford’s Vista Lab
see https://purl.stanford.edu/ng782rw8378 for details on dataset.
The data for the first subject and first session are downloaded locally (by default into the users home directory) under:
.dipy/stanford_hardi/
Anatomical data (anat) and Diffusion-weighted imaging data (dwi) are
then extracted, formatted to be BIDS compliant, and placed in the AFQ
data directory (by default in the users home directory) under:
AFQ_data/stanford_hardi/
This data represents the required preprocessed diffusion data necessary for initializing the GroupAFQ object (which we will do next)
The clear_previous_afq is used to remove any previous runs of the afq object stored in the AFQ_data/stanford_hardi/ BIDS directory. Set it to None if you want to use the results of previous runs. Setting it to “track” as here will only clear derivatives that depend on the tractography stage (i.e., bundle delination and tract profile calculation), as well as the tractography itself, to save time on recomputation. If you want to only clear derivatives that depend on bundle delineation, and keep the tractography, you can set clear_previous_afq to “recog” instead.
afd.organize_stanford_data(clear_previous_afq="track")
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[2], line 1
----> 1 afd.organize_stanford_data(clear_previous_afq="track")
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/AFQ/data/fetch.py:1819, in organize_stanford_data(path, clear_previous_afq)
1817 # fetches data for first subject and session
1818 logger.info("fetching Stanford HARDI data")
-> 1819 dpd.fetch_stanford_hardi()
1821 if path is None:
1822 if not op.exists(afq_home):
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/data/fetcher.py:494, in _make_fetcher.<locals>.fetcher(include_optional)
491 continue
492 files[str(n)] = (baseurl + f, md5_list[i] if md5_list is not None else None)
--> 494 fetch_data(files, folder, data_size=data_size, use_headers=use_headers)
496 if msg is not None:
497 logger.info(msg)
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:201, in warning_for_keywords.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
194 # Check if the current version is within the warning range
195 if (
196 version.parse(from_version)
197 <= version.parse(current_version)
198 <= version.parse(until_version)
199 ):
200 # Convert positional to keyword arguments and issue a warning
--> 201 return convert_positional_to_keyword(func, args, kwargs)
203 # If the version is greater than the until_version,
204 # pass the arguments as they are
205 elif version.parse(current_version) > version.parse(until_version):
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/testing/decorators.py:192, in warning_for_keywords.<locals>.decorator.<locals>.wrapper.<locals>.convert_positional_to_keyword(func, args, kwargs)
182 warnings.warn(
183 f"Pass {positionally_passed_kwonly_args} as keyword args. "
184 f"From version {until_version} passing these as positional "
(...) 187 stacklevel=3,
188 )
190 return func(*positional_args, **corrected_kwargs)
--> 192 return func(*args, **kwargs)
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/data/fetcher.py:397, in fetch_data(files, folder, data_size, use_headers, raise_on_error)
395 logger.info(f"From: {url}")
396 try:
--> 397 _get_file_data(fullpath, url, use_headers=use_headers, stored_md5=md5)
398 successful_downloads += 1
399 except (FetcherError, Exception) as e:
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/site-packages/dipy/data/fetcher.py:262, in _get_file_data(fname, url, use_headers, timeout, max_retries, stored_md5)
260 with open(fname, "wb") as data:
261 if response_size is None:
--> 262 copyfileobj(opener, data)
263 else:
264 copyfileobj_withprogress(opener, data, response_size)
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/shutil.py:203, in copyfileobj(fsrc, fdst, length)
201 fsrc_read = fsrc.read
202 fdst_write = fdst.write
--> 203 while buf := fsrc_read(length):
204 fdst_write(buf)
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/http/client.py:478, in HTTPResponse.read(self, amt)
475 return b""
477 if self.chunked:
--> 478 return self._read_chunked(amt)
480 if amt is not None and amt >= 0:
481 if self.length is not None and amt > self.length:
482 # clip the read to the "end of response"
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/http/client.py:608, in HTTPResponse._read_chunked(self, amt)
605 self.chunk_left = chunk_left - amt
606 break
--> 608 value.append(self._safe_read(chunk_left))
609 if amt is not None:
610 amt -= chunk_left
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/http/client.py:648, in HTTPResponse._safe_read(self, amt)
641 """Read the number of bytes requested.
642
643 This function should be used when <amt> bytes "should" be present for
644 reading. If the bytes are truly not available (due to EOF), then the
645 IncompleteRead exception can be used to detect the problem.
646 """
647 cursize = min(amt, _MIN_READ_BUF_SIZE)
--> 648 data = self.fp.read(cursize)
649 if len(data) >= amt:
650 return data
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/socket.py:719, in SocketIO.readinto(self, b)
717 raise OSError("cannot read from timed out object")
718 try:
--> 719 return self._sock.recv_into(b)
720 except timeout:
721 self._timeout_occurred = True
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/ssl.py:1304, in SSLSocket.recv_into(self, buffer, nbytes, flags)
1300 if flags != 0:
1301 raise ValueError(
1302 "non-zero flags not allowed in calls to recv_into() on %s" %
1303 self.__class__)
-> 1304 return self.read(nbytes, buffer)
1305 else:
1306 return super().recv_into(buffer, nbytes, flags)
File /opt/hostedtoolcache/Python/3.13.13/x64/lib/python3.13/ssl.py:1138, in SSLSocket.read(self, len, buffer)
1136 try:
1137 if buffer is not None:
-> 1138 return self._sslobj.read(len, buffer)
1139 else:
1140 return self._sslobj.read(len)
KeyboardInterrupt:
Generate left thalamus ROI from freesurfer segmentation file#
Load the segmentation file that was generated by Freesurfer for the specific subject.
Identify the left thalamus within the file, which has the label number 41
Create a Nifti image representing the left thalamus ROI:
Assign a value of 1 to the voxels that Freesurfer has labeled as 41 (i.e., the left thalamus).
Assign a value of 0 to all other voxels. This binary mask format is the expected input for pyAFQ when dealing with subject space ROIs. If it’s already in binary format, there is no need to do this step.
freesurfer_subject_folder = op.join(
afd.afq_home, "stanford_hardi",
"derivatives", "freesurfer",
"sub-01", "ses-01",
"anat")
seg_file = nib.load(op.join(
freesurfer_subject_folder, "sub-01_ses-01_seg.nii.gz"))
left_thal = seg_file.get_fdata() == 41
nib.save(
nib.Nifti1Image(
left_thal.astype(np.float32),
seg_file.affine),
op.join(
freesurfer_subject_folder,
"sub-01_ses-01_desc-leftThal_mask.nii.gz"))
# Fetch LV1 ROI
# which was already generated using the process above
afd.fetch_stanford_hardi_lv1()
Set tractography parameters (optional)#
We make this tracking_params which we will pass to the GroupAFQ object
which specifies that we want 10,000 seeds randomly distributed
only within the endpoint ROIs and not throughout the white matter.
This is controlled by passing
"seed_mask": RoiImage() in the tracking_params dict.
We only do this to make this example faster and consume less space.
tracking_params = dict(n_seeds=10000,
random_seeds=True,
rng_seed=42,
seed_mask=RoiImage(use_endpoints=True))
Define custom BundleDict object#
In a typical BundleDict object, ROIs are passed as paths to Nifti files.
Here, we define ROIs as dictionaries instead, containing BIDS filters.
Then pyAFQ can find the respective ROI for each subject and session.
bundles = abd.BundleDict({
"L_OR": {
"start": {
"scope": "freesurfer",
"suffix": "mask",
"desc": "leftThal"},
"end": {
"scope": "freesurfer",
"suffix": "anat",
"desc": "LV1"
},
"cross_midline": False,
"space": "subject"
}})
Initialize a GroupAFQ object:#
Creates a GroupAFQ object, that encapsulates tractometry, passing in our custom bundle info. Then we run the pipeline and generate a visualization of the bundle we found.
myafq = GroupAFQ(
bids_path=op.join(afd.afq_home, 'stanford_hardi'),
dwi_preproc_pipeline='vistasoft',
t1_preproc_pipeline='freesurfer',
tracking_params=tracking_params,
bundle_info=bundles)
bundle_html = myafq.export("indiv_bundles_figures")
plotly.io.show(bundle_html["01"]["L_OR"])