System and method for template-based automated detection of anatomy

Abstract

To automatically detect or identify boundaries of an anatomy.SOLUTION: A system and method for detecting an anatomy includes, for each training subject of a plurality of training subjects, generating an initial anatomical template based on a corresponding MR image and a first training subject of the plurality of training subjects. A computing device may map the MR images of other training subjects onto the template space by applying a global transformation followed by a local transformation, may average the mapped MR images with the initial anatomical template to generate a final anatomical template, and may delineate the boundaries of the anatomy of interest in the final anatomical template. The computing device may fine-tune the boundaries using an edge detection algorithm. The final anatomical template can be used to automatically identify the boundaries of the anatomy of interest in an untrained subject.SELECTED DRAWING: Figure 1

Description

Technical Field

[0001] Cross - Reference to Related Applications This application claims the benefit and priority of U.S. Provisional Patent Application No. 63 / 170,229, filed on April 2, 2021, entitled "SYSTEMS AND METHODS FOR AUTOMATIC TEMPLATE - BASED DETECTION OF ANATOMICAL STRUCTURES", the content of which is hereby incorporated by reference in its entirety.

[0002] The present disclosure generally relates to the field of detecting or identifying the boundaries of anatomical structures. Specifically, the present disclosure relates to methods and systems for generating a standardized template and using such a template to automatically detect or identify the boundaries of anatomical structures.

Background Art

[0003] One disease that requires the identification of the boundaries of anatomical structures is Parkinson's disease (PD), which is a chronic progressive neurodegenerative disorder that affects approximately 1% of individuals over the age of 60. PD is pathologically characterized by early neurodegeneration of the substantia nigra (NM) in the substantia nigra pars compacta (SNpc) and an increase in iron deposition in the substantia nigra (SN). Degeneration of the SN is characteristic of the progression of many neurodegenerative diseases. In atypical Parkinsonian disorders, including progressive supranuclear palsy (PSP) and multiple system atrophy (MSA), extensive neuronal loss in the SNpc also occurs, but in these disorders, various sub - regions of the SN are affected.

[0004] The SN is composed of two anatomically and functionally distinct regions, the substantia nigra pars reticulata (SNpr) and the substantia nigra pars compacta (SNpc). While the SNpc contains a dense distribution of NM that houses dopaminergic neurons, the iron content tends to be higher in the SNpr. However, clusters of SNpc dopaminergic neurons (known as nigrosomes) are deeply embedded within the SNpr, and thus, the boundary between the SNpr and SNpc is difficult to delineate, particularly in the caudal region of the SN. The regional selectivity of PD is relatively specific to a 50% - 70% loss of pigmented neurons in the ventrolateral layer of the SNpc (at symptom onset). Regarding the SN, iron deposition and volume changes in the red nucleus (RN) and subthalamic nucleus (STN) have been reported to be associated with disease status and progression rate, and also function as important targets for deep brain stimulation (DBS) therapy in PD patients.

Summary of the Invention

[0005] According to at least one aspect, a magnetic resonance imaging (MRI) system comprises an MRI scanner configured to acquire magnetic resonance (MR) data, at least one processor, and a memory having computer code instructions stored thereon. When executed by the at least one processor, the computer code instructions cause the at least one processor to acquire, for each of a plurality of training subjects via the MRI scanner, a corresponding MR image with contrast exemplifying one or more anatomical structures of interest. The at least one processor can generate an initial anatomical template image based on the first MR data of a first training subject among the plurality of training subjects. The initial anatomical template image can define a template space. For each of the plurality of training subjects other than the first training subject, the at least one processor can (i) apply a global transformation to the MR image of the training subject to generate a first morphed version of the MR image representing a first estimate of the MR data of the training subject in the template space, and (ii) apply a local transformation to the first morphed version of the MR image of the training subject to generate a second morphed version of the MR image representing a second estimate of the MR data of the training subject in the template space. The at least one processor can average the initial anatomical template image and the second morphed versions of the MR images of the plurality of training subjects to generate a final anatomical template image. The at least one processor can delineate the boundaries of one or more anatomical structures of interest in the final anatomical template and use the final anatomical template to identify the boundaries of one or more anatomical structures of interest of other non-training subjects.

[0006] According to at least one aspect, the method can include acquiring, via an MRI scanner, for each of a plurality of training subjects, a corresponding MR image with a contrast that exemplifies one or more anatomical structures of interest. The method can include a computing device generating an initial anatomical template image based on first MR data of a first training subject among the plurality of training subjects. The initial anatomical template image can define a template space. For each of the plurality of training subjects other than the first training subject, the computing device can (i) apply a global transformation to the MR image of the training subject to generate a first morphed version of the MR image representing a first estimate of the MR data of the training subject in the template space, and (ii) apply a local transformation to the first morphed version of the MR image of the training subject to generate a second morphed version of the MR image representing a second estimate of the MR data of the training subject in the template space. The method can include averaging the initial anatomical template image with the second morphed versions of the MR images of the plurality of training subjects other than the first training subject to generate a final anatomical template image. The method can include depicting boundaries of one or more anatomical structures of interest in the final anatomical template and using the final anatomical template to identify boundaries of one or more anatomical structures of interest of other non-training subjects.

[0007] According to at least one aspect, a non-transitory computer-readable medium can include computer code instructions stored on the non-transitory computer-readable medium. When executed by a processor, the computer code instructions can cause the processor to obtain, for each of a plurality of training subjects via an MRI scanner, corresponding MR images with contrast exemplifying one or more anatomical structures of interest. The processor can generate an initial anatomical template image based on the first MR data of a first training subject among the plurality of training subjects. The initial anatomical template image can define a template space. For each of the plurality of training subjects other than the first training subject, the processor can (i) apply a global transformation to the MR image of the training subject to generate a first morphed version of the MR image representing a first estimate of the MR data of the training subject in the template space, and (ii) apply a local transformation to the first morphed version of the MR image of the training subject to generate a second morphed version of the MR image representing a second estimate of the MR data of the training subject in the template space. The processor can average the initial anatomical template image with the second morphed versions of the MR images of the plurality of training subjects other than the first training subject to generate a final anatomical template image. The processor can delineate boundaries of one or more anatomical structures of interest in the final anatomical template and use the final anatomical template to identify boundaries of one or more anatomical structures of interest of other non-training subjects.

Brief Description of the Drawings

[0008]

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7A

Figure 7B

Figure 8A

Figure 8B

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 14

Figure 15

Figure 16

Figure 17

Figure 18

[0009] For some diseases such as Parkinson's disease (PD), information regarding the disease state and / or information useful in some treatments can be inferred from the volume or other geometric properties of some anatomical structures. In the case of PD, regarding the substantia nigra (SN), it is known that iron depositions and volume changes in the red nucleus (RN) and the subthalamic nucleus (STN) are associated with the disease status and progression rate, and also function as important targets for deep brain stimulation (DBS) treatment in PD patients. Therefore, an accurate and comprehensive in vivo depiction of the SN and sub-regions of the SN, RN, and STN may be useful for fully investigating changes in the composition of iron and neuromelanin (NM) in PD and other movement disorders affecting the midbrain. Therefore, an accurate in vivo depiction of the SN, sub-regions of the SN, and other midbrain structures such as the red nucleus (RN) and the subthalamic nucleus (STN) may be useful for fully investigating changes in iron and NM in PD.

