System and method for template-based automated detection of anatomical structures

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 of related applications This application claims priority and interest to U.S. Provisional Patent Application No. 63 / 170,229, filed on 2 April 2021, entitled “SYSTEMS AND METHODS FOR AUTOMATIC TEMPLATE-BASED DETECTION OF ANATOMICAL STRUCTURES,” the contents of which are incorporated herein by reference in their entirety.

[0002] This disclosure generally relates to the field of detecting or identifying the boundaries of anatomical structures. Specifically, this disclosure relates to methods and systems for generating standardized templates and for automatically detecting or identifying the boundaries of anatomical structures using such templates. [Background technology]

[0003] One disease that requires the identification of anatomical boundaries is Parkinson's disease (PD), a chronic progressive neurodegenerative disorder affecting approximately 1% of individuals over 60 years of age. PD is pathologically characterized by early neurodegeneration of the substantia nigra (NM) in the substantia nigra compacta (SNpc) and increased iron deposition in the substantia nigra (SN). Degeneration of the SN is characteristic of the progression of numerous neurodegenerative diseases. While widespread neuronal loss of the SNpc also occurs in atypical Parkinsonian disorders, including progressive supranuclear palsy (PSP) and multiple system atrophy (MSA), various subregions of the SN are affected in these disorders.

[0004] The SN (syndrome) consists of two anatomically and functionally distinct regions: the SN reticular region (SNpr) and the SNpc (syndrome reticular region). The SNpc contains a dense distribution of NMs (non-mammary neurons) containing dopaminergic neurons, while iron content tends to be higher in the SNpr. However, clusters of dopaminergic neurons (known as neglosomes) in the SNpc are deeply embedded within the SNpr, and therefore, the boundary between the SNpr and SNpc is difficult to describe, particularly in the caudal region of the SN. The regional selectivity of PD is relatively specific to the loss of 50%–70% of pigment neurons in the ventrolateral layer of the SNpc (at the time of symptom onset). For 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, and also serve as important targets for deep brain stimulation (DBS) therapy in PD patients. [Overview of the project]

[0005] According to at least one embodiment, a magnetic resonance imaging (MRI) system comprises an MRI scanner configured to acquire magnetic resonance (MR) data, at least one processor, and a memory on which computer code instructions are stored. 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, a corresponding MR image with contrast illustrating one or more anatomical structures of interest via the MRI scanner. The at least one processor can generate an initial anatomical template image based on first MR data of a first of the plurality of training subjects. The initial anatomical template image can define a template space. For each of the multiple trained subjects other than the first trained subject, at least one processor may (i) apply a global transformation to the trained subject's MR image to generate a first morphed version of the MR image representing a first estimate of the trained subject's MR data in template space, and (ii) apply a local transformation to the first morphed version of the trained subject's MR image to generate a second morphed version of the MR image representing a second estimate of the trained subject's MR data in template space. At least one processor may average the initial anatomical template image and the second morphed versions of the multiple trained subjects' MR images to generate a final anatomical template image. At least one processor may draw the boundaries of one or more anatomical structures of interest on the final anatomical template and use the final anatomical template to identify the boundaries of one or more anatomical structures of interest on other untrained subjects.

[0006] In at least one embodiment, the method may include obtaining, via an MRI scanner, corresponding MR images with contrast illustrating one or more anatomical structures of interest for each of a plurality of training subjects. The method may 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 may define a template space. For each of the plurality of training subjects other than the first training subject, the computing device may (i) apply a global transformation to the training subject's MR image to generate a first morphed version of the MR image representing a first estimate of the training subject's MR data in the template space, and (ii) apply a local transformation to the first morphed version of the training subject's MR image to generate a second morphed version of the MR image representing a second estimate of the training subject's MR data in the template space. The method may 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 may include drawing the boundaries of one or more anatomical structures of interest on a final anatomical template, and using the final anatomical template to identify the boundaries of one or more anatomical structures of interest on other untrained subjects.

[0007] According to at least one embodiment, a non-temporary computer-readable medium may include computer code instructions stored on the non-temporary computer-readable medium. When executed by a processor, the computer code instructions can cause the processor to acquire corresponding MR images with contrast illustrating one or more anatomical structures of interest via an MRI scanner for each of a plurality of training subjects. The processor may generate an initial anatomical template image based on first MR data of a first of the plurality of training subjects. The initial anatomical template image may define a template space. For each of the plurality of training subjects other than the first training subject, the processor may (i) apply a global transformation to the training subject's MR image to generate a first morphed version of the MR image representing a first estimate of the training subject's MR data in the template space, and (ii) apply a local transformation to the first morphed version of the training subject's MR image to generate a second morphed version of the MR image representing a second estimate of the training subject's MR data in the template space. The processor can generate a final anatomical template image by averaging the initial anatomical template image with a second morphed version of MR images from multiple trained subjects other than the first trained subject. The processor can then draw the boundaries of one or more anatomical structures of interest on the final anatomical template and use the final anatomical template to identify the boundaries of one or more anatomical structures of interest on other untrained subjects. [Brief explanation of the drawing]

