System and method for template-based automatic detection of anatomical structures - Patents.com

Abstract

The system and method for detecting an anatomical structure includes generating an initial anatomical template for each of a plurality of training subjects based on a corresponding MR image and a first training subject of the plurality of training subjects. The computing device can map the MR images of the other training subjects onto the template space by applying a global transformation followed by a local transformation, can average the mapped MR images with the initial anatomical template to generate a final anatomical template, and can draw a boundary of the anatomical structure of interest on the final anatomical template. The computing device can fine-tune the boundary using an edge detection algorithm. The final anatomical template can be used to automatically (e.g., without human intervention) identify the boundary of the anatomical structure of interest in the non-training subjects.

Description

[Technical field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims priority to and the benefit of U.S. Provisional Patent Application No. 63 / 170,229, filed April 2, 2021, and entitled “SYSTEMS AND METHODS FOR AUTOMATIC TEMPLATE-BASED DETECTION OF ANATOMICAL STRUCTURES,” the contents of which are incorporated by reference in their entirety herein.

[0002] The present disclosure relates generally to the field of detecting or identifying anatomical boundary. In particular, the present disclosure relates to a method and system for generating standardized templates and using such templates to automatically detect or identify anatomical boundary. [Background technology]

[0003] One disease that requires the identification of anatomical boundaries is Parkinson's disease (PD), 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 increased iron deposition in the substantia nigra (SN). Degeneration of the SN is a hallmark of the progression of many neurodegenerative diseases. In atypical parkinsonian disorders, including progressive supranuclear palsy (PSP) and multiple system atrophy (MSA), widespread neuronal loss of the SNpc also occurs, although different subregions of the SN are affected in these disorders.

[0004] The SN is composed of two anatomically and functionally distinct regions, the SN pars reticularis (SNpr) and the SNpc. The SNpc contains a dense distribution of NMs containing dopaminergic neurons, while 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 therefore the border between the SNpr and SNpc is difficult to delineate, especially in the caudal region of the SN. The regional selectivity of PD is relatively specific with a 50%-70% loss (at the onset of symptoms) of pigmented neurons in the ventral lateral layer of the SNpc. For the SN, iron deposition and volumetric changes in the red nucleus (RN) and subthalamic nucleus (STN) have been reported to be associated with disease status and progression rate, and also serve as important targets for deep brain stimulation (DBS) treatment in PD patients. Summary of the Invention

[0005] According to at least one aspect, a magnetic resonance imaging (MRI) system includes an MRI scanner configured to acquire magnetic resonance (MR) data, at least one processor, and a memory having computer code instructions stored thereon. The computer code instructions, when executed by the at least one processor, cause the at least one processor to acquire, via the MRI scanner, for each training subject of a plurality of training subjects, a corresponding MR image with contrast illustrating one or more anatomical structures of interest. The at least one processor can generate an initial anatomical template image based on a first MR data of a first training subject of the plurality of training subjects. The initial anatomical template image can define a template space. For each training subject of the plurality of training subjects other than the first training subject, the at least one processor may (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 may 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 may 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 in the other non-training subjects.

[0006] According to at least one aspect, the method can include acquiring, via an MRI scanner, for each training subject of the plurality of training subjects, a corresponding MR image with contrast illustrating one or more anatomical structures of interest. The method can include a computing device generating an initial anatomical template image based on a first MR data of a first training subject of the plurality of training subjects. The initial anatomical template image can define a template space. For each training subject 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 version of the MR image of the plurality of training subjects other than the first training subject to generate a final anatomical template image. The method may include delineating boundaries of one or more anatomical structures of interest in a final anatomical template, and using the final anatomical template to identify boundaries of one or more anatomical structures of interest in another untrained subject.

[0007] According to at least one embodiment, a non-transitory computer-readable medium may include computer code instructions stored on the non-transitory computer-readable medium. The computer code instructions, when executed by a processor, may cause the processor to acquire, for each training subject of the plurality of training subjects, a corresponding MR image with contrast illustrating one or more anatomical structures of interest via an MRI scanner. The processor may generate an initial anatomical template image based on a first MR data of a first training subject of the plurality of training subjects. The initial anatomical template image may define a template space. For each training subject of the plurality of training subjects other than the first training subject, the processor may (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 may average the initial anatomical template image with second morphed versions of the MR images of multiple training subjects other than the first training subject to generate a final anatomical template image. The processor may 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 in the other non-training subjects. [Brief description of the drawings]

[0008] [Figure 1] 1 is a flow chart illustrating a method for generating an anatomical template in accordance with the inventive concepts of the present disclosure. [Diagram 2] 1 is a flow chart illustrating a method for automatic template-based identification of anatomical structure characteristics in accordance with the inventive concepts of the present disclosure. [Diagram 3]1 shows a diagram illustrating a system for detecting or identifying boundaries of anatomical structures in accordance with the inventive concepts of the present disclosure. [Figure 4] Magnetic resonance (MR) images depicting the NM perspective in the NM-MRI template space and the SN, RN, and STN perspectives in the QSM template space are shown. [Diagram 5] 13 shows images illustrating the different steps for mapping NM boundaries from the template space to the original space. [Figure 6] 1 shows images illustrating the various steps for mapping the SN, RN and STN boundaries from the template space to the original space. [Figure 7A] Figure 1 shows the Dice similarity coefficient plotted against the volume fraction values ​​of neuromelanin (NM) based on various scenarios. [Figure 7B] Figure 1 shows the Dice similarity coefficient plotted against the volume fraction values ​​of neuromelanin (NM) based on various scenarios. [Figure 8A] Figure 1 shows Dice similarity coefficients plotted against the volume fraction values ​​of the substantia nigra (SN) based on different scenarios. [Figure 8B] Figure 1 shows Dice similarity coefficients plotted against the volume fraction values ​​of the substantia nigra (SN) based on different scenarios. [Figure 9] Figure 1 shows the Dice similarity coefficient plotted against the volume ratio values ​​of the red nucleus (RN) based on various scenarios. [Figure 10] Figure 1 shows the Dice similarity coefficient plotted against the volume ratio values ​​of the subthalamic nucleus (STN) based on various scenarios. [Figure 11] 1 shows plots illustrating the correlation between SN, RN, and STN iron content resulting from manual and template segmentation of 30 training cases (top plots) and 57 validation cases (bottom plots). [Figure 12] 1 shows a plot depicting the agreement between manual and template neuromelanin (NM) background measurements for 30 healthy controls. [Figure 13]1 shows images illustrating background boundaries versus NM boundaries identified using a dynamic programming algorithm (DPA). [Figure 14] 13 shows an image illustrating the background boundary and the NM DPA boundary after transformation from template space to original space. [Figure 15] 13 shows images illustrating the reconstructed boundary for different values ​​of the DPA parameter α. [Figure 16] We present simulation results for using DPA to identify the boundary of a rectangular region with two different intensities. [Figure 17] We present simulation results for using DPA to identify the boundary of a crescent-shaped region with two different intensities. [Figure 18] We present simulation results for using DPA to identify the boundary of a cashew-shaped region with two different intensities. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0009] For some diseases, such as Parkinson's disease (PD), information about the state of the disease and / or information useful in some treatments can be estimated from the volume or other geometric properties of some anatomical structures. In the case of PD, iron deposition and volume changes in the substantia nigra (SN), red nucleus (RN) and subthalamic nucleus (STN) are known to be associated with disease status and progression rate, and also serve as important targets for deep brain stimulation (DBS) treatment in PD patients. Therefore, accurate and comprehensive in vivo delineation of the SN and SN subregions, RN, and STN may be useful to fully investigate the changes in iron and neuromelanin (NM) composition in PD and other movement disorders affecting the midbrain. Therefore, accurate in vivo delineation of the SN, SN subregions, and other midbrain structures, such as the red nucleus (RN) and subthalamic nucleus (STN), may be useful to fully investigate the changes in iron and NM in PD.

[0010] To date, many studies still use manual or semi-automated approaches to demarcate deep brain gray matter. However, manual segmentation is time-consuming, especially when large amounts of data need to be evaluated. Also, unless the assessor is well trained, manual delineation has low reproducibility across individuals or across sites. Some approaches for in vivo delineation can include the use of templates to map iron and NM content. In such approaches, creating a standardized template can significantly impact (a) recognizing changes in iron and NM distribution, (b) automatically calculating the volumes associated with such distribution, (c) quantifying changes in iron content and NM signal, and (d) the reliability of these measurements. Anatomical templates of the SN using conventional structural MR sequences have been used previously in several studies. T1-weighted (T1W) and T2-weighted (T2W) images can be employed to create a high-resolution probabilistic in vivo subcortical nuclei atlas, which can segment the SN into the substantia nigra pars compacta (SNpc) and the substantia nigra pars reticulata (SNpr). Such atlases, when based on images from a young adult population (e.g., mean ± SD age: 28.9 ± 3.6 years), may not be suitable for studies with elderly subjects. Moreover, both T1W and T2W contrast images cannot be easily used to delineate the highly meshed boundary between the SNpc and SNpr in elderly subjects. Based on the anatomical connections of SN subregions to different parts of the brain, diffusion-based tractography can be used to segment the SN into SNpc and SNpr, or to subdivide the SN / ventral tegmental area (VTA) into dorsomedial and ventrolateral subregions. However, the estimated structural connectivity has well-known biases and is highly dependent on data acquisition and fiber tracking algorithms, and diffusion imaging suffers from low-resolution data acquisition. Hence, direct visualization and segmentation of the SN and SN subregions using high-resolution imaging would be a more accurate option, especially when it comes to detecting subtle pathological changes in the SN.

[0011] One approach to determine the boundaries of the midbrain nuclei and to examine pathological changes in PD patients is to use T2 * This can be based on the use of gradient-weighted echo (GRE) imaging and susceptibility-weighted imaging (SWI). Meanwhile, the development of quantitative susceptibility mapping (QSM) allows 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. In addition, for preoperative targeting of DBS, QSM can be used to obtain a more accurate estimate of the iron content of the brain than traditional T2 * It is superior to W imaging. However, QSM alone cannot separate the SNpc from the SNpr because both contain iron. This limitation can be overcome by using neuromelanin-sensitive MRI (NM-MRI), which has been developed in the past few years. The SNpc and ventral tegmental area (VTA) are composed of dopaminergic neurons that contain mainly NM, whereas the SNpr is not. Hence, the hyperintense signals seen in NM-MRI in the midbrain are spatially associated with the regions of the SNpc and VTA (which has been verified by postmortem histological studies). Hence, the overlap between the volume of the NM (SNpc plus VTA) and the volume of the iron-containing SN (SNpc plus SNpr) is considered to represent the SNpc.

[0012] So far, existing approaches do not accurately and reliably segment the SNpc using QSM and / or NM-MRI imaging. SN atlases may be obtained using QSM only. Alternatively, NM templates based on NM-MRI imaging, created for example using manual drawing, automatic segmentation, or artificial intelligence, may be employed. The purpose of using templates would be to facilitate the identification of boundaries. However, template mapping is not perfect, and the ideal space for marking boundaries is in the original pristine data space. Whether in template space (standard fixed brains or volunteers to which all other brains are mapped) or original space (original brain imaging data for a given individual), simple thresholding methods have drawbacks. Setting the threshold to relatively high or relatively low values ​​may lead to dramatic changes in the estimates of 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 may make it difficult to use some algorithms such as region growing. Although SNs have a high iron content compared to the surrounding regions, iron is not distributed uniformly and there is iron reduction in nigrosome 1 (N1) territories. Gaps in structures such as N1 territories can also increase the difficulty of automatically segmenting regions of interest (ROIs) in SNs and NMs.

[0013] In the current 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 calculate 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) the boundaries of each structure in the template space can be calculated, (ii) these boundaries can be mapped back to the original space, and then (iii) the boundaries in the original space can be fine-tuned using a dynamic programming algorithm (DPA) to match the details of each individual's NM and SN characteristics. Experimental results for SN NM and iron content of SN, STN, and RN all 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, and 1.05, respectively. These high-quality results indicate that tissue attributes such as volume, iron content, and neuromelanin content can be measured with accuracy. NM and iron sensitive sequences may also be used, in which case co-alignment between the two may be required.

[0014] Referring to FIG. 1, a flow chart illustrating a method 100 for generating an anatomical template according to the inventive concept of the present disclosure is shown. The method 100 may include a computing device or system acquiring corresponding MR images for each training subject 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 of the plurality of training subjects (step 104). The method 100 may include a computing device or system mapping 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 may include averaging the mapped MR images and the initial anatomical template to generate a final anatomical template (step 108). The method 100 may include drawing a boundary of an anatomical structure of interest on the final anatomical template (step 110). The method 100 may include identifying a boundary of an anatomical structure of interest in a non-training subject using the final anatomical template (step 112).

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

[0016] MR images of the training subjects may be acquired via an MRI scanner. For example, in the conducted study, MR imaging was performed on a 3T Ingenia scanner (Philips Healthcare, Netherlands) using a 15-channel head array coil. The imaging parameters for a 3D gradient echo SWI sequence with a differential 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, 0.67 × 1.34 mm. 2 Interpolated spatial in-plane resolution = 0.67 × 0.67 mm 2 , a sensing factor of 2, elliptical sampling of k-space, and a total scan time equal to 4 minutes and 47 seconds. Magnetization transfer (MT) resonance state radio frequency pulses were used with a nominal flip angle of 90°, zero frequency offset, and a set of three block pulses, each of 1.914 milliseconds (ms) duration. The minimum allowed TR was used based on safety considerations for the particular absorption rate. Due to this long repetition time of 62 ms, seven echoes were collected.

[0017] In some implementations, the computing device or system can use the first echo of the MTC-SWI magnitude image (TE=7.5 ms) to delineate NM content since that echo provides the primary MT contrast. The computing device or system can use the second echo (TE=15 ms) or a combination of QSM images from two or more echoes for QSM reconstruction to assess iron deposition in the SN. The computing device or system can create a susceptibility map by (i) segmenting the brain using a brain extraction tool, BET, (ii) unfolding the original phase data using a 3D phase unfolding algorithm (3DSRNCP), (iii) removing unwanted background fields using Advanced Harmonic Artifact Reduction (SHARP), and (iv) using a truncated k-space partitioning (TKD)-based inverse filtering technique with an iterative approach to reconstruct the final QSM map.

[0018] Manual region of interest (ROI) segmentation can be employed in tracing the boundaries of anatomical structures of interest. To measure NM and iron content, NM, SN, RN, and STN regions of interest (ROIs) can be manually traced by a single rater on 4x magnified MTC magnitude and QSM maps using SPIN software (SpinTech, Inc., Bingham Farms, MI, USA). NM-based SN boundaries can be traced starting from the last caudal slice of 4-5 slices until the NM is no longer visible when 2 mm thick slices are used (or across a total thickness of 8-10 mm). Iron-based SN boundaries can be traced starting from one slice below the most cranial slice where the subthalamic nucleus is visible and continued for 4-6 consecutive slices up to the most caudal slice or across a total thickness of 8-12 mm. RN ROI can be outlined starting from the last caudal slice and continued for 3-4 slices as cranial or across a total thickness of 6-8 mm. STN ROIs can be traced skull-like for two slices or over a total thickness of about 4 mm. For all ROIs, a dynamic programming algorithm (DPA) can be used to determine the final boundaries to reduce subjective bias. All of these boundaries can then be reviewed one by one by a second reviewer and modified accordingly in agreement with the first reviewer.

[0019] It should be understood that the above-mentioned 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 a slice thickness of 1 mm is used, these slice numbers will all be doubled. Furthermore, although the description focuses mainly on the application to PD and anatomical structures of the midbrain, the methods and techniques described herein may also be applied to other applications or other anatomical regions, such as the deep cerebellar nuclei or the striate and caudate nuclei, or any other region of interest.

[0020] In some implementations, the computing device or system can cause the MR scanner to perform MR data acquisition, for example, 7.5 ms echo time MTC data acquisition, for each training subject of the multiple training subjects, starting from the original whole brain 64 slice data. The MR scanner can construct a corresponding three-dimensional (3D) image (or multiple corresponding two-dimensional (2D) images) for each training subject based on the corresponding acquired MR data. In acquiring the MR data of the multiple training subjects, the MR scanner can be configured with parameters selected to generate an MR image with / having contrast illustrating 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 image) and / or QSM data (QSM image) for each training subject of the multiple training subjects.