[0010] So far, many studies still define the deep gray matter of the brain using manual or semi-automatic approaches. However, manual segmentation is time-consuming, especially when it is necessary to evaluate a large amount of data. Also, unless the evaluator is well-trained, manual drawing has low reproducibility reliability between individuals or between regions. Some approaches for in vivo imaging can include the use of templates for mapping iron and NM content. In such approaches, creating a standardized template can (a) recognize changes in the distribution of iron and NM, (b) automatically calculate the volume associated with such a distribution, (c) quantify the changes in iron content and NM signal, and (d) have a significant impact on the reliability of these measurements. Anatomical templates of the SN using conventional structural MR sequences have been previously used in some studies. High-resolution probabilistic in vivo subcortical nucleus atlases can be created using T1-weighted (T1W) and T2-weighted (T2W) images, and these images can now segment the SN into the substantia nigra pars compacta (SNpc) and the substantia nigra reticulata (SNpr). Such atlases may not be suitable for studies targeting the elderly when based on images from a young adult population (e.g., age: 28.9 ± 3.6 years, mean ± standard deviation). Furthermore, both T1W and T2W contrast images cannot be easily used to depict the highly intermeshed boundary between the SNpc and the SNpr in the elderly. Based on the anatomical connections of the SN subregions to different parts of the brain, diffusion-based tractography can be used to segment the SN into the SNpc and the SNpr, or to subdivide the SN / ventral tegmental area (VTA) into dorsal medial and ventral lateral subregions. However, the presumed structural connectivity has well-known biases and highly depends on data acquisition and fiber tracking algorithms, and diffusion imaging is troubled by low-resolution data acquisition. Therefore, using high-resolution imaging for direct visualization and segmentation of the SN and its subregions would be a more accurate option, especially with regard to the detection of subtle pathological changes in the SN.

[0011] One approach for determining the boundaries of the midbrain nuclei and examining the pathological changes in PD patients can be obtained based on the use of T2* weighted gradient echo (GRE) imaging and susceptibility weighted imaging (SWI), in which iron-containing regions appear as low intensity. On the other hand, the development of quantitative susceptibility mapping (QSM) enables the quantification of iron stored in ferritin and hemosiderin. QSM-based techniques have shown that tissue magnetic susceptibility correlates well with brain iron in PD patients. Additionally, in preoperative target guidance for DBS, QSM is superior to traditional T2* weighted imaging. However, QSM alone cannot separate the SNpc from the SNpr because both contain iron. This limitation can be resolved by using neuromelanin-sensitive MRI (NM-MRI), which has been developed over the past few years. The SNpc and ventral tegmental area (VTA) are mainly composed of dopaminergic neurons containing NM, while the SNpr is not. Therefore, the high-intensity signal seen on NM-MRI in the midbrain is spatially associated with the regions of the SNpc and VTA (this has been verified by postmortem histological studies). Therefore, the overlap between the volume of NM (the sum of the SNpc and VTA) and the volume of the iron-containing SN (the sum of the SNpc and SNpr) is thought to represent the SNpc. * weighted gradient echo (GRE) imaging and susceptibility weighted imaging (SWI), in which iron-containing regions appear as low intensity. On the other hand, the development of quantitative susceptibility mapping (QSM) enables the quantification of iron stored in ferritin and hemosiderin. QSM-based techniques have shown that tissue magnetic susceptibility correlates well with brain iron in PD patients. Additionally, in preoperative target guidance for DBS, QSM is superior to traditional T2* * weighted imaging. However, QSM alone cannot separate the SNpc from the SNpr because both contain iron. This limitation can be resolved by using neuromelanin-sensitive MRI (NM-MRI), which has been developed over the past few years. The SNpc and ventral tegmental area (VTA) are mainly composed of dopaminergic neurons containing NM, while the SNpr is not. Therefore, the high-intensity signal seen on NM-MRI in the midbrain is spatially associated with the regions of the SNpc and VTA (this has been verified by postmortem histological studies). Therefore, the overlap between the volume of NM (the sum of the SNpc and VTA) and the volume of the iron-containing SN (the sum of the SNpc and SNpr) is thought to represent the SNpc.

[0012] So far, existing approaches do not accurately and reliably segment the SNpc using QSM and / or NM-MRI imaging. The SN atlas can be obtained using only QSM. Alternatively, an NM template based on NM-MRI imaging, created using, for example, manual delineation, automatic segmentation, or artificial intelligence, can be employed. The purpose of using the template would be to facilitate boundary identification. However, template mapping is not perfect, and the ideal space for marking the boundary is in the original pristine data space. Whether it is the template space (a standard fixed brain or volunteer into which all other brains are mapped) or the original space (the original brain imaging data for a given individual), simple thresholding methods have drawbacks. Setting the threshold to a relatively high or relatively low value can lead to dramatic changes in the estimated NM or iron content, especially in PD patients with severe NM degeneration and iron deposition in the SN. Also, the variable contrast of the images can make it difficult to use some algorithms such as region growing. Despite the SN having a high iron content compared to the surrounding regions, the iron is not uniformly distributed, and there is a reduction in iron in the nigrosome 1 (N1) territory. Structural gaps such as the N1 territory can also increase the difficulty of automatically segmenting the regions of interest (ROIs) of the SN and NM.

[0013] In the present disclosure, a system and method for detecting anatomical structures can include using high-resolution imaging from a single multi-echo NM MRI sequence to create both NM and iron templates and calculating the boundaries of each structure in the template space. Specifically, using both NM and iron templates from a single high-resolution NM MRI sequence, (i) calculate the boundaries of each structure in the template space, (ii) map these boundaries back to the original space, and then (iii) use a dynamic programming algorithm (DPA) to fine-tune the boundaries in the original space to match the details of the NM and SN characteristics of each individual. All experimental results regarding the NM of the SN and the iron content of the SN, STN, and RN show strong agreement in terms of Dice values and volume ratios, with the former being 0.85, 0.87, 0.75, and 0.92 respectively, and the latter being 0.99, 0.95, 0.89, 1.05 respectively. These high-quality results indicate that it is possible to measure tissue attributes such as volume, iron content, and neuromelanin content with accuracy. Multiple sequences sensitive to NM and iron can also be used, in which case co-registration between the two may be required.

[0014] Referring to FIG. 1, a flowchart is shown that illustrates a method 100 for generating an anatomical template in accordance with the concepts of the present invention of the present disclosure. The method 100 can include a computing device or system obtaining a corresponding MR image for each of a plurality of training subjects (step 102), and generating an initial anatomical template based on a selected MR image of a first training subject among the plurality of training subjects (step 104). The method 100 can include a computing device or system mapping the MR images of other training subjects onto a template space defined by the initial anatomical template by applying a global transformation followed by a local transformation (step 106). The method 100 can include averaging the mapped MR images and the initial anatomical template to generate a final anatomical template (step 108). The method 100 can include drawing the boundaries of the anatomical structures of interest on the final anatomical template (step 110). The method 100 can include using the final anatomical template to identify the boundaries of the anatomical structures of interest in a non-training subject (step 112).

[0015] The method 100 can include a computing device or system obtaining a corresponding MR image for each of a plurality of training subjects (step 102). The training subjects can be selected based on several predefined criteria such as age, health status, or the pathology and / or medical history of each subject. The studies described herein were approved by the local ethics committee, and all subjects signed an informed consent. Eighty-seven healthy subjects (mean age: 63.4 ± 6.2 years, range: 45 - 81 years, 53 females) were recruited from the region by advertisement. Exclusion criteria for potential training subjects included (a) structural abnormalities such as tumors, subdural hematomas, or contusions due to previous head trauma, (b) a history of stroke, poisoning, neurological or psychiatric disorders, and (c) one or more major vascular diseases with large-volume white matter lesions (e.g., Fazekas grade III).

[0016] MR images of the training subjects can be acquired via an MRI scanner. For example, in the study conducted, MR imaging was performed on a 3T Ingenia scanner (Philips Healthcare, Netherlands) using a 15-channel head array coil. The imaging parameters of a 3D gradient echo SWI sequence using an operating magnetization transfer contrast (MTC) pulse were: time to echo (TE) = 7.5 ms, ΔTE = 7.5 ms with a total of 7 echoes, repetition time (TR) = 62 ms, flip angle = 30°, pixel bandwidth = 174 Hz / pixel, matrix size = 384×144, slice thickness = 2 mm, number of slices = 64, and in-plane resolution interpolated to 0.67×1.34 mm 2 in-plane resolution interpolated to 0.67×0.67 mm 2 including a sensitivity factor of 2, elliptical sampling of k-space, and a total scan time equal to 4 minutes and 47 seconds. The magnetization transfer (MT) resonance state high-frequency pulse used a set of 3-block pulses each with a nominal flip angle of 90°, zero frequency offset, and a duration of 1.914 milliseconds (ms). The minimum allowable TR was used based on considerations of specific absorption rate safety. Due to this long repetition time of 62 ms, 7 echoes were collected.

