EP4616221A1

EP4616221A1 - Longitudinal magnetic resonance imaging method

Info

Publication number: EP4616221A1
Application number: EP23809976.6A
Authority: EP
Inventors: Jan Sijbers; Arjan DEN DEKKER
Original assignee: Universiteit Antwerpen; Interuniversitair Microelektronica Centrum vzw IMEC
Current assignee: Universiteit Antwerpen; Interuniversitair Microelektronica Centrum vzw IMEC
Priority date: 2022-11-10
Filing date: 2023-11-10
Publication date: 2025-09-17
Also published as: WO2024100240A1

Abstract

The present disclosure relates to the field of magnetic resonance imaging, specifically addressing a method for processing magnetic resonance image to identifying structural changes between magnetic resonance images acquired at different points of time. Additionally, the disclosure describes a system configured for the same. An aspect of the disclosure relates to a method comprising the steps of: acquiring a first image dataset through MRI scanning of the target volume at a first point in time; acquiring a second image dataset through MRI scanning of the target volume at a second point in time; reconstructing a magnitude image from the first image dataset, thereby obtaining a reference image; transforming the reference image to correspond with the second image dataset, thereby obtaining a transformed image comprising transformed image data, based on a set of deformation parameters; and optimizing a cost function configured to adjust the values of the set of deformation parameters until a stopping condition is satisfied, thereby obtaining a deformation vector field describing the structural changes in the target volume from the first time point to the second time point.

Description