[0021] The method 100 may include a computing device or system generating an initial anatomical template based on a selected MR image of a first training subject of the multiple training subjects (step 104). The computing device or system (or its user) may select a first 3D MR image (or a corresponding 2D image) of the first training subject from among the MR images of the multiple training subjects. The selection may be random or may follow predefined preferences or criteria. In generating the initial anatomical template, the computing device or system may expand the first 3D image (or the corresponding 2D image) in a plane for all acquired slices of the first training subject, for example, by a factor of 2x, 4x or more depending on the desired final resolution. This step may also include an interpolation of the slice selection direction to obtain an isotropic resolution in the template space. The template space may be defined at a desired resolution or at a higher resolution than the resolution of the acquired MR data. In some implementations, the computing device or system may 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 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] The method 100 may include a computing device or system mapping the MR images of the other training subjects onto the template space defined by the initial anatomical template (step 106) by applying a global transformation followed by a local transformation. The computing device or system may apply the global transformation and the local transformation to each of the images of the other training subjects (other than the first training subject). The computing device or system may apply the global transformation over a series of slices covering the region of interest in each MR image of the other training subjects other than the first training subject after expanding each MR image in the plane by a magnification factor (e.g., a factor of 2, 4 or more) for all slices. For example, for the midbrain, if the slice thickness is 2 mm, the series of slices covering the region of interest may be a set of 50 central slices. The global transformation may include a rigid transformation followed by an affine transformation applying b-spline interpolation, for example, using the Insight Segmentation and Registration Toolkit freeware (ITK). The computing device or system may apply the same global transformation to map the QSM data for each of the other training subjects onto the initial QSM template.

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