[0017] In some implementations, a computing device or system can depict NM content using a first echo of an MTC-SWI magnitude image (TE = 7.5 ms) since the echo provides a major MT contrast. The computing device or system can evaluate iron deposition in the SN using a second echo (TE = 15 ms), or a combination of QSM images from two or more echoes, for QSM reconstruction. The computing device or system can create a susceptibility map by (i) segmenting the brain using a brain extraction tool, BET, (ii) unwrapping the original phase data using a 3D phase unwrapping algorithm (3DSRNCP), (iii) removing unwanted background fields using spectral harmonic artifact reduction (SHARP), and (iv) using a truncated k-space division (TKD)-based inverse filtering technique with an iterative approach for reconstructing the final QSM map.

[0018] When tracing the boundaries of the anatomical structures of interest, manual region of interest (ROI) segmentation can be employed. To measure the NM and iron content, the regions of interest (ROIs) of the NM, SN, RN, and STN can be manually traced by a single evaluator on the size of the MTC magnified 4 times and the QSM map using SPIN software (SpinTech, Inc., Bingham Farms, MI, USA). The NM-based SN boundary can be traced starting from the last tail slice of 4-5 slices until the NM disappears when 2-mm-thick slices are used (or over a total thickness of 8-10 mm). The iron-based SN boundary can be traced starting from one slice below the most cranial slice where the subthalamic nucleus is visible and continuing for 4-6 consecutive slices up to the most caudal slice or over a total thickness of 8-12 mm. The RN ROI can be outlined starting from the last tail slice and continuing for 3-4 slices like the cranium or over a total thickness of 6-8 mm. The STN ROI can be traced for two slices like the cranium or over a total thickness of about 4 mm. For all ROIs, a final boundary can be determined using a dynamic programming algorithm (DPA) to reduce subjective bias. Then, all these boundaries can be reviewed one by one by a second evaluator and appropriately corrected in agreement with the first reviewer.

[0019] It should be understood that the above-described data acquisition parameters and / or data processing scenarios represent exemplary implementations provided for illustrative purposes, not for limiting purposes. For example, the slice numbers provided above are specific to a given resolution with a slice thickness equal to 2 mm. If 1-mm-thick slices were used, all these slice numbers would double. Further, although the description mainly focuses on the application to PD and anatomical structures of the midbrain, the methods and techniques described herein can also be applied to other uses or other anatomical regions such as the deep cerebellar nuclei or the globus pallidus and caudate nucleus, or any other region of interest.

[0020] In some implementations, a computing device or system can cause an MR scanner to perform MR data acquisition, e.g., starting from original whole-brain 64-slice data and performing 7.5 ms echo time MTC data acquisition, for each of a plurality of training subjects. The MR scanner can construct a corresponding three-dimensional (3D) image (or a plurality of corresponding two-dimensional (2D) images) based on the corresponding acquired MR data for each training subject. When acquiring MR data for a plurality of training subjects, the MR scanner can be configured with parameters selected to generate MR images with contrast that exemplifies one or more anatomical structures of interest (e.g., regions of the SN, RN, and / or STN). The MR scanner can acquire NM-MRI data (e.g., NM-MRI images) and / or QSM data (QSM images) for each of the plurality of training subjects.

[0021] Method 100 can include a computing device or system generating an initial anatomical template based on a selected MR image of a first training subject among a plurality of training subjects (step 104). The computing device or system (or its user) can select a first 3D MR image (or corresponding 2D image) of the first training subject from among the MR images of the plurality of training subjects. The selection can be random or can follow predefined preferences or criteria. When generating the initial anatomical template, the computing device or system can magnify the first 3D image (or corresponding 2D image) in the plane by a factor of 2, 4, or more, for all of the acquired slices of the first training subject, e.g., according to a desired final resolution. This step can also include interpolation in the slice selection direction to obtain isotropic resolution in the template space. The template space can be defined at a resolution higher than the desired resolution or the resolution of the acquired MR data. In some implementations, the computing device or system can generate an initial NM-MRI template and an initial QSM template based on the NM data and QSM data of the first training subject, respectively. As will be discussed in more detail below, the computing device or system can use the initial anatomical template to map MR images of other training subjects into the template space.

[0022] Method 100 can include a computing device or system mapping MR images of other training subjects onto a template space defined by an initial anatomical template by applying a global transformation followed by a local transformation (step 106). The computing device or system can apply the global transformation and the local transformation to each of the images of other training subjects (other than the first training subject). The computing device or system can magnify each MR image of the MR images in-plane with a magnification factor (e.g., a factor of 2, 4, or more) for all slices and then apply the global transformation over a series of slices covering the region of interest within each MR image of other training subjects other than the first training subject. For example, for the midbrain, if the slice thickness is 2 mm, the series of slices covering the region of interest can be a set of 50 central slices. The global transformation can include a rigid transformation and, for example, an affine transformation following the application of B-spline interpolation using Insight Segmentation and Registration Toolkit freeware (ITK). The computing device or system can apply the same global transformation to map QSM data for each of the other training subjects to the initial QSM template.

[0023] The computing device or system can apply the global transformation to each MR image of other training subjects (e.g., QSM images and / or NM-MRI images) to generate a corresponding first morphed version of the MR image representing the corresponding first estimate of the initial template in the template space. That is, the global transformation is used to match each MR image of other training subjects (other than the first training subject) with the initial anatomical template in the template space. However, the global transformation typically does not provide an exact match with the anatomical structures in the initial anatomical template.

[0024] To improve the matching between the transformed MR image and the initial anatomical template in the template space, a computing device or system can apply a local transformation to a first morphed version of the MR images of other training subjects (other than the first training subject). The computing device or system can crop the first morphed version of the MR image in the plane so as to cover the midbrain and into approximately 16 slices (for example, when the slice thickness is 2 mm) to ensure coverage of the midbrain territory. The computing device or system can apply a local transformation to the cropped image volume to generate a second morphed version of the MR images of other training subjects. For each MR image of a training subject (other than the training subject), the corresponding second morphed version represents a better match with the initial anatomical template in the template space.

[0025] As an exemplary implementation, a total of 26 training subjects can be used. The computing device or system can select the MR data of one of these training subjects (referred to herein as the first training subject) to generate an initial anatomical template. Then, the computing device or system can apply a global transformation and a local transformation to the MR data of each of the other 25 training subjects, and after applying the global transformation and the local transformation, map the MR data of each of the other 25 training subjects to the initial anatomical template in the template space.

[0026] Method 100 can include averaging the mapped MR image and the initial anatomical template to generate a final anatomical template (step 108). A computing device or system can average the mapped MR image with the initial anatomical template to generate a final anatomical template. The result is an averaged template defined for, for example, 16 slices that include the midbrain territory. The averaged anatomical template has more representative boundaries of the anatomical structure of interest compared to an initial anatomical template generated based on the MR data of a single training subject, or compared to a template generated by averaging a first morphing version and the initial anatomical template. The computing device or system can linearly interpolate the averaged anatomical template in the slice selection direction, for example, to create a template having an isotropic resolution of 0.167 (192 slices in total).

[0027] Method 100 can include delineating the boundaries of the anatomical structure of interest on the final anatomical template (step 110). Tracing of the boundaries of the anatomical structure (or region) of interest can be performed manually, for example, by a physician or a radiologist. The final anatomical template generated according to steps 102-108 of method 100 shows enhanced contrast of the anatomical structure of interest and thus enables a more accurate tracing of the boundaries of such a structure.