[0008] [Figure 1] This flowchart illustrates a method for generating an anatomical template according to the concept of the present invention as disclosed herein. [Figure 2] This flowchart illustrates a method for automated template-based identification of anatomical structure characteristics according to the concept of the present invention as disclosed herein. [Figure 3]The diagram illustrates a system for detecting or identifying the boundaries of anatomical structures according to the concept of the present invention as disclosed herein. [Figure 4] This image shows magnetic resonance (MR) images depicting the NM viewpoint in the NM-MRI template space, and the SN, RN, and STN viewpoints in the QSM template space. [Figure 5] The image illustrates various steps involved in mapping the NM boundary from the template space to the original space. [Figure 6] The image illustrates various steps for mapping SN, RN, and STN boundaries from the template space to the original space. [Figure 7A] The Dice similarity coefficients plotted against the volume ratio values ​​of neuromelanin (NM) based on various scenarios are shown. [Figure 7B] The Dice similarity coefficients plotted against the volume ratio values ​​of neuromelanin (NM) based on various scenarios are shown. [Figure 8A] The die similarity coefficients plotted against the volume ratio values ​​of black matter (SN) based on various scenarios are shown. [Figure 8B] The die similarity coefficients plotted against the volume ratio values ​​of black matter (SN) based on various scenarios are shown. [Figure 9] The dice similarity coefficients plotted against the volume ratio values ​​of the red nucleus (RN) are shown based on various scenarios. [Figure 10] The Dice similarity coefficients plotted against the volume ratio values ​​of the subthalamic nucleus (STN) based on various scenarios are shown. [Figure 11] This plot illustrates the correlations between iron content in SN, RN, and STN obtained from manual and template-based segmentation of 30 training cases (upper plot) and 57 validation cases (lower plot). [Figure 12] A plot illustrating the agreement between manual and template neuromelanin (NM) background measurements for 30 healthy controls is shown. [Figure 13]Figures show images illustrating background boundaries against NM boundaries identified using a dynamic programming algorithm (DPA). [Figure 14] Figures show images illustrating background boundaries after conversion from template space to original space, and NM DPA boundaries. [Figure 15] Figures show images illustrating boundaries reconstructed at various values of DPA parameter α. [Figure 16] Shows simulation results for identifying boundaries of rectangular regions having two different intensities using DPA. [Figure 17] Shows simulation results for identifying boundaries of crescent-shaped regions having two different intensities using DPA. [Figure 18] Shows simulation results for identifying boundaries of cashew-shaped regions having two different intensities using DPA. DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

[0009] For several diseases such as Parkinson's disease (PD), information relating to the disease state and / or information useful in some treatments may be estimated from the volume or other geometric characteristics of several anatomical structures. In the case of PD, for the substantia nigra (SN), iron deposition and volume changes in the 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) therapy in PD patients. Accordingly, accurate and comprehensive in vivo delineation of SN and sub-regions of SN, RN, and STN would be useful to fully investigate changes in the composition of iron and neuromelanin (NM) in PD and other movement disorders affecting the midbrain. Accordingly, accurate in vivo delineation of SN, sub-regions of SN, and other midbrain structures such as the red nucleus (RN) and subthalamic nucleus (STN) would be useful to fully investigate changes in iron and NM in PD.

[0010] To date, many studies have still used manual or semi-automatic approaches to segment the deep gray matter of the brain. However, manual segmentation is time-consuming, especially when large amounts of data need to be evaluated. Furthermore, manual drawing has low reliability of reproducibility between individuals or regions unless the evaluators are sufficiently trained. Several approaches for in vivo drawing may include the use of templates to map iron and NM content. In such approaches, creating a standardized template can significantly impact (a) the recognition of changes in iron and NM distribution, (b) the automatic calculation of volumes associated with such distributions, (c) the quantification of changes in iron content and NM signal, and (d) the reliability of these measurements. Anatomical templates of the substantia nigra using conventional structural MR sequences have been used in some studies previously. High-resolution stochastic in vivo subcortical nucleus atlases can be created by employing T1-weighted (T1W) and T2-weighted (T2W) images, and these images can now segment the substantia nigra pars compacta (SNpc) and substantia nigra reticularis (SNpr). Such atlases, when based on images from a young adult population (e.g., mean ± standard deviation age: 28.9 ± 3.6 years), may not be suitable for studies involving older adults. Furthermore, both T1W and T2W contrast images cannot be readily used to depict the highly interlocked boundaries between SNpc and SNpr in older subjects. Based on the anatomical connectivity 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 region (VTA) into dorsomedial and ventrolateral subregions. However, the estimated structural connectivity is known to have biases and is highly dependent on data acquisition and fiber tracing algorithms, and diffusion imaging suffers from low-resolution data acquisition. Therefore, direct visualization and segmentation of the SN and its subregions using high-resolution imaging would be a more accurate option, particularly with regard to the detection of 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 look at T2, where iron-containing regions appear at low intensity. * This can be achieved using weighted gradient echo (GRE) imaging and sensitivity-weighted imaging (SWI). On the other hand, the development of quantitative sensitivity mapping (QSM) allows for the quantification of iron stored in ferritin and hemosiderin. QSM-based techniques have shown that tissue magnetosensibility correlates well with brain iron in PD patients. In addition, in preoperative target guidance for DBS, QSM is used in traditional T2 * QSM is superior to W imaging. However, QSM alone cannot separate SNpc from SNpr because both SNpc and SNpr contain iron. This limitation can be overcome by using neuromelanin-sensitive MRI (NM-MRI), which has been developed in recent years. SNpc and the ventral tegmental region (VTA) are composed primarily of dopaminergic neurons containing NM, while SNpr is not. For this reason, high-intensity signals seen on NM-MRI in the midbrain are spatially associated with the regions of SNpc and VTA (this has been verified by postmortem histological studies). Therefore, the overlap between the volume of NM (SNpc and VTA combined) and the volume of iron-containing SN (SNpc and SNpr combined) is thought to represent SNpc.