[0024] To improve the matching between the transformed MR images and the initial anatomical template in the template space, the computing device or system can apply a local transformation to a first morphed version of the MR images of the other training subjects (other than the first training subject). The computing device or system can crop the first morphed version of the MR images in a plane covering the midbrain and to about 16 slices (e.g., if 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 the other training subjects. For each MR image of the training subjects (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 example implementation, a total of 26 training subjects can be used. The computing device or system can select MR data of one of these training subjects (referred to herein as the first training subject) to generate an initial anatomical template. The computing device or system can then apply global and local transformations to the MR data of each of the other 25 training subjects to map the MR data of each of the other 25 training subjects to the initial anatomical template in the template space after applying the global and local transformations.

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

[0027] Method 100 may include drawing (step 110) the boundaries of the anatomical structures of interest on the final anatomical template. Tracing of the boundaries of the anatomical structures (or regions) of interest may be performed manually, for example by a physician or radiologist. The final anatomical template generated according to steps 102-108 of method 100 exhibits enhanced contrast of the anatomical structures of interest, thus enabling more accurate tracing of the boundaries of such structures.

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

[0029] The method 100 may further include using a boundary detection algorithm to refine the boundary of the anatomical structure of interest in the final anatomical template. The boundary detection algorithm may be implemented as a dynamic programming algorithm (DPA). A computing device or system may run the DPA for boundary detection to determine the template boundary. The DPA may 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] The method 100 may include a computing device or system using the final anatomical template to identify boundaries of an anatomical structure of interest in the untrained subject (step 112). The use of the final anatomical template to identify boundaries of an anatomical structure of interest in the untrained subject is discussed in detail below with respect to FIG.

[0031] Referring now to FIG. 2, a flow chart illustrating a method 200 of automatic template-based identification of characteristics of an anatomical structure in accordance with the inventive concepts of the present disclosure is shown. The method 200 may include acquiring an MR image (or MR data) of an untrained subject (step 202). The method 200 may include mapping the MR image of the untrained subject to a final anatomical template in a template space by applying a global transformation followed by a local transformation (step 204). The method 200 may 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). The method 200 may also include applying an inverse transform to the projected boundaries of the anatomical structures of interest to generate estimates of the boundaries in the MR image of the untrained subject (step 208).

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