[0028] In the final anatomical template, the mean value can be regarded as indicating a probability map for finding the boundaries. In some implementations, two or more final anatomical templates can be generated. For example, one QSM template and one NM-MRI template can be generated based on NM data and QSM data, respectively. In some implementations, the boundaries of NM, the red nucleus (RN), the substantia nigra (SN), and the subthalamic nucleus (STN) can all be drawn manually. For the neuromelanin (NM) data, slices 44 to 98 from a total of 192 interpolated slices can be used to draw the boundaries, while QSM data slices 44 to 126 can be used. The actual selection of the slice numbers depends on the resolution and degree of interpolation used.

[0029] Method 100 can further include using a boundary detection algorithm to fine-tune the boundaries of the anatomical structures of interest in the final anatomical template. The boundary detection algorithm can be implemented as a dynamic programming algorithm (DPA). A computing device or system can run (or execute) the DPA for boundary detection to determine the template boundaries. The DPA can use a cost function that depends on the local radius of curvature and the signal gradient. Further details of this algorithm are provided below.

[0030] Method 100 can include a computing device or system using the final anatomical template to identify the boundaries of the anatomical structures of interest in an untrained subject (step 112). The use of the final anatomical template to identify the boundaries of the anatomical structures of interest in an untrained subject is discussed in detail below with respect to Figure 2.

[0031] Referring now to FIG. 2, a flowchart is shown that illustrates a method 200 for the automatic template-based identification of characteristics of an anatomical structure in accordance with the concepts of the present invention of the present disclosure. Method 200 can include obtaining an MR image (or MR data) of an untrained subject (step 202). Method 200 can include mapping the MR image of the untrained subject to a final anatomical template within a template space by applying a global transformation and a subsequent local transformation (step 204). Method 200 can include projecting boundaries of one or more anatomical structures of interest from the final anatomical template onto the mapped MR image of the untrained subject (step 206). Method 200 can also include applying an inverse transformation to the projected boundaries of the anatomical structures of interest to generate an estimated value of the boundaries in the MR image of the untrained subject (step 208).

[0032] Method 200 enables the automatic identification of the boundaries of the anatomical structures of interest of an untrained subject using the final anatomical template. The method 100 for generating the anatomical template can involve human intervention (e.g., manually drawing the boundaries in step 110), but the method 200 for the automatic identification of the boundaries of the anatomical structures of interest for an untrained subject can be automatically performed / executed without any human intervention.

[0033] Method 200 can include the computing device obtaining an MR image (or MR data) of an untrained subject (step 202). The MR scanner can obtain MR data (e.g., a 3D MR image or a plurality of 2D MR images) of the untrained subject in a manner similar to that discussed above with respect to step 102 of method 100 in FIG. 1. The untrained subject can be a person (e.g., a patient or a research subject) who does not belong to the plurality of trained subjects used in method 100 to generate the final anatomical template. The MR scanner can obtain NM data and / or QSM data of the untrained subject.

[0034] Method 200 can include mapping an MR image of a non-trained subject to a final anatomical template within a template space (step 204) by applying a global transformation and a subsequent local transformation. The global transformation and the local transformation can be similar to those described with respect to step 106 of method 100 in FIG. 1. A computing device or system can first apply a global transformation to an MR image of a non-trained subject (e.g., an NM image and a QSM image) to generate a first morphed version of the MR image of the non-trained subject.

[0035] After applying the global transformation, the computing device or system can crop the first morphed version of the MR image of the non-trained subject. For example, the computing device or system can label the RNs in two central slices of the template data and map them back to the original whole brain. If the RNs are still visible only in two slices, the lowest slice is set to slice 10. If the RNs are in three slices, the middle slice is set to slice 10. This provides a means for optimally centering the data prior to the final transformation back to the template space. The computing device or system can magnify the cropped image volume by a factor of four (or another magnification factor) in the plane, and then apply a local transformation to the magnified cropped image volume to obtain a second morphed version of the MR image of the non-trained subject with an isotropic image (the transformation enables the image to be isotropic with interpolation in the slice selection direction). The second morphed version of the MR image of the non-trained subject represents a relatively good match to the final anatomical template in the template space. The interpolated in-plane resolution provides a more accurate means for estimating the final object size using the DPA algorithm.

[0036] Method 200 can include projecting the boundaries of one or more anatomical structures of interest from a final anatomical template onto a mapped MR image of a non-trained subject (step 206). Once the MR image of the non-trained subject is matched or mapped to the final anatomical template (by applying a global transformation and a local transformation), the computing device or system can project the boundaries of the anatomical structures of interest drawn on the final anatomical template onto a second morphed version of the MR image of the non-trained subject.

[0037] Method 200 can include applying an inverse transformation to the projected boundaries of the anatomical structures of interest to generate an estimate of the boundaries in the MR image of the non-trained subject (step 208). The inverse transformation is performed to map the template boundaries of the anatomical structures of interest back onto the midbrain (e.g., on the original image space) of the MR image of the non-trained subject. However, since everyone is different and the template mapping is not perfect, there is no guarantee that the projected boundaries will fit well.

[0038] Method 200 can further include a computing device or system using a boundary detection algorithm (or DPA) to fine-tune the boundary of an anatomical structure of interest in an original acquired MR image of a non-trained subject. The computing device or system can apply a threshold to the region inside the converted boundary for QSM data to remove negative values and can use image thinning to determine a corresponding centerline. To best select an initial starting point for DPA, both Otsu histogram analysis and a threshold-based approach can be used to determine whether the original boundary extends extremely outside the structure of interest. For NM data, the computing device or system can determine a background intensity and a constant equal to four times the background standard deviation averaged across all 25 trained subjects (or trained subjects other than the first trained subject) added to create a threshold below which signals are set to zero. For QSM data, the starting threshold can be set to zero. When the Otsu threshold gives a value that removes pixels inside the structure (or region) of interest, the template-converted boundary can be correspondingly shrunk. For QSM data, if the Otsu threshold exceeds 30 ppb, the threshold can be set to 30 ppb. Finally, the computing device or system can use the same DPA used in the template space (as discussed above with respect to FIG. 1) to correct the resulting boundary again. This DPA can prevent leakage of SN to RN together with the RN boundary and can distinguish STN from SN with the help of the template boundary.

[0039] Referring to FIG. 3, a diagram illustrating a system 300 for detecting or identifying the boundaries of anatomical structures according to the concepts of the present invention of the present disclosure. The system 300 can include an MR scanner 302 for acquiring MR data, and a computing device 304 for processing the acquired MR data to detect or identify the boundaries of anatomical structures. The computing device 304 can include a processor 306 and a memory 308. The memory 308 can store computer code instructions for execution by the processor 306. When executed by the processor 306, the computer code instructions can cause the computing device 304 to execute the method 100 and / or the method 200 discussed above. The computing device 304 can further include a display device 310 for displaying MR images and outputting the results of the method 100 and / or the method 200, or other data.

[0040] In some implementations, the computing device 304 can be an electronic device separate from the MR scanner 302. In such implementations, the computing device 304 may or may not be communicatively coupled to the MR scanner 302. For example, the MR data acquired by the MR scanner 302 may be transferred to the computing device 304 via a flash memory, a compact disc (CD), or other storage device. When the computing device 304 is communicatively coupled to the MR scanner 302, the MR data acquired by the MR scanner 302 can be transferred to the computing device 304 via any communication link between the MR scanner 302 and the computing device 304.

[0041] In some implementations, computing device 304 can be integrated within MR scanner 302. For example, processor 306, memory 308, and / or display device 310 can be integrated within MR scanner 302. In such implementations, MR scanner 302 can acquire MR data and perform method 100 and / or method 200.

[0042] In summary, for each of SN, RN, and STN, a boundary can be manually drawn first in the template space, and system 300 or computing device 304 can run DPA to fine-tune the boundary. System 300 or computing device 304 can map these boundaries to the original space, run DPA again to provide the final boundary, and make this a fully automated process. After the manual drawing is created, system 300 or computing device 304 can run DPA that fully automates the final boundary determination. Using the final boundary, system 300 or computing device 304 can calculate volume, signal intensity, sensitivity, and / or overlapping fractions.