[0012] To date, existing approaches do not accurately and reliably segment the SNpc using QSM and / or NM-MRI imaging. SN atlases may be obtained using QSM alone. Alternatively, NM templates based on NM-MRI imaging, created for example by manual drawing, automatic segmentation, or using artificial intelligence, may be employed. The purpose of using a template would be to facilitate boundary identification. However, template mapping is not perfect, and the ideal space for marking boundaries lies in the original raw data space. Whether in template space (a standard fixed brain or volunteer to which all other brains are mapped) or 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 estimates of NM or iron content, especially in PD patients with severe NM degeneration and iron deposition in the SN. Also, variable image contrast can make it difficult to use some algorithms such as region growing. Although the SN has a higher iron content compared to surrounding areas, iron is not uniformly distributed, and there is a reduction of iron in the nigrosome 1 (N1) territory. Gaps in structures such as the N1 territory can also increase the difficulty of automatically segmenting the regions of interest (ROI) of the SN and NM.

[0013] In the present disclosure, a system and method for detecting anatomical structures may include creating both NM and iron templates using high-resolution imaging from a single multi-echo NM MRI sequence and calculating the boundaries of each structure in template space. Specifically, using both NM and iron templates from a single high-resolution NM MRI sequence, (i) the boundaries of each structure in template space can be calculated, (ii) these boundaries can be mapped back into 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 feature details of each individual's NM and SN. Experimental results for NM in SN and iron content in SN, STN, and RN all showed strong agreement with respect to 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 demonstrate 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 may also be used, in which case joint repositioning between the two may be required.

[0014] Referring to Figure 1, a flowchart illustrating a method 100 for generating an anatomical template in accordance with the concept of the present invention as disclosed herein is shown. Method 100 may include a computing device or system acquiring a corresponding MR image for each of a group of training subjects (step 102), and generating an initial anatomical template based on a selected MR image of a first of the group of training subjects (step 104). Method 100 may also 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). Method 100 may also include averaging the mapped MR images and the initial anatomical template to generate a final anatomical template (step 108). Method 100 may also include drawing the boundaries of anatomical structures of interest on the final anatomical template (step 110). Method 100 may also include using the final anatomical template to identify the boundaries of anatomical structures of interest in untrained subjects (step 112).

[0015] Method 100 may include a computing device or system acquiring a corresponding MR image for each of multiple training subjects (step 102). Training subjects may be selected based on several predefined criteria, such as age, health status, or the condition and / or medical history of each subject. The study described herein was approved by the local ethics committee, and all subjects signed informed consent. Eighty-seven healthy subjects (mean age: 63.4 ± 6.2 years, range: 45–81 years, 53 women) were recruited locally by advertisement. Exclusion criteria for potential training subjects included (a) structural abnormalities such as tumors, subdural hematomas, or contusions from 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 training subjects can be acquired via an MRI scanner. For example, in the study conducted, MR imaging was performed using a 3T Ingenia scanner (Philips Healthcare, Netherlands) with a 15-channel head array coil. The imaging parameters for the 3D gradient echo SWI sequence using operating magnetization transfer contrast (MTC) pulses 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 resolution = 0.67 × 0.67 mm 2 This included a sensing coefficient of 2, elliptic sampling in k-space, and a total scan time equal to 4 minutes and 47 seconds. The magnetization transfer (MT) resonance state high-frequency pulses used were a nominal flip angle of 90°, a zero-frequency offset, and a set of three block pulses, each with a duration of 1.914 milliseconds (ms). The minimum acceptable TR was used based on considerations of safety for a specific absorption rate. Due to this long repetition time of 62 ms, seven echoes were collected.

[0017] In some implementations, a computing device or system can depict NM content using a first echo of an MTC-SWI size image (TE=7.5ms) since the echo provides primary MT contrast. For QSM reconstruction, the computing device or system can assess iron deposition in the SN using a second echo (TE=15ms) or a combination of QSM images from two or more echoes. A sensitivity map can be created 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 high-harmonic artifact reduction (SHARP); and (iv) using a sectioned k-space segmentation (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 when tracing the boundaries of anatomical structures of interest. To measure NM and iron content, ROIs of NM, SN, RN, and STN can be manually traced by a single evaluator on a 4x magnified MTC size and QSM map 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 over a total thickness of 8-10 mm). Iron-based SN boundaries can be traced starting from the first slice below the most cranial slice where the subthalamic nucleus is visible, and continued with 4-6 consecutive slices to the most caudal slice, or over a total thickness of 8-12 mm. RN ROIs can be outlined starting from the last caudal slice and continued with 3-4 slices, like the skull, or over a total thickness of 6-8 mm. STN ROIs can be traced across two slices, like a skull, or over a total thickness of approximately 4 mm. For all ROIs, a dynamic programming algorithm (DPA) can be used to determine the final boundaries to mitigate subjective bias. These boundaries can then all be reviewed one by one by a second evaluator and modified accordingly in agreement with the first reviewer.

[0019] It should be understood that the data acquisition parameters and / or data processing scenarios described above represent exemplary implementations provided for illustrative purposes only, not for limitation. For example, the slice numbers provided above are specific to a given resolution with a slice thickness equal to 2 mm. If slices with a thickness of 1 mm were used, these slice numbers would all be doubled. Furthermore, although the description focuses primarily on applications to the PD and anatomical structures of the midbrain, the methods and techniques described herein may also be applicable to other anatomical regions, such as the deep cerebellar nuclei or the striatal and caudate nuclei, 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 for each of multiple training subjects, for example, starting with the original whole-brain 64-slice data and performing MTC data acquisition with an echo time of 7.5 ms. For each training subject, the MR scanner can construct a corresponding three-dimensional (3D) image (or multiple corresponding two-dimensional (2D) images) based on the corresponding acquired MR data. When acquiring MR data for multiple training subjects, the MR scanner can be configured with parameters selected to generate MR images with / having contrast illustrating one or more anatomical structures of interest (e.g., SN, RN, and / or STN regions). For each of multiple training subjects, the MR scanner can acquire NM-MRI data (e.g., NM-MRI images) and / or QSM data (QSM images).