LONGITUDINAL MAGNETIC RESONANCE IMAGING METHOD FIELD OF THE INVENTION The present disclosure relates to the field of magnetic resonance imaging specifically addressing a method for processing magnetic resonance image to identifying structural changes between magnetic resonance images acquired at different points of time. Additionally, the disclosure describes a system configured for the same. BACKGROUND In longitudinal MRI studies, repeated MRI scans are performed in specific clinical scenarios, such as patient follow-up (e.g., tumor monitoring), therapy response assessment, and large-scale longitudinal studies. These studies constitute essential tools for tracking changes in pathology and evaluating treatment efficacy in various diseases. Patients are scanned with the same imaging protocol at different time points, typically every few days, weeks or months. Conventional longitudinal MRI analysis involves multiple steps. First, image data is acquired in k-space (Fourier space) at different time points, and then images are individually reconstructed from the k-space data. Second, the reconstructed images from different time points are geometrically aligned with each other and normalized in terms of intensity. Finally, local structural deformations are detected, often through voxelwise or ROI-wise comparisons. Unfortunately, conventional longitudinal MRI workflows are time-inefficient because a significant amount of redundant information is acquired with each subsequent MRI image during the longitudinal follow-up process. Conventional workflow relies on acquiring fully sampled MRI images, aiming for high-quality images at each time frame, which is considered necessary for image registration and subsequent detection of structural changes. However, acquisition time is directly proportional to the number of samples acquired in k-space, leading to lengthy acquisition times due to oversampling. Various methods are available for identifying the changes between matching images, depending on the specific application and the required level of accuracy. For instance, Nadeem Saad et al. in Medical Physics. vol.47, no.1 (2019) describes a method employing the Iterative Closest Point algorithm iteratively refines the transformation aligning two sets of image points until the minimum distance between corresponding points is achieved. In the context of structural change detection, one set of points corresponds to the points in the first image, and the other set corresponds to the points in the second image that need to be registered with the first image. However, similar to other conventional image processing approaches, this method requires the complete reconstruction of two MRI images, which involves acquiring two full MRI scans at two different points in time. Another example of a conventional image processing approach for determining structural changes is documented in US 2021/166391 A1, which encounters the same limitations as described above. To shorten acquisition time, techniques based on compressed sensing (CS) theory can be considered, exploiting sparsity in some spatial transform domain and reconstructing images from substantially less data (through k-space subsampling). Additionally, since changes between two subsequent images in a longitudinal MRI study are generally small (except for global transformations), temporal redundancy can be further leveraged to shorten image acquisition time. CS-based methods can accelerate longitudinal MRI scanning by extending spatial regularization with temporal regularization. However, these methods often require careful tuning of multiple regularization parameters to balance sparsity constraints (both spatial and temporal) and data consistency, making them less practical. Alternatively, deep learning methods have been explored to automatically detect and utilize the distribution of expected images in both spatial and temporal domains. However, data-driven methods heavily rely on extensive training with a large number of datasets, which can be particularly challenging to obtain for longitudinal MRI applications. Finally, conventional longitudinal MRI protocols also face two additional problems: (1) they carry an inherent risk of error propagation, as the current detection of changes strongly depends on the quality of individually reconstructed images, and (2) current k-space (sub)sampling schemes are optimized for image quality rather than change detection. Accordingly, there is a need to address the limitations of conventional longitudinal MRI protocols by providing a more efficient longitudinal imaging protocol, without risking error propagation throughout successive data acquisitions. SUMMARY OF THE INVENTION As described above, there is a need to remedy the limitations of conventional longitudinal MRI protocols. Longitudinal MRI is an important diagnostic imaging tool for evaluating the effects of treatment and monitoring disease progression. However, MRI, and particularly longitudinal MRI, is known to be time consuming. To accelerate imaging, deep learning methods have been considered to exploit sparsity, both on single image as on image sequence level. State-of-the-art methods, however, are generally focused on image reconstruction, and consider analysis (e.g., alignment, change detection) as a postprocessing step. In the present disclosure, a novel framework is described to estimate longitudinal image changes directly from a reference image and subsequently acquired, preferably subsampled, image data. In contrast to state-of-the-art longitudinal imaging, the herein disclosed method avoids the conventional multistep process of image reconstruction of subsequent images, image alignment, and deformation vector field computation. Instead, the set of follow-up images, along with motion and deformation vector fields that describe their relation to the reference image, are estimated in one go. In this way, the longitudinal MRI protocol can be improved, for instance, by omitting the steps related to image registration and subsequent change detection present in conventional (e.g. multi-step) process. Moreover, the acquisition time for each successive MRI scan at a different point in time can be decreased by (heavily) undersampling the k-space, thereby omitting the acquisition of redundant information, without sacrificing the effectiveness of change detection. A first overview of various aspects of the technology according to the present disclosure is given hereinbelow, after which specific embodiments will be described in more detail. This first overview is meant to aid the reader in understanding the technological concepts more quickly, but it is not meant to identify the most important or essential features thereof, nor is it meant to limit the scope of the present disclosure. An aspect of the present disclosure relates to a (computer-implemented) method for identifying structural changes in a target volume of an object or a subject using a Magnetic Resonance Imaging (MRI) device, comprising the steps of: - acquiring a first image dataset through MRI scanning of the target volume at a first point in time; wherein the first image dataset comprises sufficient data to reconstruct an MRI image representing the target volume at the first time point; - acquiring a second image dataset through MRI scanning of the target volume at a second point in time; wherein the second image dataset comprises sufficient data to capture structural changes in the target volume from the first time point to the second time point; - reconstructing a magnitude image from the first image dataset, thereby obtaining a reference image; - transforming the reference image to correspond with the second image dataset based on a set of deformation parameters describing the deformation from the reference image to the transformed image, and optionally a set of phase parameters describing a possible phase mismatch between the reference image and the transformed image, thereby obtaining a transformed image comprising transformed image data; - optimizing a cost function configured to adjust the values of the deformation parameters, when receiving the set of deformation parameters and the optional phase parameters as input, until a stopping condition is satisfied, thereby obtaining a deformation vector field (DVF) describing the structural changes in the target volume from the first time point to the second time point; wherein the cost function comprises a data consistency term configured to quantify the degree of similarity between the transformed image data and the second image data, and a regularization term configured to impose a constraint on the optimization process based on the similarities between the reference image data and the transformed image data. Another aspect of the present disclosure relates to a (computer-implemented) method for longitudinal Magnetic Resonance Imaging of a target volume of a subject; comprising the steps of - acquiring a first image dataset, preferably a first k-space dataset, by MRI scanning of the target volume at a first point in time; wherein the first image dataset, preferably the first k-space dataset comprises sufficient k-space datapoints to reconstruct a first magnitude image of the whole target volume; - acquiring a second image dataset, preferably a second k-space dataset, by MRI scanning of the target volume at a second point in time; wherein the second image dataset, preferably the second k-space dataset comprises a portion of the k-space datapoints that must be acquired to reconstruct a second magnitude image of the whole target volume; - estimating a deformation vector field (DVF) that describes a geometrical deformation from the first magnitude image to the second magnitude image; wherein the DVF is estimated by reconstructing a reference image corresponding to the first magnitude image, said reference image comprising a magnitude parameter based on the first image dataset, preferably k-space dataset, and mapping said reference image such that the magnitude parameter adheres to the second image dataset, preferably k- space dataset. In some embodiments the data consistency term includes calculating the squared difference between the transformed image data and corresponding part of the second image data, whereby the transformed image data corresponds to the reference image data after applying the set of deformation and optionally the set of phase parameter. In some embodiments the set of phase parameters are determined independently from the set of deformation parameters. In some embodiments the set of phase parameters are determined from a central portion of the first image data and the second image data. In some embodiments wherein the cost function is optimized by applying an optimization algorithm; preferably wherein the optimization algorithm comprises coordinate descent. In some embodiments the regularization term comprises applying a total variation minimization and/or using trained neural networks. In some embodiments the DVF is estimated by minimising a loss function that depends on the magnitude parameter of the reference image and the second image dataset, preferably k-space dataset. In some embodiments the estimated DVF is subjected to regularisation; preferably by applying a total variation minimization and/or using trained neural networks. In some embodiments the second point in time is later than the first point in time; preferably days, months and/or years. In some embodiments the DVF is estimated without reconstructing the second magnitude image based on the second dataset. In some embodiments the method comprises calculating an image warping operator based on the DVF; and generating the second magnitude image by reconstructing the first magnitude image based on the first image dataset, preferably k-space dataset, and applying said image warping operator onto said first magnitude image; preferably by applying the geometrical deformation described by the DVF onto at least a portion of said first magnitude image. In some embodiments the first and/or second image datasets, preferably k-space datasets, further comprise a nuisance parameter, that includes at least one of a phase parameter, an intensity parameter and/or an alignment parameter; and wherein the DVF is estimated along with the corresponding nuisance parameter. In some embodiments the image warping operator comprises a rotation, scaling, and/or translation parameter, preferably determined based on one or more nuisance parameters; and whereby applying said image warping operator onto the first magnitude image rotates at least a portion of said first magnitude image based on the rotation parameter, scales at least a portion of said first magnitude image based on the scaling parameter, and/or translates at least a portion of said first magnitude image based on the translation parameter; preferably before applying the geometrical deformation described by the DVF. In some embodiments the second image dataset, preferably k-space dataset, is acquired by subsampling the MRI scan of the target volume. In some embodiments the first image dataset, preferably k-space dataset comprises sufficient k-space datapoints to reconstruct a first magnitude image with a predetermined image quality parameter, and/or wherein the second image dataset, preferably k-space dataset, comprises insufficient k-space datapoints to reconstruct a second magnitude image with the predetermined image quality parameter; wherein said image quality parameter includes at least one of a resolution, signal-to-noise ratio (SNR) and/or a field of view (FOV). In some embodiments the method comprises identifying at least a portion of the second magnitude image that is similar to the first magnitude image; and reconstructing the first magnitude image based on the first image dataset, preferably k-space dataset, and said similar portion of the second magnitude image. In some embodiments the method comprises acquiring a third image dataset, preferably k-space dataset, by MRI scanning of the target region at a third point in time; wherein the third image dataset, preferably k-space dataset, comprises a portion of the k-space data that must be acquired to reconstruct a third magnitude image of the whole target volume; and determining a second DVF that describes a geometrical deformation from the second magnitude image to the third magnitude image. In some embodiments the second DVF is adapted such that there is a smooth transition from the first DVF to the second DVF. In some embodiments the method comprises the step of identifying one or more deformable neuroanatomy parameters based on the DVF. In some embodiments the target volume corresponds to a least a portion of an object, such as a manufactured or assembled object; wherein the structural change represents a defect or error within the object. In some embodiments the target volume corresponds to a least a portion of an object, such as a utilised or modified object; wherein the structural change represents a degradation or modification to the object. Another aspect of the present disclosure relates to a computer program product for implementing, when executed on a processor, a method of the present disclosure, when provided with data as input, preferably image data, preferably k-space data, acquired by an MRI device. Another aspect of the present disclosure relates to a system comprising an MRI device and a processor, wherein said MRI device is adapted for acquiring MRI data by MRI scanning of a target volume, and wherein said processor is adapted for receiving said MRI data as input and performing the steps of the method of the present disclosure. DESCRIPTION OF THE FIGURES The following description of the figures relate to specific aspects of the present invention that are merely exemplary in nature and are not meant to restrict the scope of the of the present invention, its application or uses. Throughout the drawings, the following symbols and abbreviations are used: first magnitude image (^_^); second magnitude image (^_^); deformation vector field (^); image warping operator based on the deformation vector field ^(^); Deformation (D); Ground Truth Deformation Vector Field (GT DVF); reference Ground Truth Deformation Vector Field (ref GT DVF); zero filled data - inverse Discrete Fourier transform (Z-IDFT); Magnetic Resonance Imaging (MRI); Temporal Compressed Sensing MRI (TCS-MRI); “Delta-MRI” corresponds to a preferred embodiment of a method of the present disclosure – as discussed in Example 1 of the present disclosure. Figure 1 Fshows a flow diagram of the (computer-implemented) method for identifying one or more structural changes in a target volume according to an embodiment of the present disclosure by determining a deformation vector field (^). Figure 2 shows a flow diagram of the (computer-implemented) method for generating a transformed image corresponding to a second image set acquired at time point 2 from a reference image by applying a warping operator ^(^). Figures 3 to 10 are discussed in Example 1 of the present disclosure. Figure 3A shows a cross section of magnitude image ^_^. Figure 3B shows a cross section of magnitude image ^_^ aligned to magnitude image ^_^. Figure 3C shows a cross section of magnitude image ^_^. Figure 4A shows the GT DVF x-component. Figure 4B shows the GT DVF y-component. Figure 4C shows the GT DVF z-component. Figure 5A shows a direct reconstruction of ^_^ using 1% subsampled k-space data. Figure 5B shows a reconstruction of ^_^ obtained by TCS-MRI based on 1% subsampled k-space data. Figure 5C shows a reconstruction of ^_^ obtained by Delta-MRI based on 1% subsampled k-space data. Figure 6A shows a direct reconstruction of ^_^ using 5% subsampled k-space data. Figure 6B shows a reconstruction of ^_^ obtained by TCS-MRI based on 5% subsampled k-space data. Figure 6C shows a reconstruction of ^_^ obtained by Delta-MRI based on 5% subsampled k-space data. Figure 7A shows the ref GT DVF x-component. Figure 7B shows the ref GT DVF y-component. Figure 7C shows the ref GT DVF z-component. Figure 8A shows the DVF x-component estimated by Delta-MRI based on 5% subsampled k-space data. Figure 8B shows the DVF y-component estimated by Delta-MRI based on 5% subsampled k-space data. Figure 8C shows the DVF z-component estimated by Delta-MRI based on 5% subsampled k-space data. FIG.9A shows a sampling scheme of a 1% normally distributed subsampling percentage. FIG.9B shows a sampling scheme of a 5% normally distributed subsampling percentage. FIG. 10 shows the normalized reconstruction error as a function of the subsampling percentage for Z- IDFT ①, TCS-MRI ② and Delta-MRI ③. DETAILED DESCRIPTION In the following detailed description, the technology underlying the present invention will be presented in various technological aspects. It should be noted that these aspects, as outlined herein and illustrated in the figures, can be arranged, substituted, combined, and designed in numerous configurations. All such configurations are explicitly considered and form a part of this disclosure. This description aims to facilitate the reader's comprehension of the technological concepts but is not intended to limit the scope of the present invention, which is solely defined by the claims. The present invention aims to improve the efficiency of a longitudinal Magnetic Resonance Imaging (MRI) workflow. Longitudinal MRI is an important diagnostic imaging tool for evaluating structural changes, for example physiological or anatomical changes, in specific clinical or nonclinical scenarios by MRI scanning an object or a subject with the same imaging protocol at different points in time, for example, every several days, weeks or months. As used herein, the subject may refer to a human subject, for example, a healthy individual or a patient suffering from a specific pathology. Alternatively, the subject can be nonhuman, for example, an animal subject. As used herein, the object may refer to a non-living object, for example, a manufactured object or an object assembled from a plurality of parts in a wide range of industrial applications, for example, for defect detection or metrology. Imaging through Magnetic Resonance (MR) techniques is well-known and widely applied in medical imaging. In essence, a typical MRI technique acquires imaging data by manipulating the magnetic spins in a selected target volume of the subject under examination using an MRI device and processing the measured responses from the magnetic spins. An MRI device may include hardware to generate different magnetic fields for imaging, such as a static magnetic field along the z-direction to polarize the magnetic spins, gradient fields along the x, y, or z directions to select a portion of the target volume for imaging, and a radiofrequency magnetic field to manipulate the spins. Various advanced MRI techniques are known in the field, for example, to improve sampling efficiency, but the present disclosure is not limited to a particular method. The MR signal is an induced current generated by the precession of the net magnetization after stimulation by the radiofrequency pulse. An MRI detects the MR signal from two channels that measure the same precessing magnetization from two different orthogonal perspectives. The MR signal can be represented as a vector with real and imaginary components recorded by the corresponding real and imaginary channels, from which the magnitude and phase can be derived. The measured MRI signal is recorded in a k-space data matrix containing the raw MRI data in Fourier Space. Creating a single MR image commonly involves collecting a series of data frames, referred to as "acquisitions." In each acquisition, an RF excitation produces new transverse magnetization, which is then sampled along a specific trajectory in k-space. Due to various physical and physiological constraints, most MRI imaging methods use a sequence of acquisitions, with each one sampling part of k-space. The data from this sequence of acquisitions are then used to reconstruct an image. Therefore, as used here, "MRI data" refers to the measurement data acquired from one or more MRI scans by an MRI device, which can be stored as k-space data in the frequency domain and/or as imaging data in the imaging domain. Accordingly, a "parameter" with reference to the MRI data, such as a "magnitude parameter" or "phase parameter," refers to a subset of the MRI data associated with the referenced parameter, from the k-space data and/or image data, depending on the context. Consequently, a complete data representation in k-space or in image space requires displaying either the real and imaginary components or an equivalent display of magnitude and/or phase of the complex- valued signal. However, in most clinical scans, the phase is discarded, and only the magnitude image is utilized for diagnosis. For example, a "magnitude image" can be used to display background tissue with spin-density-like contrast. The reconstructed image may consist of a plurality of image elements, such as pixels or voxels, each characterized by both a magnitude and a phase component. Therefore, a "component" concerning reconstructed image data, like a "magnitude component" or "phase component," signifies a subset of the image elements associated with the referenced component within the image data, depending on the context. The number of k-space samples directly influences spatial resolution. Therefore, a higher sampling of the MRI scan acquires more k-space data, allowing for the reconstruction of a higher quality image. Image quality may be defined based on relevant parameters, such as resolution, signal-to-noise ratio (SNR), and field of view (FOV). Tomographic techniques are known for reconstructing an MRI image based on the recorded MRI data. As previously described, conventional longitudinal MRI analysis typically follows a 'multi-step' approach. Initially, data is acquired in k-space (Fourier space) through MRI scans of a target volume at various time points, and individual images are then reconstructed from this k-space data. Subsequently, the reconstructed images from different time points are aligned geometrically and normalized in terms of intensity. Finally, local structural deformations are identified. Traditional longitudinal MRI analysis prioritizes high-quality imaging of the target volume at each time frame, as this is deemed essential for image registration and subsequent change detection. Consequently, this approach relies on acquiring fully sampled k-space data to reconstruct a high-quality MRI image of the entire target volume. However, acquiring fully sampled k-space data is time-consuming, as acquisition time is directly proportional to the number of samples obtained. This leads to longer acquisition times due to oversampling, as a substantial amount of redundant information is acquired with each subsequent MRI image. Typically, only a portion of the entire target volume will exhibit anatomical deformation, contributing to longer acquisition times in conventional longitudinal imaging. The innovation presented in this technology is based on the recognition that, in longitudinal imaging, the primary focus is on the anatomical and structural changes in the target volume of a subject over time. In contrast to conventional longitudinal acquisition schemes that emphasize high-quality imaging, this technology aims to optimize the quantification of these changes within a given acquisition time. When change detection and analysis are the main objectives, there's no need to target high-quality images at every time frame. In this disclosure, a novel framework is described that estimates longitudinal image changes directly from, preferably subsampled, k-space datasets acquired at different time points. This framework is designed to emphasize change quantification by estimating deformable neuroanatomy parameters of interest from the longitudinal image data, making the most of temporal redundancy in successive scans concerning a reference image. This approach improves the efficiency of a longitudinal MRI protocol by eliminating the steps related to image registration and subsequent change detection found in conventional 'multi-step' approaches. Furthermore, the acquisition time for each successive MRI scan at a different time point can be reduced by heavily subsampling the MRI scanning area, thus avoiding the acquisition of redundant information without compromising the effectiveness of change detection. This technology can be considered 'general-purpose' (medical) imaging and can be adapted for imaging longitudinal changes in various targets. The 'target volume' of a subject can refer to any part of the subject's body for which MRI data can be obtained through MRI scanning, such as one or more organs or anatomical features. It is important to note that the disclosed technology is not limited to a specific anatomy or pathology, provided that MRI data can be acquired for the target volume or its parts. In principle, this technology can be directly applied to non-medical imaging of any object where longitudinal changes are expected between successive images. For example, quality inspection of manufactured or assembled. While this technology has broad applicability, we provide specific examples related to tracking changes in local tissue growth of a subject, which is a common MRI application, especially in scenarios like tumor monitoring. However, it can also find utility in various other medical applications, such as post-surgical patient follow-up, therapy response assessment, large-scale longitudinal studies, and more. In essence, this technology allows for a novel approach in which changes between MRI data acquired by successive MRI scans of a target volume at different time points can be directly estimated. Specifically, it is possible to estimate a set of follow-up images along with motion and deformation vector fields (DVFs) describing the relationship between the reference image and the follow-up images. In this innovative approach, change detection is performed by directly estimating the DVF from the originally acquired (raw) image data, which helps avoid redundancy. An overview of various aspects of the technology of the present disclosure is given hereinbelow, after which specific embodiments will be described in more detail. This overview is meant to aid the reader in understanding the technological concepts more quickly, but it is not meant to identify the most important or essential features thereof, nor is it meant to limit the scope of the present disclosure. When describing specific embodiments, reference is made to the accompanying drawings, which are provided solely to aid in the understanding of the described embodiment. An aspect of the present disclosure relates to a (computer implemented) method for identifying one or more structural changes in a target volume of an object or a subject using a Magnetic Resonance Imaging (MRI) system, comprising the steps of: - acquiring a first image dataset through MRI scanning of the target volume at a first time point; wherein the first image dataset comprises sufficient image data for a reconstruction of an MRI image representing the target volume at the first time point; - acquiring a second image dataset through MRI scanning of the target volume at a second time point; wherein the second image dataset comprises sufficient image data to capture the structural changes in the target volume from the first time point to the second time point; - reconstructing an MRI image from the first image dataset, preferably a magnitude image, thereby obtaining a reference image consisting of reconstructed image data; - transforming the reference image to correspond with the image data of the second image dataset, preferably transforming the reference magnitude image to correspond with the image data of the second image data, thereby obtaining a transformed image consisting of transformed image data; - defining a set of deformation parameters describing the structural changes from the reference image to the transformed image, preferably from the reference magnitude image to the transformed magnitude image, and optionally a set of phase parameters describing a possible phase mismatch between the reference image and the transformed image; and, - optimizing a cost function configured to adjust the values of the deformation parameters, when receiving the set of deformation parameters and the optional set of phase parameters as input, until a stopping condition is satisfied, thereby obtaining a deformation vector field (DVF) describing the structural changes in the target volume between the first and the second time points; wherein the cost function comprises a data consistency term configured to quantify the degree of similarity between the transformed image data of the transformed image and the second image data of the second time point, and a regularization term configured to impose constraints on the optimization process based on the similarities between the reference image and the transformed image, preferably the reference magnitude image and the transformed magnitude image. As described above, the present invention differs from a conventional image processing approach in that the DVF can be determined ‘directly’ from (a limited amount of) imaging data by leveraging the temporal redundancy existing between an MRI scan of an object or subject at a first point in time (serving as a reference image) and an MRI scan of the same object or subject at a subsequent point in time. Temporal redundancy implies that the imaging data in a sequence of successive images exhibits structural similarities. For example, successive medical images of the same anatomical region for a single patient typically show organs of largely similar size and relative positions. Alternatively, successive nonmedical images of the same region for a single object or a plurality of similar objects typically show features of largely similar size and relative positions. The determination of the DVF can be achieved by transforming the reference image (reconstructed from the ‘full’ MRI scan acquired at a first point in time) to correspond to the image data of a second MRI scan (acquired at a second point in time). In this way a transformed image is generated that represents a ‘prediction’ of how an image reconstructed from the image data of the second MRI scan is expected to look without performing this reconstruction. The transformed image is generated by defining a set of deformation parameters describing the transformation from the reference image (corresponding to the first image data) to the transformed image (representing an estimation of the second image data), and optionally a set of phase parameters describing a possible phase mismatch between these images based on the datasets. The transformed image consists of transformed image data that is estimated or predicted based on the deformation between the reference magnitude image and the transformed magnitude image. This prediction is made using the defined set of deformation parameters. In other words, it represents what the image data at the second time point would look like if the deformation between the two time points is accurately modelled based on the parameters and deformation information. In embodiments, the transformation from the reference image to the transformed image can be performed using image data from a ‘partial’ MRI scan of the same object at a subsequent point in time. Specifically, it can be sufficient to capture enough image data to describe the structural changes in the target region from the first time point to the second time point, but it does not require the same amount of image data as needed for a ‘complete’ reconstruction of an MRI image. This means that the scanning speed during the successive scan can be significantly increased, for example, through subsampling the MRI scan. Further, the values of these sets of parameters can then be optimised by implementing a cost function configured to receive these parameters as input. In an embodiment the optimisation can be performed using an optimisation algorithm, optionally iteratively, until the values are obtained that directly align the first image dataset with the second image dataset, thereby obtaining the DVF without the need to reconstruct a second magnitude image from the second image dataset. To ensure proper optimisation, the cost function can include a data consistency term configured to quantify the similarity between the ‘predicted data’ of the transformed image representing the second time point (combination of the image data acquired through MRI scanning at the first time point and the set of deformation parameters) and the ‘actual data’ corresponding to the second time point (the image data obtained through MRI scanning at the second time point), and a regularization term configured to introduce a form of regularization to the optimization process that leverages the similarity between the respective magnitude images from the first and second time points. The data consistency term is included to ensure that the estimated data more closely matches the actual measured data such that an accurate and effective image registration can be realized. In an embodiment, the data consistency term quantifies how well the estimated data (obtained using the current set of parameters) matches the actual measured data. It can measure the degree of similarity or consistency between the two datasets and assign a value based on a predefined matching accuracy. Since the data consistency term is a component of the cost function, minimizing this term indicates that the estimated data fits the measured data as closely as possible. This encourages the optimization process to adjust the parameters in a way that makes the estimated data more similar to the actual data. In an embodiment, the data consistency term can also serve as an indicator of convergence during the optimization process. As the optimization iteratively adjusts the parameters, the data consistency term should ideally decrease. When it reaches a minimum or a sufficiently low value, it indicates that the estimated data closely matches the measured data, and the optimization process is converging. The optimization process can be stopped once a stopping condition is satisfied. The stopping condition may be advantageously implemented to force the algorithm to terminate the refinement process once a desired outcome has been produced. The stopping condition can be any meaningful condition, such as number of iterations, quality of solutions, statistical values, and so on. In an embodiment, the regularization term can introduce a form of regularization or constraint to the optimization process to promote smoothness or stability in the estimated deformation field. For example, the regularization term can discourage abrupt or erratic changes in the deformation field. Smoothness in the DVF is relevant because real-world deformations, such as those occurring in biological tissues, tend to be continuous and smooth. Alternatively or in combination, the regularization term can improve the robustness of the deformation model by reducing sensitivity to noise or outliers in the data. It can help prevent overfitting, which occurs when the deformation field excessively adjusts to data noise rather than capturing meaningful changes. In embodiment, the impact of the regularization term can be controlled by a weighting parameter that determines its influence on the optimization process. A higher weight can place more emphasis on the regularization term, which can lead to smoother deformations, while a lower weight allows more flexibility in the deformation field. The choice of the regularization weight helps strike a balance between fitting the data well and achieving a smooth deformation. The optimization of the cost function involves finding the set of parameters that minimizes the cost function. This process can be carried out using numerical optimization techniques. An exemplary approach may include the following steps: - Initialization: Initial values are assigned to the parameters provided as input which represent an initial ‘estimate’. These set of parameters can include the set of deformation parameters and the optional phase parameters. These parameters are the parameters that will be adjusted while running the cost function. - Evaluation: The values of the cost function can be evaluated based on the (initial) values of the provided parameter. The evaluation provides a measure of how well these (initial) values of the parameters currently align or fit the data. - Optimization: An optimization algorithm can guide the optimization process by providing information about the direction and rate of change of the cost function, as well as imposing constraints on possible changes to the values of the parameters. The choice of optimization method and the way it is configured can impact the efficiency and effectiveness of the cost function optimization. For example, the optimization algorithm can include coordinate descent method. The goal of the optimization is to find the set of parameters that align the estimated and actual datasets, thereby providing the optimal fit while minimizing the cost function. - Parameter update: The values of the provided parameters can be adapted based on the optimization algorithm. The goal is to adjust the parameters in a way that reduces the cost function at the next evaluation step. - Termination: The optimization process can continue until a stopping condition is met. This criterion can be a specified number of iterations, reaching a predefined level of accuracy, or other criteria depending on the specific optimization method. In an embodiment, the optimization of the cost function can be an iterative process, specifically, by iteratively adapting the values of the provided parameters until termination. With each iteration the cost function can be advantageously re-evaluated to determine if the cost function is lowered, and the optimization is proceeding in the right direction. The optimization algorithm aims to find the minimum of the cost function.A preferred aspect of the present disclosure relates to a (computer implemented) method for longitudinal Magnetic Resonance Imaging (MRI) of a target volume of or an object or a subject, comprising the steps of: - acquiring a first image dataset, preferably comprising k-space datapoints, by MRI scanning of the target volume at a first point in time; wherein the first image dataset comprises sufficient k-space datapoints to at least reconstruct a first magnitude image of the whole target volume; - acquiring a second image, preferably comprising k-space datapoints, by MRI scanning of the target volume at a second point in time; wherein the second image dataset comprises a portion of the k-space datapoints that must be acquired to at least reconstruct a second magnitude image of the whole target volume; and - estimating a deformation vector field (DVF) that describes a geometrical deformation from the first magnitude image to the second magnitude image. Another aspect of the present disclosure relates to a computer program product for implementing, when executed on a processor, a method in accordance with any embodiment described herein when provided with image data as input, preferably from a medical imaging device. Another aspect of the present disclosure relates to a system comprising an MRI imaging device and a processor wherein said imaging device is adapted for acquiring a plurality of MRI images, and wherein said processor is adapted for receiving said MRI images as MRI image data and performing a method in accordance with any of the herein described embodiments when provided with a plurality of MRI image data. Another aspect of the present disclosure relates to a system comprising a medical imaging device and a processor wherein said medical imaging device is adapted for acquiring a plurality of medical images, and wherein said processor is adapted for receiving said medical images as medical image data and performing a method in accordance with any of the herein described embodiments when provided with a plurality of medical image data. An aspect of the present invention will be discussed in more detail with reference to Figure 1, which schematically illustrates a (computer implemented) method for longitudinal MRI imaging. The illustrated method can generally be performed by or under the control of a processor of a computing unit, for instance, the computing unit of a system described in the present disclosure. The method can be partially or fully automated according to some embodiments thereof. To aid the reader in the understanding of the present disclosure, the method will be explained as being performed by a processor of a medical imaging system. The skilled person will, however, appreciate that the method may be performed by any other processor if appropriately configured, for instance, of a personal computer. At step 11, the method may comprise the acquiring of a first MRI dataset, for example comprising a plurality of k-space datapoints, through MRI scanning of a target volume of an object or subject. The MRI data may, for instance, be acquired by an MRI device. In any of the embodiments described herein, the image data may be stored in a memory of a computing system, remotely or locally (for example, a memory of a database, a server, or any other memory). For example, the processor may download the image data from the memory. The image data is preferably stored in a format that is suitable for processing by the processor or can be converted thereto, for example, for the reconstruction of an image. At step 12, the method may comprise the acquiring of at least a second MRI dataset, comprising a plurality of k-space datapoints, through MRI scanning of the (same) target volume of the (same) object or subject at a different point in time. For ease of explanation, the method will be described for two image datasets only, specifically a first image dataset acquired at time 1 and a second image dataset acquired at time 2. Nevertheless, it should be noted that embodiments of the present method can be expanded to include further datasets, for example, a third image dataset acquired at time 3, a fourth image dataset acquired at time 4, and so on, as will be discussed later. Advantageously, the one or more further points in time are later than the first point in time, for example days, months and/or years after the first point in time. Alternatively or in combination, the method may be performed on datasets acquired at an earlier point in time to retrospectively improve the image quality, for example, days, months and/or years before the first point in time. In some embodiments the first image dataset comprises sufficient image data to reconstruct a magnitude image of the whole target volume with a predetermined image quality parameter. Advantageously the first image dataset comprises sufficient k-space datapoints to reconstruct a complex image of the whole target volume, as will be described further below. As such, the first image dataset is preferably acquired by sufficiently sampling an MRI scan of the target region at a first point in time. In this way, the first image dataset may form the basis for the reconstruction of one or more (high quality) reference images, for example, of a portion of the target volume that does not have an (anatomical) deformation. In some embodiments the second image dataset comprises a portion of the image data that must be acquired to reconstruct a magnitude image of the whole target volume with a predetermined image quality parameter. Alternatively, the second image dataset comprises insufficient k-space datapoints to reconstruct a magnitude image with the predetermined image quality parameter. As such, the second image dataset is preferably acquired by subsampling the MRI scan of the target volume. In this way, the acquisition of redundant information can be avoided, speeding up the scanning process of successive MRI scan. As used herein, “subsampling” refers to a method of acquiring a subset of the image data required to reconstruct a complex image of the whole target volume. Typically, subsampling is performed by speeding up scan time at the expense of lower spatial resolution. For example, subsampling the k-space during scanning time can be performed using Cartesian trajectories, such as rectilinear sampling. While Cartesian trajectories are the most common in clinical settings, some other trajectories can be used, for example, sampling along radial lines, and/or sampling along spiral trajectories. Hence, it should be understood that the present technology method is not limited to a specific subsampling technique, as it can be readily adjusted for different subsampling patterns. At step 13, , a reference image may be reconstructed based the first image dataset, that corresponds to a magnitude image of the first image data. Accordingly, the reference image may consist of a plurality of image elements, such as pixels or voxels, each image element being characterized by a magnitude and/or a phase component. The magnitude component based on the first image dataset, more specifically, based on the amplitude of the sampled MRI signal, is distinct from the phase component. Additionally, the reference image may comprise one or more nuisance parameters that can be, at least to an extent, unspecified but must be accounted for mapping the magnitude component which is of interest. For example, the nuisance parameters may include at least one of a phase parameter, an intensity parameter data and/or an alignment data parameter. In any of the above-described embodiments, a (high) quality reference image can be reconstructed on sufficiently sampled MRI data to improve the accuracy of the following image processing steps. At step 14, using the reference image, the goal is then to find the deformation vector field (DVF), hereinunder also expressed as “^”, by mapping the reference image (corresponding to the image data acquired at time point 1 – in the image domain) such that it adheres to the image data acquired at time point 2, more specifically by defining a set of deformation parameters describing the deformation from the reference image to the transformed image, optionally along with possible nuisance parameters from the acquired datasets. As used herein, ‘mapping’ indicates that the reference image is adapted to correspond with an image dataset other than the image dataset that was used to reconstruct the reference image. The objective of the longitudinal imaging is, therefore, to generate an output that can be, for instance, displayed to a user of the method or the system, for example, a healthcare professionals. In an embodiment the output may comprise an image of the deformation vector field. Alternatively or in combination, the output may comprise an image that is generated by applying the deformation vector field onto a reference image, as discussed later. Advantageously, the output is provided on a dedicated user interface to improve an interpretation by a user. Hereinunder various embodiments are described for estimating the DVF that can be implemented as part of the above-described aspects of the herein disclosed method. However, the skilled person understands that various (image) parameter mapping techniques may be contemplated, for instance, that improve the mapping accuracy and/or computational efficiency, and although examples are provided, the present technology is not limited to the herein discussed embodiments. Moreover, when applied onto medical image data specifically, the DVF may thus represent a physiological change and/or a pathophysiological change in a target volume. However, as previously discussed, the present method can also be implemented for various other (for example, non-medical) applications for imaging changes in a target volume (for example, of an object) between successive MRI scans. In some embodiments, the DVF can be estimated by reconstructing a reference image corresponding to the magnitude image of the first image dataset, and mapping it such that it adheres to the second image dataset, thereby obtaining a transformed magnitude image representing the second image dataset. Preferably, said reference image comprises at least a magnitude component based on the first image dataset, and the DVF is estimated by mapping said reference image such that at least the magnitude parameter adheres to the second image dataset. The transformed magnitude image may consist of ‘estimated’ image data representing an estimation or prediction of the second image dataset – this can at an initial stage be different from the ‘actual’ second image data acquired through MRI scanning. Advantageously, the DVF is estimated by mapping the magnitude parameter along with the corresponding one or more nuisance parameters. In some embodiments, the DVF can be estimated by formulating a loss function that receives the set of deformation parameters as input. In this way, the DVF can be estimated by minimising a loss function that depends on the magnitude component from the reference image, which corresponds to the magnitude image of the first image dataset that has been transformed to adhere to the image data of the second dataset. Advantageously, the DVF is estimated by minimising a loss function that depends on the magnitude parameter along with the corresponding one or more nuisance parameters. In an embodiment, the optimization of the cost function can be performed by implementing an optimization algorithm known in the art. The choice of optimization algorithm can impact the efficiency and effectiveness of the cost function optimization depending on the desired outcome. For example, the choice of optimization algorithm can be dependent on the optimization of the calculation speed or accuracy, or a combination of speed and accuracy. In a preferred embodiment the optimization algorithm may include the Coordinate Descent that is a variable cost function configured to (iteratively) optimize one variable (coordinate) at a time while holding the others fixed. It is particularly useful when dealing with high-dimensional problems or when the cost function is not easily differentiable. Depending on the complexity, this can involve solving sub-problems using other optimization techniques. In any of the above embodiments, the estimated DVF is subjected to regularisation to improve the model quality. Preferably, the regularisation can be performed by a applying a total variation minimization. Alternatively or in combination, the regularisation can be performed using trained neural networks as known in the art, for example, by feeding the data as input to a (trained) convolutional 2D/3D neural network. Within the context of medical imaging, the regularization term can be used to encourage the model to generate images that are visually consistent with the previous scans. This can help in denoising images, improving image registration (alignment of images), and tracking changes in a subject’s condition over time. In some embodiments the method comprises acquiring a further image dataset, for example third or fourth, by MRI scanning of the target region at a further point in time; wherein the further image dataset comprises a portion of the image data that must be acquired to reconstruct a further magnitude image of the whole target volume; and determining a further DVF that describes a geometrical deformation from the former magnitude image to the following magnitude image, for example, the second magnitude image to the third magnitude image, the third magnitude image to the fourth magnitude image, and so on. Advantageously, the DVF is adapted such that there is a smooth transition from a former, preferably preceding, DVF to the following, preferably succeeding, DVF. As mentioned earlier, as part of a further implementation of the herein disclosed technology, the method may comprise the step of calculating an image warping operator based on the DVF, hereinunder expressed as “^(^)”. This aspect is shown, for example, in Figure 2. Specifically, in an embodiment the magnitude image corresponding for the second, preferably subsampled, image dataset can be generated by reconstructing the first magnitude image based on the first image dataset and applying the image warping operator ^(^) onto said first magnitude image, such that the geometrical deformation described by the DVF is applied onto at least a portion of said first magnitude image. In this way, the magnitude image of the second image dataset may be represented as a warped version of the magnitude image of the first image dataset. For example, with further reference to Figure 2, at step 21, data representing a complex image can be acquired at time point 1: ^_^ + ^_^^ = ^_^^^{^^^}; wherein ^_^ corresponds to a magnitude image of the data acquired at time point 1. Although not shown, subsampled Fourier data is acquired at time point 2: ^_^ = ^^^^_^^^{^^^}^. Since ^_^, which is the magnitude image corresponding to the data acquired at time point 2, is assumed to be a deformation of ^_^, the corresponding complex image at time point 2 can be determined as follows: ^_^ + ^_^^ = ^_^^^{^^^} . The DVF (^) represents the geometrical deformation from ^_^ to ^_^. Hence, at step 22, by applying a warping operator ^(^), which is based on the DVF (^), onto ^_^, an estimation of the corresponding complex image ^_^ can be generated, which represents the data of time point 2. In some embodiment, the image warping operator may comprise a rotation parameter such that said image warping operator onto the first magnitude image rotates at least a portion of said first magnitude image based on said rotation parameter. In some embodiment, the image warping operator may comprise a scaling parameter such that said image warping operator onto the first magnitude image scales at least a portion of said first magnitude image based on said scaling parameter. In some embodiment, the image warping operator may comprise a translation parameter such that said image warping operator onto the first magnitude image translates at least a portion of said first magnitude image based on said translation parameter. In some embodiments, the selected views are projection views and the transformation includes transforming the projection views to k-space. In some embodiments the imaging system is a magnetic resonance imaging system and each view samples a line in k- space. In some embodiments, the imaging system is an x-ray CT system and each view is a projection in Radon space. In any of the above embodiments, applying the image warping operator onto the first magnitude image may comprise rotating at least a portion of said first magnitude image based on the rotation parameter, scaling at least a portion of said first magnitude image based on the scaling parameter, and/or translating at least a portion of said first magnitude image based on the translation parameter. Preferably, the at least rotation parameter, a scaling parameter and/or a translation parameter are applied before the geometrical deformation described by the DVF. In some embodiments the method may comprise identifying at least a portion of the second magnitude image that is similar to the first magnitude image; and reconstructing the first magnitude image based on the first image dataset and said similar portion of the second magnitude image. In this way, the similarities between the successive images can be exploited to improve the image quality of the image reconstruction, for example, the SNR. Alternatively or in combination, this could even be performed retrospectively to improve the quality of the first magnitude image. Hereinunder a number of embodiments is described for estimating the DVF based on the complexity of the relevant transformations. This section is meant to aid the reader in understanding the implementation of the herein disclosed method more easily, but it is not meant to identify the most important or essential features thereof, nor is it meant to limit the scope of the present disclosure. In any of the below embodiment, let ^_^ = be two complex valued, preferably 3D, MRI images with ^ voxels in the image domain, of which the elements can be expressed in polar form the magnitude and phase of ^_^ , respectively. Furthermore, let ^_^, ^_^ ∈ ^^{^} be the k-space (Fourier) representations of ^_^, ^_^ respectively: ^_^ = ^^_^ , with ^ ∈ ℂ^{^×^} the discrete Fourier transform operator. Finally, let ^ ∈ ℝ^{^×^} be the (unknown) deformation vector field (DVF), comprised of one, preferably 3D, vector per voxel, describing the geometrical deformation from ^_^ to ^_^, and let ^( ∙, ^) denote the image warping operator along the DVF (^) that (approximately) yields ^_^ when applied to ^_^: ^⁽^_^, ^⁾ = ^_^. In this way, the DVF (^) can enables the construction of a forward model that transforms the magnitude image ^_^ and the phase image ^_^ into the k-space representation ^_^. In some embodiments, the DVF ^ can be estimated directly from k-space measurements of ^_^and ^_^. ^ These measurements, herein denoted by ^^^{^} _^^^^^ ∈ , may be noise disturbed and can be modelled as: ^^{^} _^ = ^_^^^_^ + ^_^ , with ^ = 1, 2, where ^_^ ∈ ℂ^{^^} is an additive noise contribution, modeled as a zero-mean, complex-valued Gaussian random variable, and ^_^ ∈ {0,1}^{^^ × ^} with ^_^ ≤ is a subsampling operator that selects the k-space points acquired for the ^-th image. In an embodiment wherein the k-space data is fully sampled, denoted by ^^{^} _^, the operator ^_^ corresponds with the identity matrix ^ ∈ ℝ^{^×^}. In some embodiments, it may be that sufficient data ^^{^} _^ is acquired to allow a good quality reconstruction ^̂_^of the reference magnitude image ^_^. If fully sampled data is indeed available, such an estimate can be calculated as ^^{^^} the inverse discrete Fourier transform and | ∙ | the point-wise modulus operator. Furthermore, it can be assumed that the phase ^_^ is slowly varying and that the low the low frequencies in k-space are sufficiently sampled to obtain a good estimate of ^_^ by calculating ^_^ = ∠^^^{^^} ^^{^} _^^, where ^^_^ is zero-filled at the non-sampled indices and ∠⁽∙⁾ denotes the phase operator. Under these assumptions, ^ can be directly estimated by minimizing a regularized least-squares functional without requiring the reconstruction of ^_^: ^ = arg min + ^ with λ^(^) a regularization term with weight λ that imposes prior knowledge about the DVF ^. In the above embodiment, the regularization term ^(^) may be chosen equal to ‖∇ ^‖_^, with ∇ the gradient operator. This embodiment reflects the assumption that ^_^ only differs from ^_^ by a local deformation. Hence, ^ can be expected to be mostly zero, except in the region of change, where it is assumed to be smooth. This decreases the amount of data points needed for its estimation, compared to a direct reconstruction of ^_^ from the data ^_^. In another embodiment, where ^_^ also differs from ^_^ by a rigid transformation, on top of a local deformation, this transformation can be included in the DVF ^, but then ^ will no longer be mostly zero. Alternatively, the rigid transformation can be isolated from the DVF, and represented by rotation and translation parameters ^, ^ ∈ ℝ^{^} ^ , as follows: ^^, ^^{^}, ^̂^ = arg min + ^, ^, ^ where ^( ∙, ^, ^, ^) is an image warping operator that rotates the image according to ^, translates the image according to ^, and then applies the local deformation described by ^. In any of the above embodiments, the optimization problem can be solved using an (iterative) computational optimization technique. Preferably, a coordinate descent method can be applied, for example, cyclic block coordinate descent where the first block consists of the parameters ^ and ^ and the second block consist of ^. Reference throughout this specification to “one embodiment” or “an embodiment” means that a particular feature, structure or characteristic described in connection with the embodiment is included in at least one embodiment of the present disclosure. Thus, appearances of the phrases “in one embodiment” or “in an embodiment” in various places throughout this specification are not necessarily all referring to the same embodiment. As used herein, the terms “comprising”, “comprises” and “comprised of” as used herein are synonymous with “including”, “includes” or “containing”, “contains”, and are inclusive or open-ended and do not exclude additional, non-recited members, elements or method steps. The terms “comprising”, “comprises” and “comprised of” when referring to recited members, elements or method steps also include embodiments which “consist of” said recited members, elements or method steps. The singular forms “a”, “an”, and “the” include both singular and plural referents unless the context clearly dictates otherwise. As used herein, the term “substantially” refers to the complete or nearly complete extent or degree of an action, characteristic, property, state, structure, item, or result. For example, an object that is “substantially” enclosed would mean that the object is either completely enclosed or nearly completely enclosed. The exact allowable degree of deviation from absolute completeness may in some cases depend on the specific context. However, generally speaking the nearness of completion will be so as to have the same overall result as if absolute and total completion were obtained. The use of “substantially” is equally applicable when used in a negative connotation to refer to the complete or near complete lack of an action, characteristic, property, state, structure, item, or result. As used herein, the term “about” is used to provide flexibility to a numerical value or range endpoint by providing that a given value may be “a little above” or “a little below” said value or endpoint, depending on the specific context. Unless otherwise stated, use of the term “about” in accordance with a specific number or numerical range should also be understood to provide support for such numerical terms or range without the term “about”. For example, the recitation of “about 30” should be construed as not only providing support for values a little above and a little below 30, but also for the actual numerical value of 30 as well. The recitation of numerical ranges by endpoints includes all numbers and fractions subsumed within the respective ranges, as well as the recited endpoints. Furthermore, the terms first, second, third and the like in the description and in the claims, are used for distinguishing between similar elements and not necessarily for describing a sequential or chronological order, unless specified. It is to be understood that the terms so used are interchangeable under appropriate circumstances and that the embodiments of the disclosure described herein are capable of operation in other sequences than described or illustrated herein. Reference in this specification may be made to devices, structures, systems, or methods that provide “improved” performance (e.g. increased or decreased results, depending on the context). It is to be understood that unless otherwise stated, such “improvement” is a measure of a benefit obtained based on a comparison to devices, structures, systems or methods in the prior art. Furthermore, it is to be understood that the degree of improved performance may vary between disclosed embodiments and that no equality or consistency in the amount, degree, or realization of improved performance is to be assumed as universally applicable. In addition, embodiments of the present disclosure may include hardware, software, and electronic components or modules that, for purposes of discussion, may be illustrated and described as if the majority of the components were implemented solely in hardware. However, one of ordinary skill in the art, and based on a reading of this detailed description, would recognize that, in at least one embodiment, the electronic based aspects of the present disclosure may be implemented in software (e.g., instructions stored on non-transitory computer-readable medium) executable by one or more processing units, such as a microprocessor and/or application specific integrated circuits. As such, it should be noted that a plurality of hardware and software-based devices, as well as a plurality of different structural components may be utilized to implement the technology of the present disclosure. For example, “servers” and “computing devices” described in the specification can include one or more processing units, one or more computer-readable medium modules, one or more input/output interfaces, and various connections connecting the components. EXAMPLES An exemplary implementation of the technology according to the present disclosure is discussed hereinbelow. The provision of examples is meant to aid the reader in understanding the technological concepts more easily, but it is not meant to identify the most important or essential features thereof, nor is it meant to limit the scope of the present disclosure. Example 1 The herein disclosed technology was validated based on an experimentally acquired dataset and benchmarked against a longitudinal compressed sensing-based method of the art. In what follows, an implementation of disclosed technology as a preferred embodiment thereof will be referred to as “Delta- MRI”. The method selected for benchmarking is Temporal Compressed Sensing MRI (TCS-MRI), which is a technique known in the art that allows accurate reconstruction from sparsely sampled image data. Experimental data: High-resolution fully-sampled image data was acquired (using a Cartesian sampling scheme) in a 12 weeks old male C57BL/6 wild-type mouse on a 9.4T Biospec 94/20 USR horizontal MR system (Bruker Biospin MRI, Ettlingen, Germany), with a mouse head 2×2 array cryo-coil (Bruker Biospin MRI, Ettlingen, Germany), using a 3D T2-weighted Turbo-RARE sequence with FOV: 20×15×10 mm³, image acquisition matrix: 256×192×128, 0.078 mm isotropic voxels, TE/TR: 67.2/1800 ms, RARE factor: 12. The fully sampled multi-coil image data was then averaged across the coils to obtain ^^{^} _^. Simulation set-up: Simulation experiments were setup to validate the Delta-MRI method, using the following stepwise procedure. First, a complex valued image ^_^ = ^^^{^^} ^^{^} _^^ was reconstructed, with corresponding magnitude ^_^ and ^_^. From ^_^, a predefined rigid transform (^, ^) and a predefined DVF ^, a magnitude image ^_^ was constructed using an inverse image warp, such that ^( ^_^, ^, ^, ^) = ^_^. FIG.3A, FIG.3B, and FIG.3C show cross sections of ^_^, its fully aligned version ^( ^_^, ^, ^, ^) and ^_^, respectively. Additionally, the area with the most notable deformation (D) is indicated with a dashed arrow in FIG.3C – this is added purely to guide the reader towards the region of interest but is not part of the magnitude image ^_^. The ground truth parameter values of ^ and ^ are given by (2.9°, 4.0°, 5.7°) and (-6, -5, -4.5) pixels, respectively. FIG.4 shows cross sections of the Ground Truth (GT) DVF ^, more specifically, FIG.4A, FIG.4B, and FIG.4C show the GT DVF “x”-component, the GT DVF “y”-component, and the GT DVF “z”- component, respectively. In what follows, the initial estimates of ^, ^ and ^ will be denoted as ^_^, ^_^ and ^_^, respectively. However, as will be described below, Delta-MRI will be benchmarked against TCS-MRI, which does depend on ^_^ and assumes it to be similar to ^_^. Therefore, a complex valued image ^_^ with magnitude ^_^ was constructed, by setting = ^_^. This complex image ^_^ was polluted with complex, Gaussian distributed white noise with standard deviation equal to 4% of the average foreground value of ^_^. This noise disturbed image is denoted by ^_^, and its magnitude and phase by ^̂_^ and ^^_^, respectively. Next, a subsampling operator ^_^ was applied that under-samples ^^_^ by selecting a given percentage of the k-space data points of the image data. In particular, a small cubic central area was complemented by k-space data points that were randomly drawn from a central Gaussian distribution with a standard deviation that corresponds with a quarter of the maximum frequency. FIG.9 shows illustrative examples of the used sampling schemes, more specifically, FIG.9A and FIG.9B show sampling schemes with normally distributed subsampling percentages equal to 1% and 5%, respectively. Finally, of ^, ^ and ^ were estimated from the under-sampled k-space data ^^{^} _^ using the estimator ^^, ^^{^}, ^̂^ = arg The above simulation experiments were repeated for different subsampling rates up to 20%. Implementation details: The optimization task of the above described estimator was carried out using a cyclic Block Coordinate Descent method approach. In the first block, ^ and ^ were optimized with the Limited-memory Broyden–Fletcher–Goldfarb–Shanno method (L-BFGS-B) of SciPy 1.9.3 (Python), with default stopping criteria and bounds ^ ∈ [-20,20]³, ^ ∈ [-0.3 rad, 0.3 rad]³. In the second block, ^ was optimized with a Barzilai-Borwein gradient method with 2000 iterations to ensure convergence. Because the deformation is local, the initial guess ^_^ = 0 was close enough to give an accurate estimation of ^ and ^. Therefore, each block was optimized only once. It was verified that adding more cycles did not change the results significantly. Benchmarking: The simulation experiments were benchmarked against a Temporal Compressed Sensing MRI method, which reconstructs the complex image ^_^ in the temporal and wavelet domain, by solving ^ ^ ^_^ = arg min ^^ d_^ − ^_^^^_^ ^_^ + λ_^ ^‖diag⁽^⁾(^_^ − ^_^)^‖ _^ + λ_^ ^‖Ψ x_^ ^‖ _^^, where Ψ is the Daubechies- x_^ 4 wavelet transform and ^ is a weight vector that controls the demand for similarity between ^_^ and ^_^, enforcing sparsity only in regions where ^_^ and ^_^ are similar. In the present experiments, ^ was fixed to the support of the ground truth DVF. It should be noted that TCS-MRI assumes that the rigid transformation aligning ^_^ and ^_^ is known. Delta- MRI doesn’t make this assumption, and instead estimates the rigid transformation, along with ^, from the available, subsampled data. In each experiment, the rigid transformation estimated by the experiments based on Delta-MRI was used to align the data before using TCS-MRI. Furthermore, the regularization λ_^ and λ_^ were set to λ_^·5.85·10⁻⁹ and λ_^·3.63·10⁻¹ ^ weights ⁰, respectively, where λ_^ = ^ d^_^ ^_^ maintains the regularization strength when varying the sampling scheme. These values were obtained by empirical tuning to optimize reconstruction quality at 1% normal sampling. To evaluate how much information was present in d^{^} _^, and how much information was obtained from ^_^, the results of TCS-MRI and Delta-MRI were compared against a reconstruction of ^_^ obtained from d^_^ alone by applying an inverse Discrete Fourier transform to the zero filled data d^{^} _^, denoted as Z-IDFT. Reconstructions: FIG.5 shows the reconstruction results obtained from Gaussian subsampled k-space data with subsampling percentages equal to 1%, more specifically, FIG.5A shows (a cross-section of) the IDFT reconstruction of ^_^, while FIG.5B and FIG.5C show the reconstructions of ^_^ obtained by TCS-MRI and Delta-MRI, respectively. Based on a visual inspection, it is clear that, despite subsampling, the image generated by Delta-MRI in FIG.5C approximates the fully sampled k-space data shown in FIG.3C. Similarly, FIG.6 shows the reconstruction results obtained from Gaussian subsampled k-space data with subsampling percentages equal to 5%, more specifically, FIG.6A shows (a cross-section of) the IDFT reconstruction of ^_^, while FIG.6B and FIG.6C show the reconstructions of ^_^ obtained by TCS-MRI and Delta-MRI, respectively. Based on a visual inspection, it is clear that, despite heavy subsampling, the image generated by Delta-MRI in FIG.6C approximates the fully sampled k-space data shown in FIG.3C. FIG.7 shows cross sections of the reference Ground Truth Deformation Vector Field (ref GT DVF) that are estimated from perfectly aligned and noiseless ^_^ to ^_^, more specifically, FIG.7A, FIG.7B, and FIG.7C show the ref GT DVF “x”-component, the ref GT DVF “y”-component, and the ref GT DVF “z”-component, respectively. The ref GT DVF has the same effect on the image as the GT DVF. Additionally, included as an illustrative example, FIG.8 shows the Deformation Vector Field (DVF) estimated by Delta-MRI using 5% subsampled k-space data, more specifically, FIG.8A, FIG.8B, and FIG.8C show estimates of the DVF “x”-component, the DVF “y”-component, and the DVF “z”-component, respectively. Results: As a quantitative performance criterion, the normalized reconstruction error of the reconstructed image ^̂_^ was calculated as: ^(^̂_^) = ‖^̂_^ − ^_^‖_^ / ‖^̂_^ − ^_^‖_^, where ^̂_^ = ^^ ^̂_^, ^^{^}, ^̂, ^^ for Delta-MRI. FIG.10 shows the normalized reconstruction error as a function of the subsampling percentage for IDFT ① TCS-MRI ② and Delta-MRI ③. It is shown that the reconstruction error for Delta-MRI is significantly lower than IDFT and the benchmarking method TCS-MRI for subsampling percentages up to 20%. Conclusion: Experiments demonstrate that at the applied sampling rates (1-20%), Delta-MRI significantly outperforms the selected benchmarking method of the art in terms of the normalized reconstruction error, as observed visually in FIGs.5-6 and confirmed quantitatively in FIG.10.

