Techniques for deformable registration of medical images

Abstract

Techniques for deformable registration of medical images. The invention relates to a technique for providing a deformable registration result of a second medical image to a first medical image. A first medical image comprising an anatomical structure is received. A second medical image comprising the anatomical structure is received. A nested iteration is performed comprising iterations over spatial resolution levels with increasing spatial resolution and iterations at each spatial resolution level with respect to a set of deformation parameters of an analytical deformable registration function for registering a deformation of the second medical image to the first medical image. The nested iteration is terminated at a predetermined fine spatial resolution level. A deformable registration result of registering the second medical image to the first medical image is provided comprising the set of deformation parameters in the deformable registration function after termination of the optimization.

Description

Technical Field

[0001] This invention relates to a technique for providing deformable registration results from a second medical image to a first medical image, specifically including a method, a computing device, a computer program product, and a computer-readable storage medium. Background Technology

[0002] Deformable registration is a key facilitator for image-guided interventions, treatment monitoring, and longitudinal studies. However, traditional methods fail to produce accurate and reliable results, hindering their application in real-world clinical workflows. While non-learning iterative methods offer more reliable and accurate registration results than learning methods for images of high quality, their computationally intensive nature is often prohibitive. On the other hand, although learning methods exhibit greater robustness to poor image quality, they typically lack the ability to handle large motions and model highly variable deformations. Furthermore, due to small training sets, these learning methods cannot handle motions in unseen test data where the distribution differs from the training data distribution. One major reason for the inability to handle large-variable motions is the insufficient number of model parameters per output dimension (e.g., the dimension of the non-parametric transformation (i.e., the displacement field) is three times the number of image voxels), resulting in insufficient model capacity to learn large distributions of possible local motions. However, the number of model parameters is limited by available GPU memory and input image size (i.e., image resolution). For example, a learning model can be trained on a small amount of images (allowing for a larger number of model parameters per output dimension), but the ability to solve motion problems involving detailed structures is lost (due to severe downsampling of the image size).

[0003] Various learning-based and non-learning-based methods have been proposed for deformable registration. A key attractive feature of learning-based methods is their insensitivity to some degree of image degradation and computational speed, as they involve only a single forward pass (single iteration) compared to non-learning methods. However, their poor registration performance on unseen data and their inability to operate at native image resolution limit their clinical applications.

[0004] To circumvent the problem of insufficient number of parameters per output dimension in deep learning (DL) methods for deformable registration, recent work has focused on iterative optimization of neural fields for test time [1]. Since they are iterative methods, their runtime tends to be longer than that of DL methods (which involve only a single forward propagation / iteration). While neural fields can be trained with a single coordinate point for other applications, for deformable registration, the application of penalties (e.g., smoothness, rigid motion of bones, incompressibility, etc.) requires deformability knowledge in the local neighborhood. Therefore, GPU / CPU memory still limits the size of neural networks and the number of input images. However, since the network only needs to solve the motion of a single image pair (compared to DL methods trained on multiple image pairs), neural field methods require fewer parameters per output dimension than DL methods. Because of this, and because neural field methods optimize the motion of each individual test data pair (compared to DL trained on training data rather than test data), neural fields show improved performance compared to DL methods [1]. However, both the DL method and the neural field method are memory-limited, which prevents their application in solving the motion of detailed structures such as lymph nodes, seminal vesicles, larynx, glottis, cochlea, and optic nerve.

[0005] Traditional methods [2-6] lack the ability to address large motions and abrupt changes in adjacent motions (in terms of amplitude and / or direction). They also cannot model motion based on tissue properties (like all neural field and learning methods), i.e., smooth, coherent motions that cause skeletal deformation and otherwise discontinuous sliding interfaces. Summary of the Invention

[0006] Therefore, one object of the present invention is to provide a solution for performing registration from one medical image to another, which is efficient in terms of computational resources such as memory and improves computational speed. Alternatively or additionally, one object is to provide a solution for improving the anatomical consistency and / or accuracy of registration from one medical image to another, which is capable of modeling large motions and / or large deformations, abrupt changes in adjacent motions and / or deformations, and / or motions and / or deformations based on tissue properties.

[0007] This objective is achieved by a method, computing device, computer program (and / or computer program product), and computer-readable storage medium for providing deformable registration results from a second medical image to a first medical image, as described in the appended independent claims. Advantageous aspects, features, embodiments, and advantages are described in the dependent claims and the following description.

[0008] In the following description, the invention will first be described with respect to the claimed method. Features, advantages, or alternative embodiments mentioned with respect to the method may be assigned to other claimed objects (e.g., computer programs or devices (especially computing devices) or systems or computer program products), and vice versa. In other words, a system, apparatus, or device may be modified using features described or claimed in the context of the method, and vice versa. In this case, the functional features of the method are embodied by the structural units of the apparatus or device or system, and vice versa. The method may refer to a software implementation, and the apparatus may refer to a hardware implementation (e.g., having a spatial physical structure) or its virtualization. Generally, in computer science, software implementations and corresponding hardware implementations (e.g., as embedded systems) are equivalent. Thus, for example, method steps for “storing” data may be performed using storage units and corresponding instructions for writing data into the storage device. To avoid redundancy, although the apparatus may also be used in alternative embodiments described with reference to the method, these embodiments are no longer explicitly described with respect to the apparatus. In principle, the corresponding apparatus or device claims are configured to perform the claimed method.

[0009] Regarding the method, a method (particularly computer-implemented) is provided for providing deformable registration results from a second medical image to a first medical image. The method includes the step of receiving a first medical image, which includes at least one anatomical structure. The method further includes the step of receiving a second medical image, which includes the at least one anatomical structure (particularly at least one identical or similar anatomical structure). The method further includes the step of performing nested iterations. The nested iterations include iterations at spatial resolution levels as the spatial resolution increases. The nested iterations further include, at each spatial resolution level, iterations with respect to a set of deformation parameters for an analytical deformable registration function (also simply: transfer function, deformation function, and / or deformation) used to deformably register the received second medical image to the received first medical image.

[0010] Nested iteration may include the step of determining the image dissimilarity between a received first medical image and a received second medical image for a predetermined spatial resolution level. Image dissimilarity per unit volume (or unit area) associated with the predetermined spatial resolution level may be determined. Nested iteration may further include the step of selecting a plurality of active control points for the predetermined spatial resolution level. The active control points may be located in regions of the first and / or second medical images where the determined image dissimilarity exceeds (or is at least equal to) a predetermined image dissimilarity threshold (particularly a level-specific predetermined image dissimilarity threshold). Nested iteration may further include the step of iteratively optimizing the set of deformation parameters by iteratively optimizing an objective function for the predetermined spatial resolution level. The objective function may be a differentiable function of an analytical deformable registration function. The analytical deformable registration function may include the sum of analytical functions centered at the selected active control points.

[0011] Iteration at the spatial resolution level may include proceeding to another predetermined spatial resolution level, which has a finer spatial resolution and / or a smaller unit volume (or unit area) than the previous predetermined spatial resolution level, and repeating the steps of determining image dissimilarity, selecting multiple active control points, and iteratively optimizing the set of deformation parameters at this other predetermined spatial resolution level. Iteration at the spatial resolution level may further include terminating the nested iteration if or when, for a predetermined spatial resolution level having the finest spatial resolution and / or a predetermined minimum unit volume (or unit area), the set of deformation parameters is iteratively optimized (and / or has reached its optimum according to the nested iteration).

[0012] The method further includes the step of providing a deformable registration result that registers a second medical image to a first medical image. This deformable registration result includes the use of the set of deformation parameters after optimization termination in a deformable registration function.

[0013] The technology of this invention improves the accuracy and speed of deformably registering one (e.g., second) medical image to another (e.g., first) medical image. Alternatively or additionally, computational resources (particularly processing and memory resources) can be saved, allowing deformable registration to be performed within short computation times (particularly real-time) and / or on any computing device (e.g., workstations in medical facilities, such as those in operating rooms and / or physician offices), particularly utilizing limited working memory space. This improves support for treatment planning (e.g., surgery and / or radiotherapy), optical guidance during treatment (e.g., surgery and / or radiotherapy), and / or post-treatment control (e.g., surgery and / or radiotherapy), and / or treatment monitoring. Alternatively or additionally, the comparison of two medical images for longitudinal and / or lateral studies can be improved.

[0014] Deformable registration using the techniques of this invention can be particularly accurate and / or reliable. Alternatively or additionally, optimization stability can be improved. Alternatively or additionally, the dependence of registration quality on image quality can be reduced, thereby allowing, in particular, high-resolution levels of deformation (and / or motion) of previously unseen data with detailed structures (e.g., the detailed structure of blood vessels compared to organs such as the liver or kidneys).

[0015] By iteratively optimizing deformation parameters (also known as transformation parameters) defined based on a number of selected active control points, the optimization problem can be simplified to computations related to multiple relevant points. Alternatively or additionally, by using a number of selected active control points, the iterative optimization of the deformation parameters can be independently parallelized and / or performed for each active control point (or at least a subset of the active control points). Parallelization of optimization can be combined in particular with shared memory to avoid duplication of storing deformation parameter values.

[0016] Alternatively or additionally, optimizing the deformation parameters iteratively based on the derivatives of analytical functions can avoid the problem of a small training dataset, which typically arises when performing deformable registration using machine learning (ML) and / or any other trainable method.

[0017] The objective function can be a combination of image similarity measures (or image dissimilarity measures), feature similarity measures (e.g., point similarity and / or grid similarity), and / or a penalty (and / or regularization) function that promotes physically meaningful structural alignment between the deformation of a first medical image and a second medical image (also referred to as: two input images). The penalty function (and / or regularization function) can typically be a function of the derivative of an analytical deformable registration function with respect to spatial location. Iterative optimization can involve updating deformation parameters, such as through the derivative of the objective function with respect to the deformation parameters.

[0018] According to the technique of the present invention, the objective function is not computed only at (e.g., active and / or inactive) control points. For example, the image dissimilarity metric can be computed "anywhere" (e.g., at all voxels or pixels, or at a predetermined subset of voxels or pixels) on the first medical image and on the deformation of the second medical image.

[0019] Alternatively or additionally, the objective function can be a combination of image similarity (dissimilarity) measures, point similarity measures, and penalty and / or regularization functions. Therefore, the objective function is not necessarily a function of the derivative of the deformation with respect to the deformation parameters.

[0020] The technology of this invention allows deformable registration to be performed in the presence of highly variable deformation, large “motion” and / or abrupt changes in the magnitude and / or direction of adjacent deformation (and / or motion).

[0021] Longitudinal studies can refer to tracking a patient's health status over time (e.g., in relation to at least one anatomical structure). In this case, the first and second medical images can be obtained from the same patient (e.g., a human or an animal, particularly a mammal, such as a dog or a horse). Alternatively or additionally, cross-sectional studies can refer to comparing the health status of different patients (e.g., in relation to at least one anatomical structure). In this case, the first and second medical images can be acquired from different patients.

[0022] At least one anatomical structure included in the first and second medical images may include at least one identical (and / or similar, and / or comparable) anatomical structure, which may be included in a set of anatomical structures. For example, the first and second medical images may only partially overlap, such as one (e.g., the first) medical image including the upper part of the patient's chest cavity, while the other (e.g., the second) medical image including the lower part of the patient's chest cavity, with some overlap in the imaged anatomical structures at the center of the patient's chest cavity.