[0021] Method 100 may include a computing device or system generating an initial anatomical template based on selected MR images of a first training subject among a plurality of training subjects (step 104). The computing device or system (or its user) may 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 may be random or may follow a predefined preference or criterion. When generating the initial anatomical template, the computing device or system may enlarge the first 3D image (or corresponding 2D image) in a plane for all acquired slices of the first training subject by a factor of 2x, 4x or more, depending on the desired final resolution. This step may also include interpolation of the slice selection direction to obtain isotropic resolution in the template space. The template space may be defined at a desired resolution or a resolution higher 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 will be discussed in more detail below, a computing device or system can use an initial anatomical template to map MR images of other training subjects into the template space.

[0022] Method 100 may 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 may apply the global and local transformations to each of the images of the other training subjects (other than the first training subject). After scaling each MR image in the plane with a scaling factor for all slices (e.g., a factor of 2, 4 or more), the computing device or system may apply the global transformation over a set of slices covering the region of interest in each MR image of the other training subjects (other than the first training subject). For example, for the midbrain, if the slice thickness is 2 mm, the set 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 with b-spline interpolation applied, for example, using the Insight Segmentation and Alignment Toolkit Freeware (ITK). The computing device or system may apply the same global transformation to map QSM data for each of the other training subjects to an initial QSM template.

[0023] A computing device or system can apply a global transformation to each MR image (e.g., QSM images and / or NM-MRI images) of other training subjects to generate a corresponding first morphing version of the MR image that represents 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 does not usually provide an accurate 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 template space, the computing device or system can apply local transformations to the first morphed version of the MR image 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 a plane to cover the midbrain and into approximately 16 slices (for example, with slice thicknesses of 2 mm) to ensure coverage of the midbrain territory. The computing device or system can apply local transformations to the cropped image volume to generate a second morphed version of the MR image of the 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 template space.

[0025] As an exemplary implementation, a total of 26 training subjects can be used. A computing device or system can select MR data from 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 and, after applying the global and local transformations, map the MR data of each of the other 25 training subjects to the initial anatomical template in the template space.

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

[0027] Method 100 may include drawing the boundaries of the anatomical structure of interest on the final anatomical template (step 110). Tracing the boundaries of the anatomical structure (or region) of interest can be done manually, for example, by a physician or radiologist. The final anatomical template produced according to steps 102-108 of Method 100 will show enhanced contrast of the anatomical structure of interest and thus allow for more accurate tracing of the boundaries of such structure.

[0028] In the final anatomical template, the mean values ​​can be considered as a probability map for finding boundaries. 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 NM data and 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 neuromelanin (NM) data, boundaries can be drawn using slices 44-98 from a total of 192 interpolated slices, while slices 44-126 of the QSM data can be used. The actual selection of slice numbers depends on the resolution and degree of interpolation used.

[0029] Method 100 may 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 may include a computing device or system using a final anatomical template to identify the boundaries of an anatomical structure of interest in an untrained subject (step 112). The use of the final anatomical template to identify the boundaries of an anatomical structure of interest in an untrained subject is discussed in detail below with respect to Figure 2.

[0031] Referring here to Figure 2, a flowchart illustrating a method 200 for automated template-based identification of anatomical structure characteristics in accordance with the concept of the present invention as disclosed herein is shown. Method 200 may include acquiring MR images (or MR data) of an untrained subject (step 202). Method 200 may include mapping the MR images of the untrained subject to a final anatomical template in template space by applying a global transformation followed by a local transformation (step 204). Method 200 may include projecting the boundaries of one or more anatomical structures of interest from the final anatomical template onto the mapped MR images of the untrained subject (step 206). Method 200 may also include applying an inverse transformation to the projected boundaries of the anatomical structures of interest to generate estimates of the boundaries in the MR images of the untrained subject (step 208).

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

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

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

[0035] After applying a global transformation, the computing device or system can crop a first morphed version of the MR image of an untrained subject. For example, the computing device or system can label RNs in two central slices of template data and return them to the original whole brain. If RNs are still evident in only two slices, the lowest slice is set to slice 10. If RNs are in three slices, the middle slice is set to slice 10. This provides a means of best centering the data before the final transformation back to template space. The computing device or system can enlarge the cropped image volume in a plane by a factor of 4 (or other scaling factor), and then apply a local transformation to the enlarged cropped image volume to obtain a second morphed version of the MR image of an untrained subject with an isotropic image (the transformation, along with interpolation of the slice selection direction, allows the image to be isotropic). The second morphed version of the MR image of an untrained subject represents a relatively good match to the final anatomical template in template space. The interpolated in-plane resolution provides a more accurate means of estimating the final object size using the DPA algorithm.

[0036] Method 200 may include projecting the boundaries of one or more anatomical structures of interest from a final anatomical template onto a mapped MR image of an untrained subject (step 206). Once the MR image of the untrained subject is matched or mapped to a final anatomical template (by applying global and local transformations), a computing device or system may 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 untrained subject.

[0037] Method 200 may include applying an inverse transform to the projected boundary of an anatomical structure of interest to generate an estimate of the boundary in MR images of an 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., on the original image space) in the MR images of the untrained subject. However, since everyone is different and template mapping is not perfect, there is no guarantee that the projected boundary will fit well.