Claims

CLAIMS 1. Computer-implemented method for identifying structural changes in a target volume of an object or a subject using a Magnetic Resonance Imaging (MRI) device, comprising the steps of: - acquiring a first image dataset through MRI scanning of the target volume at a first point in time; wherein the first image dataset comprises sufficient data to reconstruct an MRI image representing the target volume at the first time point; - acquiring a second image dataset through MRI scanning of the target volume at a second point in time; wherein the second image dataset comprises sufficient data to capture structural changes in the target volume from the first time point to the second time point; - reconstructing a magnitude image from the first image dataset, thereby obtaining a reference image; - transforming the reference image to correspond with the second image dataset, thereby obtaining a transformed image comprising transformed image data, based on a set of deformation parameters describing the deformation from the reference image to the transformed image, and optionally a set of phase parameters describing a possible phase mismatch between the reference image and the transformed image; - optimizing a cost function configured to adjust the values of the deformation parameters, when receiving the set of deformation parameters and the optional phase parameters as input, until a stopping condition is satisfied, thereby obtaining a deformation vector field (DVF) describing the structural changes in the target volume from the first time point to the second time point; wherein the cost function comprises a data consistency term configured to quantify the degree of similarity between the transformed image data and the second image data, and a regularization term configured to impose a constraint on the optimization process based on the similarities between the reference image data and the transformed image data.

