In the revision of the projection-based thickness (PBT) framework (_cat_volpbt), major parts of the code were reorganized and simplified. In addition, new ideas were implemented to improve the accuracy in general and the reconstruction of fine structures in particular: (i) the error-prone Eikonal thickness was replaced by a more robust approach, (ii) the simple thickness reconstruction of gyri was extended to reduce intersections of the outer surface, (iii) approximation of thickness for non-cortical voxels instead of complex local smoothing for continuous values, and (iv) to reduce blurring of fine sulci additional functions were introduced.
The supersimple parameter can be used to focus on major steps and switch off additional parts that have minor global but strong local effects.
The Figure shows the difference in several example brains, where the new pipeline supports smoother results (without the complex smoothing part of the original PBT function) and preserves fine anatomical details in sulci (e.g. callosal sulcus above the CC, but also small medial sulci of the swallow), but also more pronounced gyral structures.
Details
(1) Multi-level distance estimation
The PBT mapping function _cat_volpbtp.c (Dahnke et al., 2012) depends strongly on the accuracy of the distance maps. Since the simple voxel-based distance function _catvbdist.c uses only binary maps, it is quite sensitive to subtle intensity differences, smoothing, myelination, interpolation, small distances below 2 mm, and quantifying asymmetric structures. Although _cat_voleidist handles most of these problems, it is quite slow and relatively noisy in real data.
Another issue is the correct handling of myelinated GM values, which leads to an underestimation of cortical thickness by misinterpreting myelinated GM values as WM. The estimation is also quite sensitive to small local changes in intensity, i.e., if the values of the myelinated GM cross the threshold of 2.5 in the label map (Yp0 with 1=CSF, 2=GM, and 3=WM), several voxels there change their intensity, resulting in thickness changes of 1-2 mm for quite similar thresholds (e.g., 2.5 ± eps).
While trying to maintain _cat_volpbt and its subfunctions, the idea arises to measure the distance to a tissue class not only by a single estimate to the average threshold but by averaging several paired measurements to different thresholds. Using pairs allows compensation of the distance biases for lower/higher thresholds. For instance, the simplest case uses 1 pair at 0.25 and 0.75 and the tissue distance underestimation for the 0.25 cases is compensated by the distance overestimation of the 0.75 level. The 0.5 level itself is critical in the case of partial volume classes used by the AMAP.
Using multiple levels improves the accuracy logarithmically with a slight increase of computational needs (but faster than using higher resolutions or later corrections) with good values between 4 and 16 (opt.levels = 8).
However, more levels need to consider interpolation and smoothing effects that can cause overestimation of distances and therefore thickness with increased topology defects, especially in thinner regions. I.e., when we have an object with value 1 and background 0 and PVE in between, we will have paired levels for 0.5±r with 0<r<.5 in general but practical values between 0.2 and 0.4 (opt.range = 0.3).
To further stabilize the mapping, also the PVE values beyond the boundary were used with the same limit. However, values beyond this boundary (e.g. the CSF boundary for the WM distance map) have to be corrected for the additional distance to the second boundary (by a second distance map) to avoid overestimation of thickness.
As we use paired estimation (one with underestimation the other one with overestimation) with the same offset, it could be assumed that no further correction is required. However, a general offset correction by the PVE value itself and/or by half of the average intensity difference supports slight improvements (opt.correctoffset = 2).
In a simplified 1D view, the object boundary can be seen as a hill to which we want to measure the (air) distance (Figure 2A). If the sloop is uniform and there are no plateaus, the mean distance is a good approximation, but otherwise, the average of several measurements is essential. The plateau case is a side effect of the AMAP segmentation with a subclass for CSF-GM and GM-WM, which is one reason why the T1-based surface processing was previously chosen (another is the finer definition of fine CSF grades to detect asymmetries, which is practically less relevant so far).
However, the number of measurements is practically limited by (i) processing time (where _catvbdist is much faster than _cat_voleidist), (ii) image/segmentation quality, and (iii) possible interpolation artifacts with overshooting (Figure 2B). The effect of additional measurements is a logarithmic function, so the first iterations are the most effective, where we observed that 2 to 4 pairs represent a good compromise between speed and quality, before other issues (such as interpolation overshoot and other artifacts) become more prominent.
[ Figure 2 - how to measure the average distance to a hill, especially if there is a plateau (A) overshooting issue (B) ]
(2) Gyri reconstruction:
To further stabilize the result, we also use the inverse reconstruction of gyri instead of sulci by processing the inverse label map (4 - Yp0) and using the minimum thickness value (this typically results in lower gyral values). In addition, we correct for the hypothetical shift of the boundary (e.g., for a threshold of 0.25, a distance underestimation of 0.25 voxels is expected (again, this has less impact). Since PBT maps the maximum local value along the anterior to the WM, we used an extended range of voxels that are beyond the opposite tissue boundary (e.g., the CSF-GM boundary for the WM distance), but corrected for the additional distance.
Besides the reconstruction of blurred sulci also the reconstruction of blurred thin gyri can improve the resulting surface.
This can be achieved by (i) a simple independent estimation that finally uses only the minimum thickness or (ii) by a complex dependent estimation that first reconstructs the sulci, followed by the reconstruction of the gyri, and a final estimation step (opt.gyrusrecon = 2). However, the reconstruction of gyri is quite sensitive to blood vessels and requires removing extreme WM distance values that are detected by using a very smooth preliminary thickness map (opt.distcleanup = 1).
(3) Thickness smoothing:
Instead of smoothing the thickness map and a small extension (for stable downsampling), we now approximate the thickness values for the whole image using the _cat_volapprox function, which also includes a small amount of smoothing.
(4) Keeping details and sharpening:
For neocortical surfaces, a mesh resolution of about 1 mm is sufficient and supports fast processing at once. As the down-sampling of surface meshes caused sometimes topology issues and even MATLAB crashes, the default setting reconstructs the surface in the original (internal) resolution (e.g. 1 mm) rather than the interpolated PBT space (0.5 mm). However, the limited precision can result in topology defects those corrections can result in massive errors (removing of hole sulci/gyri). Two filters, keepdetails and sharpening, are used to preserve the blurring of fine structures by emphasizing extreme values in case of smoothing.
(-) Topology correction (inactive):
Since the surface-based topology correction requires re-meshing the surface, which requires further extensive re-meshing and reshaping, we tried to get by with the pure voxel-based topology correction (CS2 vs. CS3 pipeline). Although the surfaces were better in non-defective regions, the results of the voxel-based topology correction were not satisfactory. An additional soft WM topology correction, which corrects only small problems (typically by closing) was implemented in _cat_volpbtsimple (opt.WMTC) but did not improve the results satisfactorily.
Preliminary results
The Google table include preliminary tests of different parameter settings for 8 selected cases (Collins, HR075, 4397, OASIS031 (WMHs), BUSS02 (Child), BWPT (Phantom), NISALS_UTR (PVEs), ADHD200NYC (LowQuality/Motion)).
PBT simple surface reconstruction
In the revision of the projection-based thickness (PBT) framework (_cat_volpbt), major parts of the code were reorganized and simplified. In addition, new ideas were implemented to improve the accuracy in general and the reconstruction of fine structures in particular: (i) the error-prone Eikonal thickness was replaced by a more robust approach, (ii) the simple thickness reconstruction of gyri was extended to reduce intersections of the outer surface, (iii) approximation of thickness for non-cortical voxels instead of complex local smoothing for continuous values, and (iv) to reduce blurring of fine sulci additional functions were introduced.
The supersimple parameter can be used to focus on major steps and switch off additional parts that have minor global but strong local effects.
Details
(1) Multi-level distance estimation
The PBT mapping function _cat_volpbtp.c (Dahnke et al., 2012) depends strongly on the accuracy of the distance maps. Since the simple voxel-based distance function _catvbdist.c uses only binary maps, it is quite sensitive to subtle intensity differences, smoothing, myelination, interpolation, small distances below 2 mm, and quantifying asymmetric structures. Although _cat_voleidist handles most of these problems, it is quite slow and relatively noisy in real data.
Another issue is the correct handling of myelinated GM values, which leads to an underestimation of cortical thickness by misinterpreting myelinated GM values as WM. The estimation is also quite sensitive to small local changes in intensity, i.e., if the values of the myelinated GM cross the threshold of 2.5 in the label map (Yp0 with 1=CSF, 2=GM, and 3=WM), several voxels there change their intensity, resulting in thickness changes of 1-2 mm for quite similar thresholds (e.g., 2.5 ± eps).
While trying to maintain _cat_volpbt and its subfunctions, the idea arises to measure the distance to a tissue class not only by a single estimate to the average threshold but by averaging several paired measurements to different thresholds. Using pairs allows compensation of the distance biases for lower/higher thresholds. For instance, the simplest case uses 1 pair at 0.25 and 0.75 and the tissue distance underestimation for the 0.25 cases is compensated by the distance overestimation of the 0.75 level. The 0.5 level itself is critical in the case of partial volume classes used by the AMAP.
Using multiple levels improves the accuracy logarithmically with a slight increase of computational needs (but faster than using higher resolutions or later corrections) with good values between 4 and 16 (opt.levels = 8). However, more levels need to consider interpolation and smoothing effects that can cause overestimation of distances and therefore thickness with increased topology defects, especially in thinner regions. I.e., when we have an object with value 1 and background 0 and PVE in between, we will have paired levels for 0.5±r with 0<r<.5 in general but practical values between 0.2 and 0.4 (opt.range = 0.3). To further stabilize the mapping, also the PVE values beyond the boundary were used with the same limit. However, values beyond this boundary (e.g. the CSF boundary for the WM distance map) have to be corrected for the additional distance to the second boundary (by a second distance map) to avoid overestimation of thickness.
As we use paired estimation (one with underestimation the other one with overestimation) with the same offset, it could be assumed that no further correction is required. However, a general offset correction by the PVE value itself and/or by half of the average intensity difference supports slight improvements (opt.correctoffset = 2).
In a simplified 1D view, the object boundary can be seen as a hill to which we want to measure the (air) distance (Figure 2A). If the sloop is uniform and there are no plateaus, the mean distance is a good approximation, but otherwise, the average of several measurements is essential. The plateau case is a side effect of the AMAP segmentation with a subclass for CSF-GM and GM-WM, which is one reason why the T1-based surface processing was previously chosen (another is the finer definition of fine CSF grades to detect asymmetries, which is practically less relevant so far).
However, the number of measurements is practically limited by (i) processing time (where _catvbdist is much faster than _cat_voleidist), (ii) image/segmentation quality, and (iii) possible interpolation artifacts with overshooting (Figure 2B). The effect of additional measurements is a logarithmic function, so the first iterations are the most effective, where we observed that 2 to 4 pairs represent a good compromise between speed and quality, before other issues (such as interpolation overshoot and other artifacts) become more prominent.
[ Figure 2 - how to measure the average distance to a hill, especially if there is a plateau (A) overshooting issue (B) ]
(2) Gyri reconstruction:
To further stabilize the result, we also use the inverse reconstruction of gyri instead of sulci by processing the inverse label map (4 - Yp0) and using the minimum thickness value (this typically results in lower gyral values). In addition, we correct for the hypothetical shift of the boundary (e.g., for a threshold of 0.25, a distance underestimation of 0.25 voxels is expected (again, this has less impact). Since PBT maps the maximum local value along the anterior to the WM, we used an extended range of voxels that are beyond the opposite tissue boundary (e.g., the CSF-GM boundary for the WM distance), but corrected for the additional distance. Besides the reconstruction of blurred sulci also the reconstruction of blurred thin gyri can improve the resulting surface. This can be achieved by (i) a simple independent estimation that finally uses only the minimum thickness or (ii) by a complex dependent estimation that first reconstructs the sulci, followed by the reconstruction of the gyri, and a final estimation step (opt.gyrusrecon = 2). However, the reconstruction of gyri is quite sensitive to blood vessels and requires removing extreme WM distance values that are detected by using a very smooth preliminary thickness map (opt.distcleanup = 1).
(3) Thickness smoothing:
Instead of smoothing the thickness map and a small extension (for stable downsampling), we now approximate the thickness values for the whole image using the _cat_volapprox function, which also includes a small amount of smoothing.
(4) Keeping details and sharpening:
For neocortical surfaces, a mesh resolution of about 1 mm is sufficient and supports fast processing at once. As the down-sampling of surface meshes caused sometimes topology issues and even MATLAB crashes, the default setting reconstructs the surface in the original (internal) resolution (e.g. 1 mm) rather than the interpolated PBT space (0.5 mm). However, the limited precision can result in topology defects those corrections can result in massive errors (removing of hole sulci/gyri). Two filters, keepdetails and sharpening, are used to preserve the blurring of fine structures by emphasizing extreme values in case of smoothing.
(-) Topology correction (inactive):
Since the surface-based topology correction requires re-meshing the surface, which requires further extensive re-meshing and reshaping, we tried to get by with the pure voxel-based topology correction (CS2 vs. CS3 pipeline). Although the surfaces were better in non-defective regions, the results of the voxel-based topology correction were not satisfactory. An additional soft WM topology correction, which corrects only small problems (typically by closing) was implemented in _cat_volpbtsimple (opt.WMTC) but did not improve the results satisfactorily.
Preliminary results
The Google table include preliminary tests of different parameter settings for 8 selected cases (Collins, HR075, 4397, OASIS031 (WMHs), BUSS02 (Child), BWPT (Phantom), NISALS_UTR (PVEs), ADHD200NYC (LowQuality/Motion)).