[0038] Method 200 may further include a computing device or system using a boundary detection algorithm (or DPA) to fine-tune the boundary of the anatomical structure of interest in the original acquired MR images of the untrained subject. For QSM data, the computing device or system may apply a threshold to the region inside the transformed boundary to remove negative values ​​and use image thinning to determine the corresponding centerline. For DPA, both Otsu histogram analysis and threshold-based approaches may be used to determine whether the original boundary extends excessively outside the structure of interest in order to best select the initial starting point. For NM data, the computing device or system may determine the background intensity and a constant equal to four times the average background standard deviation across all 25 trained subjects (or trained subjects other than the first trained subject), which is added to create a threshold below which signals are set to zero. For QSM data, the starting threshold may be set to zero. When the Otsu threshold gives a value that removes pixels inside the structure (or region) of interest, the template-transformed boundary may be reduced accordingly. 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 modify the resulting boundaries using the same DPA used in the template space (as discussed above with respect to Figure 1). This DPA, along with the RN boundary, prevents leakage of SN to the RN and, with the help of the template boundary, can distinguish STN from SN.

[0039] Referring to Figure 3, a diagram illustrating a system 300 for detecting or identifying the boundaries of anatomical structures according to the concept of the present invention as disclosed herein. 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 the boundaries of anatomical structures. 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. When the computer code instructions are executed by the processor 306, the computing device 304 can be caused to perform methods 100 and / or method 200 as considered above. The computing device 304 may further include a display device 310 for displaying MR images and outputting the results of methods 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, MR data acquired by the MR scanner 302 may be transferred to the computing device 304 via 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 into the MR scanner 302. For example, the processor 306, memory 308, and / or display device 310 can be integrated into the MR scanner 302. In such implementations, the MR scanner 302 can acquire MR data and perform methods 100 and / or 200.

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

[0043] System 300 or computing device 304 can calculate the overlap between the NM composite volume (SNpc and VTA added) and the iron containing the SN volume (SNpc and SNpr) by superimposing two ROIs from the MTC data and QSM data, respectively. System 300 or computing device 304 can normalize the overlap by iron containing the SN volume to create an overall overlap fraction (OOF) scale, which is essentially a measure of the SNpc fraction relative to the total SN volume, as shown below. JPEG0007915321000001.jpg16169

[0044] The process for generating the final anatomical template described in Figure 1 above was evaluated using MR data from various control subjects. The evaluation was based on a comparison of the performance of the final anatomical template generation process against manually drawn boundaries of the anatomical structures of interest.

[0045] A total of 87 healthy controls (HCs) were scanned. SNs, STNs, and RNs were manually traced for all 87 cases. Of these, 30 (test dataset: including 17 males and 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 with the next 57 cases (validation dataset: including 17 males and 40 females, age range 61.9 ± 5.0 years). Two measures were used to evaluate the performance of the anatomical template generation process: spatial overlap between structures associated with the manual segmentation method and the template segmentation method, and the volume ratio (VR) of structures from the template volume divided by the manual segmentation volume, as well as the Dice similarity coefficient.

[0046] Finally, all data were combined to generate quantitative information regarding structural volume, NM and iron content, and overlap between the SN and NM regions. Total iron content was calculated by summing the products of volume and the average sensitivity of the structure across all slices in which the structure was drawn. Similarly, total NM content was obtained as a result of summing the products of NM volume and NM contrast across the corresponding slices.

[0047] Referring to Figure 4, MR images are shown illustrating 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) to (c) illustrate various viewpoints of NM402 in the NM-MRI template using the validation dataset. The boundary of NM402 is drawn manually. Images (d) to (f) illustrate the viewpoints of SN404, RN406, and STN408. The boundaries of SN404, RN406, and STN408 in images (d) to (f) of Figure 4 are drawn manually.

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

[0049] Referring to Figure 6, images illustrating the various steps for mapping the boundaries of 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, labeled as row A, shows images of various slices of the constructed QSM template. The second row, labeled as row B, shows the same images with QSM DPA boundaries for SN602, RN604, and STN606. The third row, labeled as row C, shows the boundaries for SN602, RN604, and STN606 transformed into the original space. The final row, labeled as row D, shows images of various midbrain slices of the subject with the final boundaries of SN602, RN604, and STN606 obtained after applying DPA to the boundaries in row C.

[0050] For initial template training, 30 cases were processed for MTC and QSM data. The completeness of the template's automated MTC background intensity was such that the slope relating from one to the other was 0.99, and R 2 This is demonstrated by the fact that the value is 0.53 and the p-value is less than 0.001. Background values ​​are important for properly thresholding the NM and iron content signals.

[0051] The Dice similarity coefficient and volume ratio (VR) were measured for both MTC and QSM images against different thresholds. The data associated with the 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 the Dice similarity coefficient against the VR value of neuromelanin (NM) are shown based on various scenarios. 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 were generated using data from 30 training cases, while the plots in Figure 7B were generated using data from 57 validation cases. Mean Dice, volume ratio, and volume loss, associated with template data and manually drawn data, are cited for each threshold. A threshold of 1000 units minimizes volume loss and produces excellent results.

[0052] Referring to Figures 8A and 8B, plots of Dice similarity coefficients plotted against VR values ​​of SN are shown based on various scenarios. In both Figures 8A and 8B, plot A corresponds to the scenario with no applicable threshold, plot B corresponds to the threshold for sensitivity values ​​above 50 ppb, plot C corresponds to the threshold for sensitivity values ​​above 75 ppb, and plot D corresponds to the threshold for sensitivity values ​​above 100 ppb. The plots in Figure 8A were generated using data from 30 training cases, while the plots in Figure 8B were generated using data from 57 validation cases. Mean Dice, volume ratio, and volume loss, associated with template data and manually drawn data, are cited for each threshold. A threshold of 50 ppb results in minimal loss in SN volume.

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

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

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