[0043] System 300 or computing device 304 can calculate the overlap between the NM complex volume (the sum of SNpc and VTA) and the iron containing the SN volume (SNpc and SNpr) by overlaying two ROIs from the MTC data and QSM data, respectively. System 300 or computing device 304 can normalize the overlap by the iron containing the SN volume that creates the overall overlap fraction (OOF) measure, which is basically a measure of the SNpc fraction relative to the entire SN volume, as shown below. JPEG2025105622000002.jpg16169

[0044] The process of generating the final anatomical template described with respect to FIG. 1 above was evaluated using MR data of various control subjects. The evaluation was based on a comparison of the performance of the process of generating the final anatomical template against manually drawn boundaries of the anatomical structures of interest.

[0045] A total of 87 healthy controls (HC) were scanned. The SN, STN, and RN were manually traced for all 87 cases. Of these, 30 individuals (test dataset: including 17 males and 13 females, age range 66 ± 7.2 years) were used for the initial training of the template approach described above. Once all aspects of the algorithm were in place, the method was then verified in the next 57 cases (validation dataset: including 17 males and 40 females, age range 61.9 ± 5.0 years). Two metrics were used to evaluate the performance of the anatomical template generation process. These are the Dice similarity coefficients indicating the spatial overlap between the structures associated with the manual segmentation method and the template segmentation method, and the volume ratio (VR) of the structure from dividing the template volume by the volume of the manual segmentation.

[0046] Finally, all the data were combined to generate quantitative information regarding the structure volume, NM and iron content, and the overlap between the SN region and the NM region. The total iron content was calculated by summing the product of the volume and the average susceptibility of the structure over all slices where the structure was drawn. Similarly, the total NM content was obtained as a result from the sum of the products of the NM volume and the NM contrast over the corresponding slices.

[0047] Referring to FIG. 4, MR images are shown that depict the perspective of the NM in the NM-MRI template space and the perspectives of the SN, RN, and STN in the QSM template space. Specifically, images (a)-(c) illustrate various perspectives of NM402 in the NM-MRI template using the validation dataset. The boundary of NM402 is manually drawn. Images (d)-(f) illustrate the perspectives of SN404, RN406, and STN408. The boundaries of SN404, RN406, and STN408 in images (d)-(f) of FIG. 4 are manually drawn.

[0048] Referring to FIG. 5, an image is shown that exemplifies various steps for mapping the NM boundary from the template space to the original space. Specifically, boundaries from 16 slices taken for every 12th slice in the 0.167 mm isotropic template space are overlaid on the original 2 mm thick midbrain slices and shown for the NM data. Since the midbrain is only seen on these 4 slices, the boundary only appears on these 4 slices. Each column of images in FIG. 5 represents a separate slice. The topmost row of images shown as row A depicts images of various slices of the constructed neuromelanin template. The second row (from the top) shown as row B depicts the same images as row A but with the neuromelanin DPA boundary 502. The third row shown as row C depicts the overlaid template boundary 504 representing the boundary 502 transformed onto the original image 504. The bottommost row shown as row D depicts a midbrain image of the subject with the final boundary 506 of the NM obtained by fine-tuning the overlaid template boundary 504 using DPA.

[0049] Referring to FIG. 6, an image is shown that illustrates various steps for mapping the boundaries of SN, RN, and STN from the template space to the original image space for QSM data. Each column represents a different slice. The first (topmost) row shown as row A shows an image of various slices of the constructed QSM template. The second row shown as row B shows the same image with the QSM DPA boundaries for SN602, RN604, and STN606. The third row shown as row C shows the boundaries converted to the original space for SN602, RN604, and STN606. The final row shown as row D shows an image of various midbrain slices of a subject with the final boundaries of SN602, RN604, and STN606 obtained after applying DPA to the boundaries in row C.

[0050] Regarding the initial training of the template, 30 cases of MTC data and QSM data were processed. The integrity of the template automatic MTC background intensity has a slope related from one to the other of 0.99, and R 2 is demonstrated by the fact that it is 0.53 and the p-value is less than 0.001. The background value is important for properly thresholding the NM and iron content signals.

[0051] The measurement of the dice similarity coefficient and volume ratio (VR) was performed for different thresholds for both MTC images and QSM images, and the data associated with the NM and SN dice coefficients plotted against VR are shown in FIGS. 7A - 7B and FIGS. 8A - 8B, respectively. Referring to FIGS. 7A and 7B, plots of the dice similarity coefficient plotted against the VR value of neuromelanin (NM) are shown based on various scenarios. In both FIGS. 7A and 7B, plot A corresponds to the scenario where no threshold is applied, plot B corresponds to a threshold of NM contrast greater than 1000, plot C corresponds to a threshold of NM contrast greater than 1250, and plot D corresponds to a threshold of NM contrast greater than 1500. The plots in FIG. 7A are generated using data for 30 training cases, while the plots in FIG. 7B are generated using data for 57 validation cases. The average dice, volume ratio, and volume loss associated with the template data and manually drawn data are cited for each threshold. A threshold of 1000 units minimizes the volume loss and yields excellent results.

[0052] Referring to FIGS. 8A and 8B, plots of the dice similarity coefficient plotted against the VR value of SN are shown based on various scenarios. In both FIGS. 8A and 8B, plot A corresponds to the scenario where no threshold is applied, plot B corresponds to a threshold of sensitivity value greater than 50 ppb, plot C corresponds to a threshold of sensitivity value greater than 75 ppb, and plot D corresponds to a threshold of sensitivity value greater than 100 ppb. The plots in FIG. 8A are generated using data for 30 training cases, while the plots in FIG. 8B are generated using data for 57 validation cases. The average dice, volume ratio, and volume loss associated with the template data and manually drawn data are cited for each threshold. A threshold of 50 ppb results in a minimal loss in SN volume.

[0053] The higher the threshold, the denser the distribution becomes, and theoretically it will ultimately approach a single one on both axes. However, a higher threshold also causes a higher loss of volume. Therefore, there must be a trade-off between a sufficient dice and a volume ratio with volume loss. For the NM contrast using a threshold of 1000, the average volume loss of the template data for all cases is less than 10%, resulting in average dice and VR values of 0.90 and 0.95, respectively. Higher thresholds such as 1250 result in an average volume loss slightly exceeding 10% (average dice = 0.92, average VR = 0.96), and data showing an NM contrast exceeding 1500 resulted in relatively higher dice and VR (0.94 and 0.98, respectively) than lower thresholds, with an average loss of 20%. The threshold of 1000 seems to result in acceptable dice and VR values while keeping the volume loss below 10%. Similarly, for SN, as shown in FIGS. 8A - 8B, an average sensitivity threshold of 50 ppb resulted in an average volume loss of slightly less than 10%, while the average dice and VR values were 0.88 and 0.94, respectively.

[0054] Figure 9 shows a plot of the dice similarity coefficient plotted against the VR value of the red nucleus (RN) based on various scenarios. Plot A (left column) corresponds to the scenario where no threshold is applied, and plot B (right column) corresponds to the scenario where a threshold of 50 ppb is applied. The upper plots are generated using a dataset of 30 training cases, while the lower plots are generated using a dataset of 57 validation cases. The average dice, volume ratio, and volume loss associated with the template data and the manually drawn data for the selected threshold are cited within each plot. The threshold of 50 ppb limits the volume loss to approximately 10%.

[0055] Figure 10 shows plots of the Dice similarity coefficient plotted against the VR value of the subthalamic nucleus (STN) based on various scenarios. Plot A (in the left column) represents the scenario where no threshold is applied, while plot B (in the right column) represents the scenario with a 50 ppb threshold. The upper plots are generated using a dataset of 30 training cases, while the lower plots are generated using a dataset of 57 validation cases. The average Dice, volume ratio, and volume loss associated with the template data and manually drawn data for the selected threshold are cited within each plot. The 50 ppb threshold keeps the volume loss at approximately 10%.