[0033] The method 200 may include a computing device acquiring (step 202) an MR image (or MR data) of an untrained subject. The MR scanner may acquire the MR data (e.g., a 3D MR image or multiple 2D MR images) of the untrained subject in a manner similar to that discussed above with respect to step 102 of the method 100 of FIG. 1. The untrained subject may be one (e.g., a patient or a research subject) that is not one of the multiple training subjects used in the method 100 to generate the final anatomical template. The MR scanner may acquire the NM data and / or the QSM data of the untrained subject.

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

[0035] After applying the global transformation, the computing device or system can crop the first morphed version of the MR image of the untrained subject. For example, the computing device or system can label the RNs in the two central slices of the template data and return them to the original whole brain. If the RNs are still evident in only 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 to best center the data before the final transformation back to the template space. The computing device or system can expand the cropped image volume by a factor of four (or other expansion factor) in the plane and then apply a local transformation to the expanded cropped image volume to obtain a second morphed version of the MR image of the untrained subject with an isotropic image (the transformation, together with the interpolation of the slice selection direction, allows the image to be isotropic). The second morphed version of the MR image of the untrained subject represents a relatively good match of 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 may include projecting boundaries of one or more anatomical structures of interest from the final anatomical template onto the mapped MR image of the non-trained subject (step 206). Once the MR image of the non-trained subject has been matched or mapped to the final anatomical template (by applying global and local transformations), a computing device or system may project the boundaries of the anatomical structures of interest depicted in the final anatomical template onto a second morphed version of the MR image of the non-trained subject.

[0037] Method 200 may include applying an inverse transform to the projected boundary of the anatomical structure of interest to generate an estimate of the boundary in the MR image of the untrained subject (step 208). The inverse transform is performed to map the template boundary of the anatomical structure of interest back onto the midbrain (e.g., onto the original image space) in the MR image of the untrained subject. However, because everyone is different and template mapping is not perfect, there is no guarantee that the projected boundary will fit well.

[0038] The method 200 may further include a computing device or system using a boundary detection algorithm (or DPA) to refine the boundary of the anatomical structure of interest in the original acquired MR image of the untrained subject. The computing device or system may apply a threshold to the area inside the transformed boundary for the QSM data to remove negative values ​​and may use image thinning to determine the corresponding centerline. To best select an initial starting point for the DPA, both Otsu histogram analysis and a threshold-based approach may be used to determine whether the original boundary extends too far outside the structure of interest. For the NM data, the computing device or system may determine the background intensity and a constant equal to four times the background standard deviation averaged over all 25 training subjects (or training subjects other than the first training subject) added to create a threshold below which signals are set to zero. For the QSM data, the starting threshold may be set to zero. When the Otsu threshold provides a value that causes pixels inside the structure (or region) of interest to be removed, the template transformed boundary may be correspondingly reduced. 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 again refine the resulting boundaries using the same DPA used in the template space (as discussed above with respect to FIG. 1). This DPA, together with the RN boundaries, can prevent leakage of SN into RN, and with the help of the template boundaries, STN can be distinguished from SN.

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