[0056] Regarding Figures 9 and 10, before applying an arbitrary 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 threshold of 50 ppb to the QSM data results in average dice and VR of 0.95 and 1.04 for RN and 0.83 and 0.98 for STN. An average VR greater than 1 indicates that most structures found by the fully automated template / DPA approach tend to be larger than those found by the manual / DPA approach.

[0057] Referring to Figure 11, plots illustrating the correlations between iron content of SN, RN, and STN resulting from manual and template segmentation are shown. The correlations between iron content of SN, RN, and STN resulting from manual and template segmentation for 30 training cases are shown in the upper plot, while the correlations between iron content of SN, RN, and STN resulting from manual and template segmentation for 57 validation cases are shown in the lower plot. The final template results for iron content compared to manual plots are shown in various plots, using ROI associated with SN, RN, and STN, and a threshold of 50 ppb. For each of the SN, RN, and STN structures, the slope and R are shown.2 This 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] Regarding template validation, the validation dataset included 57 cases used to process QSM and MTC data. Dice similarity coefficients were plotted against both NM and SN VR values, and these results are shown in Figures 7B and 8B, respectively. For NM contrast using a threshold of 1000, the mean volume loss in the template data for all cases was approximately 5%, resulting in mean Dice and VR values ​​of 0.88 and 1.06, respectively. Similarly, for mean sensitivity of SN, a threshold of 50 ppb resulted in an mean volume loss of approximately 5%, while the mean Dice and VR values ​​were 0.89 and 0.97, respectively. Similar to the previous dataset, a threshold of 1000 for NM contrast and a threshold of 50 ppb for mean sensitivity of SN yielded sufficient results for mean Dice, VR, and volume loss across 57 cases.

[0059] The dice similarity coefficients for VR for RN and STN are shown in Figures 9 and 10, respectively. Before applying an arbitrary threshold to the QSM data, the mean dice and VR values ​​for RN and STN were 0.92 and 1.05, and 0.74 and 0.86, respectively. Applying a threshold improved these values, as was the case with the previous dataset. Using a threshold of 50 ppb for the QSM data resulted in mean 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 approximately 12% for both structures.

[0060] Figure 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 2This is 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 scale, and Dice similarity coefficient for each structure. The mean and standard deviation values ​​for the first and second datasets and the merged data are shown. [Table 1] Table 1. Mean and standard deviation (SD) values ​​of volume estimates, die similarity coefficients, and volume ratios (VR) for each structure.

[0061] In the studies described in Figures 4–12, NM and QSM images derived from a single sequence (less than 5 minutes) were used for automated 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 the return transformation from the template space is described and validated. Both dice values ​​and volume ratios showed good agreement for measures of volume and iron content between the automated 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, the MNI template is 1 mm isotropic, while the approach described herein uses interpolated 0.67 mm in-plane isotropic data, which is then interpolated to a 3D isotropic resolution of 0.167 mm. In practice, it was found that using either ANT or SpinITK with a 0.67 mm isotropic resolution for whole-brain global deformable registration did not yield morphological transformations that consistently matched the shape of the midbrain structure across the entire subject. 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, both iron and neuromelanin templates were created from a single multi-echo sequence.

[0063] Previous studies have attempted to segment SNs using structural MRI atlases based on T1 and / or T2 weighted sequences. However, SNs are small structures in the midbrain with low contrast in these images, making it difficult to accurately define 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 composed of NM-enhanced brain images for automated SN segmentation, one group achieved 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 SNs. NM contrast can be significantly improved by using MT-MRI acquisition sequences. To the inventors' knowledge, no studies have attempted to automatically segment SNs, SNpr, and SNpc using both NM-MRI and QSM.

[0064] Previous studies segmented SNs by including either younger (over 5 years and under 18 years) or older (over 60 years) subjects. The use of atlases based solely on young, healthy subjects may not be appropriate for studying patients with neurodegenerative diseases affecting midbrain structure. Age-related and disease-related morphometric changes between and within these subjects are significant and can impact the success of using template approaches when localizing SNs. The most appropriate atlas for a given study requires minimizing global or local distortion from the native subject space to the atlas space; therefore, an atlas for studying patients with neurodegenerative diseases was created using older HC subjects (mean age: 63.4 years). To mitigate bias caused by inter-subject variability, several studies have used stochastic atlases and attempted to investigate microstructural abnormalities of SNpcs using diffusion MRI in PD patients by constructing stochastic atlases of SNpcs based on NMS-MRI sequences. These researchers applied symmetric differential homeomorphism to align T1 images and 27HC SNpc masks in MNI space, thresholding the atlas with a 50% probability. However, the average Dice coefficient for atlas vs. human evaluators was less than 0.61. Thresholding can lead to dramatic changes in NM depiction and therefore reduces the reproducibility of the atlas.

[0065] Regarding the validation phase, the second dataset shows that the Dice similarity coefficients and VR values ​​for the structure are very close to those from the first dataset, thereby confirming the consistency of the results generated by the automated templating approach proposed herein. Furthermore, the fact that these values ​​are very close to a single value demonstrates the excellent performance of this approach.

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

[0067] In conclusion, the novel approach described herein for fully automated segmentation of deep gray matter in the midbrain (or anatomical structures in general) uses a template approach as a guide, but fine-tunes boundaries using the original data along with a dynamic programming algorithm. This approach makes it possible to quantify neuromelanin in the SN, as well as 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 a variety of neurodegenerative diseases without requiring manual tracing of these structures.

[0068] Regarding neuromelanin (NM) background measurement, since background values ​​are important for properly thresholding the NM signal, 30 MTC cases were processed to obtain neuromelanin (NM) background measurements. Referring to Figure 12, a plot showing the agreement between manual and templated neuromelanin (NM) background measurements for 30 healthy controls (or trained subjects) is shown. Specifically, the plot illustrates the completeness of the templated automated MTC background intensity scale over the manual plot. The agreement between the manual and templated MTC background scales is a slope of 0.99 and an R of 0.53. 2 , and a p-value of less than 0.001 is shown.