[0023] At least one anatomical structure included in the first and second medical images may be similar and / or similar in size and / or shape across the entire population (or at least across a set of comparative candidates, such as patients similar in age, developmental status, and / or sex, e.g., adult males in a cross-sectional study). Alternatively or additionally, at least one anatomical structure included in the first and second medical images may be similar and / or similar for the same patient (and / or in a longitudinal study) in the sense that modifications due to surgery (e.g., removal of the patient's appendix or part of an organ such as the lung) and / or due to lesions (e.g., tumor growth and / or tumor shrinkage during treatment). For example, both the first and second medical images may provide images of the patient's organs (e.g., liver, kidneys, lungs, and / or heart), bone structures, and / or ducts (e.g., blood vessels).

[0024] Examples of anatomical structures include, but are not limited to, organs, bone structures, and / or ducts.

[0025] A first medical image (also referred to simply as the first image) can be used as a reference image for defining a deformable registration of a second medical image (also referred to simply as the second image) in order to map at least one (e.g., common, similar, or analogous) anatomical structure from the second image onto the first image. Alternatively or additionally, the first image may also be represented as a still image. Alternatively or additionally, the second image may also be represented as a deformable or moving image.

[0026] Using a deformable registration function, at least one (e.g., common, similar, or analogous) anatomical structure and / or organ included in the two images is deformed in the second image to match (or map) its location in the first image.

[0027] The deformation may include regions (and / or areas) where the deformable registration function acts nontrivially on the second medical image (e.g., not merely as a global rotation of the second medical image). Nontrivial action may include the deformable registration function taking different values ​​at different points in the neighborhood.

[0028] Examples of functions associated with deformable registration are radially symmetric functions (also known as radial basis functions), such as the impulse function, Gaussian function, Vanderland function, Wu function, Busserman function, and / or oscillatory function (also known as Gnetin function).

[0029] One way to view the distortion between the second and first images is as motion, such as due to breathing of soft tissues, and / or rotation (especially rigid rotation), such as the skeletal structures used to map the two images onto each other. For example, the first image (e.g., a CT scan) can be acquired when the patient is inhaling, the lungs are fully expanded, and / or the liver is in a lower position within the thoracic cavity. Alternatively or additionally, the second image can be obtained when the patient is exhaling and / or the liver is in a higher position within the thoracic cavity. Such images can be obtained, for example, during radiotherapy.

[0030] In the context of this invention, the terms "motion" and "deformation" are essentially used synonymously. Alternatively or additionally, "moving image" may correspond to an image that has been deformed to register it to another image (such as a second medical image) (such as a first medical image and / or a reference image).

[0031] In some cases, prior information about structural deformation (e.g., caused by motion) and / or physiology is available. For example, imaging results acquired for diagnostic purposes (e.g., the acquisition of a first image) can be considered as pre-planned information for a second image. Prior information may include prior therapy information (e.g., surgical area and / or applied radiation dose), contrast agent dosage, contrast agent type, and / or motion trajectory.

[0032] The deformation can cover a range, for example, corresponding to up to approximately 15 cm in real life. For instance, the deformation can cover the range affected by changes in lung volume between complete inhalation and complete exhalation.

[0033] Deformation (and / or motion) may be subject to one or more physical constraints (also known as: penalty terms and / or regularization terms), such as the rigidity of bone structures and / or restrictions on soft tissue deformation, such as the relative motion of adjacent organs or restrictions on organ volume changes.

[0034] Segmentation can be based on (e.g., semantics) using physical constraints related to organization type. For example, the first and second images can be segmented separately (and / or independently) and / or include segmentation masks (also known as segmentation labels).

[0035] Each medical image can be acquired using a medical imaging modality and / or a medical scanner. The scanner can be configured to provide 3D medical images (also known as volumetric images), which can correspond to a set of slices of a 3D medical image, such as a CT scan and / or an MRT scan.

[0036] Each (e.g., first and / or second) medical image may have been acquired independently using a 3D medical imaging modality. The 3D medical imaging modality may be (or may include) computed tomography (CT), magnetic resonance imaging (MRT), positron emission tomography (PET), single-photon emission computed tomography (SPECT), magnetic particle imaging (MPI), ultrasound (including, for example, echocardiography, elastography, and / or photoacoustic imaging), and / or 3D thermography.

[0037] Deformable registration can register a second 3D medical image to a first 3D medical image. Alternatively or additionally, it is possible to register a slice of the second 3D medical image to a first 2D medical image using the techniques of the present invention. It can be assumed that the motion is substantially 2D and / or that a slice of the second 3D medical image (also referred to as: a 2D slice) is mapped to the first 2D medical image. Alternatively or additionally, the organ motion of the second 3D medical image can (in particular, still) be 3D motion, and / or the slice can be selected for optimal mapping.

[0038] 2D medical imaging modalities may include radiography (also known as X-ray imaging) and / or scintigraphy.

[0039] The techniques of this invention can be used for various medical imaging modalities. Alternatively or additionally, the techniques of this invention incorporate computationally efficient algorithms for multimodal deformable registration with physical constraints. Parameterization of terms such as contributions to the objective function can depend on the medical imaging modality (or, if the first and second medical images are acquired using different medical imaging modalities, on the medical imaging modality pair).

[0040] In some embodiments, the first and second images are acquired using the same (e.g., type) medical imaging modality and / or the same scanner. For example, two CT scans can be registered to each other, which can be achieved using the same scanner or using different scanners (e.g., with different technical specifications and / or different technical parameters). In other embodiments, the second image is acquired using a different medical imaging modality than the first image. For example, the first image may be a CT scan, while the second image may correspond to an MRT scan.

[0041] A deformable registration function can transform a set of input points (specifically, the set of input points of a medical image) into a set of output points (specifically, the set of output points in the domain of another medical image). For example, specifically by applying a deformable registration function, the set of input points of a second medical image can be transformed into the set of output points in the domain of a first medical image. Alternatively or additionally, specifically by applying the inverse of a deformable registration function, the set of input points of a first medical image can be transformed into the set of output points in the domain of a second medical image.

[0042] A deformable registration function can be an analytic function of the set of deformation parameters. Analytic functions are typically functions locally given by a convergent power series.

[0043] Deformable registration functions may include, in particular, linear terms, which consist of the sum of basis functions multiplied by torsion coefficients, affine matrices (and / or affine deformation parameters, and / or affine transformation parameters), and offset terms.

[0044] Affine deformation is typically deformation that can be fully described by an affine transformation. Such transformations generally consist of linear transformations (such as rotation, shear, stretching, and / or compression) and rigid body translations. Affine deformation can also be referred to as uniform deformation.

[0045] Basis functions may have local support around points, particularly around (especially) control points of initial activity and / or around local support around landmarks included in (e.g., first and / or second) images, and / or may be radial basis functions (RBF).

[0046] Alternatively or additionally, basis functions can have a tensor structure whose size is determined by the number of input points multiplied by the number of (in particular, active) control points.

[0047] (Especially radial) basis functions can be derived from a central and scalar value. To parameterize, the center can be represented as a control point (c) in the context of this invention. i ), the scalar value Define the width of the brace (also simply referred to as the brace and / or width). A tight brace RBF (CSRBF) can have a bell-shaped curve and / or a convex shape. Alternatively or additionally, the parameters required to define the CSRBF can be similar to those required to define a Gaussian function (especially including the center c). i and width ). Scalar value For example, you can define the width of a CSRBF.

[0048] In the context of this invention, an active control point can represent that the pre-multiplication factor of the radial basis function centered at that control point is being optimized, and the pre-multiplication factor can be represented as a torsion coefficient.

[0049] In the context of this invention, the inactivity of a control point can indicate that the optimization process of the pre-multiplication factors and / or torsion coefficients of the radial basis functions centered at that control point has been optimized and / or has reached a constant value. Alternatively or additionally, after the control point becomes inactive, the corresponding torsion coefficients can be regarded as fixed values ​​(and / or frozen parameters) within the deformable registration function.

[0050] As iterations progress to increasingly finer spatial resolution levels, initial (and / or earlier) active control points may become inactive. In some embodiments, new active control points may emerge, which may be associated with more refined (and / or smaller) deformation requirements.

[0051] In some examples, newly activated control points may coincide in location with control points that have been inactive so far. The centers of the newly activated and inactive control points may be different basis functions. For example, their widths may be different, and / or their twist coefficients may be different. In particular, to account for more refined deformation structures, newly activated control points may contribute smaller values. With the maximum value of the basis functions normalized to a fixed value (e.g., normalized to 1), more refined deformations may correspond to smaller twist coefficient values.

[0052] The marker can be a point (and / or region) that is identifiable in each (e.g., the first and / or the second) image. Alternatively or additionally, the marker can be a point that needs to be mapped to the same point in both the second and first images during deformation.

[0053] Radial basis functions (RBFs) are typically real-valued functions. Its value depends only on the distance between the input and a fixed point, which can be the origin, such that... It could also be some other fixed point. (referred to as the center (and / or, in particular, the activity control point)) so that Satisfies the property Any function These are all radial functions. The distance is typically Euclidean, although sometimes one or more other metrics can be used. They are often used as a set to form the basis for some function space of interest. .

[0054] Deformable registration functions can be spline-based, particularly B-spline-based. Splines are typically functions defined piecewise by polynomials.

[0055] B-splines, or basis splines, are typically spline functions with minimal support for a given order, smoothness, and / or domain partitioning. Any spline function of a given order can be represented as a linear combination of B-splines of that order.

[0056] Alternatively, or additionally, examples of RBF may include a center point as And the width is The impact function, .

[0057] In the following text, the arrow at the top of the vector is often omitted. Alternatively, or additionally, this can be used for multiple (N) control points. The basis functions are usually expressed as Omitting the arrow at the top of the vector results in a vector having support (and / or width). Alternative symbols for the basis functions of N control points .

[0058] Multiple active control points can vary throughout the nested iterations. Alternatively or additionally, active control points can be unevenly distributed across the first and / or second images. In particular, active control points can be densely distributed within misaligned (and / or dissimilar) regions of the first and (especially deformed) second images. Alternatively or additionally, in well-aligned (and / or similar) regions of the first and (especially deformed) second images, (especially active) control points can be sparse (and / or non-existent).

[0059] The initial number of active control points can be on the order of 1000 or more.

[0060] Multiple active control points can be a sparse set. Alternatively, or additionally, control points do not need to be located at the center of an image voxel.

[0061] For any newly selected active control point, a basis function can be chosen (e.g., in terms of the shape and / or width of the support). In one example, each basis function can be the same (e.g., the impact function) until different support widths are selected.

[0062] The basis functions can be determined before executing the method (also known as before running the application). In some examples, the user may have already selected (e.g., by means of manual input on a user interface (UI) such as a graphical user interface) which basis functions will be used, and only the user-selected basis functions will be selected throughout the optimization process and / or to determine deformable registration. Therefore, which basis functions(s) are used may depend on a one-time definition by the user before implementing the techniques of this invention.

[0063] The technique of this invention allows for the simultaneous estimation of deformation (and / or deformable registration function) and inverse deformation (and / or the inverse of deformable registration function). Alternatively or additionally, regardless of the output mapping direction (e.g., from a second medical image to a first medical image, or vice versa), the technique of this invention may involve only one set of active control points. The control points may be defined, for example, in the domain of the second medical image, and an initial transformation (and / or initial deformation) may be used to map them to the domain of the first medical image, or vice versa. In other words, for the technique of this invention, it is not important in which domain the control points to be worked are defined. The control points can be mapped between the two domains using the initial transformation (and / or initial deformation). In other examples, at the beginning of each optimization iteration, such as for the next spatial resolution level, the (and / or refined) transformation (and / or deformation) determined up to that point is used to map (particularly active) control points between the two domains.