[0040] In some implementations, the computing device 304 may 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. If the computing device 304 is communicatively coupled to the MR scanner 302, the MR data acquired by the MR scanner 302 may 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, the computing device 304 can be integrated within the MR scanner 302. For example, the processor 306, the memory 308, and / or the display device 310 can be integrated within the MR scanner 302. In such implementations, the MR scanner 302 can acquire MR data and perform the method 100 and / or the method 200.

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

[0043] The system 300 or computing device 304 can calculate the overlap between the NM complex volume (SNpc plus VTA) and the iron containing SN volume (SNpc and SNpr) by overlaying the two ROIs from the MTC data and the QSM data, respectively. The system 300 or computing device 304 can normalize the overlap by the iron containing SN volume creating an overall overlap fraction (OOF) measure, which is essentially a measure of the SNpc fraction relative to the total SN volume, as shown below. JPEG0007663710000001.jpg18153

[0044] The process of generating the final anatomical template described above with respect to Figure 1 was evaluated using MR data from 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 (HCs) were scanned. The SN, STN, and RN were manually traced for all 87 cases. Of these, 30 (test dataset: 17 males, 13 females, age range 66±7.2 years) were used for initial training of the template approach described above. Once all aspects of the algorithm were in place, the method was then validated on the next 57 cases (validation dataset: 17 males, 40 females, age range 61.9±5.0 years). Two measures were used to evaluate the performance of the anatomical template generation process. These are the Dice similarity coefficient, which indicates 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 the template volume divided by the volume of the manual segmentation.

[0046] Finally, all data were combined to generate quantitative information regarding structure volume, NM and iron content, and overlap between SN and NM regions. Total iron content was calculated by summing the product of the volume and the average sensitivity of the structure across all slices in which it was imaged. Similarly, total NM content resulted from the sum of the products of NM volume and NM contrast across the corresponding slices.

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

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

[0049] Referring to FIG. 6, images illustrating various steps for mapping the boundaries of the SN, RN, and STN from template space to the original image space for QSM data are shown. Each column represents a different slice. The first (top) row, designated as row A, shows images showing various slices of the constructed QSM template. The second row, designated as row B, shows the same image with QSM DPA boundaries for SN 602, RN 604, and STN 606. The third row, designated as row C, shows the boundaries transformed to the original space for SN 602, RN 604, and STN 606. The final row, designated as row D, shows images showing various midbrain slices of a subject with the final boundaries of SN 602, RN 604, and STN 606 obtained after applying DPA to the boundaries in row C.

[0050] For the initial training of the template, 30 cases were processed for MTC and QSM data. The completeness of the template automatic MTC background intensity was confirmed by a slope relating one to the other of 0.99 and R 2 This is demonstrated by the fact that the Δp value is 0.53 and p-value is less than 0.001. The background value is important for proper thresholding of the NM and iron content signals.

[0051] Measurements of Dice similarity coefficient and volume ratio (VR) were performed for different thresholds for both MTC and QSM images, and data associated with NM and SN Dice coefficients plotted against VR are shown in Figures 7A-7B and 8A-8B, respectively. Referring to Figures 7A and 7B, plots of Dice similarity coefficient plotted against VR values ​​for neuromelanin (NM) based on various scenarios are shown. In both Figures 7A and 7B, plot A corresponds to a scenario where no threshold is applied, plot B corresponds to a threshold of NM contrast above 1000, plot C corresponds to a threshold of NM contrast above 1250, and plot D corresponds to a threshold of NM contrast above 1500. The plots in Figure 7A are generated using data for 30 training cases, while the plots in Figure 7B are generated using data for 57 validation cases. The average Dice, volume ratio, and volume loss associated with the template data and the manually drawn data are quoted for each threshold. A threshold of 1000 units keeps volume loss to a minimum and produces excellent results.

[0052] Referring to Figures 8A and 8B, plots of Dice similarity coefficient plotted against SN VR values ​​based on various scenarios are shown. In both Figures 8A and 8B, plot A corresponds to the scenario of no threshold applied, plot B corresponds to a susceptibility value threshold of more than 50 ppb, plot C corresponds to a susceptibility value threshold of more than 75 ppb, and plot D corresponds to a susceptibility value threshold of more than 100 ppb. The plot in Figure 8A was generated using data for 30 training cases, while the plot in Figure 8B was generated using data for 57 validation cases. The average Dice, volume ratio, and volume loss associated with the template data and the manually drawn data are quoted for each threshold. A threshold of 50 ppb produces the least loss in SN volume.

[0053] The higher the threshold, the tighter the distribution will be, theoretically eventually approaching unity on both axes. However, a higher threshold also causes a higher loss of volume. Thus, there must be a trade-off between sufficient Dice and the volume ratio with volume loss. For 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. A higher threshold such as 1250 caused an average volume loss of just over 10% (average Dice=0.92, average VR=0.96), while data showing NM contrast above 1500 caused an average loss of 20%, with relatively higher Dice and VR (0.94 and 0.98, respectively) than the lower threshold. A threshold of 1000 appears to result in acceptable Dice and VR values ​​while keeping the volume loss below 10%. Similarly, as shown in Figures 8A-8B for SN, an average susceptibility threshold of 50 ppb resulted in an average volume loss of just under 10%, while the average Dice and VR values ​​were 0.88 and 0.94, respectively.

[0054] FIG. 9 shows plots of Dice similarity coefficients 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 top row of plots is generated using a dataset of 30 training cases, while the bottom row of plots is 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 quoted in each plot. The 50 ppb threshold limits the volume loss to roughly 10%.

[0055] FIG. 10 shows plots of Dice similarity coefficients plotted against VR values ​​for the subthalamic nucleus (STN) under various scenarios. Plot A (in the left column) represents the scenario where no threshold is applied, while plot B (right column) represents the 50 ppb threshold scenario. The top row of plots is generated using a dataset of 30 training cases, while the bottom row of plots is 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 thresholds are quoted within each plot. The 50 ppb threshold keeps the volume loss at approximately 10%.