[0056] Regarding Figures 9 and 10, before applying any threshold to the QSM data, the average Dice and volume ratio were 0.93 and 1.06 for the RN, and 0.76 and 0.95 for the STN. However, applying a 50 ppb threshold to the QSM data results in average Dice and VR of 0.95 and 1.04 for the RN, and 0.83 and 0.98 for the STN. An average VR greater than 1 indicates that most of the structures found by the fully automated template / DPA approach tend to be larger than the structures in the manual / DPA approach.

[0057] Referring to Figure 11, plots are shown that illustrate the correlations between the iron contents of the SN, RN, and STN resulting from manual segmentation and template segmentation. The correlation relationships between the iron contents of the SN, RN, and STN resulting from manual segmentation and template segmentation for 30 training cases are shown in the upper plots, while the correlation relationships between the iron contents of the SN, RN, and STN resulting from manual segmentation and template segmentation for 57 validation cases are shown in the lower plots. The ROIs associated with the SN, RN, and STN, and the final results of the template for the manually drawn iron content using a 50 ppb threshold are shown in various plots. For each of the SN, RN, and STN structures, the slope and R2 is shown for each of the plots in FIG. 11. The p-values are less than 0.001 for SN and RN and equal to 0.006 for STN.

[0058] For template verification, the validation dataset included 57 cases used in the processing of QSM data and MTC data. The Dice similarity coefficients were plotted for the VR values of both NM and SN, and these results are shown in FIGS. 7B and 8B, respectively. For the NM contrast using a threshold of 1000, the average volume loss of the template data for all cases was about 5%, yielding average Dice and VR values of 0.88 and 1.06, respectively. Similarly, for the average sensitivity of SN, a threshold of 50 ppb yielded an average volume loss of approximately 5%, while the average Dice and VR values were 0.89 and 0.97. As with the previous dataset, a threshold of 1000 for the NM contrast and a threshold of 50 ppb for the average sensitivity of SN yielded sufficient results regarding the average Dice, VR, and volume loss across 57 cases.

[0059] The Dice similarity coefficients for VR for RN and STN are shown in FIGS. 9 and 10, respectively. Before applying any threshold to the QSM data, the average Dice values and VR values for RN and STN were 0.92 and 1.05, and 0.74 and 0.86, respectively. Similar to the results from the previous dataset, applying a threshold improved these values. Using a threshold of 50 ppb for the QSM data yielded average Dice and VR of 0.95 and 1.03, and 0.81 and 0.90 for RN and STN, respectively, and an average template volume loss of about 12% for both structures.

[0060] FIG. 11 shows the correlation between the iron content of SN, RN, and STN for template data and manual data. For each structure, the corresponding slope and R 2is shown in FIG. 10. The p-values are less than 0.001 for SN and RN and equal to 0.008 for STN. Table 1 below summarizes the results associated with the estimated template volume, VR measure, and Dice similarity coefficient for each structure. Mean and standard deviation values for the first and second datasets and the merged data are shown.

Table 1

[0061] In the studies described with respect to FIGS. 4 - 12, NM images and QSM images derived from a single sequence (less than 5 minutes) were used for the automatic segmentation of SN, STN, and RN in the original space without the need to manually draw ROIs for non-template-based subjects. A multi-contrast atlas combined with DPA for boundary detection in both the template space and the original image after back-transformation from the template space is described and validated. Both Dice values and volume ratios showed excellent agreement for measures of volume and iron content between the automatic template approach and manual drawing.

[0062] Most existing templates do not include NM, SN, STN, and RN and do not have the resolution presented in this work. For example, while the MNI template is 1mm isotropic, the approach described herein uses interpolated in-plane isotropic data of 0.67mm and further interpolates it to an isotropic resolution of 0.167mm in 3D. In fact, it has been found that using whole-brain global deformable registration with an isotropic resolution of 0.67mm using either ANT or SpinITK does not result in a morphological transformation that can consistently and well match the shape of the midbrain structure across all subjects. Therefore, to solve this problem and lead to significantly improved results, a local registration (also referred to herein as local transformation) step was added. From this local transformation approach, templates for both iron and neuromelanin were created from a single multi-echo sequence.

[0063] Previous studies have attempted to segment the SN using structural MRI atlases based on T1 and / or T2 weighted sequences. However, the SN is a small structure in the midbrain with low contrast in these images, making it difficult to accurately define the SN boundary. To overcome this limitation, some other studies have attempted to obtain an SN atlas using either NM-MRI or QSM. By using a dynamic atlas composed of NM-enhanced brain images for automatic segmentation of the SN, one group showed a Dice value of less than 0.75. The NM-sensitive T1 weighted fast spin echo sequence applied in their study is not optimal for clearly depicting the SN. The contrast of NM can be significantly improved by using an MT-MRI acquisition sequence. To the inventors' knowledge, there have been no studies that have attempted to automatically segment the SN, SNpr, and SNpc using both NM-MRI and QSM.

[0064] Previous studies have segmented the SN by including either younger (over 5 years and under 18 years) or older (over 60 years) subjects. The use of atlases based only on young healthy subjects may not be appropriate for studying patients with neurodegenerative diseases that affect the midbrain structure. Age-related and disease-related morphometric changes between and within these subjects are important and can affect the success of using a template approach when localizing the SN. The most appropriate atlas for a given study needs to minimize global or local distortion from native subject space to atlas space. Therefore, an atlas for studying patients with neurodegenerative diseases was created using older HC subjects (mean age: 63.4 years). To reduce bias caused by inter-subject variability, some studies have used probabilistic atlases and attempted to construct a probabilistic atlas of the SNpc based on an NMS-MRI sequence to investigate microstructural abnormalities of the SNpc using diffusion MRI in PD patients. They applied symmetric diffeomorphic registration to align T1 images and SNpc masks of 27 HCs onto the MNI space and thresholded the atlas at a probability of 50%. However, the average Dice coefficient for the atlas versus human raters was less than 0.61. Thresholds can lead to dramatic changes in the depiction of the NM and thus reduce the reproducibility of the atlas.

[0065] Regarding the validation stage, the second dataset shows that the Dice similarity coefficient and VR values are very close to those from the first dataset for the structure, thereby confirming the consistency of the results generated by the automatic template processing approach proposed herein. Also, the fact that these values are very close to a single one indicates the excellent performance of this approach.

[0066] In the approach described herein, to further improve the validity of the NM atlas, DPA (implementing an edge detection algorithm) is adopted after creating a probabilistic atlas of the SNpc.

[0067] In conclusion, the new approach described herein for the fully automated segmentation of deep gray matter in the midbrain (or generally anatomical structures) uses the template approach as a guide but uses the original data along with a dynamic programming algorithm to fine-tune the boundaries. Using this approach, it is possible to quantify neuromelanin in the SN and the volume and iron content for all of the SN, STN, and RN. This approach should enable the study of changes in these imaging biomarkers for various neurodegenerative diseases without the need for manual tracing of these structures.

[0068] Regarding neuromelanin (NM) background measurements, since the background value is important for properly thresholding the NM signal, 30 MTC cases were processed to obtain neuromelanin (NM) background measurements. Referring to FIG. 12, a plot showing the agreement between the manual and template neuromelanin (NM) background measurements for 30 healthy controls (or training subjects) is shown. Specifically, the plot illustrates the completeness of the template automated MTC background intensity scale with respect to the manual drawing. The agreement between the manual and template MTC background scales shows a slope of 0.99, an R 2 , and a p-value of less than 0.001.