[0064] Nested iterations may include a first type of iteration for optimizing deformation parameters at a predetermined spatial resolution level (also known as: multi-resolution level) for the first and second images. At the first (and / or coarsest-grained) spatial resolution level, the deformation parameters may be initialized by a deformable registration function that acts as an identity function (and / or global scaling, in the case that the size of the second image is different from that of the first image).

[0065] For example, for a deformable registration function with affine transformation, offset, and twist coefficients (also known as twist coefficients and / or twist factors) as deformation parameters, the initial values ​​can correspond to constant matrices as affine transformation matrix (also simply affine matrix), zero offset, and vanishing twist coefficients.

[0066] Nested iterations can include a second type of iteration by moving from the coarsest spatial resolution level (also known as: the coarsest granularity) to a finer spatial resolution level (also known as: the finer granularity) until a predetermined finest spatial resolution level (also known as: the finest granularity, e.g., corresponding to voxels or pixels of the first image and / or the second image) is reached.

[0067] Spatial resolution levels (e.g., consecutive spatial resolution levels) can include a multiplication factor of the spatial resolution (e.g., a factor of 2 per image orientation, and a factor of 2 for 3D images). 3 =8, and / or for 2D images, the factor is 2. 2 =4), for example, the number of voxels per unit volume (or the number of pixels per unit area) per predetermined spatial resolution level.

[0068] The optimization of deformation parameters can be initialized at each subsequent finer spatial resolution level using the optimization results from the previous coarser spatial resolution level.

[0069] For previously active control points that become inactive, the final spatial resolution-level optimization result from the twist coefficients (and / or widths) of the basis functions can correspond to the final (and / or frozen) result.

[0070] New active control points can be introduced, which allows for optimization of new twist coefficients (and / or widths).

[0071] If no active control points exist, and / or if there is no longer an area (e.g., per unit) of image dissimilarity exceeding a (particularly level-specific) predefined image dissimilarity threshold, then the iterative optimization of the set of deformation parameters at the predetermined spatial resolution level can end (also known as: stop and / or terminate).

[0072] The level specificity of a predefined image dissimilarity threshold may mean that the requirement for image similarity increases as the spatial resolution level becomes finer. Alternatively, or additionally, the predefined image dissimilarity threshold can be decreased with increasing spatial resolution level. Thus, image regions that are sufficiently similar at a coarse spatial resolution level (e.g., after optimizing deformation parameters for multiple selected active control points) may have caused some previously active control points to become inactive. By decreasing the predefined image dissimilarity threshold at a finer spatial resolution level, it may be necessary to select new active control points to refine the deformable registration function.

[0073] Alternatively or additionally, iterative optimization may terminate if the number of iterations reaches a predefined maximum number of iterations, and / or if one or more stopping criteria have been met. One or more stopping criteria may include minimizing the change in the magnitude of the derivative of the objective function (and / or the change from one step to the next within an iteration being below a predetermined threshold). Alternatively or additionally, the objective function value may reach a certain threshold (e.g., for minimizing, it may approach 0). Further alternatively or additionally, minimizing the change in the objective function value per iteration (and / or the change being below a predetermined threshold) may signify convergence.

[0074] A unit volume (or unit area) at any given level of spatial resolution (and / or granularity) can be provided based on the spatial expansion and / or the number of voxels (or pixels). The finest (and / or final) level can correspond to a minimum spatial expansion of one (or a few) voxels (or pixels) and / or unit volume (or unit area). Alternatively or additionally, the coarsest (especially the initial) level can correspond to a large number of voxels (or pixels) and / or a large spatial expansion of unit volume (or unit area).

[0075] Active control points can be adaptively selected, and / or can be selected based on nested iterations. In one embodiment, multiple active control points are selected (e.g., independently) at each spatial resolution level. In another embodiment, multiple active control points are selected (e.g., independently) for each iteration of deformation parameter optimization at a predetermined spatial resolution level.

[0076] Multiple selected active control points can be iteratively reduced. Alternatively or additionally, active control points can be defined by local values ​​of image anisotropy that are at least equal to (or exceed) a predefined image anisotropy threshold. When performing nested iterations (e.g., at a predetermined spatial resolution level), local values ​​of image anisotropy can decrease as deformation parameters are optimized, causing previously active control points to become inactive. Inactive control points can be discarded and / or selected (also known as removed) from additional steps of the nested iterations (e.g., optimizing additional deformation parameters by considering derivatives). By discarding inactive control points, significant processing power can be saved. It is important to note that discarding from future optimizations is equivalent to retaining the latest optimized values ​​of deformation parameters (e.g., the torsion coefficients related to the basis function width) at the currently inactive control points in shared memory. Alternatively or additionally, storing control points and deformation parameters in shared memory can reduce the storage space required to perform optimizations, especially when using parallel computation for different control points. Different control points can, for example, be spaced apart by a distance greater than the sum of their supports. In this case, deformable registration at one control point is independent of deformable registration at other control points.

[0077] A predefined image dissimilarity threshold can be dependent on the spatial resolution level. For example, the predefined image dissimilarity threshold can decrease as the spatial resolution level increases. Alternatively or additionally, multiple active control points can be selected (e.g., independently) for each predetermined spatial resolution level. For example, an active control point selected based on a first (and / or initial, and / or maximum) predefined image dissimilarity threshold may become inactive at the end of iterative optimization at the first (and / or coarsest) spatial resolution level. At the next (particularly finer) spatial resolution level, a second (and / or lower) predefined image dissimilarity threshold can cause previously active control points to be reactivated, and / or result in newly selected active control points.

[0078] The first and second medical images can be acquired using independently selected medical imaging modalities.

[0079] By enabling the registration of images from different medical imaging modalities, it is possible to select the most appropriate modality for their function (e.g., CT or MRT image acquisition for treatment planning, and real-time image acquisition during treatment using CT, X-ray imaging, or ultrasound). Alternatively or additionally, the elimination of the need to acquire the first and second images using the same medical imaging modality facilitates cross-sectional studies of different patients and / or longitudinal studies of a single patient, where medical images are acquired by scanners available at the time. This in particular allows for comparison of newer images acquired using the most modern scanners with older images acquired using scanners that were established and / or available at the time.

[0080] For example, a CT scan can be acquired in seconds. The time between acquiring a CT scan for treatment preparation and the start of treatment can vary considerably, ranging from hours to days and / or weeks. Therefore, the time required for image acquisition is usually not a limiting factor, but the time from image availability to clinical decision-making is critical. If information from previous imaging needs to be considered, it may be necessary to make this information reliably available in “daily anatomy,” and depending on the clinical task, this may need to happen quickly.

[0081] Comparing older hardware to newer hardware is not necessarily the primary objective. The quality of images acquired throughout the clinical pathway can depend on acquisition techniques and spatial and / or temporal constraints. For example, MRT acquisition may be performed prior to treatment, but in the treatment room, due to space constraints and / or because the patient cannot move, only orthogonal X-rays (e.g., for acquiring the first medical image) may be available. Therefore, it may be advantageous to attempt to bring better 3D information (e.g., as a second medical image) into the treatment room prior to treatment by using (especially deformable) registration.

[0082] To enable comparison and / or registration of images (particularly images from different medical imaging modalities), a modality-independent neighborhood descriptor (MIND) can be used. For each voxel (or pixel) in an image, the MIND can comprise a real-valued vector describing the local structure around that voxel (or pixel). This vector may, for example, include the gray level of each voxel (or pixel) and / or the voxel location (or pixel location). The MIND can be determined independently for the images to be (particularly deformably) registered (e.g., first and second). The MIND is described, for example, by Mattias P. Heinrich et al. [8], which is incorporated herein by reference.

[0083] Alternatively or additionally, for comparing and / or registering images (particularly images from different medical imaging modalities), segmentation of each image can be used. For example, the result of segmentation can be a multi-mask comprising (e.g., binary) segmentation masks for different organs. (e.g., first and second) image (e.g., deformable) registration can be initialized based on mapping matching segmentation masks onto each other.

[0084] The first and second medical images can be acquired using the same medical imaging modality.

[0085] For images that are similar in terms of image parameters (e.g., frequency, opacity, image contrast, and / or color) and / or image resolution, it may be simpler to determine the dissimilarity and / or perform a registration procedure.

[0086] The iterative optimization of this set of deformation parameters can be performed using parallel computing. Parallel computing can include distributing the iterative optimization across parallel computing resources based on selected active control points.

[0087] A deformable registration function may comprise the sum of basis functions with localization and / or tight support. Alternatively or additionally, basis functions may be located at different (e.g., active) control points and will not (or at least not significantly) affect the deformable registration function values ​​at other (e.g., active) control points. Optimization of deformation parameters (particularly the twist coefficients (also known as weighting coefficients) associated with the basis functions) can therefore be advantageously performed independently of each active control point. Alternatively or additionally, independent optimization at different active control points can be implemented in parallel on different processing units (e.g., central processing unit CPU, graphics processing unit GPU, and / or different cores of CPU and / or GPU). Thus, the optimization process can make optimal use of existing computing hardware, and / or the optimization process can be accelerated, for example, by determining multiple twist coefficients simultaneously.

[0088] Specifically, at a predetermined spatial resolution level, multiple active control points can be adaptively selected at each iterative optimization step and / or after a predetermined number of iterative optimization steps. Alternatively, or additionally, multiple active control points can be adaptively selected for each predetermined spatial resolution level.

[0089] After each iterative optimization step, or after a predetermined number of iterative optimization steps (e.g., after every second or fifth iterative optimization step), the image dissimilarity between the deformed second and first medical images can be re-evaluated (and / or re-determined). If the image dissimilarity in the region where the control point is located drops below a predefined image dissimilarity threshold, the previously active control point can be set to inactive. This may be particularly due to the distortion coefficients (and / or widths, also known as support widths) of the corresponding basis functions being changed (and / or optimized).

[0090] In iterative optimization at a predetermined spatial resolution level, the number of active control points can be reduced as the number of iterative optimization steps performed increases.

[0091] Alternatively or additionally, multiple active control points can be reset and / or newly selected when moving from one (e.g., the previous) predetermined spatial resolution level (e.g., with coarse grain) to the next predetermined spatial resolution level (e.g., with a finer grain). In particular, the new selection can be performed based on a lower predefined image dissimilarity threshold for the next predetermined spatial resolution level compared to the previous predetermined spatial resolution level.

[0092] Therefore, multiple active control points can be selected as needed to further improve the use of computing resources (e.g., processing resources and / or memory resources). Alternatively or additionally, whenever a control point is set to inactive, the associated deformation parameters (e.g., torsion coefficient and / or support width) do not need to be further optimized, thereby freeing up processing resources and / or eliminating the need for any further updates in shared memory.

[0093] Alternatively, or additionally, the control points for the supports can be obtained by querying the basis functions based on the spatial tree and / or graph structure. Only the basis functions need to be evaluated simultaneously for those control points. Therefore, the optimization process at a single point can include more than one torsion coefficient and / or more than one support width. This is especially likely to occur if the distance between two control points is less than the sum of their corresponding support widths.

[0094] Selecting multiple active control points can be based on point sampling. Optionally, performing iterations at (especially increasingly finer) spatial resolution levels can include repeated point sampling.