[0056] 9 and 10, before applying any threshold to the QSM data, the average Dice and Volume Ratios were 0.93 and 1.06 for RN and 0.76 and 0.95 for STN. However, applying a 50 ppb threshold to the QSM data results in average Dice and VRs of 0.95 and 1.04 for RN and 0.83 and 0.98 for STN. The average VRs greater than 1 indicate that most of the structures found by the fully automated template / DPA approach tend to be larger than those from the manual / DPA approach.

[0057] Referring to FIG. 11, plots illustrating the correlation between the iron content of the SN, RN, and STN resulting from manual segmentation and template segmentation are shown. The correlation between the iron content of the SN, RN, and STN resulting from manual segmentation and template segmentation for the 30 training cases is shown in the upper plot, while the correlation between the iron content of the SN, RN, and STN resulting from manual segmentation and template segmentation for the 57 validation cases is shown in the lower plot. The ROIs associated with the SN, RN, and STN, as well as the final results of the template versus manual drawing for iron content using a threshold of 50 ppb, are shown in various plots. The slope and R2 is shown in each of the plots in Figure 11. The p-values ​​are less than 0.001 for SN and RN, and equal to 0.006 for STN.

[0058] For template validation, the validation data set included 57 cases that were used to process the QSM and MTC data. The Dice similarity coefficient was plotted against the VR values ​​for both NM and SN, and these results are shown in Figures 7B and 8B, respectively. For 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 nearly 5%, while the average Dice and VR values ​​were 0.89 and 0.97, respectively. As with the previous data set, a threshold of 1000 for NM contrast and a threshold of 50 ppb for the average sensitivity of SN yielded satisfactory results in terms of average Dice, VR, and volume loss across the 57 cases.

[0059] The Dice similarity coefficients versus VR for RN and STN are shown in Figures 9 and 10, respectively. Before applying any threshold to the QSM data, the average Dice 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 50 ppb threshold on the QSM data resulted in average Dice and VR values ​​of 0.95 and 1.03, and 0.81 and 0.90 for RN and STN, respectively, and a template average 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 Figure 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. The mean and standard deviation values ​​for the first and second data sets and the merged data are shown. [Table 1] Table 1. Mean and standard deviation (SD) values ​​of volume estimates, Dice similarity coefficients, and volume ratios (VRs) for each structure.

[0061] In the study described for Figures 4-12, NM and QSM images derived from a single sequence (less than 5 min) were used for automatic segmentation of the SN, STN, and RN in original space without the need to manually draw ROIs for non-template-based subjects. A multi-contrast atlas combined with DPA for border detection in both template space and in the original image after back transformation from 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 of the existing templates do not include NM, SN, STN, and RN and do not contain the resolution presented in this work. For example, the MNI template is 1 mm isotropic, while the approach described herein used 0.67 mm in-plane isotropic data that was interpolated and further interpolated to 0.167 mm isotropic resolution in 3D. In practice, we found that using whole-brain global deformable registration with 0.67 mm isotropic resolution using either ANT or SpinITK, we were unable to obtain a morphological transformation that could consistently and well match the shape of midbrain structures across subjects. Therefore, to solve this problem and lead to significantly improved results, we added a local registration (also referred to as local transformation here) step. From this local transformation approach, both iron and neuromelanin templates 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 precisely define the SN boundaries. To overcome this limitation, several other studies have attempted to obtain SN atlases using either NM-MRI or QSM. By using a dynamic atlas constructed from NM-enhanced brain images for automatic segmentation of the SN, one group showed a Dice value below 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 the NM can be significantly improved by using MT-MRI acquisition sequences. To the inventors' knowledge, no study has 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 (>5 and <18 years) or older (>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 midbrain structures. These inter- and intra-subject age-related and disease-related morphometric changes may be significant and affect the success of using template approaches when localizing the SN. The most appropriate atlas for a given study is one that requires the least amount of global or local distortion from native subject space to atlas space, therefore, we used older HC subjects (mean age: 63.4 years) to create an atlas for studying patients with neurodegenerative diseases. To mitigate the bias caused by inter-subject variability, several studies have used probabilistic atlases and constructed a probabilistic atlas of the SNpc based on NMS-MRI sequences to investigate microstructural abnormalities of the SNpc using diffusion MRI in PD patients. They applied a symmetric diffeomorphic registration to align the T1 images and the 27HC SNpc mask onto the MNI space and thresholded the atlas with a probability of 50%. However, the average Dice coefficient for the atlas versus human rater was less than 0.61. Thresholding can lead to dramatic changes in the depiction of NMs, thus reducing the reproducibility of the atlas.

[0065] Regarding the validation phase, the second dataset shows that the Dice similarity coefficient and VR values ​​for structures are very close to those from the first dataset, which confirms the consistency of the results produced by the automated template processing approach proposed herein, and the fact that these values ​​are very close to unity indicates the good performance of the approach.

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

[0067] In conclusion, the new approach described here for fully automated segmentation of deep gray matter in the midbrain (or anatomical structures in general) uses a template approach as a guide, but uses the original data together with a dynamic programming algorithm to fine-tune the boundaries. With 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 make it possible to study changes in these imaging biomarkers for a variety of neurodegenerative diseases without the need for manual tracing of these structures.

[0068] Regarding neuromelanin (NM) background measurements, because background values ​​are important for proper thresholding of NM signals, 30 MTC cases were processed to obtain neuromelanin (NM) background measurements. Referring to FIG. 12, a plot showing agreement between manual and template neuromelanin (NM) background measurements for 30 healthy controls (or training subjects) is shown. Specifically, the plot illustrates the integrity of the template automated MTC background intensity measure versus the manual drawing. Agreement between the manual and template MTC background measures was measured with a slope of 0.99 and an R of 0.53. 2 , and p-values ​​less than 0.001 are indicated.