2. The method according to any one of the preceding claims, wherein the data consistency term includes calculating the squared difference between the transformed image data and corresponding portion of the second image data.

3. The method according to any one of the preceding claims, wherein the set of phase parameters are determined independently from the set of deformation parameters.

4. The method according to any one of the preceding claims, wherein the set of phase parameters are determined from a central portion of the first image data and the second image data.

5. The method according to any one of the preceding claims, wherein the cost function is optimized by applying an optimization algorithm; preferably by applying a coordinate descent method.

6 The method according to any one of the preceding claims, wherein the regularization term comprises applying a total variation minimization and/or a trained neural network.

7. The method according to any one of the preceding claims, further comprising determining a set of a nuisance parameters, that includes at least one of a phase parameter, an intensity parameter and/or an alignment parameter; wherein the DVF is determined along with the nuisance parameters.

8. The method according to any one of the preceding claims, wherein the method further comprises the steps of generating an image warping operator based on the DVF; and generating at least a portion of a second magnitude image representing the second image dataset by applying the image warping operator onto at least a corresponding portion of the reference image.

9. The method according to claim 8, wherein the image warping operator comprises at least one of a rotation parameter, a scaling parameter, and/or a translation parameter; whereby applying the image warping operator onto the reference image rotates at least a corresponding portion of the reference image based on the rotation parameter, scales at least a corresponding portion of the reference image based on the scaling parameter, and/or translates at least a corresponding portion of the reference image based on the translation parameter.