[0069] With respect to 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 may be as follows: 1. The initial boundary can be drawn on both NM and QSM template data. 2. Next, the background can be drawn onto the template image, as shown in Figure 13. In Figure 13, both NM702 and the background area 704 are shown. 3. Next, the boundary and background can be transformed back to the original space, as shown in Figure 14. In Figure 14, the NM region 802 and background region 804 are shown after transformation from the template space back to the original space. By using a background value increased by 1000, it becomes 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 + 4σ (where σ is the standard deviation of the background region of interest, which was found to be approximately 250 units in this study) can be used as a threshold for the MTC data to remove all points below this threshold and determine the NM refined boundary. For QSM data, pixels with sensitivity of less than zero ppb can be removed before DPA is run for all SN, STN, and RN. The initial boundary can be smoothed using a 5.3×3 Gaussian filter. 6. Next, the center line can be determined using the thinning method. 7. Boundary refinement before DPA: Since there may be some variation in intensity 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 a filtered local image of 40x40 pixels square around each single centerline point, and all points below 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 nearest remaining point (referred to as the Otsu boundary). This step helps to speed up the convergence of the DPA. 8. Implementation of DPA: For each point along the centerline, DPA can be executed in the corresponding search box. To enable curved shapes, the center of the centerline can be used to determine the initial rays. DPA can then be applied to the next set of rays by shifting one pixel along the centerline associated with the initial boundary. The centerline can be updated for each DPA iteration. Once the end point is reached, these rays can be swept through 180°. The algorithm can then travel back along the centerline until the centerline sweeps another 180° to reach the opposite end point. Finally, the center point travels back along the centerline to the starting point, and the boundary can be closed. This process can be repeated five times to explore both the inner and outer sides of the boundary to obtain optimal results. For NM, STN, and SN, after this step, a further five iterations can be continued using only inner exploration.

[0070] The search radius is limited to 4 pixels both inside and outside the boundary of the enlarged space, and the cost function used in this DPA includes a derivative term and a radius of curvature term that avoid leakage of the structure of interest into adjacent objects. For points outside the boundary, negative differential values were 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, JPEG0007915321000003.jpg1389In the formula, JPEG0007915321000004.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 derivative inside the image, and R avgis the average radius over the previous three radii. The constant α represents the relative weighting of the derivative term and the radius term, and 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 by selecting a large value of α. Since all α values ​​from 0.05 to 0.15 played a similar role in faithfully finding the edges, a value of α = 0.1 was selected for its conservative nature.

[0071] Referring to Figure 15, images illustrating the reconstructed boundaries with different values ​​of the DPA parameter α are shown. The selection of α can dramatically affect the final DPA boundary. In the upper left image, the NM boundary 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 a smaller value of α leads to less smoothing from the radius constraint, while a larger value of α leads to greater smoothing. Higher values ​​of α tend to result in a rounder structure that contracts more towards a circle. Lower values ​​lead to boundaries with very zigzag edges. Finally, candidate points for the new boundary can be selected by maximizing the cost function described above for each ray.

[0072] Several exemplary 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 five times. Image (d) is similar to image (b), except that noise with a 7:1 contrast-to-noise ratio (CNR) has been added based on the difference in signal intensity between the central and border areas within the second rectangular boundary. Image (e) shows the final boundary found after five iterations of DPA that match the correct area of ​​1400 pixels.

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

[0074] Figure 17 shows simulation results for identifying the boundary of a crescent-shaped region with two different intensities using DPA. Image (a) shows the crescent-shaped region with two different intensities, and image (b) shows that the same region was the 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 5 iterations of DPA. Images (d) and (e) show the centerline and boundary updates after 10 and 15 iterations of DPA, respectively. Image (f) shows the centerline and boundary 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 centerline and boundary after only 5 iterations of DPA when using the adaptive Otsu thresholding approach.

[0075] A crescent-shaped region was selected in Figure 17 to simulate a more challenging case of detecting the boundary of a curved object. Figure 17 shows the effect of running different numbers of iterations. The mean (image (f)) after running a total of 35 iterations was found to be 990.5 with a standard deviation of 1.5, while using the Otsu approach in the presence of a 10:1 SNR, the mean was 984 with a standard deviation of 2.0.

[0076] Figure 18 shows simulation results for identifying the boundary of a cashew-shaped region with two different intensities using DPA. Image (a) shows the cashew-shaped region with two different intensities, and image (b) shows the cashew-shaped region with an initial boundary drawn around the central area. Image (c) shows the cashew-shaped region with the corresponding centerline and boundary 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 result of the boundary and centerline after 5 iterations of DPA when using the adaptive Otsu thresholding approach and in the presence of noise is shown in image (e).

[0077] Finally, we selected a cashew-shaped region in Figure 18 to mimic the SNR and evaluated the method in the absence of sharp edges. After running the first 35 iterations, we found the mean 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 another 30 iterations, we found the mean to be 1097 with a standard deviation of 0.0.

[0078] Those skilled in the art will understand that the processes described in this disclosure can be implemented using computer code instructions that can be executed by a processor. These computer code instructions can be stored in a non-temporary or tangible computer-readable medium, such as memory. The memory may 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 comprising at least one processor and / or memory storing executable code instructions. The code instructions, when executed by at least one processor, can result in the execution of any of the processes or operations described in this disclosure. The apparatus may be, for example, an MRI scanner or a computing device.

Claims