[0069] Regarding 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. The initial boundary can be drawn on the template data for both NM and QSM. 2. The background can then be drawn into the template image as shown in Figure 13. In Figure 13, both the NM 702 and background regions 704 are shown. 3. The boundary and background can then be transformed back to the original space as shown in Figure 14, where NM region 802 and background region 804 are shown after transformation from template space to the original space. Using a background value plus 1000 makes it easier to refine the DPA boundary to provide faster convergence and to use a smaller and safer search radius to prevent the algorithm from leaking outside the original boundary. 4. The background mean value + 4σ (where σ is the standard deviation of the background region of interest, which was found to be roughly 250 units in this study) can be used as a threshold for the MTC data to remove all points below this threshold to determine the NM refined boundary. For the QSM data, pixels with sensitivity below zero ppb can be removed before DPA is run for all SN, STN, and RN. 5. A 3x3 Gaussian filter can be used to smooth the initial boundary. 6. The centerline can then be determined using thinning techniques. 7. Boundary refinement before DPA: A global threshold may not be valid for the entire structure in the MTC data, since there may be some variation in intensity around the structure. 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 40x40 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 modified to the nearest remaining point (called the Otsu boundary). This step helps to speed up the convergence of DPA. 8. DPA Implementation: For each point along the centerline, DPA can be performed in the corresponding search box. To allow for curved shapes, the center of the centerline can be used to determine the initial ray. DPA can then 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. Once an end point is reached, these rays can be swept 180°. The algorithm can then move back along the centerline until the centerline has swept another 180° to reach the opposite end point. Finally, the centerpoint can move back along the centerline back to the starting point to close the boundary. This process can be repeated five times to search both inside and outside the boundary to obtain the best result. For NM, STN, and SN, this step can be followed by five more iterations using only the inside search.

[0070] The search radius was limited to 4 pixels both inside and outside the boundary of the expanded space, but the cost function used in this DPA includes derivative and radius of curvature terms that avoid leakage of the structure of interest into adjacent objects. For points outside the boundary, the negative derivative was set to zero when searching outward from the centerline. This restriction prevented the boundary from leaking into nearby bright objects. The cost function is given below: JPEG0007663710000003.jpg1389In formula, 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. max represents the maximum derivative of the interior of the image, and R avgis the average radius over the previous three radii. The constant α represents the relative weighting of the derivative and radius terms, 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 α can be set. If there are sharp edges, these can be smoothed out by a large choice of α. Values ​​of α between 0.05 and 0.15 all played a similar role in faithfully finding edges, so a value of α=0.1 was chosen to be conservative.

[0071] Referring to FIG. 15, images illustrating reconstructed boundaries with different values ​​of the DPA parameter α are shown. The choice of α can dramatically affect the final DPA boundary. In the top left image, the NM boundary is obtained using α=0.05. In the top right image, the NM boundary is determined using α=0.10. In the bottom left image, the NM boundary is obtained using α=0.15. In the bottom right image, the NM boundary is obtained using α=0.20. Note that smaller values ​​of α lead to less smoothing from the radius constraint, while larger values ​​of α lead to more smoothing. Higher values ​​of α tend to make the structure rounder and shrink more towards a circle. Low values ​​lead to boundaries with very zigzag edges. Finally, candidate points for the new boundary can be selected by maximizing the above mentioned cost function for each ray.

[0072] Some example simulations for various shapes are shown in Figures 16-18. Referring to Figure 16, simulation results for identifying the boundary of a rectangular region with two different intensities using DPA are shown. Image (a) shows a rectangular region with two different intensities. Image (b) shows a rectangular region with a boundary drawn around the central area, and image (c) further shows the found centerline and updated boundary after applying DPA for five iterations. Image (d) is similar to image (b) except that noise with a contrast-to-noise ratio (CNR) of 7:1 has been 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 five iterations of DPA that matches the correct area of ​​1400 pixels.

[0073] For the simulation corresponding to Figure 16, the central region was set to 100 units, the outer frame region within 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 with a mean value of zero and a standard deviation of 10 units was added to the image. A first example of a rectangle with sharply defined corners and 1400 pixels is shown in Figure 16. Despite the presence of noise, the boundary was still found perfectly.

[0074] FIG. 17 shows simulation results for identifying the boundary of a crescent-shaped region with two different intensities using DPA. Image (a) shows a crescent-shaped region with two different intensities, and image (b) shows the same region was boundary drawn around the central area of ​​the crescent-shaped region. Image (c) shows the crescent-shaped region with the corresponding centerline and updated boundary after five iterations of DPA. Images (d) and (e) show the centerline and boundary update after 10 and 15 iterations of DPA, respectively. Image (f) shows the centerline and boundary after only five iterations when using the adaptive Otsu threshold approach. Image (g) is similar to image (c) with noise added, with a CNR of 7:1, and image (h) shows the centerline and boundary after only five iterations of DPA when using the adaptive Otsu threshold approach.

[0075] The crescent shaped region in Figure 17 was chosen to mimic the more challenging case of detecting the boundary of a curved object. Figure 17 shows the effect of running different numbers of iterations. After running a total of 35 iterations, the average value (image (f)) was found to be 990.5 with a standard deviation of 1.5, whereas using the Otsu approach in the presence of an SNR of 10:1, the average value was 984 with a standard deviation of 2.0.

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

[0077] Finally, the cashew-shaped region in Figure 18 was chosen to mimic SN and evaluate the method in the absence of sharp edges. After running the first 35 iterations, the mean value was found to be 1086.5 with a standard deviation of 2.5. Using the Otsu approach in the presence of a 10:1 SNR and running 30 more iterations, the mean value was found to be 1097 with a standard deviation of 0.0.

[0078] Those skilled in the art should understand that the process described in the present 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 a random access memory (RAM), a read only memory (ROM), a cache memory, a disk memory, any other memory, or any other computer readable medium. The process described in the present disclosure can be implemented by an apparatus including at least one processor and / or a memory having executable code instructions stored therein. The code instructions, when executed by the at least one processor, can cause the execution of any of the processes or operations described in the present disclosure. The apparatus can be, for example, an MRI scanner or a computing device.