[0069] Regarding the DPA, the final template boundaries of all structures of interest can be determined using a dynamic programming algorithm (DPA) for boundary detection in both the template space and the original space. Exemplary detailed steps of the DPA algorithm can be as follows. 1. An initial boundary can be drawn on the template data for both NM and QSM. 2. Then, as shown in FIG. 13, the background can be drawn on the template image. FIG. 13 shows both the NM 702 and the background region 704. 3. Next, the boundary and background can be transformed into the original space as shown in FIG. 14. In FIG. 14, the NM region 802 and the background region 804 are shown after transformation from the template space to the original space. By using a value obtained by adding 1000 to the background value, the DPA boundary can be refined to provide faster convergence and it becomes easier to use a smaller and safer search radius to prevent the algorithm from leaking outside the original boundary. 4. The background average value + 4σ (where σ is the standard deviation of the background region of interest found to be approximately 250 units in this study) can be used as a threshold for the MTC data, and all points lower than this threshold can be removed to determine the NM refined boundary. For QSM data, pixels with a sensitivity less than 0 ppb can be removed before DPA is run for all of SN, STN, and RN. 5. The initial boundary can be smoothed using a 3×3 Gaussian filter. 6. Then, the centerline can be determined using a thinning method. 7. Boundary refinement before DPA: Since there may be some intensity variations around the structure, the global threshold may not be valid for the entire structure in the MTC data. Therefore, Otsu's method [Reference: 1979 IEEE “A Threshold Selection Method from Gray-Level Histograms”, Nobuyuki Otsu] can be applied to the filtered local image of a 40×40 pixel square around each single centerline point, and all points lower than the Otsu threshold can be removed. If the initial boundary is outside the boundary obtained by the Otsu threshold, the boundary can be corrected to the closest remaining point (referred to as the Otsu boundary). This step helps to speed up the convergence of DPA. 8. Implementation of DPA: For each point along the centerline, DPA can be executed in the corresponding search box. To enable a curved shape, the center of the centerline can be used to determine the initial ray. Then, DPA can be applied to the next set of rays by shifting one pixel along the centerline associated with the initial boundary. For each DPA iteration, the centerline can be updated. When the end point is reached, these rays can be swept through 180°. Then, the algorithm can return along the centerline until the centerline sweeps through another 180° and reaches the opposite end point. Finally, the center point can return along the centerline back to the starting point, closing the boundary. This process can be repeated 5 times to explore both the inside and outside of the boundary to obtain the best results. For NM, STN, and SN, after this step, only the inside search can be used to continue with another 5 iterations.

[0070] The search radius was limited to 4 pixels both inside and outside the boundary of the enlarged space, but the cost function used in this DPA includes a differential coefficient term and a radius of curvature term that avoid leakage to adjacent objects of the structure of interest. For points outside the boundary, when searching outward from the centerline, negative differential values were set to zero. This limitation prevented leakage into bright objects near the boundary. The cost function is given by JPEG2025105622000004.jpg1389where JPEG2025105622000005.jpg1540and the gradient and radius along the m-th ray at the r-th point are denoted by G(r,m) and R(r,m), respectively. The term G max represents the maximum differential coefficient inside the image, and R avgis the average radius over the previous three radii. The constant α represents the relative weighting of the differential coefficient term and the radius term, which can be set to a value between 0 and 1. The closer the shape of the structure of interest is to a circle, the higher the α can be set. If there are sharp edges, they can be smoothed by a large selection of α. Since all values of α from 0.05 to 0.15 played a similar role in faithfully finding the edges, the value of α = 0.1 was selected to be conservative.

[0071] Referring to FIG. 15, images are shown that illustrate the boundaries reconstructed with different values of the DPA parameter α. The selection of α can dramatically affect the final DPA boundary. The NM boundary in the upper left image was obtained using α = 0.05. In the upper right image, the NM boundary was determined using α = 0.10. In the lower left image, the NM boundary was obtained using α = 0.15. In the lower right image, the NM boundary was obtained using α = 0.20. Note that as the value of α decreases, it leads to less smoothing from the radius constraint, while as the value of α increases, it leads to more smoothing. The higher the value of α, the more rounded the structure becomes and the more it tends to shrink towards a circle. Low values lead to boundaries with very jagged edges. Finally, candidate points for the new boundary can be selected by maximizing the cost function described above for each ray.

[0072] For various shapes, some exemplary simulations are shown in FIGS. 16 - 18. Referring to FIG. 16, simulation results are shown for using DPA to identify the boundary of rectangular regions having two different intensities. Image (a) shows rectangular regions having two different intensities. Image (b) shows the rectangular region with a boundary drawn around the central area, and image (c) further shows the found center line and the updated boundary after applying DPA 5 times. Image (d) is similar to image (b), except that noise with a contrast - to - noise ratio (CNR) of 7:1 is added based on the difference in signal intensity between the central region and the frame region within the second rectangular boundary. Image (e) shows the final boundary found after 5 iterations of DPA that matches the correct area of 1400 pixels.

[0073] For the simulation corresponding to FIG. 16, the central region was set to 100 units, the outer frame region that is inside the second boundary but outside the first boundary was set to 30 units, and the background outside the second boundary was set to zero. Gaussian noise having a zero mean value and a standard deviation of 10 units was added to the image. The first example of a rectangle having sharply - defined corners and 1400 pixels is shown in FIG. 16. Despite the presence of noise, the boundary was still found completely.

[0074] Figure 17 shows simulation results for identifying the boundary of crescent-shaped regions with two different intensities using DPA. Image (a) shows the crescent-shaped regions with two different intensities, and image (b) shows that the same regions were the boundaries drawn around the central area of the crescent-shaped regions. Image (c) shows the crescent-shaped regions with the corresponding centerlines and updated boundaries after 5 iterations of DPA. Images (d) and (e) show the updates of the centerlines and boundaries after 10 iterations and 15 iterations of DPA, respectively. Image (f) shows the centerlines and boundaries after only 5 iterations when using the adaptive Otsu thresholding approach. Image (g) is similar to image (c) with added noise where the CNR is 7:1, and image (h) shows the centerlines and boundaries after only 5 iterations of DPA when using the adaptive Otsu thresholding approach.

[0075] The crescent-shaped regions in Figure 17 were selected to mimic a more difficult case of detecting the boundary of a curved object. Figure 17 shows the effect of running different numbers of iterations. It was found that the average value (image (f)) after running a total of 35 iterations was 990.5 with a standard deviation of 1.5, while when using the Otsu approach in the presence of a 10:1 SNR, the average value was 984 and the standard deviation was 2.0.

[0076] Figure 18 shows simulation results for identifying the boundary of cashew-shaped regions with two different intensities using DPA. Image (a) shows the cashew-shaped regions with two different intensities, and image (b) shows the cashew-shaped regions with the initial boundaries drawn around the central area. Image (c) shows the cashew-shaped regions with the corresponding centerlines and boundaries updated after 5 iterations of DPA when using the adaptive Otsu thresholding approach. Image (d) is similar to image (b) but with added noise and a CNR of 7:1. The final results of the boundaries and centerlines after 5 iterations of DPA when using the adaptive Otsu thresholding approach and in the presence of noise are shown in image (e).

[0077] Finally, the cashew-shaped region of FIG. 18 was selected, and SN was emulated to evaluate the method when there are no sharp edges. After running the first 35 iterations, the mean value was found to be 1086.5 and the standard deviation = 2.5. Using the Otsu approach in the presence of a 10:1 SNR and running an additional 30 iterations, the mean value was found to be 1097 and the standard deviation was 0.0.

[0078] One of ordinary skill in the art should understand that the processes described in this disclosure can be implemented using computer code instructions executable by a processor. The computer code instructions can be stored in a non-transitory or tangible computer-readable medium such as a memory. The memory can be random access memory (RAM), read only memory (ROM), cache memory, disk memory, any other memory, or any other computer-readable medium. The processes described in this disclosure can be implemented by an apparatus including at least one processor and / or a memory storing executable code instructions. The code instructions can cause the execution of any of the processes or operations described in this disclosure when executed by at least one processor. The apparatus can be, for example, an MRI scanner or a computing device.

Claims

