Cross-sequence segmentation
Author: Xiao Han (see also: XiaoNotes)
This wiki pages describes the research effort for improving the subcortical segmentation accuracy for data acquired across different scanner platforms or imaging sequences.
The subcortical segmentation method within FreeSurfer is an atlas-driven method. It uses a probabilistic brain atlas to encode prior intensity models for individual brain structures at every atlas location, as well as prior models for the spatial distribution of these brain structures. As a result, the accuracy of the method often degrades when processing data acquired using scanner plaforms and/or imaging sequences that differ from the data used for the atlas training, due to changes in underlying image contrast. To resolve this problem and reduce the sensitivity of the method to changes in scanner platforms or acquisition parameters, we have developed an intensity renormalization procedure that automatically adjusts the prior intensity model of the atlas, in order to better fit it to the input data.
The intensity renormalization is performed for individual structures. The basic idea is to try to get an estimation of the intensity distribution of each structure in the input image, and then renormalize (a linear mapping) the atlas intensity model so that the global intensity distribution for each structure in the atlas matches that in the given image.
The intensity renormalization is currently implemented inside mri_ca_register (and mri_ca_label) and consists of a multi-linear image registration and a histogram matching step. It starts after the global linear registration of mri_em_register, which grossly aligns the atlas to the input image. This global linear registration has only 12 degrees of freedom, and as a result, the registration accuracy is much less affected by cross-sequence image contrast changes.
In the first step, a local linear registration is being sought for each individual struture. The registration result then defines a (rough) delineation of the particular structure on the input image, and its intensity distribution can be estimated. Since this linear registration is computed for each structure independently (although all use the global linear registration as the starting point), the final registration result for all the structures combined is no longer a single linear transform, thus comes the name "multi-linear" registration. The structure-wise linear registration works as follows: for an individual structure, a template (sub-)image is first constructed from the atlas using the maximum likelihood rule. That is, the label with the maximum prior probability at every atlas location is first found. If the label matches the structure being considered, the intensity value at the current location is set to the mean intensity of the particular label in the atlas. Otherwise, the intensity is set to zero. In practice, we build a maximum likely volume from the atlas first, where each voxel has an intensity value equal to the mean intensity of the maximum likely lable at that location. Then the template for a particualr structure can be built by masking out all other structures (i.e., reset their intensity values to zero). In practice, we dilate the initial mask first in order to keep a small neighborhood of non-zero values around the given structure, which aimes to provide sufficient contextual information to help the registration.
After the structure template is built, it need be registered to the input image. We first need to define a registration criterion. The criterion we use for the local registration is a variant of the standard SSE measure (i.e., the sum-of-squared-error or sum-of-squared-intensity-differences). From our problem setting, it is clear that the matching criterion has to take into account of possible image contrast differences. As a result, we modified the original SSE criterion to allow a (unknown) linear intensity mapping between the source and template images. The criterion is similar to what is used in the following paper:
S. Periaswamy and H. Farid, "Elastic Registration in the Presence of Intensity Variations", IEEE T. Med. Imag., 22(7):865-874, 2003
Of cause, other criteria such as the normalized correlation ratio can also be used, but the optimization process will become more complicated and time-consuming.
After the matching criterion or objective function is defined, various optimization algorithms can be applied to find the optimal linear registration (12 parameters) that optimizes (minimizes in our case) the objective function. One choice is the Powell's method, which is posed to find the global optimal solution. Unfortunately, we found that this method is not robust enough for the particular problem. We thus settled on a gradient-descent optimization approach, which finds a local optimal solution around the initial global linear mapping.
After the optimal registration is computed for a particular structure, its label in the atlas can then be transformed to the input image to define a mask for that structure. At this stage, the structure mask or segmentation cannot be very accurate, and thus the intensity normalization should rely on global features of the intensity distribution. In particular, two intensity histograms are built for the structure under consideration: one from the masked region in the input image, one from the atlas (using the intensity value of the maximum-likelihood volume for that structure). Next, a linear intensity transformation (a scaling plus an offset) is sought that maximizes the cross-correlation between the two histograms (currently, we force the offset to be zero and only optimize for the unknown scaling factor). The intensity transformation is then applied to renormalize the intensity model of the given structure at every atlas location. Note that the intensity for each structure at every atlas location is modeled as a Gaussian. The atlas intensity renormalization thus corresponds to a scaling of the mean and standard deviation of the Gaussian model.
After the intensity renormalization is finished for every structure, the new atlas created is then applied in the sequential steps of the subcortical segmentation procedure as in the original method. In a newer version, we also use the multi-linear registration results to re-initialize the non-linear atlas morphing instead of using the original global linear registration. The advantage of this new initilization approach is not yet conclusive, however.
We note that the key component for the atlas intensity renormalization procedure is the local structure-wise linear registration. This type of part-to-all registration is much harder than the whole image registration, and thus robust and accurate results are hard to obtain. Based on this consideration, we deliberately reduce the degree-of-freedom in our definition of the intensity renormalization problem. For example, as mentioned earlier, in the histogram matching, we only search for a scaling factor instead of both a scaling and an offset. We found that the histogram matching based scaling factor computation is not sensitive to registration errors.
Based on the same consideration, the local linear registration is computed only for a subset of structures (small structures like vessels are not considered; see the gca code for details). Based on the results from these structures, we can compute an average normalization factor for each of three major tissure classes (WM, GM, and CSF), and they will be assigned to the relevant remaining structures.
Another detail not mentioned earlier involves some constraints we added to detect and avoid obviously wrong registration or histogram matching results. For example, a registration result is dropped if the determinant of the registration matrix is too small or too big (to avoid unreasonable expansion or contraction of the structure). Also, a histogram matching result is dropped if the overlap in the final histograms is low. When this happens, the intensity normalization factor for the failed struture is set using the average normalization factor computed from other structures. For exmaple, the average GM normalization factor will be assigned to hippocampus if the initial linear registration and histogram matching for hippocampus failed.
Future directions:
More reliable structure-to-image registration method can be sought. One approach is to incorporate prior models of admissible structure registration, as a resort to improve the robustness of the registration method. The 1997 NeuroImage paper by Ashburner et al is a good reference for such type of approaches.