Claims

1. 1. A magnetic resonance imaging (MRI) system comprising: a magnetic resonance (MRI) scanner configured to acquire MR data; At least one processor; a memory having computer code instructions stored thereon, the computer code instructions, when executed by the at least one processor, causing the at least one processor to: acquiring, via the MRI scanner, corresponding MR images with contrast illustrating one or more anatomical structures of interest for each of the plurality of training subjects; generating 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, applying a global transformation to the MR images of the training subject to generate a first morphed version of the MR images in the template space that represents a first estimate of the MR data of the training subject; and applying a local transformation to the first morphed version of the MR image of the training subject to generate, in the template space, a second morphed version of the MR image that represents a second estimate of the MR data of the training subject; 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; delineating boundaries of the one or more anatomical structures of interest in the final anatomical template; and determining the boundaries of the one or more anatomical structures of interest for another untrained subject using the final anatomical template.

2. In determining the boundaries of the one or more anatomical structures of interest using the final anatomical template, the at least one processor: The MRI system of claim 1 , configured to refine the boundaries of the one or more anatomical structures of interest in the final anatomical template using a boundary detection algorithm.

3. In determining the boundaries of the one or more anatomical structures of interest using the final anatomical template, the at least one processor: acquiring, via the MRI scanner, MR images of an untrained subject with contrast illustrating the one or more anatomical structures of interest; applying the global transformation to the MR image of the untrained subject to generate a first morphed version of the MR image of the untrained subject; applying the local transformation to the first morphed version of the MR image of the untrained subject to generate a second morphed version of the MR image of the untrained 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 untrained subject to determine an estimate of the boundaries of the one or more anatomical structures of interest in the template space for the untrained subject; and applying an inverse transformation to the estimates of the boundaries of the one or more anatomical structures of interest in the template space to determine a first estimate of the boundaries of the one or more anatomical structures of interest in the original MR image of the untrained subject.

4. The at least one processor:

4. The MRI system of claim 3, configured to use a boundary detection algorithm to fine-tune the first estimate of the boundary of the one or more anatomical structures of interest in the original MR image of the untrained subject to generate a second estimate of the boundary of the one or more anatomical structures of interest.

5. The MRI system of claim 3 , wherein the inverse transformation comprises a concatenation of the inverse of the local transformation followed by the inverse of the global transformation.

6. The MRI system of claim 1 , wherein applying the local transformation to the first morphed version of the MR image comprises applying the local transformation to a cropping region of the first morphed version of the MR image.

7. The MRI system of claim 1 , wherein the corresponding MR image comprises a quantitative susceptibility map (QSM) or a neuromelanin image.

8. 4. The MRI system of 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 an original MR image of the untrained subject.

9. 4. The MRI system of claim 3, wherein the at least one processor is configured to determine an average intensity 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 an original MR image of the untrained subject.

10. 4. The MRI system of 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 original MR image of the untrained subject.

11. 1. A method comprising: acquiring, via a magnetic resonance imaging (MRI) scanner, corresponding MR images with contrast illustrating one or more anatomical structures of interest for each of the plurality of training subjects; generating, by a computing device, 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, applying a global transformation to the MR images of the training subject to generate a first morphed version of the MR images in the template space that represents a first estimate of the MR data of the training subject; and applying a local transformation to the first morphed version of the MR image of the training subject to generate, in the template space, a second morphed version of the MR image that represents a second estimate of the MR data of the training subject; 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; delineating boundaries of the one or more anatomical structures of interest in the final anatomical template; and using the final anatomical template to determine the boundaries of the one or more anatomical structures of interest for other untrained subjects.

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

13. determining the boundaries of the one or more anatomical structures of interest in the MR image using a final anatomical template; acquiring, via the MRI scanner, MR images of an untrained subject with contrast illustrating the one or more anatomical structures of interest; applying the global transformation to the MR image of the untrained subject to generate a first morphed version of the MR image of the untrained subject; applying the local transformation to the first morphed version of the MR image of the untrained subject to generate a second morphed version of the MR image of the untrained 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 untrained subject to determine an estimate of the boundaries of the one or more anatomical structures of interest in the template space for the untrained subject; and applying an inverse transform to the estimate of the boundary of the one or more structures of interest in the template space to determine a first estimate of the boundary of the one or more anatomical structures of interest in the original MR image of the untrained subject.

14. 14. The method of claim 13, further comprising using a boundary detection algorithm to fine-tune the first estimate of the boundary of the one or more anatomical structures of interest in the original MR image of the untrained subject to generate a second estimate of the boundary of the one or more anatomical structures of interest in the second morphed version of the MR image of the untrained subject.

15. The method of 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. The method of claim 11 , wherein applying the local transformation to the first morphed version of the MR image comprises applying the local transformation to a cropping region of the first morphed version of the MR image.

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

18. The method of claim 13, further comprising determining 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 a mean intensity 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 14. The method of 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 morphed version of the MR image of the untrained subject.

20. A non-transitory computer readable medium comprising computer code instructions stored on the non-transitory computer readable medium, the computer code instructions, when executed by a processor, causing the processor to: acquiring, via a magnetic resonance imaging (MRI) scanner, corresponding MR images with contrast illustrating one or more anatomical structures of interest for each of the plurality of training subjects; generating 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, applying a global transformation to the MR images of the training subject to generate a first morphed version of the MR images in the template space that represents a first estimate of the MR data of the training subject; and applying a local transformation to the first morphed version of the MR image of the training subject to generate, in the template space, a second morphed version of the MR image that represents a second estimate of the MR data of the training subject; 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; delineating boundaries of the one or more anatomical structures of interest in the final anatomical template; and determining the boundaries of the one or more anatomical structures of interest for other, untrained subjects using the final anatomical template.

Images