**Claim 1** A magnetic resonance imaging (MRI) system comprising: an MRI scanner configured to acquire magnetic resonance (MR) data; at least one processor; a memory having computer code instructions stored thereon, wherein when the computer code instructions are executed by the at least one processor, the at least one processor is caused to: acquire, via the MRI scanner, for each training subject of a plurality of training subjects, a corresponding MR image with contrast exemplifying one or more anatomical structures of interest; generate an initial anatomical template image based on first MR data of a first training subject of the plurality of training subjects, the initial anatomical template image defining a template space; for each training subject of the plurality of training subjects other than the first training subject: apply a global transformation to the MR image of the training subject to generate a first morphed version of the MR image representing a first estimate of the MR data of the training subject in the template space; and apply a local transformation to the first morphed version of the MR image of the training subject to generate a second morphed version of the MR image representing a second estimate of the MR data of the training subject in the template space; average the initial anatomical template image with the second morphed versions of the MR images of the plurality of training subjects other than the first training subject to generate a final anatomical template image; depict the boundaries of the one or more anatomical structures of interest on the final anatomical template; and use the final anatomical template to identify the boundaries of the one or more anatomical structures of interest for other non-training subjects. A magnetic resonance imaging (MRI) system. **Claim 2** When determining the boundaries of the one or more anatomical structures of interest in the MR image, the at least one processor is configured to: use a boundary detection algorithm to fine-tune the boundaries of the one or more anatomical structures of interest in the final anatomical template. The MRI system according to claim 1. **Claim 3** When determining the boundary of the one or more anatomical structures of interest within the MR image, the at least one processor acquires, via the MRI scanner, an MR image of a non-trained subject with a contrast that exemplifies the one or more anatomical structures of interest; applies the global transformation to the MR image of the non-trained subject to generate a first morphing version of the MR image of the non-trained subject; applies the local transformation to the first morphing version of the MR image of the non-trained subject to generate a second morphing version of the MR image of the non-trained subject; projects the boundary from the final anatomical template of the one or more anatomical structures of interest onto the second morphing version of the MR image of the non-trained subject to determine an estimated value of the boundary of the one or more anatomical structures of interest within the template space for the non-trained subject; applies an inverse transformation to the estimation of the boundary within the template space of the one or more anatomical structures of interest to determine a first estimated value of the boundary of the one or more anatomical structures of interest within the MR image of the original non-trained subject, the MRI system according to claim 1, configured to perform the above.

4. The at least one processor is configured to use a boundary detection algorithm to fine-tune the first estimated value of the boundary of the one or more anatomical structures of interest within the original MR image of the non-trained subject and generate a second estimated value of the boundary of the one or more anatomical structures of interest within the second MR image of the non-trained subject, the MRI system according to claim 3.

5. The inverse transformation according to claim 3 includes a concatenation of the inverse of the local transformation followed by the inverse of the global transformation.

6. Applying the local transformation to the first morphing version of the MR image includes applying the local transformation to a cropping region of the first morphing version of the MR image, the MRI system according to claim 1.

7. The corresponding MR image includes a quantitative susceptibility map (QSM) or a neuromelanin image, the MRI system according to claim 1.

8. The MRI system according to claim 3, wherein the at least one processor is configured to determine a volume of at least one of the one or more anatomical structures of interest based on the boundaries of the one or more anatomical structures of interest in the original MR image of the untrained subject.

9. The MRI system according to claim 3, wherein the at least one processor is configured to determine an average intensity value of at least one of the one or more anatomical structures of interest based on the boundaries of the one or more anatomical structures of interest in the original MR image of the untrained subject.

10. The MRI system according to claim 3, wherein the at least one processor is configured to determine a total signal of at least one of the one or more anatomical structures of interest based on the boundaries of the one or more anatomical structures of interest in the second MR image.

11. A method comprising: acquiring, via an MRI scanner, for each of a plurality of trained subjects, a corresponding MR image with a contrast that exemplifies one or more anatomical structures of interest; generating, by a computing device, an initial anatomical template image based on first MR data of a first trained subject among the plurality of trained subjects, wherein the initial anatomical template image defines a template space; for each of the plurality of trained subjects other than the first trained subject, applying a global transformation to the MR image of the trained subject to generate a first morphed version of the MR image representing a first estimate of the MR data of the trained subject in the template space, and applying a local transformation to the first morphed version of the MR image of the trained subject to generate a second morphed version of the MR image representing a second estimate of the MR data of the trained subject in the template space; averaging the initial anatomical template image with the second morphed versions of the MR images of the plurality of trained subjects other than the first trained subject to generate a final anatomical template image; depicting the boundaries of the one or more anatomical structures of interest on the final anatomical template. A method comprising: using the final anatomical template to identify the boundaries of the one or more anatomical structures of interest for other non-trained subjects.

12. Determining the boundaries of the one or more anatomical structures of interest in the MR image, The method according to claim 11, comprising using a boundary detection algorithm to fine-tune the boundaries of the one or more anatomical structures of interest in the final anatomical template.

13. Determining the boundaries of the one or more anatomical structures of interest in the MR image, Acquiring an MR image of a non-trained subject with contrast exemplifying the one or more anatomical structures of interest via the MRI scanner; Applying the global transformation to the MR image of the non-trained subject to generate a first morphed version of the MR image of the non-trained subject; Applying the local transformation to the first morphed version of the MR image of the non-trained subject to generate a second morphed version of the MR image of the non-trained subject; Projecting the boundaries from the final anatomical template of the one or more anatomical structures of interest onto the second morphed version of the MR image of the non-trained subject to determine an estimated value of the boundaries of the one or more anatomical structures of interest in the template space for the non-trained subject; The method according to claim 11, comprising applying an inverse transformation to the estimated value of the boundaries in the template space of the one or more anatomical structures of interest to determine a first estimated value of the boundaries of the one or more anatomical structures of interest in the original MR image of the non-trained subject.

14. The method according to claim 13, further comprising using a boundary detection algorithm to fine-tune the first estimated value of the boundaries of the one or more anatomical structures of interest in the original MR image of the non-trained subject to generate a second estimated value of the boundaries of the one or more anatomical structures of interest in the second MR image of the non-trained subject.

15. The method according to claim 13, wherein the inverse transformation comprises a concatenation of the inverse of the local transformation followed by the inverse of the global transformation.

16. Applying the local transformation to the first morphing version of the MR image includes applying the local transformation to a cropping region of the first morphing version of the MR image. The method according to claim 11.

17. The method according to claim 11, wherein the corresponding MR image includes a quantitative susceptibility map (QSM) or a neuromelanin image.

18. The method according to claim 13, wherein the at least one processor is configured to determine a volume of at least one of the one or more anatomical structures of interest based on the boundaries of the one or more anatomical structures of interest in the original MR image of the untrained subject.

19. Determining an average intensity value of at least one of the one or more anatomical structures of interest based on the boundaries of the one or more anatomical structures of interest in the original MR image of the untrained subject, or The method according to claim 13, further comprising determining a total signal of at least one of the one or more anatomical structures of interest based on the boundaries of the one or more anatomical structures of interest in the second MR image.

20. A non-transitory computer-readable medium comprising computer code instructions stored on the non-transitory computer-readable medium, which, when executed by a processor, cause the processor to Acquire, via an MRI scanner, for each of a plurality of trained subjects, a corresponding MR image with a contrast exemplifying one or more anatomical structures of interest. Generate an initial anatomical template image based on the first MR data of a first trained subject among the plurality of trained subjects, wherein the initial anatomical template image defines a template space. For each of the plurality of trained subjects other than the first trained subject, Apply a global transformation to the MR image of the trained subject to generate a first morphing version of the MR image representing a first estimate of the MR data of the trained subject in the template space, and Apply a local transformation to the first morphing version of the MR image of the trained subject to generate a second morphing version of the MR image representing a second estimate of the MR data of the trained subject in the template space. Averaging the initial anatomical template image with the second morphing version of the MR images of the plurality of training subjects other than the first training subject to generate a final anatomical template image; Depicting the boundaries of the one or more anatomical structures of interest on the final anatomical template; A non-transitory computer-readable medium for causing the boundaries of the one or more anatomical structures of interest for other non-training subjects to be identified using the final anatomical template.

Images