[0095] Point sampling can be a method of determining an appropriate tree for sampling within a sampling or cruise system. Alternatively or additionally, point sampling can be a method of distributing active control points across (e.g., connected) regions that exceed (and / or at least reach) a predefined image dissimilarity threshold.

[0096] The number of points that meet the control point criteria (e.g., those exhibiting high dissimilarity values) may be greater than the number of control points specified by the user. In this case, it is necessary to sample active control points from all control point candidates. For example, random sampling or sampling using the farthest point sampling method can be used.

[0097] Deformable registration functions may include spline functions, particularly B-spline functions. Alternatively or additionally, deformable registration functions may include the sum of basis functions, particularly radial basis functions (RBF), and more particularly compactly supported RBF (CSRBF)).

[0098] Deformable registration functions can be locally represented as a sum of polynomials (such as spline functions) or (particularly) radial basis functions. Advantageously, deformable registration functions can be extended to a sum of basis functions supported only in the local neighborhood of the center and / or (e.g., active) control points. For example, deformable registration functions can include B-spline functions and / or can be scalable in terms of CSRBF. This facilitates parallel computation for (particularly simultaneously) determining deformation parameters associated with different active control points.

[0099] This set of deformation parameters may include an affine transformation matrix, offset, a set of torsion coefficients, and / or a set of associated support widths.

[0100] This set of twist coefficients can also be represented as a twist matrix. Alternatively or additionally, this set of twist coefficients can correspond to the coefficients of a twist function that is scalable in terms of RBF. Any spatial (and / or non-general) dependencies can be encoded in affine transformation terms that can be identical for all points (and / or volumes, areas, and / or regions) of the second image, and can be encoded in the sum of RBFs, each of which depends on the input point, and in the case of CSRBF, advantageously each of which is non-zero only in a local neighborhood around the control point (e.g., by being parameterized by the support width). Thus, parallelization of parameter optimization is facilitated, particularly parallelization of the twist coefficients and / or support width.

[0101] In some embodiments, the basis functions to be used (specifically the CSRBF to be used) and a set of (e.g., discrete) support widths can be selected before performing deformation parameter optimization. Alternatively or additionally, the support width is not necessarily a parameter to be optimized.

[0102] The set of deformation parameters (e.g., including the support width, even if the technique according to the invention is pre-selected and / or not optimized) and / or the selected multiple control points (in particular including active control points and / or inactive control points) can be stored in shared memory.

[0103] Alternatively or otherwise, the input and output point sets of the deformable registration function can be stored in shared memory.

[0104] By using shared memory, each quantity only needs to be stored once. Alternatively or additionally, any optimizations and / or updates performed in a parallel computation become (particularly immediately) accessible to all other optimization processes.

[0105] In shared memory, input and output points can be stored. Input points can be points used in determining (e.g., calculating) terms in an objective function (particularly point similarity and / or penalty functions). Input points can be used in each iteration to compute the objective function and update deformation parameters (and / or to optimize the update function). Alternatively or additionally, input points can be used to determine image dissimilarity, for example, at each level of spatial iteration. Input points can, for example, be uniformly distributed across one of the medical images (e.g., a second medical image or a first medical image).

[0106] Input points are typically different from control points. Control points can be part of a transformation function, specifically serving as the center points of basis functions. Alternatively or additionally, input points can be transformed by a transformation function (which can be defined solely based on control points, since the control points define the transformation) to produce output points. Output points can be used together with the corresponding input points to determine (specifically, to compute) the objective function and transformation parameter updates.

[0107] According to the technique of the invention, input points and control points are advantageously stored in shared memory throughout the deformation (and / or registration) process. According to the technique of the invention, output points are also stored in shared memory, but their values ​​can change, such as changing with each iteration.

[0108] Multiple active control points can be selected for one image, such as a second medical image. Alternatively or additionally, multiple active control points for another image (e.g., a first medical image) can be determined by applying a deformable registration function (and / or its inverse). The mapping of (particularly active) control points between the first and second images can be performed, particularly at an initial (and / or coarsest) resolution level, using initial values ​​of the deformation parameters. In some embodiments, the mapping of (particularly active) control points between the first and second images can be performed at any given finer spatial resolution level, using optimized values ​​of the deformation parameters from the previous (and / or the next coarser) spatial resolution level.

[0109] The supports for the radial basis functions can correspond to scalar values, thus defining the deformation function together with the control points. Further definitions of the deformation function can include affine transformations, offsets, and / or torsion coefficients.

[0110] The objective function may include a measure of point cloud similarity and / or a measure of image similarity.

[0111] Point cloud similarity can be replaced by continuous, especially analytic, functions to replace discrete point similarity.

[0112] The objective function can be the sum of two or more terms. For example, point cloud similarity and / or image similarity can be measured by various analytical functions such as mean squared error (MSE), mean distance (MD, also denoted as L2 distance) and / or another function such as Gaussian mixture of discrete point sets, deformable registration function and / or its inverse.

[0113] Dissimilarity can be a negative (and / or inverse) of similarity. For example, a zero value for MSE or MD can correspond to maximum similarity and / or minimum dissimilarity.

[0114] The inverse of a deformable registration function can deform a first image into a second image. The inverse of a deformable registration function (and / or its deformation parameters) can be determined entirely by the deformable registration function (and / or the corresponding deformation parameters, for example, for a preset basis function and / or multiple sets of support widths).

[0115] Therefore, the objective function for optimization can correspond to maximizing the similarity between images (and / or points on images). Alternatively or additionally, the objective function for optimization can correspond to minimizing the dissimilarity between images (and / or points on images).

[0116] Optimizing the objective function may include imposing a set of physical constraints, which includes at least one physical constraint. Optionally, imposing the set of physical constraints may include including corresponding physical constraint terms in the objective function. Alternatively or additionally, imposing the set of physical constraints may include optimization using corresponding Lagrange multipliers.