1. It is a system, One or more processors, A memory containing computer code instructions, When the computer code instruction is executed by one or more processors, the one or more processors: To obtain MR images of a subject with contrast illustrating one or more anatomical structures of interest, Applying a global transformation to the MR image of the subject to generate a first morphing version of the MR image representing a first estimate of the subject's MR data in the template space. Applying a local transformation to the first morphing version of the MR image of the subject to generate a second morphing version representing a second estimate of the MR data in the template space, Projecting the boundaries of one or more anatomical structures of interest from an anatomical template image onto the second morphing version of the subject's MR image, and generating projected boundaries that represent estimates of the boundaries of one or more anatomical structures of interest in the subject's template space. A system that performs the following: applying an inverse transform to the projection boundary to determine a first boundary estimate of the corresponding boundary of one or more anatomical structures of interest in the MR image of the subject.

2. The one or more processors described above The system according to claim 1, configured to use a boundary detection algorithm to fine-tune the first boundary estimate of the corresponding boundaries of the one or more anatomical structures of interest in the MR images of the subject, and to generate a second boundary estimate of the corresponding boundaries of the one or more anatomical structures of interest.

3. The system according to claim 1, wherein the inverse transform includes the concatenation of the inverse of the local transform and the subsequent inverse of the global transform.

4. In order to apply the local transformation to the first morphing version of the MR image, one or more processors The system according to claim 1, configured to apply the local transformation to the cropping region of the first morphing version of the MR image.

5. In order to apply the local transformation to the first morphing version of the MR image, one or more processors The region of the first morphing version of the MR image is cropped, The cropping region is expanded, The system according to claim 1, configured to apply the local transformation to the enlarged cropping region.

6. The system according to claim 1, wherein the MR image of the subject includes a quantitative sensitivity map (QSM) or a neuromelanin image.

7. The one or more processors described above The system according to claim 1, configured to determine the volume of at least one of the one or more anatomical structures of interest based on the corresponding boundaries of the one or more anatomical structures of interest in the MR image of the subject.

8. The one or more processors described above The system according to claim 1, configured to determine the mean intensity of at least one of the one or more anatomical structures of interest based on the corresponding boundaries of the one or more anatomical structures of interest in the MR image of the subject.

9. The one or more processors described above The system according to claim 1, configured to determine the sum signal of at least one of the one or more anatomical structures of interest based on the corresponding boundaries of the one or more anatomical structures of interest in the MR image of the subject.

10. The system according to claim 1, wherein the anatomical template image is generated using multiple MR images of multiple training subjects with contrast illustrating one or more anatomical structures of interest.

11. It is a method, A computer system with one or more processors acquires MR images of a subject with contrast illustrating one or more anatomical structures of interest, The computer system applies a global transformation to the MR image of the subject to generate a first morphing version of the MR image that represents a first estimate of the subject's MR data in the template space. The computer system applies a local transformation to the first morphing version of the subject's MR image to generate a second morphing version representing a second estimate of the MR data in the template space. The computer system projects the boundaries of one or more anatomical structures of interest from an anatomical template image onto the second morphing version of the subject's MR image, and generates projected boundaries that represent estimates of the boundaries of one or more anatomical structures of interest in the subject's template space. A method comprising: applying an inverse transform to the projection boundary using the computer system to determine a first boundary estimate of the corresponding boundary of one or more anatomical structures of interest in the MR image of the subject.

12. The method according to claim 11, comprising using a boundary detection algorithm to fine-tune the first boundary estimate of the corresponding boundaries of the one or more anatomical structures of interest in the MR images of the subject, and generating a second boundary estimate of the corresponding boundaries of the one or more anatomical structures of interest.

13. The method according to claim 11, wherein the inverse transform includes the concatenation of the inverse of the local transform and the subsequent inverse of the global transform.

14. Applying the local transformation to the first morphing version of the MR image is The method according to claim 11, comprising applying the local transformation to the cropping region of the first morphing version of the MR image.

15. Applying the local transformation to the first morphing version of the MR image is Cropping the region of the first morphing version of the MR image, To expand the aforementioned cropping region, Applying the local transformation to the expanded cropping region, The method according to claim 11, including the method described in claim 11.

16. The method according to claim 11, wherein the MR image of the subject includes a quantitative sensitivity map (QSM) or a neuromelanin image.

17. The method according to claim 11, comprising determining the volume of at least one of the one or more anatomical structures of interest based on the corresponding boundaries of the one or more anatomical structures of interest in the MR image of the subject.

18. The method according to claim 11, comprising determining the mean intensity of at least one of the one or more anatomical structures of interest based on the corresponding boundaries of the one or more anatomical structures of interest in the MR image of the subject.

19. The method according to claim 11, comprising determining the sum signal of at least one of the one or more anatomical structures of interest based on the corresponding boundaries of the one or more anatomical structures of interest in the MR image of the subject.

20. A non-temporary computer-readable medium, which includes computer code instructions stored on the non-temporary computer-readable medium, and when the computer code instructions are executed by one or more processors, the one or more processors To obtain MR images of a subject with contrast illustrating one or more anatomical structures of interest, Applying a global transformation to the MR image of the subject to generate a first morphing version of the MR image representing a first estimate of the subject's MR data in the template space. Applying a local transformation to the first morphing version of the MR image of the subject to generate a second morphing version representing a second estimate of the MR data in the template space, Projecting the boundaries of one or more anatomical structures of interest from an anatomical template image onto the second morphing version of the subject's MR image, and generating projected boundaries that represent estimates of the boundaries of one or more anatomical structures of interest in the subject's template space. A non-temporary computer-readable medium that performs the following: applying an inverse transform to the projection boundary to determine a first boundary estimate of the corresponding boundary of one or more anatomical structures of interest in the MR image of the subject.

Images