10. The method according to any one of the preceding claims, wherein the MRI scanning comprises acquiring a plurality of views of the target volume, wherein each view comprises a sampling of a line in k- space.

11. The method according to claim 10, wherein the second image dataset is acquired through a subsampling of the k-space.

12. The method according to any one of the preceding claims, wherein the second image dataset lacks sufficient image data to reconstruct a second magnitude image that meets a predetermined minimum threshold of an image quality parameter; wherein the image quality parameter includes a at least one of a resolution, a signal-to-noise ratio (SNR), and/or a field of view (FOV).

13. The method according to any one of the preceding claims, wherein the method further comprises the steps of acquiring a third image dataset through MRI scanning of the target volume at a third point in time; wherein the third image dataset comprises sufficient data to capture structural changes in the target volume from the second time point to the third time point; and determining a second DVF describing the structural changes in the target volume from the second time point to the third time point.

14. The method according to claim 13, wherein the method further comprises the step of transforming the transformed image to correspond with the third image dataset, thereby obtaining a second transformed image comprising second transformed image data, based on a set of deformation parameters describing the deformation from the transformed image to the second transformed image, and optionally a set of phase parameters describing a possible phase mismatch between the transformed image and the second transformed image.

15. The method according to any one of the preceding claims, wherein the cost function is adapted to improve the transition from the first DVF to the second DVF.