[0117] Physical constraints (also known as penalties and / or regularization terms) can ensure that rigid structures, such as bone, do not deform in physically impractical ways. Alternatively or additionally, constraints can be imposed on deformable registration regarding the relative positions of different tissue types and / or different organs and / or their compressibility (and / or volume changes, e.g., due to respiration and / or phases in the patient's cardiac cycle). This can thereby improve the accuracy and / or reliability of deformable registration.

[0118] Each physical constraint can be added to the unconstrained objective function as one or more terms (e.g., two, one for the deformable registration function and the other for its inverse). Alternatively or additionally, the Lagrange multiplier method can be used to modify the objective function.

[0119] The optimization objective function may be symmetric because it depends on the deformable registration function and its inverse. Alternatively or additionally, the optimization (and / or nested iterations) may be independent, for example, at least in terms of which image is interpreted as a still image and which image is interpreted as a moving image, and / or in terms of which of the two images is deformed.

[0120] This set of physical constraints (especially those with at least one physical constraint) may include affinity or bending energy constraints, marker matching, sliding motion constraints, skeletal rigidity constraints, tissue incompressibility constraints, and / or topology preservation constraints.

[0121] Physical constraints can improve well-posedness, optimize stability, and / or the physical rationality of deformation results.

[0122] Affinity or bending energy constraints can also be expressed as smoothness constraints.

[0123] Landmark matching can utilize identifiable structures (and / or landmarks) in two (e.g., first and second) images. This can improve the geometric accuracy of deformable registration. For example, landmarks can be provided through segmentation of the first and second images. For instance, in cases where the patient is known to have an implant (e.g., a stent), the implant can be used as a landmark to register the second image to the first image, or vice versa. Alternative examples of landmarks include anatomical structures such as heart valves.

[0124] Sliding motion constraints can be used to model organ interfaces using surface normal vectors at surface points. For example, the lungs and / or bladder can move relative to adjacent organs according to predetermined rules.

[0125] Skeletal stiffness constraints can impose rigid motion and / or rotation matrices that are orthogonal to the skeletal structure within the image. For example, it can be assumed that the patient's spine is not moving.

[0126] Tissue incompressibility constraints can impose limitations on compression and / or expansion based on tissue type (e.g., soft tissue). Tissue incompressibility constraints can depend on the type of tissue (and / or the organ to which the tissue belongs).

[0127] Topology preservation constraints can also be expressed as reverse consistency constraints.

[0128] The dissimilarity of an image or point cloud can be determined by the mean distance (MD), mean squared error (MSE), gradient information, and / or the mutual information of a combination thereof.

[0129] MD can be minimized to obtain the maximum similarity (e.g., image and / or point cloud).

[0130] Alternatively or additionally, MSE can be minimized to obtain maximum (e.g., image and / or point cloud) similarity. MSE may be appropriate if the first and second images were acquired using the same medical imaging modality.

[0131] Gradient information can be maximized to obtain maximum (e.g., image and / or point cloud) similarity. Gradient information can be independent of the medical imaging modality. In particular, gradient information can be determined if the medical imaging modalities of the first and second images are selected independently.

[0132] The mutual information of the combination can be maximized to obtain the maximum image similarity.

[0133] Alternatively or additionally, point cloud anisotropy can be based on a mixture of parameters, particularly (e.g., density) weights, widths, and / or the center of radial functions (e.g., Gaussian functions). If the uncertainty at a point is high, the weights (e.g., density weights) may be low, and vice versa.

[0134] The technique of this invention can utilize the first derivative at any point on the second image (and / or the first image). The basis functions, and especially their (particularly first-order) derivatives, can each take a constant value at predetermined points. Therefore, it is only necessary to calculate the basis function value and the value of each derivative once for each point on the image. The set of values ​​for each point can then be stored in shared memory and used for optimizing deformation parameters.

[0135] In alternative examples, the second derivative (and / or Hessian derivative) of the deformable registration function and / or the objective function can be used.

[0136] Receiving a first medical image and receiving a second medical image may include receiving segmentation masks, particularly multiple masks, for the first medical image and for the second medical image, respectively.

[0137] Performing nested iterations may include the steps of receiving an indication of a downsampling factor and / or an initial predetermined spatial resolution level.

[0138] Alternatively or additionally, performing nested iterations may further include the step of downsampling the first medical image and / or the second medical image to achieve an initial predetermined spatial resolution level. Alternatively or additionally, performing nested iterations may further include the step of upsampling the first medical image and / or the second medical image to achieve an additional predetermined spatial resolution level.

[0139] To perform nested iterations, an indication of the downsampling factor can be received (e.g., 2 for numbers n=3, 4, or 5). nFor each spatial direction, the downsampling factor can be different (and / or independent). The downsampling factor can be the maximum downsampling factor to achieve the initial predetermined spatial resolution level.

[0140] For each finer spatial resolution level, execution can be performed according to a predetermined factor (e.g., for a spatial orientation of 2). 1 Upsampling can be performed independently for each spatial direction.

[0141] Thus, each spatial resolution level can be determined, for example, until the finest (and / or final) spatial resolution of the first medical image and / or the second medical image is reached.

[0142] This method can be used to support image-guided interventions, for treatment monitoring, for image comparison, and / or for radiotherapy planning.

[0143] This method can be executed by a computing device.

[0144] This method (especially in computer implementation) can be executed by multiple processors (e.g., 10 CPUs or 10 GPUs), each communicating with shared memory. The processors and memory can be implemented as hardware within a computing device. These processors can operate independently and / or in parallel on different threads.

[0145] Regarding the device, a computing apparatus is provided for providing deformable registration results from a second medical image to a first medical image. The computing apparatus includes a first input interface configured to receive a first medical image, which includes at least one anatomical structure. The computing apparatus further includes a second input interface configured to receive a second medical image, which also includes at least one anatomical structure. The computing apparatus further includes a nested iteration module configured to perform nested iterations. The nested iterations include iterations at spatial resolution levels as the spatial resolution increases. The nested iterations further include iterations at each spatial resolution level with respect to a set of deformation parameters for an analytical deformable registration function used to deformably register the received second medical image to the received first medical image.

[0146] Nested iteration may include: determining, for a predetermined spatial resolution level, the image dissimilarity between a received first medical image and a received second medical image. Image dissimilarity per unit volume (or unit area) associated with the predetermined spatial resolution level may be determined. Nested iteration may further include: selecting multiple active control points for the predetermined spatial resolution level. These active control points may be located in regions of the first and / or second medical images where the determined image dissimilarity exceeds (or is at least equal to) a predefined image dissimilarity threshold (particularly a level-specific predefined image dissimilarity threshold). Nested iteration may further include: iteratively optimizing the set of deformation parameters by iteratively optimizing an objective function for the predetermined spatial resolution level. The objective function may be a differentiable function of an analytical deformable registration function. The analytical deformable registration function may include the sum of analytical functions centered at the selected active control points.

[0147] Iteration at spatial resolution levels may include: proceeding to another predetermined spatial resolution level, which has a finer spatial resolution and / or a smaller unit volume (or unit area) than the previous predetermined spatial resolution level; and repeating the steps of determining image dissimilarity, selecting multiple active control points, and iteratively optimizing the set of deformation parameters at this other predetermined spatial resolution level. Iteration at spatial resolution levels may further include: terminating the nested iteration if (or when) the set of deformation parameters has been iteratively optimized for a predetermined spatial resolution level having the finest spatial resolution and / or a predetermined minimum unit volume (or unit area).

[0148] The computing device further includes an output interface configured to provide a deformable registration result that registers a second medical image to a first medical image. This deformable registration result includes the utilization of the set of deformation parameters in a deformable registration function after optimization termination.

[0149] Nested iterative modules can include multiple processors suitable for parallel computing.

[0150] The computing device can be configured to perform any of the steps disclosed in the context of the method, and / or may include any of the features disclosed in the context of the method.

[0151] In another aspect, a computer program product is provided, comprising program units that, when loaded into the memory of a computing device, cause the computing device to perform steps of a method according to this aspect for providing deformable registration results from a second medical image to a first medical image.

[0152] In another aspect, a computer-readable medium is provided having program units stored thereon that can be read and executed by a computing device, such that when the program units are executed by the computing device, the steps of the method according to this aspect of the method for providing deformable registration results from a second medical image to a first medical image are performed.

[0153] The nature, features, and advantages of the invention described above, as well as the ways in which they are implemented, will become clearer and more understandable in light of the following description and embodiments, which will be described in more detail in the context of the accompanying drawings.

[0154] The following description does not limit the invention to the included embodiments. In different figures, the same components or parts may be labeled with the same reference numerals. Generally, the figures are not drawn to scale.

[0155] It should be understood that the preferred embodiments of the present invention may also be any combination of the dependent claims or the foregoing embodiments with their respective independent claims.

[0156] These and other aspects of the invention will become apparent and will be explained with reference to the embodiments described below. Attached Figure Description

[0157] Figure 1 This is a flowchart of a method for providing deformable registration results from a second medical image to a first medical image according to a preferred embodiment of the present invention; Figure 2 This is an overview of the structure and architecture for providing deformable registration results from a second medical image to a first medical image according to a preferred embodiment of the present invention; Figure 3 This schematically illustrates a computing architecture that performs parallel computations using shared memory. The computing device could be... Figure 2 The equipment; Figure 4 The diagram illustrates the storage of Figure 2 and / or Figure 3 Data in the shared memory of the computing device, including input and output points, control points, and deformation parameters used to determine image anisotropy; Figure 5 (a) Figure 5 (b) Figure 5 (c) and Figure 5 (d) illustrates an example of the adaptive control point distribution overlaid on the MIND map and the first medical image for two different spatial resolution levels; and Figure 6 The technique according to the invention is shown, in particular using Figure 1 An example comparing the organization-related deformation of the physical constraints of the method with traditional techniques lacking consistency constraints.

[0158] Any reference marks in the claims should not be construed as limiting the scope. Detailed Implementation

[0159] Figure 1 An exemplary flowchart, particularly a computer-implemented method, is schematically illustrated for providing deformable registration results from a second medical image to a first medical image. This method is typically referred to by reference numeral 100.

[0160] Method 100 includes a step S102 of receiving a first medical image. The first medical image includes at least one anatomical structure. Method 100 further includes a step S104 of receiving a second medical image. The second medical image includes (e.g., the same and / or similar) at least one anatomical structure. Method 100 further includes a step S106 of performing nested iterations. The nested iterations include iterations at spatial resolution levels as the spatial resolution increases. The nested iterations further include iterations at each spatial resolution level with respect to a set of deformation parameters of an analytically deformable registration function used to deformably register the received second medical image of S104 to the received first medical image of S102.

[0161] Figure 1 The nested iteration S106 in the example includes a sub-step S106-D, for a predetermined spatial resolution level, determining the image dissimilarity between the received first medical image of S102 and the received second medical image of S104. The image dissimilarity per unit volume (or unit area) is determined in relation to the predetermined spatial resolution level. The nested iteration S106 further includes a sub-step S106-S, for a predetermined spatial resolution level, selecting a plurality of active control points. These active control points are located in regions of the first and / or second medical images where the image dissimilarity determined in S106-D exceeds (or is at least equal to) a predefined image dissimilarity threshold (particularly level-specific). The nested iteration S106 further includes a sub-step S106-O, for a predetermined spatial resolution level, iteratively optimizing the set of deformation parameters by iteratively optimizing an objective function. The objective function is a differentiable function of an analytical deformable registration function. The analytical deformable registration function comprises the sum of analytical functions centered on the selected active control points of S106-S.

[0162] Figure 1The iteration at the spatial resolution level includes a sub-step S106-P proceeding to another predetermined spatial resolution level, which has a finer spatial resolution and / or a smaller unit volume (or unit area) than the previous predetermined spatial resolution level, and repeating the steps of determining the image dissimilarity in S106-D, selecting multiple active control points in S106-S, and iteratively optimizing the set of deformation parameters in S106-O at this other predetermined spatial resolution level. The iteration at the spatial resolution level further includes a sub-step S106-T: if the set of deformation parameters has already been iteratively optimized in S106-O for the predetermined spatial resolution level with the finest spatial resolution and / or with the predetermined minimum unit volume (or unit area) (and in particular, no further optimization is performed when repeating steps S106-D, S106-S, and / or S106-O), then the nested iteration terminates.

[0163] Method 100 further includes step S108 of providing a deformable registration result that registers the second medical image to the first medical image. The deformable registration result includes utilizing the set of deformation parameters in the deformable registration function after optimization termination S106-T.

[0164] Multiple active control points S106-S can be adaptively selected during each iterative optimization step and / or after a predetermined number of iterative optimization steps, particularly at a predetermined spatial resolution level. Alternatively or additionally, multiple active control points S106-S can be adaptively selected at each predetermined spatial resolution level.

[0165] Selecting multiple active control points S106-S can be based on point sampling. Optionally, performing iterations at (especially finer) spatial resolution levels includes repeated point sampling.

[0166] Performing the nested iteration S106 may further include a sub-step S106-I of receiving an indication of a downsampling factor and / or an initial predetermined spatial resolution level. Performing the nested iteration S106 may further include a sub-step S106-DS of downsampling the first medical image and / or the second medical image to achieve the initial predetermined spatial resolution level. Performing the nested iteration S106 may further include a sub-step S106-US of upsampling the first medical image and / or the second medical image to achieve an additional predetermined spatial resolution level.

[0167] Figure 2 An exemplary architecture of a computing device for providing deformable registration results from a second medical image to a first medical image is schematically illustrated. The computing device is generally referred to by reference numeral 200.

[0168] The computing device 200 includes a first input interface 202 configured to receive a first medical image. The first medical image includes at least one anatomical structure. The computing device 200 further includes a second input interface 204 configured to receive a second medical image. The second medical image includes (e.g., the same or similar) at least one anatomical structure. The computing device 200 further includes a nested iteration module 206 configured to perform nested iterations. The nested iterations include iterations at spatial resolution levels as the spatial resolution increases. The nested iterations further include iterations at each spatial resolution level with respect to a set of deformation parameters of an analytically deformable registration function used to deformably register the received second medical image to the received first medical image.

[0169] The nested iterative module 206 may include an image dissimilarity determination submodule 206-D, configured to determine the image dissimilarity between a received first medical image and a received second medical image. This determination is made per unit volume (or unit area) in relation to a predetermined spatial resolution level.

[0170] The nested iteration module 206 may further include an active control point selection submodule 206-S, which is configured to select a plurality of active control points. These active control points may be located in the following regions of the first medical image and / or the second medical image, where the determined image dissimilarity exceeds (or is at least equal to) a predefined image dissimilarity threshold, particularly specific to the level.

[0171] The nested iterative module 206 may further include an iterative optimization submodule 206-0, configured to iteratively optimize the set of deformation parameters by iteratively optimizing an objective function. The objective function may be a differentiable function of an analytical deformable registration function. The analytical deformable registration function may include the sum of analytical functions centered on the selected active control points.

[0172] Iteration at the spatial resolution level may include proceeding to another predetermined spatial resolution level, which has a finer spatial resolution and / or a smaller unit volume (or unit area) than the previous predetermined spatial resolution level, and repeating the steps of determining image anisotropy, selecting multiple active control points, and iteratively optimizing the set of deformation parameters at this other predetermined spatial resolution level. An additional processing submodule 206-P may be configured to proceed to this other predetermined spatial resolution level.

[0173] Iteration at the spatial resolution level may further include terminating nested iterations if the set of deformation parameters has been iteratively optimized for a predetermined spatial resolution level having the finest spatial resolution and / or a predetermined minimum unit volume (or unit area). The termination submodule 206-T can be configured to terminate nested iterations.

[0174] The computing device 200 further includes an output interface 208 configured to provide a deformable registration result of registering a second medical image to a first medical image. The deformable registration result includes the use of this set of deformation parameters in a deformable registration function after optimization has been terminated.

[0175] The nested iteration module 206 may further include a downsampling factor indication receiving submodule 206-I, which is configured to receive an indication of the downsampling factor and / or an initial predetermined spatial resolution level.

[0176] The nested iteration module 206 may further include a downsampling submodule 206-DS, which is configured to downsample the first medical image and / or the second medical image to achieve an initial predetermined spatial resolution level.

[0177] The nested iteration module 206 may further include an upsampling submodule 206-US, which is configured to perform upsampling of the first medical image and / or the second medical image to achieve an additional predetermined spatial resolution level.

[0178] Alternatively, any module can be represented as a cell, and vice versa.

[0179] The computing device 200 may include an input / output (I / O) interface 212. The I / O interface may be embodied as a first input interface 202, a second input interface 204, an output interface 208, and / or an optional downsampling factor indication receiving submodule 206-I.

[0180] The computing device 200 may include a set of parallel processors 206, which may (e.g., each) be embodied as nested iterative modules 206 (and / or any of their optional submodules 206-D; 206-S; 206-O; 206-P; 206-T; 206-DS; 206-US).

[0181] The computing device 200 may include a memory 210. The memory 210 may be shared memory.

[0182] The computing device 200 can be configured to execute method 100.

[0183] The system may include a computing device 200 and a database, such as a Picture Archiving and Communication System (PACS) database, from which a first medical image and a second medical image may be received. Alternatively or additionally, the system may include a scanner from which the first medical image and the second medical image may be received.

[0184] The system can be configured to execute method 100.

[0185] The techniques of this invention (e.g., including method 100 and / or computing device 200) provide a rapid implementation of non-learning deformable registration for optimizing parameter transformations under various physical constraints. To allow for the resolution of large and / or varying motions while minimizing memory usage, the techniques of this invention estimate deformation based on a set of active control points, which can be determined in-run based on the quality of image alignment in their local neighborhoods. Alternatively or additionally, to accelerate optimization, derivative-based optimization with analytical gradients is performed on all terms in the objective function. This contrasts with traditional learning methods whose performance is currently limited by CPU / GPU memory.

[0186] The techniques of this invention (e.g., including method 100 and / or computing device 200) focus on accelerating non-learning iterative methods. Previous work has introduced methods to accelerate nonparametric deformable registration methods (i.e., displacement field estimation) that involve at most smoothness constraints (e.g., diffusion regularization and differential homeomorphism regularization) [2-6]. These methods lack the ability to handle large motions and abrupt changes in adjacent motions (in terms of amplitude and / or direction). They also cannot model motion based on tissue properties (like all neural field and learning methods), i.e., smooth, coherent motions of sliding interfaces that cause deformation of the skeleton and discontinuities in other respects. Unlike prior work, the techniques of this invention can accelerate constraint-parameter-based methods, taking advantage of the properties of parametric deformation, the parallelism of the problem, and / or analytical computation. The techniques of this invention can be implemented in both CPUs and GPUs, for example using OpenMP for the former and CUDA for the latter.

[0187] The techniques of the present invention (e.g., including method 100 and / or computing device 200) can alternatively be represented as parallel computation utilizing analytical derivatives and deformable registration of constraint parameters based on misaligned adaptive control points.

[0188] The technique of this invention can accelerate parametric deformable registration. This method is applicable to spline-based deformations (such as those using B-splines and / or radial basis functions), i.e., deformations of the following form: (1), in It is an affine matrix. Indicates offset, and It has N control points and distortion coefficient Related support The basis functions. Equation (1) covers a variety of deformation models, including but not limited to those involving a single torsion coefficient matrix and those with global support (i.e., Transformation of basis functions of ) and only involving deformation or The proposed multiprocessing method can be implemented in a CPU (using, for example, OpenMP) or a GPU (using, for example, CUDA).

[0189] In equation (1), the input point and output point Each is written as a d-dimensional vector (d=3 for 3D images and d=2 for 2D images), and the sum of the basis functions at the control points is written as a product of matrices, i.e., in, yes The i-th row, for those with the same center point Different radial basis functions centered on Implicit in support index and distortion coefficient In the sum of .

[0190] It can be seen that equation (1) is independent of the input point. (and / or for any input point) (Both are valid). This means that different processing and / or compute nodes and / or cores (e.g., cores in CPU and / or GPU units) can be used to process each input point individually. As an example, Figure 3 This diagram schematically illustrates a basic computing architecture found in typical modern computers. For illustrative purposes, Figure 3 The exemplary architecture shown in the diagram illustrates three processors 306-1; 306-2; 306-3, each having four computing cores 308 that are data-connected to shared memory 210.

[0191] exist Figure 3In the example, for a conventional serial program where each point is processed sequentially, executing a single process using only a single core 308 each time is highly inefficient. Conversely, according to the technique of the present invention, computation can be accelerated using hardware 206;306-1;306-2;306-3 by distributing the evaluation of each point (e.g., x) in equation (1) to different processes (and / or cores 308) (referred to as multiprocessing). Alternatively or additionally, multiprocessing may include different control points. And / or (especially tightly supported) basis function assignment computation. An example is provided in Algorithm III.

[0192] The simplest way to use equation (1) to transform the set of input points in a multiprocessing scheme is to apply it to all... Deformation parameters , , and Copy the data to all processes involved in transforming the set of input points. This typically means that each process must allocate a large memory space to store all the data. of , , and This copying is highly redundant and inefficient, and may be infeasible if any matrix is ​​too large. Furthermore, because all of , and These are all parameters to be optimized in deformable registration optimization, so their values ​​need to be updated iteratively. Each process has a copy of the parameter update values, which, by convention, increases computation time. Therefore, in addition to increasing memory usage of redundant data, having a separate copy of the information for each process increases runtime.

[0193] In some embodiments, the support width is a user-defined value entered (e.g., with the help of a UI, such as a GUI) before deformable registration is performed.

[0194] alternative sites , , , and It can be stored in shared memory because equation (1) is only read (and / or used). , , , and (But not modifying and / or writing). Stored in shared memory. , , , and The technology of this invention benefits from the transformation of the input point set. Required reduction in memory footprint, where index This indicates "input". This processing method follows a shared-memory multiprocessing architecture supported by various multiprocessing libraries (e.g., OpenMP, CUDA), such as... Figure 4 As illustrated in the example.

[0195] Figure 4 This schematically illustrates the use of equation (1) to connect multiple input points. Simultaneously transformed into an output point Shared memory multiprocessing. Note the point set. and index and Let "input" and "output" be simply represented, respectively. For the deformation function used to deformably register the second medical image to the first medical image, the set of input points... The points correspond to those within the original first medical image. Alternatively, or additionally, the output set of points. The point corresponding to the second medical image is intended to maximize similarity to the first medical image, or vice versa.

[0196] exist Figure 4 In this context, 1…K represents each… Flattened index, such as .therefore, The i-th row and j-th column correspond to the flatten index. Furthermore, unlike the j-th input vector used throughout the entire document... symbols ,exist Figure 4 In China, for Figure 4 The output vector and control points, coordinate entries of the three-dimensional vector (and / or d=3). They have similar notation. Algorithm I shows methods for, for example... Figure 4 The shared-memory multiprocessing architecture described in the text integrates multiple input points. Simultaneously converted to an output point The following is an example pseudocode. Although Example Algorithm I assumes that multiprocessing is implemented via OpenMP (i.e., using #pragma omp parallelism), other multiprocessing libraries such as CUDA can also be used.

[0197] Algorithm I: Simultaneous Spline-Based Transformation of Multiple Input Points

[0198] In addition to performing equation (1) in parallel on multiple input points. (and / or (This is for each point) and each control point Entries with scalar values In addition to the vector, according to the technique of the present invention, each scalar value entry It can be computed in parallel because each It is independent. Alternatively, or otherwise, assuming there are N computational units, computation is performed on a single computational unit. This would require N sequential computations, or it could be computed on each machine (especially independently). Each element The complete vector can be obtained by performing the calculation only once on each machine. .

[0199] In addition, for each support (and / or support width) of Independent of other supports, it can be independently of, as shown in Exemplary Algorithm II All support and All control points Simultaneously calculate a set of basis functions .

[0200] Algorithm II: Multiprocessing of Basis Functions

[0201] In addition, for tight support RBF (CSRBF). It is sparse, because if At the control point In addition to the support, The evaluation is zero. Aside from the speed and memory usage benefits of sparse matrix operations (especially compared to dense matrix operations), not all... All of them need to be evaluated.

[0202] Exemplary Algorithm III combines spatial trees / graphs to query its support cover. The control points, and only these selected control points are evaluated. .

[0203] Algorithm III: Multiprocessing of Compactly Supported Basis Functions

[0204] Note that for each input point Each control point And each (especially the fixed) support , It is fixed and / or constant, regardless of changes in optimization iterations, multi-resolution levels, or deformation parameter values. Therefore, It can be computed once and stored for use in subsequent optimization iterations and the next multi-resolution level.

[0205] Note that in Algorithm III, Indicates size is Multiply by a tensor of N, where N is the number of control points. It is a scalar value entry.

[0206] In some embodiments, the support is fixed and set by the user before performing the deformable registration technique of the present invention. According to other examples, adaptive and / or automatically determined supports are also possible.

[0207] Combining Algorithm I and Algorithm III may result in Algorithm IV. If RBF is globally supported, Algorithm III in lines #5-6 can be replaced with Algorithm II in lines #3-4, or if If it is pre-computed and stored elsewhere, then Algorithm III in #5-6 can be replaced by a read operation.

[0208] Algorithm IV: Simultaneous Spline-Based Transformation of Multiple Input Points

[0209] Note that in Algorithm IV, yes ( A vector. In other words, in Algorithm III It is a scalar, and in Algorithm IV yes vector.

[0210] Similar to other registration methods, the transformation in equation (1) is estimated by optimizing the objective function involving data terms and regularization terms. To further accelerate optimization and reduce memory usage, the technique of this invention uses a derivative-based optimizer based on analytical derivatives (e.g., compared to gradient-free optimization methods or gradient-based optimization methods using automatic gradient computation). In particular, the first derivative is used for optimization.

[0211] Optimizations based on the second derivative (e.g., involving Hessen) are also possible, but the computational complexity may outweigh the optimization speed and stability.

[0212] Below, examples of data terms and regularization terms that can be used to optimize equation (1) according to the technique of the present invention are provided. Other data terms and regularizations not described below can also be integrated into the technique of the present invention. Hereinafter, dissimilarity measures used as data terms and penalties used to regularize deformations (and / or as physical constraints) are described. Analytical derivatives of these terms with respect to the parameters of equation (1) are provided, including A, t, and , .

[0213] An active control point scheme is described that can help reduce memory requirements while maintaining or improving alignment accuracy.

[0214] Given a set of points from a moving image domain (and / or a second medical image). and a set of points from a fixed image domain (and / or a first medical image, for example, serving as a reference image). Instead of using a distance function that requires the distance between discrete points corresponding to known points (such as L...), p The norm is used to measure the distance between them. Each set of points can be modeled as a mixture and can be represented as follows: L 2 To measure the dissimilarity between the two transformed mixtures. (2), in, It is Points in the map are mapped to The transformation, and Will Mapped to .For example, It can be a deformable registration function of equation (1). That is the opposite.

[0215] L represents the mixture 2 Distance, which will be further described in the context of point cloud anisotropy below. Alternatives to or appended to point cloud anisotropy can minimize the transformed moving image. (e.g., second medical image) and fixed image Image dissimilarity (such as mean squared error) between (e.g., the first medical image): (3) Further alternatives or alternatives can maximize gradient information [7], .

[0216] It is worth noting that equation (4) is a modification of the gradient information proposed in [7] to improve the stability of the optimization.

[0217] Alternatively to or supplemented by image dissimilarity and point cloud dissimilarity, the technique of this invention can be combined with marker point matching constraints, which are defined as: (5), in and These are the second medical images (and / or ) and first medical images (and / or The set of known markers in ().

[0218] Alternatively or otherwise, various penalties (and / or physical constraints) can be imposed to constrain the solution space, thereby improving well-posedness, optimization stability, and physically plausible deformation results.

[0219] The penalty can include the following per-slide structure Calculated sliding motion constraints, (6), in Represents surface points Surface normal at the location, yes Displacement at that point yes about The derivative of, and yes and surface points The minimum distance between them. The sliding motion constraint on the surface can be derived in a similar way. In equation (6), It is the set of surface points (and / or vectors) of an organ that can slide relative to other organs.

[0220] Alternative sites or other options, for Rigid motion of the skeletal structure can be applied using the following formula (also known as: rigid skeletal constraint). , And can be for Similar constraints can be derived. To reduce computational overhead, they can be defined using only one Frobenius norm term from equation (4) (e.g., the first or last term). .

[0221] Alternatively, or otherwise, tissue incompressibility can be applied via the following formula: (8).

[0222] Alternatively, or additionally, smoothness can be applied via affinity or bending energy penalty as follows: (9), in It is the basis function matrix of the control points. It can be derived in a similar way for... of and .

[0223] U is typically calculated between the input and control points. The matrix in equation (9) is calculated in essentially the same way. It simply replaces the input points with control points. In other words, K is at the control point. Calculated between [the two].

[0224] Alternatively or otherwise, topology preservation can be imposed via a reverse consistency constraint as defined below. (10), in as well as Similar to To reduce computational overhead, this can be achieved using a single term from equation (7) (e.g., the first or last term). .

[0225] Therefore, the overall objective function can be read as: .

[0226] As previously mentioned, other data items for measuring image dissimilarity and / or other regularizations can be readily integrated into the technology of this invention. Hereinafter, parameters (including) for transformations (and / or deformation functions) of the form defined in equation (1) are provided. , The analytic derivatives of each loss and penalty term of (A and t).

[0227] Point cloud heterogeneity As an alternative to discrete dissimilarity, discrete point set It can be represented as a mixture of continuously distributed elements (e.g., Gaussian mixture). , in It is a density weight that can be used to reduce the weights of points with high uncertainty, and It has a mean Covariance Gaussian density, .

[0228] Transformation The set of moving points after (For example, representing a second medical image) and a set of fixed points The dissimilarity between mixtures (e.g., representing the first medical image) can be calculated as L. 2 distance, , in It is a set of mixed parameters, including .focus on The derivative, because It can be similarly derived, given: (12).

[0229] The first term is independent of the transformation parameters (i.e., its derivative with respect to the transformation parameters is zero). Therefore, only the last two terms contribute to the derivative of (12). For Gaussian mixtures, the following Gaussian mixture identity can be used. .

[0230] The second term in equation (12) can be written as: (13) in Represents a set The j-th point in, and Represents a set The i-th point in the equation, and assuming the normalization term It is absorbed into the variance. Similarly, the last term in equation (12) becomes: (14).

[0231] Equation (13) regarding The derivative is: (15), in .

[0232] The derivative of equation (13) with respect to A is: (16).

[0233] The derivative of equation (13) with respect to t is: (17).

[0234] Similarly, equation (14) regarding The derivative is: (18) in, .

[0235] The derivative of equation (14) with respect to A is: (19).

[0236] Notice .

[0237] Mean square error Mean square error (MSE) measurement in equation (3) (and / or the first medical image) and the transformed (and / or a second medical image) The dissimilarity between them, and and the transformed ( The dissimilarity between ) is assumed and They were acquired from the same imaging modality. Because and The derivative of the MSE is similar to and The derivative of the MSE between them is taken, so we only focus on the derivative of the MSE between them. and The MSE derivative between (i.e., the first term in equation (3)) leads to the following result. Regarding The derivative of the mn-th element is: .

[0238] Similarly, the derivative with respect to the mn-th element of A is: , Furthermore, the derivative with respect to the m-th element of t is: .

[0239] Here, image area, image grid, and / or image domain, and / or Can be with same.

[0240] gradient information Similar to equation (3), the gradient information (GI) in equation (4) includes two terms. Unlike equation (3), GI is measured... (and / or the first medical image) and the transformed (and / or a second medical image) The similarity between them, and and the transformed ( The similarity between them, without assuming and The data were acquired from the same imaging modality. Since the derivatives of the first and second terms can be derived in a similar manner, only the derivative of the first term is provided.

[0241] .set up and .

[0242] GI about The derivative of the mn-th element is: (20) in , and .

[0243] Similarly, derivative Following equation (20), where , and .

[0244] Finally, the derivative Similarly, it follows equation (20), where and .

[0245] Marker point matching constraints The marker matching constraint defined in equation (5) includes two terms, which respectively apply the estimated transformation. and Geometric accuracy. Because each term about the parameter can be derived in a similar way. and The derivative of , so only information about The derivative of equation (5) with respect to equation (1). The derivative is: Similarly, the derivative with respect to A is: , And the derivative with respect to t is: .

[0246] Sliding motion constraints The sliding motion constraints defined in Equation (6) can be used to model organ interfaces that slide relative to each other, such as the lung sliding relative to the ribs and / or the liver sliding relative to the spine and nearby soft tissue organs (e.g., the kidneys). In the following sections, only Equation (6) is derived with respect to the parameters. The derivative of the sliding motion constraint with respect to the parameters can be derived in a similar manner. The derivative of equation (6) with respect to equation (1). The derivative is: (twenty one), For the sake of simplicity, dots Displacement at point Represented as And located at that point The normal vector of the structure at that location Represented as Furthermore, for the sake of brevity, Represented as .

[0247] are basis functions about The derivative of .

[0248] Algorithms II and III can be used to speed up equation (21). and The calculation. Moreover, with Same, It is a constant and can be pre-computed, stored, and used in subsequent optimizations. Similarly, the derivative of equation (6) with respect to A is: , And about t The derivative is: .

[0249] Rigid skeletal constraints As shown in equation (7), rigid motion of the bone structure can be applied using the orthogonality of the rotation matrix. In other words, for a structure to become rigid... The spatial derivative of the transformation function (also called the transformation) in equation (1) must be rotated at the point of motion. Equation (7) is derived below. The derivative of equation (1) with respect to the parameters can be derived similarly. about The derivative and constraints with respect to The derivative of . In equation (7) about The derivative of the mn-th element is: , in Dimension and For the same matrix, all elements are zero except for the mn-th element, which has a value of 1. Similarly, Dimension and For the same matrix, only the mn-th element has a value of 1. are basis functions about The derivative of, and similar to Its computation can be facilitated by spatial trees / graphs and can be parallelized per support and per point (see Algorithms II and III).

[0250] The derivative of A with respect to the mn-th element is: .

[0251] Note that the penalty (and / or bone stiffness constraint) is independent of t, so its derivative with respect to t is a zero vector.

[0252] Organizational incompressibility constraint In addition to equation (8), the organizational incompressibility constraint can be defined as Although this form supports symmetric penalties for both compression and expansion, when When approaching zero and when When it degenerates and / or singularizes, it can lead to unstable optimization because its derivative involves... The reverse.

[0253] Similar to equation (7), equation (8) is independent of t. Therefore, the derivative of equation (8) with respect to t is the zero vector.

[0254] Equation (8) regarding The derivative of the mn-th element is: .

[0255] The derivative of equation (8) with respect to the mn-th element of A is: .

[0256] Affinity The affinity or bending energy in equation (9) is used to apply the smoothness of the deformation. It penalizes any local deformation that is not affine, implying the name. As shown in equation (9), the penalty depends only on And about The derivative is: .

[0257] Topology preservation The reverse consistency constraint in equation (10) is used to impose topology preservation. Because The derivative and The derivation of the derivative is similar, so only the former with respect to the parameters is provided. The derivative of the former can be calculated in the same way with respect to the parameters. The derivative of . about The derivative of the mn-th element is: , in and yes The affine and torsion coefficients, and yes control points and The basis functions.

[0258] The derivative of A with respect to the mn-th element is: .

[0259] The derivative of the m-th element with respect to t is: .

[0260] Activity Control Point Scheme The location of control points is important for the accuracy of deformable registration because displacement (e.g., ...) Control points are defined only at their locations. This means that, to accurately address distortion, more control points should be distributed around areas exhibiting misalignment, and control points should not (or are unnecessary) be defined in well-aligned areas. It is well known that at the start of registration, most areas in an image are misaligned, and as registration progresses, some misalignments in certain areas will be improved or mitigated. Leveraging this knowledge, the technique of this invention uses a scheme to select active control points for each optimization iteration and / or each multi-resolution level based on a point-by-point image dissimilarity metric (MSE). Examples of point-by-point dissimilarity metrics, besides MSE, include L between modality-independent neighborhood descriptors (MINDs). 2 Distance [8] (e.g.) Figure 5 (Exemplary depictions in (a) and (c)) and negative values ​​of mean absolute error and / or GI.

[0261] Figure 5 An example of adaptive control point distribution is shown. The first (and / or coarsest) multi-resolution level can involve a relatively dense distribution of control points based on a threshold of image dissimilarity measures (in... Figure 5 In particular, marked as light spots), such as L between modally independent neighborhood descriptors. 2Distance[8]. Figure 5 (a) shows the L overlaid between modality-independent neighborhood descriptors (MINDs) at the coarsest spatial resolution level. 2 The distance to the active control point (as a light spot). Figure 5 (b) shows the active control points (as light spots) overlaid on a fixed image (e.g., the first medical image) at the coarsest spatial resolution level. Figure 5 (c) shows the L overlaid between modality-independent neighborhood descriptors (MINDs) at the next finer spatial resolution level. 2 The number of activity control points has decreased in terms of distance. Figure 5 (d) shows a reduced number of active control points overlaid on a fixed image (e.g., the first medical image) at the next finer spatial resolution level.

[0262] Figure 5 (a) and Figure 5 (b) shows active control points at the same coarsest spatial resolution level, respectively overlaid on the MIND map and the fixed image (and / or the first medical image). Therefore, they show the same dense set of control points.

[0263] Figure 5 (c) and Figure 5 (d) shows the active control points at the next finer spatial resolution level on the MIND map and the fixed image (and / or the first medical image), respectively. Therefore, compared with the MIND map... Figure 5 (a) In comparison, Figure 5 (c) shows a sparser set of active control points, and for coverage on a fixed image, Figure 5 (d) and Figure 5 (b) Similar.

[0264] In addition to this disclosure Figure 5 (a) to Figure 5 Aside from minor image artifacts in (d) artifacts, control points of the same level appear identical on both the MIND map and the fixed image. Figure 5 (a) and Figure 5 (b) covers the same control points at the coarsest level, which is different from the previous level. Figure 5 (c) and Figure 5 (d) The same control points covered at the next finer level.

[0265] After loading the first and second medical images (also referred to as: the image pair to be registered; e.g., according to steps S102 and S104) and before running registration optimization (e.g., according to step S106-O), the algorithm calculates the image anisotropy mapping (e.g., according to step S106-D), with... Figure 5 (a) and Figure 5 The same as those exemplified in (c). Based on simple thresholding and subsequent farthest point sampling or any other more complex point sampling method (e.g., according to steps S106-S), the algorithm selects control points from locations exhibiting high image dissimilarity, such as... Figure 5 The example is marked by a light spot. Figure 5 This shows the mapping overlaid on MIND dissimilarity (in Figure 5 (a) and Figure 5 (c) and fixed images (such as the first medical image; in Figure 5 (b) and Figure 5 (d) The light control points on the image. Initially, (e.g., the first and second) images are mostly misaligned, resulting in... Figure 5 (a) and Figure 5 (b) shows a large number of control points. In later optimization iterations and / or at multi-resolution levels, active control points can be further sampled from the previously densely defined control points based on the same image dissimilarity metric. For example... Figure 5 (c) and Figure 5 As shown in (d), relative to Figure 5 (a) and Figure 5 (b) shows that the number of control points (especially active control points) has been reduced.

[0266] By allowing the set of active control points to be a subset of the previously densely defined control points, all derivative calculations defined in the previous section involve a smaller memory footprint. Let... It is an index A set of indices for active control points. In other words, The following terms in the above derivative calculations can be modified into smaller matrices: •if If it has been calculated in advance, then from arrive Otherwise, Algorithms II and III can be evaluated for a smaller set of control points; •from arrive That is, indexing to the displacement associated with the active control point; •from arrive That is, indexing to the displacement associated with the active control point.

[0267] Adaptive selection of active control points can be enabled for each optimization iteration and / or each multi-resolution level. The former may involve expensive and frequent computations of image anisotropy mapping. Therefore, the latter may be a more attractive option.

[0268] Traditional non-learning iterative methods are slow and cannot handle large motions (and / or large deformations), model sudden changes in motion, or model tissue-related motions. Although deep learning (DL) methods involve only a single forward propagation, their ability to handle deformations has so far been inferior to non-learning methods. Neural field methods have demonstrated better performance than DL, but like DL, their performance is limited by the memory of the computational unit (e.g., CPU and / or GPU memory). Because DL and neural field methods cannot operate at native image resolution (or less coarse resolution), they struggle to handle motions (and / or deformations) of detailed structures compared to non-learning methods.

[0269] Since non-learning methods remain superior in terms of registration accuracy and memory usage, the technique of this invention leverages non-learning and parallelizes its computation. Unlike previous work on parallelization of non-learning deformable registration, the deformation model of this invention is parameterized, capable of capturing varying degrees of deformation and modeling abrupt changes in motion amplitude and direction (see, for example, equation (1) for multiple tightly supported RBFs). The technique of this invention advantageously incorporates an adaptive control point distribution scheme to improve computational efficiency and registration accuracy. The technique of this invention can be additionally equipped with prior knowledge to model sliding motion between organ interfaces, model tissue-related motions (e.g., rigid motion of skeletal structures and / or motion of incompressible tissues), and preserve topology (e.g., through local differential homeomorphisms). According to the technique of this invention, these data items and penalized analytical gradients are used to further accelerate optimization while reducing memory consumption.

[0270] The algorithm's input, parameters, and computation speed are likely characteristics of this invention. Its output displacement field can further help identify whether the algorithm has been used, as it estimates tissue-related deformations and provides clear visualizations of sliding motion, rigid skeletal motion, and almost independent motions of organs within body cavities, such as... Figure 6 As illustrated in the example.

[0271] Figure 6 Exemplary examples illustrate tissue-related deformation. In the upper left, conventional unconstrained methods fail to maintain vertebral rigidity, while the inventive technique in the upper right maintains skeletal rigidity while allowing deformation of surrounding soft tissues. In the lower left, conventional unconstrained methods fail to model the motion of the liver and spleen due to respiratory motion in dual-energy CT imaging, while the inventive technique in the lower right effectively captures the sliding motion of the liver and spleen.

[0272] Where not explicitly described, the various embodiments or aspects and features thereof described in conjunction with the accompanying drawings may be combined or interchanged with each other without limiting or expanding the scope of the described invention, provided that such combination or interchange is meaningful and in the sense of the invention. Advantages described with respect to specific embodiments of the invention or with respect to specific figures are also advantages of other embodiments of the invention wherever applicable.

[0273] By referring to the prior art cited herein: [1] X. Li et al., Neural Graphics Primitives-based Deformable ImageRegistration for On-the-fly Motion Extraction, 2024 [2] X. Gu et al., Implementation and evaluation of various demonsdeformable image registration algorithms on a GPU, Phys. Med. Biol., 2009 [3] K.PunithaA, GPU-Accelerated Deformable Image RegistrationAlgorithm With Applications to Right Ventricular Segmentation, IEEE Access, 2017 [4] S. Ekstrom et al., Faster dense deformable image registration by utilizing both CPU and GPU, Journal of Medical Imaging, 8(1): 2021 [5] P. Muyan-Ozcelik et al., Fast Deformable Registration on the GPU: ACUDA Implementation of Demons, IEEE ICCSA, 2008 [6] H. Mousazadeh et al., GPU implementation of a deformable 3D image registration algorithm, 2011 # finite element method [7] J. P. Pluim et al., Image Registration by Maximization of Combined Mutual Information and Gradient Information, MICCAI 2000, pp. 452–461 [8] M. P. Heinrich et al., Modality independent neighbourhood descriptor for multi-modal deformable registration, Med. Imag. Anal., 16(7):1423-35, 2012

Claims

1. A computer-implemented method (100) for providing deformable registration results from a second medical image to a first medical image, comprising the following method steps: - Receive (S102) a first medical image, the first medical image including at least one anatomical structure; - Receive (S104) a second medical image, wherein, The second medical image includes the at least one anatomical structure; - Perform nested iteration (S106), wherein the nested iteration includes iteration at spatial resolution levels as spatial resolution increases, and wherein the nested iteration further includes iteration at each spatial resolution level with respect to a set of deformation parameters of an analytically deformable registration function used to deformably register the received (S104) second medical image to the received (S102) first medical image, the nested iteration including: for a predetermined spatial resolution level, o Determine (S106-D) the image dissimilarity between the first medical image received (S102) and the second medical image received (S104), wherein the image dissimilarity per unit volume or per unit area is determined in relation to the predetermined spatial resolution level; o Select (S106-S) multiple active control points, wherein the active control points are located in the following regions of the first medical image and / or the second medical image, in which the determined (S106-D) image dissimilarity exceeds or is at least equal to a predefined image dissimilarity threshold, particularly specific to the level; The set of deformation parameters (S106-S) is iteratively optimized by optimizing the objective function, wherein the objective function is a differentiable function of the analytical deformable registration function, and the analytical deformable registration function comprises the sum of analytical functions centered on the selected active control points (S106-S); and - The iterations at the spatial resolution level include: The process proceeds (S106-P) to another predetermined spatial resolution level, which has a finer spatial resolution and / or a smaller unit volume or unit area than the previous predetermined spatial resolution level, and repeats the steps of determining (S106-D) image dissimilarity, selecting (S106-S) multiple active control points, and iteratively optimizing (S106-O) the set of deformation parameters at the other predetermined spatial resolution level; and If, for a predetermined spatial resolution level with the finest spatial resolution and / or a predetermined minimum size for a unit volume or unit area, the set of deformation parameters has been iteratively optimized (S106-O), then the nested iteration (S106-T) is terminated; and - Provide (S108) a deformable registration result that registers the second medical image to the first medical image, wherein the deformable registration result includes the use of the set of deformation parameters in the deformable registration function after termination (S106-T) optimization.

2. The method (100) according to claim 1, wherein, The first medical image and the second medical image are acquired using independently selected medical imaging modalities.

3. The method (100) according to any one of the preceding claims, wherein, The first medical image and the second medical image were acquired using the same medical imaging modality.

4. The method (100) according to any one of the preceding claims, wherein, The step of performing iterative optimization (S106-O) of the set of deformation parameters using parallel computing includes distributing the iterative optimization (S106-O) across parallel computing resources according to the selected (S106-S) active control points.

5. The method (100) according to any one of the preceding claims, wherein, The multiple active control points (S106-S) are adaptively selected in the following manner: - After each iteration of optimization steps or after a predetermined number of iterations, especially at a predetermined spatial resolution level; and / or - Per predetermined spatial resolution level.

6. The method (100) according to any one of the preceding claims, wherein, The selection of the multiple active control points (S106-S) is based on point sampling, optionally... Specifically, performing the iteration at a particularly finer spatial resolution level includes repeated point sampling.

7. The method (100) according to any one of the preceding claims, wherein, The deformable registration function includes: - Spline functions, especially B-spline functions; and / or - The sum of basis functions, especially the sum of radial basis functions (RBFs), and even more particularly the sum of compactly supported basis functions (CSRBFs).

8. The method (100) according to any one of the preceding claims, wherein, This set of deformation parameters includes the affine transformation matrix, offset, a set of torsion coefficients, and / or a set of basis function support widths.

9. The method (100) according to any one of the preceding claims, wherein, The set of deformation parameters and / or the selected multiple control points, especially the selected multiple active control points, are stored in shared memory.

10. The method (100) according to any one of the preceding claims, wherein, The objective function includes: - A measure of point cloud similarity; and / or - A measure of image similarity.

11. The method (100) according to any one of the preceding claims, wherein, Optimizing the objective function includes applying a set of physical constraints, which may include at least one physical constraint, optionally Applying this set of physical constraints includes including corresponding physical constraint terms in the objective function; and / or Applying this set of physical constraints includes optimization using the corresponding Lagrange multipliers.

12. The method (100) according to claim 11, wherein, The set of physical constraints, including at least one of the physical constraints, includes: Affinity or bending energy constraint; Marker point matching; Sliding motion constraints; Rigid skeletal constraints; Organizational incompressibility constraints; and / or Topology preservation constraints.

13. The method (100) according to any one of the preceding claims, wherein, The dissimilarity of an image or point cloud is determined by the following: - Average distance MD; - Mean Squared Error (MSE); - Gradient information; and / or - Mutual information of combinations.

14. The method (100) according to any one of the preceding claims, wherein, Receiving (S102) the first medical image and receiving (S104) the second medical image includes receiving segmentation masks, particularly multiple masks, associated with the first medical image and associated with the second medical image, respectively.

15. The method (100) according to any one of the preceding claims, wherein, Executing (S106) nested iteration further includes the following steps: - Receive (S106-I) an indication of the downsampling factor and / or the initial predetermined spatial resolution level; - Downsample the first medical image and / or the second medical image (S106-DS) to achieve the initial predetermined spatial resolution level; - Perform upsampling on the first medical image and / or the second medical image (S106-US) to achieve the additional predetermined spatial resolution level.

16. Use of the method (100) according to any one of the preceding claims for supporting image-guided interventions, for treatment monitoring, for image comparison and / or for radiotherapy planning.

17. A computing device (200) for providing deformable registration results from a second medical image to a first medical image, comprising: - A first input interface (202) is configured to receive a first medical image, the first medical image including at least one anatomical structure; - A second input interface (204) is configured to receive a second medical image, wherein the second medical image includes the at least one anatomical structure; - A nested iteration module (206) configured to perform nested iterations, wherein the nested iterations include iterations at spatial resolution levels as spatial resolution increases, and wherein the nested iterations further include iterations at each spatial resolution level with respect to a set of deformation parameters of an analytically deformable registration function used to deformably register a received second medical image to a received first medical image, wherein for a predetermined spatial resolution level, the nested iterations include: o Determine the image dissimilarity between the received first medical image and the received second medical image, wherein the image dissimilarity per unit volume or per unit area is determined in relation to the predetermined spatial resolution level; o Select multiple active control points, wherein the active control points are located in the following regions of the first medical image and / or the second medical image, in which the determined image dissimilarity exceeds or is at least equal to a predefined image dissimilarity threshold, particularly specific to the level; The set of deformation parameters is iteratively optimized by optimizing the objective function, wherein the objective function is a differentiable function of the analytical deformable registration function, and the analytical deformable registration function comprises the sum of analytical functions centered on the selected active control points (S106-S); and - The iterations at the spatial resolution level include: The process proceeds to another predetermined spatial resolution level, which has a finer spatial resolution and / or a smaller unit volume or unit area than the previous predetermined spatial resolution level, and repeats the steps of determining image dissimilarity, selecting multiple active control points, and iteratively optimizing the set of deformation parameters at the other predetermined spatial resolution level; and If the set of deformation parameters has been iteratively optimized for a predetermined spatial resolution level with the finest spatial resolution and / or a predetermined minimum size for a unit volume or unit area, then the nested iteration is terminated; and - Output interface (208) configured to provide a deformable registration result of registering the second medical image to the first medical image, wherein the deformable registration result includes the use of the set of deformation parameters in the deformable registration function after the termination optimization.

Images