16. The method according to any one of the preceding claims, wherein the second point in time is later than the first point in time; preferably days, months and/or years later.

17. The method according to any one of the preceding claims, wherein the second point in time is earlier than the first point in time; preferably days, months and/or years earlier.

18. The method according to any one of the preceding claims, wherein the DVF is determined without performing a step of reconstructing a second magnitude image based on the second image dataset.

19. The method according to any one of the preceding claims, wherein the target volume corresponds to a least a portion of a body of the subject; wherein the structural change represents an anatomical or pathological change in the body of the subject.

20. The method according to any one of the preceding claims, wherein the target volume corresponds to a least a portion of an object; wherein the structural change represents a defect or assembly error within the object, and/or a degradation or modification to the object.

21. Computer program product for implementing, when executed on a processor, a method in accordance with any one of the preceding claims when provided with data as input, preferably wherein the data comprises image data acquired by an MRI device.

22. System comprising an MRI device and a processor, wherein the MRI device is configured for acquiring MRI data by performing an MRI scan, and wherein the processor is configured for receiving the MRI data as input and performing the steps of the method in accordance with any one of the preceding claims.

23. System for identifying structural changes in a target volume of an object or a subject, comprising an MRI device and a processor, wherein the MRI device is configured for acquiring imaging data through MRI scanning, and wherein the processor is configured to receive the imaging data from the MRI device and perform the steps of: - receiving a first image dataset acquired through MRI scanning of the target volume at a first point in time; wherein the first image dataset comprises sufficient data to reconstruct an MRI image representing the target volume at the first time point; - receiving a second image dataset acquired through MRI scanning of the target volume at a second point in time; wherein the second image dataset comprises sufficient data to capture structural changes in the target volume from the first time point to the second time point; - reconstructing a magnitude image from the first image dataset, thereby obtaining a reference image; - transforming the reference image to correspond with the second image dataset based on a set of deformation parameters describing the deformation from the reference image to the transformed image, and optionally a set of phase parameters describing a possible phase mismatch between the reference image and the transformed image, thereby obtaining a transformed image comprising transformed image data; - optimizing a cost function configured to adjust the values of the deformation parameters, when receiving the set of deformation parameters and the optional phase parameters as input, until a stopping condition is satisfied, thereby obtaining a deformation vector field (DVF) describing the structural changes in the target volume from the first time point to the second time point; wherein the cost function comprises a data consistency term configured to quantify the degree of similarity between the transformed image data and the second image data, and a regularization term configured to impose a constraint on the optimization process based on the similarities between the reference image data and the transformed image data.

24. Use of a system according to claim 23 for identifying structural changes in a target volume of a subject, wherein the structural change represents an anatomical or pathological change in a body of the subject.

25. Use of a system according to claim 23 for identifying structural changes in a target volume of an object, wherein the structural change represents a defect or error within the object, and/or a degradation or modification to the object.