Techniques for de-noising labeled medical imagery

Abstract

Techniques for de-noising label-differenced medical images includes operating a medical scanner to collect a control image of a region of interest in a subject. The medical scanner is also operated to collect a labeled image of the region of interest. Multiple patches within the control image are selected. Each patch has a similarity better than a threshold compared to a co-located patch in the labeled image. A control low-rank variability model is derived for the patches within the control image to produce a de-noised control image. A labeled low-rank variability model is derived for the patches within the labeled image. A de-noised labeled image in the region of interest is produced based on the labeled low-rank variability model. A difference image in the region of interest is output based on a difference between the de-noised labeled image and the de-noised control image in the region of interest.

Description

TECHNIQUES FOR DE-NOISING LABELED MEDICAL IMAGERY CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims benefit of Provisional Appln. 63 / 534,206, filed August 23, 2023, the entire contents of which are hereby incorporated by reference as if fully set forth herein, under 35 U.S.C. §119(e). STATEMENT OF GOVERNMENTAL INTEREST

[0002] This invention was made with government support under Grant Numbers EB031080 and AG060054 awarded by the National Institutes of Health. The government has certain rights in the invention. BACKGROUND

[0003] Cross-sectional imaging is an imaging technique which produces a large series of two- dimensional (2D) images of a subject, e.g., a human subject. Examples of cross-sectional imaging techniques include computerized tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET), SPECT scanning, ultrasonography (US), among others. A set of cross-sectional images for a single patient, e.g., for different axially located cross-sections or for the same cross section at different times can be considered three-dimensional (3D) image data, and even four-dimensional (4D) image data for combinations of axial and temporal cross sectional images.

[0004] Various analytical approaches can be applied to the cross-sectional images to detect and highlight portions of the patient's anatomy of interest. For example, the cross-sectional images can be processed by labeling which causes certain materials to register in the scanning device with higher amplitude than in an unlabeled image. Labeling includes the use of x-ray dense materials, radioactive emitters, proton spin precession excitation, among others. In these different image modalities, the labeled image, the unlabeled image, and the difference between the labeled and unlabeled images can all provide the desired information. However, in some circumstances, these images including the labeled, unlabeled, and difference images can have low signal to noise ratio (SNR) that can interfere with performing the desired analysis to obtain a desired result.SUMMARY

[0005] Techniques are provided for de-noising labeled, unlabeled and label-differenced medical images. In a first set of embodiments, a method includes operating a medical scanner to collect a control image of a region of interest in a subject. The method also includes operating the medical scanner to collect a labeled image of the region of interest in the subject. The method further includes selecting, automatically on a processor, multiple patches within the control image. Each patch has a similarity metric value better than a threshold when compared to a patch in the labeled image at a same location in the region of interest. Each patch of the plurality of patches contains fewer than half a number of voxels in the region of interest in the control image. The method still further includes deriving, automatically on the processor, a control low-rank variability model for the multiple patches within the control image; and producing a de-noised control image in the region of interest based on the control low-rank variability model. The method even further includes deriving, automatically on the processor, a labeled low-rank variability model for the multiple patches within the labeled image; and producing a de-noised labeled image in the region of interest based on the labeled low-rank variability model. Yet further still, the method includes outputting, automatically from the processor to a display device, a difference image in the region of interest based on a difference between the de-noised labeled image in the region of interest and the de-noised control image in the region of interest.

[0006] In some embodiments of the first set, each patch has a first dimension size in a range from 4 voxels to 32 voxels and a perpendicular second dimension size in a range from 4 voxels to 32 voxels.

[0007] In some embodiments of the first set, each patch has a square shape.

[0008] In some embodiments of the first set, the similarity metric is a low vector distance between vectors representing each patch. In some of these embodiments the similarity metric is a low Euclidean distance and the threshold is selected so that a number of patches is at least 8.

[0009] In some embodiments of the first set, the low-rank variability model for each of the control image and the label image is based on a subset of eigenvectors having the highest eigenvalues for a low-rank matrix having a low-rank first dimension comprising the voxels in each patch and a low-rank second dimension comprising the plurality of patches. In some of these embodiments, the eigenvalues are based on a non-convex regularization of the low-rank matrix. In some of these embodiments, the non-convex regularization of the low-rank matrix isbased on a log-determinant of a square covariance matrix for the low-rank matrix. Thus, the log- determinant of covariance is used to replace low-rank regularization.

[0010] In some embodiments of the first set, the medical scanner is a Magnetic Resonance Imager (MRI) scanner, the labeled image is an arterial spin labeling (ASL) image, and the difference image is a perfusion image. In some of these embodiments, the region of interest transects a brain of the subject and the perfusion image indicates cerebral blood flow (CBF).

[0011] In other sets of embodiments, a non-transient computer-readable medium or an apparatus or a system is configured to perform one or more steps of the above methods.

[0012] Still other aspects, features, and advantages are readily apparent from the following detailed description, simply by illustrating a number of particular embodiments and implementations, including the best mode contemplated for carrying out the invention. Other embodiments are also capable of other and different features and advantages, and its several details can be modified in various obvious respects, all without departing from the spirit and scope of the invention. Accordingly, the drawings and description are to be regarded as illustrative in nature, and not as restrictive. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Embodiments are illustrated by way of example, and not by way of limitation, in the figures of the accompanying drawings in which like reference numerals refer to similar elements and in which:

[0014] FIG. 1A is a block diagram that illustrates scan elements in a 2D scan, such as one scanned image from a CT scanner;

[0015] FIG. 1B is a block diagram that illustrates scan elements in a 3D scan, such as stacked multiple scanned images from a CT imager or true 3D scan elements from volumetric CT imagers or ultrasound;

[0016] FIG. 2A is a block diagram that illustrates an example of a labeled scan difference system, according to an embodiment;

[0017] FIG. 2B is a block diagram that illustrates examples of low-rank data structures for the labeled scan difference system, according to an embodiment;

[0018] FIG. 3 is a flow diagram that illustrates an example of a method to determine a de-noised difference between a labeled medical scan and an un-labeled (control) medical scan of the same region of interest inside a subject, according to an embodiment;

[0019] FIG.4 is a plot that illustrates an example of the upper convex nature of the log function, used according to the example embodiment;

[0020] FIG. 5A through FIG.5C are labeled / control (L / C) difference images without de noising and with two prior art de-noising methods, according to an example image modality;

[0021] FIG. 6 is a block diagram showing the application of the method of FIG. 3 to MRI ASL perfusion, according to an example embodiment;

[0022] FIG. 7A through FIG. 7D are labeled / control (L / C) difference images as in FIG. 5A through FIG. 5C plus an example of an image with de-noising using the method of FIG. 6, according to the example embodiment;

[0023] FIG. 8A and FIG. 8B are images that illustrate examples of L / C difference images using, respectively, without and with collaborative data selection, respectively, according to example embodiments;

[0024] FIG. 8C is a plot that illustrates the improvement in image quality (reduction of noise), according to the example embodiment;

[0025] FIG. 9 is a block diagram that illustrates a computer system upon which an embodiment of the invention may be implemented; and

[0026] FIG. 10 illustrates a chip set upon which an embodiment of the invention may be implemented. DETAILED DESCRIPTION

[0027] A method and apparatus are described for de-noising label-differenced medical images. In the following description, for the purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the present invention. It will be apparent, however, to one skilled in the art that the present invention may be practiced without these specific details. In other instances, well-known structures and devices are shown in block diagram form in order to avoid unnecessarily obscuring the present invention.

[0028] Some example embodiments of the invention are described below in the context of arterial spin labeling (ASL) in magnetic resonance imaging (MRI) scanning for the brain as a region of interest in a human subject. However, the invention is not limited to this context, and can be used with any pair of labeled and control (L / C) images

[0029] For example, in some embodiments, the de-noising method also incorporates inter-slice dependency of ASL images. In still other embodiments, the de-noising method is configured tode-noise other imaging modalities and signal modulations having similar subjects in which multiple images or multiple slices are acquired at different time points or different receiving coils, such as MRI, fMRI. These modalities variously include diffusion MRI (e.g., receiving image pairs acquired at different diffusion encoding directions or different encoding gradient strengths), functional MRI (e.g., receiving image pairs acquired at different time points), cardiac MRI (e.g., receiving image pairs acquired at different time points), chemical exchange saturation transfer MRI (e.g., receiving multiple MR images that have been acquired at different off- resonance frequencies), MR fingerprinting (e.g., receiving multiple images acquired at different time), multiple TE imaging, inversion recovery imaging for T1 and T2 estimation, and contrast enhanced MRI, alone or in combination. In still other embodiment, the de-noising method can be configured to denoise other imaging modalities such as positron emission tomography (PET) data, computing tomography (CT) data, and ultrasound data, alone or in combination with each other and or in combination with the above MRI modalities.

[0030] In various embodiments, the de-noising method is implemented using a single processor or using parallel processor computing techniques, or some combination. For example, each slice of an image (e.g., in a signal) can be handled separately in parallel.

[0031] Thus, according to several embodiments, a method is provided for de-noising images to improve image quality based on a single pair of labeled / control (L / C images). This is in contrast, for example, with the mean perfusion image obtained from multiple L / C pairs and the output of several other state-of-art denoising methods. 1. Structures

[0032] FIG. 1A is a block diagram that illustrates scan elements in a 2D scan 110, such as one scanned image from a CT scanner. The two dimensions of the scan 110 are represented by the x direction arrow 102 and the y direction arrow 104. The scan 110 consists of a two-dimensional array of 2D scan elements (also called picture elements and abbreviated as pixels) 112 each with an associated position. Typically, a 2D scan element position is given by a row number in the x direction and a column number in the y direction of a rectangular array of scan elements. A value at each scan element position represents a measured or computed intensity or amplitude that represents a physical property (e.g., X-ray absorption, or resonance frequency of an MRI scanner) at a corresponding position in at least a portion of the spatial arrangement 132a, 132b of the living body. The measured property is called amplitude hereinafter and is treated as a scalarquantity. In some embodiments, two or more properties are measured together at a pixel location and multiple amplitudes are obtained that can be collected into a vector quantity, such as spectral amplitudes in MRSI. Although a particular number and arrangement of equal sized circular scan elements 112 are shown for purposes of illustration, in other embodiments, more elements in the same or different arrangement with the same or different sizes and shapes are included in a 2D scan.

[0033] FIG. 1B is a block diagram that illustrates scan elements in a 3D scan 120, such as stacked multiple scanned images from a CT imager or true 3D scan elements from volumetric CT imagers or MRI or US. The three dimensions of the scan are represented by the x direction arrow 102, the y direction arrow 104, and the z direction arrow 106. The scan 120 consists of a three-dimensional array of 3D scan elements (also called volume elements and abbreviated as voxels) 122 each with an associated position. Typically, a 3D scan element position is given by a row number in the x direction, column number in the y direction and a scanned image number (also called a scan number) in the z (axial) direction of a cubic array of scan elements or a temporal sequence of scanned slices. A value at each scan element position represents a measured or computed intensity that represents a physical property (e.g., X-ray absorption for a CT scanner, or resonance frequency of an MRI scanner) at a corresponding position in at least a portion of the spatial arrangement 132a, 132b of the living body. Although a particular number and arrangement of equal sized spherical scan elements 122 are shown for purposes of illustration, in other embodiments, more elements in the same or different arrangement with the same or different sizes and shapes are included in a 3D scan.

[0034] The term voxels is used herein to represent either 2D scan elements (pixels) or 3D scan elements (voxels), or 4D scan elements, or some combination, depending on the context.

[0035] FIG. 2A is a block diagram that illustrates an example of a labeled scan difference system 200, according to an embodiment. Although subject 290, such as a human patient, is depicted for purposes of illustration, subject 290 is not a part of system 200. The system includes an imaging scanner 210, such as a CT, MRO, PET, SPECT or ultrasound (US) scanner, or some combination, configured to produce one or more images representing properties inside human or nonhuman, animate or inanimate, biological, geological, artificial, mechanical or natural subject 290. Often, the scanner is configured to include in the image captured one or more regions of interest inside the subject, which correspond to one or more regions of interest voxels in thecaptured image. In some embodiments, the scanner is operated by a scanner controller module in a computer system 220, such as computer system 900 or chipset 1000 described below and depicted in FIG.9 and FIG. 10, respectively.

[0036] Thus, in some embodiments, the medical scanner is a MRI, labeled and unlabeled images only differ by time. In other words, they can be the repeated scans of the same time of labeled or unlabeled images. Or they can be images acquired with different parameters such as the diffusion weighting direction and or strength in diffusion tensor imaging. Or they can be images under different dynamic physiological or pharmacological modulation (for example, images at different timepoints after contrast agent injection, or images at different time when human subjects are under physical or functional operation or performing). In some embodiments, the medical scanner can be an ultrasound scanner or x-ray or CT or PET scanner and the labeled images and unlabeled images can be images acquired at different timepoints with or without using different parameters with or without contrast injection, with or without physiological or pharmacological modulations.

[0037] The system also includes computer system 220 such as computer system 900 or chipset 1000 described below and depicted in FIG. 9 and FIG. 10, respectively. The computer system includes label difference detection module 230 configured to collect one or more pairs of images for each area of interest, one image of each pair called the control image, and the other image of the pair called a labeled image or simply a label image and to output a label difference output image 229. The output may be used to measure the condition of a subject, recommend treatment for the subject, or control subsequent operation of scanner 210, or some combination.

[0038] The label image refers to an image of the same region of interest but using an imaging modality that differs in some respect, such as at different sampling times, from the image modality used to collect the control image. For example, as described above, the modalities include MRI, fMRI. diffusion MRI (, functional MRI, cardiac MRI, chemical exchange saturation transfer MRI, MR fingerprinting, multiple TE imaging, inversion recovery imaging for T1 and T2 estimation, contrast enhanced MRI, PET, SPECT, or CT, alone or in some combination.

[0039] Noise levels in L / C images are much lower than the in the corresponding difference image 229. For independent and identically distributed (i.i.d.) noise, one assumes each sample has the same probability distribution while white noise samples could follow differentprobability distributions. For example, if the variance of the i.i.d. random noise in the L image is σ 2 L and that in the C image is σC2, then the noise variance of the difference C-L should be σL2+σC2, making the modeling in the difference image more challenging.

[0040] In the illustrated embodiment, the label difference detection module 230 includes a low- rank regularization module 234 used to substantially reduce noise and improve usefulness of the label difference output image 229, a process termed de-noising hereinafter, as described in more detail below. The process exploits the underlying redundancy in the image pairs to facilitate extracting highly correlated local contents in patch groups (subsets of image voxels), which is more robust against noise and better preserve actual structures compared to previous denoising methods that employ voxel or slice level regularization. To support low-rank regularization module 234, the label difference detection module 230 includes one or more low-rank data structures 232.

[0041] FIG. 2B is a block diagram that illustrates examples of low-rank data structures 250 for the labeled scan difference system, according to an embodiment. The low rank data structures 250 include a control image data structure 251a configured to hold the voxel amplitudes and positions for a control image in the region of interest 292 of subject 290 collected from scanner 210 at a first time with a first modality. The low rank data structures 250 include a label image data structure 251b configured to hold the voxel amplitudes and positions for a label image in the region of interest 292 of subject 290 collected from scanner 210 at a second, optionally different, time with a second, optionally different, modality. On each image the location and size of various patches 252 are depicted.

[0042] A patch is a spatially defined subset of contiguous voxels for which amplitude variations are especially similar across the pair of control and labeled (C / L) images. Typically a captured image includes about 500 to 1000 rows and about 500 to 1000 columns of voxels, for a total voxel count between about one quarter of a million to one million voxels. A patch is much smaller than the image so that many patches can be selected within an image to attain statistical significance. A patch is also advantageously large enough to capture the underlying structures to be maintained in the de-noising process. For example, in some embodiments each patch has a first dimension size in a range from 4 voxels to 32 voxels and a perpendicular second dimension size in a range from 4 voxels to 32 voxels. Thus, each patch includes a number of voxels in a range from about 16 to about 1000. In some embodiments, the number of voxels in eachdimension is the same and the patch shape is square. In some embodiments, a patch is curved or circular in shape.

[0043] Similarity can be defined as high covariance or low vector difference of voxel amplitudes between the control and labeled images. The two-dimensional (2D) patches can be represented by voxel amplitudes arranged in a 2D matrix or, alternatively, as arranged as a one- dimensional (1D) vector in which all voxels are listed in sequence, e.g., the first row amplitudes to which the second row of amplitudes is concatenated, to which the third row is concatenated, and so on. Vector distance can be measured using any of several metrics, including, the largest absolute value difference between vector elements (l0), the sum of the absolute values of the differences (l1), the square root of the sum of the squares of the differences (l2, also known as the Euclidean distance), among others using higher powers of the differences between vector elements, or other combinations of differences in vector elements.

[0044] The similar patches are found by sliding a window the size and shape of the patch and measuring the similarity (e.g., based on 2D or 1D covariance, or vector distance). When the similarity of the window between the control and label images exceeds a threshold, the voxels of that window are saved as a patch in the patch information data structure 254. In some embodiments, the patch information data structure 254 includes a patch size field 255 that holds data that indicates the patch size or shape or both. The patch information data structure 254 includes a patch record 256 for each patch. The patch record 256 for each found patch includes a patch location field 257 that holds data that indicates the location of the patch in the two images, e.g., using the upper left most voxel location, or the center voxel location, or any other voxel location and the size and shape of the patch, or the locations of two diagonally opposite corners. In some embodiments, the patch record 256 for each found patch includes a patch similarity field 257 that holds data that indicates the similarity metric, such as covariance or vector distance or the inverse of vector distance. Thus, image patches selected from the same locations in the L and C images should present nearly identical structure when the patches are grouped into a data matrix. Thus, at least one embodiment uses the collaborative data selection through a joint distance-based patch selection.

[0045] A matrix is low-rank if it has many fewer linearly independent rows than columns (i.e., a matrix with a rank that's much lower than the minimum number of rows and columns it has). Such matrices can be efficiently represented using rank-factorizations, which can be used toperform various computations rapidly. In various embodiments, a low rank matrix is produced from the multiple patches found in any pair of C / L images, where each patch vector is one row of the low-rank matrix and thus columns are the voxel values for that patch in the same order for each patch. In various embodiments, the low rank matrix for the control image is stored in a control low rank matrix data structure 260a and the low rank matrix for the labeled image is stored in a label low rank matrix data structure 260b.

[0046] It is common that the number of patches (about 10 to about 100) is less than the number of voxels in a patch (e.g., about 16 to about 1000) so the low-rank matrices 250a, 260b are low rank on their face. However, even if not low rank on their face, the linearly independent vectors that account for most of the variability are expected to be fewer in number than the number of columns (voxels in a patch) and thus a matrix of these vectors is low rank. One way to determine the number of independent rows in a matrix is to compute the eigenvectors of the matrix using any well-known method for doing so. The eigenvectors are orthogonal, i.e., are not linear combinations of each other, and are considered basis vectors that span the space represented by the patches. Each eigenvector has an associated eigenvalue that indicates the variance in pixel values for all patches fit by that eigenvector alone. The eigenvectors ranked by their eigenvalues are the “principal components” in principal component analysis (PCA). The set of eigenvectors with eigenvalues that explain most of the variance constitute a low-rank regularization of the variability in the patches and are expected to capture the denoised structures in the image, as described in more detail below. In some embodiments, the principal eigenvectors and their eigenvalues are stored in separate data structures, not shown.

[0047] Although data structures and fields are depicted in FIG. 2A, and subsequent diagram in FIG. 2B, as integral blocks in a particular order for purposes of illustration, in other embodiments, one or more data structures or fields, or portions thereof, are arranged in a different order, in the same or different number of data structures or databases in one or more hosts or messages, or are omitted, or one or more additional fields are included, or the data structures and messages are changed in some combination of ways. 2. Methods

[0048] FIG. 3 is a flow diagram that illustrates an example of a method 300 to determine a de- noised difference between a labeled medical scan and an un-labeled (control) medical scan of the same region of interest inside a subject, according to an embodiment. Although steps aredepicted in FIG. 3 as integral steps in a particular order for purposes of illustration, in other embodiments, one or more steps, or portions thereof, are performed in a different order, or overlapping in time, in series or in parallel, or are omitted, or one or more additional steps are added, or the method is changed in some combination of ways.

[0049] In step 301, 2D or 3D medical control scans without labeling are collected from a scanner 210 directed to a region of interest in a subject 290 and stored in data structure 251a. In step 303, 2D or 3D medical control scans with labeling are collected from a scanner 210 directed to a region of interest in a subject 290 and stored in data structure 251b.

[0050] In step 305, multiple patches are selected, each patch is co-located in the control and labeled scans and the patterns of amplitude variation in the patch are similar in the control and labeled scans. For example, 2D or 1D covariance is above a threshold or the vector distance is below some threshold. The patches can overlap.. The similarity is high because these patches have relatively low noise and large signal and thus the signal to noise ratio (SNR) is favorable. In some embodiments, step 305 includes selecting a patch size and shape. For example, each patch has a first dimension size in a range from 4 voxels to 32 voxels and a perpendicular second dimension size in a range from 4 voxels to 32 voxels. In some embodiments, the number of voxels in first dimension and second dimension are the same so that each patch has a square shape. In some embodiments, N patches with smallest distance are selected, in which the number N increases with the noise level, e.g., for heavier noise more patches (lager N) are used. In some embodiments, set a threshold τ is set and all patches with distance < τ are included, and meanwhile guarantee there are at least n similar patches. Here N, τ and n are empirically chosen to give advantageous results. Thus, in some of these embodiments, the similarity metric is a low Euclidean distance, and the threshold is selected so that the number of patches, N, is at least 8

[0051] As a further example, if the vectorized version of the ^ × ^ patch at location ^ in the label ^^^ image is ^ (^) ∈ ℝ and that in the control image is ^^(^), the joint Euclidean distance between patches ^(^) and ^(^) is calculated using Equation 1a. ^(^, ^) =^^^^(‖^^(^) − ^^(^)‖^^ + ‖^^(^) − ^^(^)‖^^). (1a) Inway, errors to typical methods. The ^ coordinates with smallest distances are included in the ^-th set of similar patches ^^, which is shared by both ^^(^) and ^^(^).

[0052] In step 307, a low rank variability model is based on the multiple patches in control scan. In some embodiments the variability model is a straightforward eigenvector and eigenvalue determination that keeps the eigenvectors with the largest eigenvalues that sum up to a certain percentage of the total variability, e.g., 50% or 75% or 90% or 95 Since the matrices formed by grouped similar patches are highly low-rank in nature while noise is not, the image signals of interest could be well separated from the noise via low-rank regularization [42-44]. However, in some embodiments it is advantageous to use the following method for regularization.

[0053] In some embodiments, the low rank variability model is computed as described here. The potentially low-rank matrices extracted from the label image, for example, are formed by concatenating the N vectorized patches using Equation 1b. ^ (^) = [ ^ ^^ (1) ( ) ( )^^×# ^^ ^^, ^^^^^2 ^, … ^^^^(^) ] ∈ ℝ (1b) Thosein regional intensities that are ^^(^ and ^ ^ are denoised separately to better)^( )preserve such information. Since the L / C images are treated independently after the data selection stage, the subscripts “$” and “%” can be omitted and denote ^(^)by ^^.

[0054] In step 309, the control scan is de-noised based on control low rank variability model. This amounts to rank minimization that throws out the contributions due to noise. Rank minimization is difficult to solve directly but can be approximated through convex or non-convex regularizations [45, 46]. A non-constant function with more than one global minima is highly likely to be nonconvex. The method according to one embodiment uses non-convex regularization to provide better performance compared to a convex approach [47-49].

[0055] In one embodiment of the method, the non-convex rank surrogate used is the log- determinant of covariance matrix calculated from the patch groups ^^

[0047] . For example, let & ∈ ℝ'×(denote a latent noise-free slice in the L or C image that is H number of voxels in height and W number of voxels in width. Let ) be the corresponding noisy observation, the slice-wise restoration is summarized in the objective function given by Equation 2. * g min^ l ^ ^6^ & = ar ∑og4^5^54 + − ) 8. (2) Here ∙‖^‖8is : by the noise level, which means that heavier noise leads to less contribution from the observation ) to the estimate &;.^5^is the centralized version of ^^extracted from &, given by Equation 3. ^5^= ^^− <^× =6, <^= ∑#>?^^(^^(^)) , (3) where =

[0056] The objective function (Eq. 2) is still difficult to solve directly because the rank penalty is applied to the patch groups while the data fidelity term is applied to the whole slice. According to one embodiment, the noisy version of ^5^is denoted by @5^, which is generated from corresponding patches from the input image ). Assuming that noise in ) is zero-meaned and that the times of each voxel appearing in the patch groups are roughly the same, the slice level regularization can be broken down to patch group level using Equation 4 ‖& − )‖^8 ≈ B‖^5^− @5^‖^8 , (4)level optimization problem using Equation 5. ^5 *^ = arg mC5in log4^5^^5^64 +^ ^E‖^5^ − @5^‖^ 8D, (5) withaggregating all the patches using Equation 6: &*= (∑^H^H^6)I^∑^JH^6^^5 *^ + <^× =6^K , (6)

[0057] Rank minimization using this method is generally mathematically equivalent to a content-adaptive sparsity regularization in the latent PCA domain. The patches stored in ^5^can be viewed as samples derived from the same distribution as ^(^), thus ^5^^5^6∈ ℝ^^×^^reflects the estimate of the corresponding covariance matrix. Suppose L^(M) is the M-th eigenvalue of matrix ^5^^5^6obtained by eigenvalue decomposition given by Equation 7. ^5^^5^6= N^O^N^6, O^= diag[… , L^(M), … ] (7) where N^6∈ ℝ^^×^^could serve as the content-adaptive PCA transform that thoroughly ^(^) to achieve near-optimal sparse representation, and L^(M) corresponds to thevariance of the PCA coefficients in the M-th transform band. Since 4^5^^5^64 = ∏RL^(M) the relation in Equation 8 holds. log4^5^^5^64 =∑Rlog L^(M) . (8)^5 *^ = arg min ∑ log‖^5− @5‖^ ^C5 RDL^(M) +^E ^ ^ 8. (9)

[0058] for Eq 9. At one end of the spectrum, assuming the total energy (i.e., sum of the eigenvalues) remains constant after the optimization to provide: ^5 *^ = arg m*C5in ∑Rlog L^(M) s. t. ∑RL^(M) = ∑R( ) D L^M , (10)observation ). Note that the logarithm function is upper convex, hence for any LV> LX> 0 such that the sum of LVand LXis constant, log LV+ log LXdecreases as the difference between LVand L increases, asin FIG. 4. FIG. 4 is a plot that illustrates an example of the convex nature of the log function, used according to the example embodiment. Such observation can be generalized to L^> L^> ⋯ > L#> 0, where the sum of these values is held constant. That is, more diversified eigenvalues lead to smaller∑Rlog L^(M) , and objective function (Eq. 10) basically encourages the energy of the signal toconcentrated in very few PCA coefficients so that optimal sparsity is attained.

[0059] At the other end, penalizing ∑Rlog L^(M) alone would shrink the eigenvalues towards zero. Sitting between these twocases, minimizing the combination of the log- determinant and the data fidelity term achieves both effects. This means that the energy is concentrated on a small number of PCA coefficients while the coefficients are adaptively reduced towards zero.

[0060] In some embodiments of step 309, the objective function is optimized using numerical methods. In one such embodiment, the denoising method uses an empirical Bayesian procedure

[0050] to solve the optimization problem represented by Eq.5. Since denoising can be viewed as the foundation of most image restoration tasks

[0051] , such numerical solutions can be extended for other applications like compressive sensing reconstruction

[0050] . To simplify the notation, thesubscript "^" can be omitted for the following equations in this section without loss in clarity.

[0061] From a Bayesian point of view, Eq.5 can be considered as a maximum a posteriori (MAP) estimation given by Equation 11. ^5*= arg mC5ax \(^5|@5) = arg mC5ax \(@5|^5)\(^5), (11) where the a priori probability is given by Equation 12. \(@5|^5)∝ expa^ ^ ^^E‖^5− @5‖8b, \(^5)∝4C5C5c4. (12) 12 isdifficult, one embodiment of the method can be configured to be bounded, such as by using a strict convex upper-bound on the log-determinant in Eq.5 via majorization-minimization

[0047] expressed in Equation 13 log|^5^56^ | = mdin# Tr(^5^56fI^) + log|f| + %, (13) whereparameters, and % is a constant. Thus, \(^5) in (12) can be replaced with its lower bound \g(^5) given by Eq. 13 to produce Equation 14. \g(^5; f) = exp a^ 6 I^# Tr(^5^5f ) + log|f|b, (14) andby Equation 15. \g(^5|@5; f)=V(i5|C5)Vg(C5;d)V(i5|C5)Vg(C5;d)kC5. (15)and the method can be configured to estimate ^5using Equation 16. ^5 *^ = f(f + Fl)I^@5, (16)of the first column in @5.

[0063] In another embodiment, determining the hyperparameter f includes minimizing the difference between \(^5) and \g(^5) when the likelihood \(@5|^5) is significant

[0050] as given by Equation 17. f = arg mdinj|\(^5) − \g(^5|f)|\(@5|^5)d^5, ( 17) which leads to an iterative solution: f*^ =J^5^56+ n^fn(f)= f − f(f + Fl)I^f, (18) where M is the number of iteration and f(*o)= @5@56 / ^.

[0064] Returning to FIG. 3, to step 307 and 309, described above,but for the labeled image.

[0065] In step 315, the de-noised control and de-noised label images are differenced to produce the de-noised label difference output 229.

[0066] In step 317, the condition of the subject is assessed based on the label difference output. For example, a patient is treated for insufficient perfusion based on the perfusion measurements made using arterial spin labeling (ASL). In step 319, the scanner 210 is further operated based on the label difference output 229. For example, the scanner is operated to obtain another time slice to track perfusion changes in response to some treatment or enhanced labeling is used to make up for an insufficient label difference measurement. In some embodiments, step 319 includes using the denoised label difference output as “ground truth” in a training instance to train a neural network to automatically produce a de-noised difference image from an input control image and input labeled image.

[0067] Thus, several embodiments provide a denoising method for denoising images using a collaboratively selected image patch set-based low-rank regularized denoising algorithm for labeled imagery. Through enforcing a low-rank property and constraining the data consistency, the denoising method improves denoising and tissue detail preservation performance, compared to several state-of-art denoising methods, as described in more detail below. In some embodiments, the method uses a non-convex low-rank approximation that is connected to a signal-adaptive sparsity in PCA domain. 3. Example Embodiments

[0068] Provided herein are non-limiting examples of the denoising method, as described herein according to many embodiments, for purposes of illustration. These embodiments use an implementation called Locally adaptive low-rank regularization with Collaborative data Selection (LACS) that was evaluated on a real-world ASL MRI dataset.

[0069] In the LACS implementation, the set of similar patches ^^can be efficiently obtained by calculating distances of patch ^^with all the patches andthe distances in an ascending order, then take the first ^ patches (smallest distances) to include the corresponding indexes inthe ^-th set ^^. The parameters ^ and ^ are chosen empirically. According to previous studies, removing heavier noise requires larger ^ and ^ [41, 51]. This example uses similar settings without tuning parameters since no ground-truth ASL image is available for parameter training [41, 51].

[0070] For context, endogenous tracer MR perfusion methods take advantage of signal loss resulting from magnetically labeled water protons (aligned spins) flowing into the imaging plane and exchanging with (random spin) tissue protons. Water protons within inflowing arterial blood are magnetically labeled (or “tagged”) by the application of a special radiofrequency pulse designed to invert spins in a thick slab proximal to the slice of interest. By measuring signal changes between tagged (labeled) images and baseline untagged (control) images, qualitative or quantitative images of cerebral blood flow can be obtained. Inflowing blood may be tagged continuously or intermittently [31, 32]. Although continuous-labeling techniques afford twice as much signal contrast compared with pulsed techniques, they produce substantially more radiofrequency pulse-induced power deposition to the subject. This safety consideration can ultimately limit slice coverage and acquisition time.

[0071] Arterial spin labeling (ASL) perfusion is thus a non-invasive and non-ionizing MRI technique for measuring cerebral blood flow (CBF) [1-5]. ASL MRI labels (with a T1 precession decay) the arterial blood water using radio-frequency pulses and acquires the perfusion weighed MR images after the labeled arterial blood perfuses into the brain tissues in a region of interest. The non-invasive and non-radioactive nature of ASL MRI makes it particularly suitable for perfusion studies in pediatric populations, healthy individuals, patients requiring serial follow-up or with certain conditions (e.g., impaired renal function) [6]. The capability of measuring regional perfusion makes it highly valuable for assessing regional brain functional markers or disease related targets. Over the past decades, ASL MRI has been increasingly used in many different clinical applications [7-11], including vascular malformation [12, 13], cerebrovascular stroke [14, 15], and neurodegenerative diseases [16-21].

[0072] The full value of ASL is still hindered by the inherently low signal-to-noise-ratio (SNR) of the perfusion data, which is mainly caused by the weakness of labeled blood signal due to the varying labeling efficiency, limited time for labeling, and the signal decay (exponentially decay in a rate of blood T1) during the post-labeling transit time. To mitigate this problem, it is standard to acquire multiple pairs of label and control (L / C) images and take the average of thecorresponding perfusion measurements. The rationale behind this approach is that the objects in image pairs acquired within a short period of time can be roughly considered static and the averaging procedure could suppress random noises. This approach will require many L / C pairs to substantially increase SNR, which is impractical for many applications because of the long scan time and the increased risk of motions. Another issue is that averaging across many pairs could blur the image because of the inevitable subtle between frame movements.

[0073] During the past couple of decades, several ASL MRI denoising methods have been proposed [22-27]. Early denoising approaches include spatial smoothing

[0028] , principal component decomposition-based noise removal

[0016] , noise component filtering [29, 30], wavelet-domain filtering

[0031] , independent component decomposition-based noise suppression

[0032] . More recent work has extended into a spatio-temporal denoising framework, including the support vector machine-based CBF quantification and denoising algorithm

[0033] , the total generalized variation regularization-based spatio-temporal filtering

[0034] , and the robust PCA based low-rank and sparse decomposition approach

[0035] . Most of the previous ASL denoising methods either apply a global model to the whole slice or attempt to recover the image at the voxel level. The former lacks attention to the local structures while the latter is prone to be affected by noise.

[0074] In general, these methods applied global models to the images or slices and did not fully exploit the structural characteristics of ASL images. Liang et al. combined nonlocal means

[0036] with dual-tree complex wavelet transform for ASL denoising

[0037] and obtained improved results, but the voxel level solution is easily affected by noise and could be unstable for many image contents. Inspired by the success of deep learning, Xie et al. designed a deep convolutional neural network for ASL denoising

[0038] but it might not work if the signal and noise distribution of the new subjects’ ASL MRI differ dramatically from the training samples. The demand of advanced ASL denoising methods is still outstanding.

[0075] The arterial spin labeling (ASL) technique labels arterial blood water near the tissue under study with radio-frequency pulses, and acquires perfusion-weighted MRI after the labeled perfuse into brain. Meanwhile, it also acquires a control image with phase modulations to the labeling pulses so that arterial spins are roughly unaffected. To remove the background signal, the perfusion is obtained by pair-wise subtraction of the spin label (L) image and the spinuntagged image control (C) image. Then the perfusion is used to calculate the quantitative cerebral blood flow (CBF) in ml / 100g / min [7].

[0076] The labeled blood signal is relatively weak due to the limited time available for labeling arterial blood and the longitudinal signal decay of the labeled blood during the post-labeling delay time. These factors all contribute to a low signal-to-noise-ratio (SNR). A typical way to mitigate this issue is to acquire multiple pairs of L / C images and obtain the mean perfusion map to cancel out some random noise and improve the SNR. This approach inevitably increases the total scan time and may not be practical for clinical use when the patients have difficulty staying in the scanner. Besides, motions caused by breath, heartbeat etc. are inevitable during the acquisition, hence averaging across the measurements could blur the edges and textures. Therefore, it is desirable to enhance the quality of perfusions while reducing the acquisition time.

[0077] To examine the effectiveness of the LACS method, twelve healthy volunteers (age 37.8 ± 14.6 years, 7 males, 5 female) were recruited. Written informed consent was received from each subject as approved by the local institutional review board. MRI studies were performed on a 3 tesla (T) whole-body system (MAGNETOM Prisma, version MR VE11E, Siemens Healthcare, Erlangen, Germany) equipped with a 64-channel head coil. A three-dimensional pseudo-continuous ASL (pCASL) sequence was used [52, 53]. The scan parameters were TR / TE=4600 / 9.3 milliseconds (ms), slice thickness=2 millimeters (mm), FOV=220 × 220mm2, matrix size=110 × 110, slice partial Fourier 5 / 8, 72 partitions, 6 in-plane spiral interleaves acquired in 6 shots (one spiral in one shot), 2-fold parallel imaging acceleration along the kz (partition encoding) direction, 10 repetitions (5 labeled images and 5 control images), labeling duration=2 seconds (s), post label delay=1.8 s. Background suppression was achieved by three frequency offset corrected inversion (FOCI) pulses. The scan time was 5 minutes (min) 53 seconds. The resolution of the scan was 2 × 2 × 2 mm3.

[0078] Evaluating the effectiveness of the LACS method included comparing it with several typical and state-of-the-art denoising methods in the area, including Marchenko-Pastur principle component analysis (MPPCA)

[0054] , tensor MPPCA (tMPPCA)

[0055] , non-local means combined with dual-tree complex wavelet transform (DTCWT)

[0037] and the standard pipeline of taking the average of multiple measurements. DTCWT, MPPCA, tMPPCA and LACS used the first pair of L / C images, while the standard pipeline takes all the available pairs as input. Previously published code, without adjusting settings, was used for the implementation of the comparisons.When the input parameter needs to be specified, the optimal setting indicated herein were used. In one example of comparing the LACS denoising method, an existing technique

[0056] for noise level estimation was used as a comparison.

[0079] FIG. 5A through FIG. 5C are labeled / control (L / C) difference images without de noising and with two prior art de-noising methods, according to an example image modality. These images illustrate ASL perfusion denoising by depicting representative perfusion weighted image slices. FIG. 5A depicts the non-processed data from a single L / C image pair in which noise creates apparent white dot and black space speckle, among other effects. FIG. 5B depicts output of a standard processing pipeline using five C / L pairs. FIG. 5C depicts output of tMPPCA

[0048] using a C / L pair.

[0080] LACS is configured to use structural characteristics of ASL images to select data for low-rank regularization. Compared to typical denoising methods, the method can include collaborative data selection for more robust low-rank modeling and locally adaptive low-rank regularization that corresponds to more efficient signal-adaptive sparse representation without explicit training. The overall framework of the LACS denoising method is shown in FIG. 6. FIG. 6 is a block diagram showing the application of the method of FIG. 3 to MRI ASL perfusion, according to an example embodiment, the LACS implementation. The leftmost panel 610 depicts collaborative data selection between the label and control images, 612a, 512b, respectively, to obtain groups of highly correlated patches 614a, 614b, respectively. The middle panel 620 depicts nonconvex surrogate of low-rank regularization 624a, 624b. respectively, applied to the matrices 622a, 622b, respectively, formed by the patch groups 614a, 614b, respectively, which basically performs content-adaptive filtering in the underlying PCA domain. The rightmost panel 630 depicts a perfusion weighted signal 634 obtained from the difference between the denoised L / C images 632a, 632b, respectively.

[0081] FIG. 7A through FIG. 7D are labeled / control (L / C) difference images as in FIG. 5A through FIG. 5C plus an example of an image with de-noising using the method of FIG. 6, according to the example embodiment. The images of FIG. 7A through FIG. 7C are identical to the images of FIG. 5A through FIG. 5C, respectively, described above. FIG. 7D depicts label difference output 229 of the LACS denoising method using a single C / L pair. Not only is the speckle missing, but the elongated light structures are more continuous, indicating the underlying structure.

[0082] Since the ground truth is not available for ASL denoising, it is not feasible to calculate conventional metrics that require high-quality ground-truth images to serve as references, such as peak signal-to-noise-ratio (PSNR). Instead, a no-reference image quality score known as blind / referenceless image spatial quality evaluator (BRISQUE)

[0057] was used to test the denoising algorithms. BRISQUE compares the evaluated image to a default model computed from images with similar distortions, including noise and blurring. A smaller BRISQUE score indicates better image quality. The BRISQUE results are summarized in Table 1.

[0083] As shown in Table I, LACS significantly reduced BRISQUE in all the cases. On average, LACS decreased BRISQUE by 20.64 compared with raw input, by 18.75 compared with the standard pipeline, by 13.80 compared with DTCWT, by 18.41 compared with MPPCA, TABLE I COMPARISON OF IMAGE QUALITY MEASURED BY BRISQUE. SMALLER BRISQUE SCORE INDICATES BETTER QUALITY. THE BEST BRISQUE SCORE IS MARKED BY BOLD. Image ID Raw Input Standard DTCWT MPPCA tMPPCA LACS 1 42.10 38.07 35.32 40.33 36.58 19.83 2 43.47 42.67 35.43 40.85 34.93 21.94 3 42.69 40.71 34.23 39.93 37.19 19.57 4 42.49 40.40 34.14 40.17 35.04 20.25 5 41.53 39.38 34.95 37.10 33.20 22.07 6 42.12 40.06 39.43 40.41 37.61 20.70 7 42.72 40.96 35.85 39.03 36.35 23.02 8 43.73 41.16 38.87 42.71 40.43 33.13 9 43.52 40.95 34.26 42.12 38.32 22.70 10 43.03 40.75 37.28 42.37 40.38 23.79 11 42.59 40.06 32.92 38.49 35.49 22.10 12 43.58 42.71 35.83 40.76 36.63 21.28 13 41.84 42.97 38.06 41.27 35.02 21.44 14 41.92 39.42 34.89 40.07 34.48 18.80 15 43.71 41.86 37.87 40.99 36.62 21.47 16 43.49 42.09 35.75 42.27 35.31 22.22 AVE 42.78 40.89 35.94 40.55 36.47 22.14 and by 14.33 compared with tMPPCA, corresponding to 48.25%, 45.85%, 38.4%, 45.4% and 39.29% of the BRISQUE scores of the compared methods respectively. Such evident improvement was in line with the visual quality, as shown in FIG. 7D. The reconstructed ASL images produced by the LACS denoising method were much clearer with much less noise and blur, compared to typical methods. Particularly, as can be seen from Fig. 7, some fine structures that other methods failed to reconstruct were successfully recovered by LACS, making thepatterns more consistent between the adjacent slices, even though LACS processed different slices independently.

[0084] In one example, the relative contributions of two steps of the LACS denoising method were determined, i.e., locally adaptive low-rank regularization and collaborative data selection. LACS is compared with the raw input and locally adaptive low-rank regularization without collaborative data selection (LA) in FIG. 8A through FIG. 8C. FIG. 8A and FIG. 8B are images that illustrate examples of L / C difference images without and with collaborative data selection, respectively, according to embodiments. Thes images depict perceptual quality of images re- constructed by LA and LACS, and demonstrate that the LACS denoising method does better recover the textures (e.g., the patterns in the bottom region) and remove noise more effectively, compared to typical methods. FIG. 8C is a plot that illustrates the improvement in image quality (reduction of noise), according to the example embodiment. FIG. 8C plots the relative contributions of LACS versus LA and demonstrates denoising improvement as measured by a decrease in BRISQUE brought by LACS. On average, LA decreased BRISQUE by 10.23 compared with the raw input, and LACS further decreased BRISQUE by 10.41 compared with LA.

[0085] A pair of ASL images allows the LACS denoising method to achieve superior denoising performance compared with the standard pipeline that requires multiple image pairs. For example, the LACS denoising method removes noise more effectively and better preserves useful image structures, compared to typical methods. Furthermore, by using a single C / L pair, the LACS denoising method can be configured to significantly reduce the total acquisition time of ASL MRI, compared to typical methods, since no repetition is required. 4. Computational Hardware Overview

[0086] FIG. 9 is a block diagram that illustrates a computer system 900 upon which an embodiment of the invention may be implemented. Computer system 900 includes a communication mechanism such as a bus 910 for passing information between other internal and external components of the computer system 900. Information is represented as physical signals of a measurable phenomenon, typically electric voltages, but including, in other embodiments, such phenomena as magnetic, electromagnetic, pressure, chemical, molecular atomic and quantum interactions. For example, north and south magnetic fields, or a zero and non-zero electric voltage, represent two states (0, 1) of a binary digit (bit). Other phenomena canrepresent digits of a higher base. A superposition of multiple simultaneous quantum states before measurement represents a quantum bit (qubit). A sequence of one or more digits constitutes digital data that is used to represent a number or code for a character. In some embodiments, information called analog data is represented by a near continuum of measurable values within a particular range. Computer system 900, or a portion thereof, constitutes a means for performing one or more steps of one or more methods described herein.

[0087] A sequence of binary digits constitutes digital data that is used to represent a number or code for a character. A bus 910 includes many parallel conductors of information so that information is transferred quickly among devices coupled to the bus 910. One or more processors 902 for processing information are coupled with the bus 910. A processor 902 performs a set of operations on information. The set of operations include bringing information in from the bus 910 and placing information on the bus 910. The set of operations also typically include comparing two or more units of information, shifting positions of units of information, and combining two or more units of information, such as by addition or multiplication. A sequence of operations to be executed by the processor 902 constitutes computer instructions.

[0088] Computer system 900 also includes a memory 904 coupled to bus 910. The memory 904, such as a random access memory (RAM) or other dynamic storage device, stores information including computer instructions. Dynamic memory allows information stored therein to be changed by the computer system 900. RAM allows a unit of information stored at a location called a memory address to be stored and retrieved independently of information at neighboring addresses. The memory 904 is also used by the processor 902 to store temporary values during execution of computer instructions. The computer system 900 also includes a read only memory (ROM) 906 or other static storage device coupled to the bus 910 for storing static information, including instructions, that is not changed by the computer system 900. Also coupled to bus 910 is a non-volatile (persistent) storage device 908, such as a magnetic disk or optical disk, for storing information, including instructions, that persists even when the computer system 900 is turned off or otherwise loses power.

[0089] Information, including instructions, is provided to the bus 910 for use by the processor from an external input device 912, such as a keyboard containing alphanumeric keys operated by a human user, or a sensor. A sensor detects conditions in its vicinity and transforms those detections into signals compatible with the signals used to represent information in computersystem 900. Other external devices coupled to bus 910, used primarily for interacting with humans, include a display device 914, such as a cathode ray tube (CRT) or a liquid crystal display (LCD), for presenting images, and a pointing device 916, such as a mouse or a trackball or cursor direction keys, for controlling a position of a small cursor image presented on the display 914 and issuing commands associated with graphical elements presented on the display 914.

[0090] In the illustrated embodiment, special purpose hardware, such as an application specific integrated circuit (IC) 920, is coupled to bus 910. The special purpose hardware is configured to perform operations not performed by processor 902 quickly enough for special purposes. Examples of application specific ICs include graphics accelerator cards for generating images for display 914, cryptographic boards for encrypting and decrypting messages sent over a network, speech recognition, and interfaces to special external devices, such as robotic arms and medical scanning equipment that repeatedly perform some complex sequence of operations that are more efficiently implemented in hardware.

[0091] Computer system 900 also includes one or more instances of a communications interface 970 coupled to bus 910. Communication interface 970 provides a two-way communication coupling to a variety of external devices that operate with their own processors, such as printers, scanners and external disks. In general the coupling is with a network link 978 that is connected to a local network 980 to which a variety of external devices with their own processors are connected. For example, communication interface 970 may be a parallel port or a serial port or a universal serial bus (USB) port on a personal computer. In some embodiments, communications interface 970 is an integrated services digital network (ISDN) card or a digital subscriber line (DSL) card or a telephone modem that provides an information communication connection to a corresponding type of telephone line. In some embodiments, a communication interface 970 is a cable modem that converts signals on bus 910 into signals for a communication connection over a coaxial cable or into optical signals for a communication connection over a fiber optic cable. As another example, communications interface 970 may be a local area network (LAN) card to provide a data communication connection to a compatible LAN, such as Ethernet. Wireless links may also be implemented. Carrier waves, such as acoustic waves and electromagnetic waves, including radio, optical and infrared waves travel through space without wires or cables. Signals include man-made variations in amplitude, frequency, phase,polarization or other physical properties of carrier waves. For wireless links, the communications interface 970 sends and receives electrical, acoustic or electromagnetic signals, including infrared and optical signals, that carry information streams, such as digital data.

[0092] The term computer-readable medium is used herein to refer to any medium that participates in providing information to processor 902, including instructions for execution. Such a medium may take many forms, including, but not limited to, non-volatile media, volatile media and transmission media. Non-volatile media include, for example, optical or magnetic disks, such as storage device 908. Volatile media include, for example, dynamic memory 904. Transmission media include, for example, coaxial cables, copper wire, fiber optic cables, and waves that travel through space without wires or cables, such as acoustic waves and electromagnetic waves, including radio, optical and infrared waves. The term computer-readable storage medium is used herein to refer to any medium that participates in providing information to processor 902, except for transmission media.

[0093] Common forms of computer-readable media include, for example, a floppy disk, a flexible disk, a hard disk, a magnetic tape, or any other magnetic medium, a compact disk ROM (CD-ROM), a digital video disk (DVD) or any other optical medium, punch cards, paper tape, or any other physical medium with patterns of holes, a RAM, a programmable ROM (PROM), an erasable PROM (EPROM), a FLASH-EPROM, or any other memory chip or cartridge, a carrier wave, or any other medium from which a computer can read. The term non-transitory computer- readable storage medium is used herein to refer to any medium that participates in providing information to processor 902, except for carrier waves and other signals.

[0094] Logic encoded in one or more tangible media includes one or both of processor instructions on a computer-readable storage media and special purpose hardware, such as ASIC *420.

[0095] Network link 978 typically provides information communication through one or more networks to other devices that use or process the information. For example, network link 978 may provide a connection through local network 980 to a host computer 982 or to equipment 984 operated by an Internet Service Provider (ISP). ISP equipment 984 in turn provides data communication services through the public, world-wide packet-switching communication network of networks now commonly referred to as the Internet 990. A computer called a server 992 connected to the Internet provides a service in response to information received over theInternet. For example, server 992 provides information representing video data for presentation at display 914.

[0096] The invention is related to the use of computer system 900 for implementing the techniques described herein. According to one embodiment of the invention, those techniques are performed by computer system 900 in response to processor 902 executing one or more sequences of one or more instructions contained in memory 904. Such instructions, also called software and program code, may be read into memory 904 from another computer-readable medium such as storage device 908. Execution of the sequences of instructions contained in memory 904 causes processor 902 to perform the method steps described herein. In alternative embodiments, hardware, such as application specific integrated circuit 920, may be used in place of or in combination with software to implement the invention. Thus, embodiments of the invention are not limited to any specific combination of hardware and software.

[0097] The signals transmitted over network link 978 and other networks through communications interface 970, carry information to and from computer system 900. Computer system 900 can send and receive information, including program code, through the networks 980, 990 among others, through network link 978 and communications interface 970. In an example using the Internet 990, a server 992 transmits program code for a particular application, requested by a message sent from computer 900, through Internet 990, ISP equipment 984, local network 980 and communications interface 970. The received code may be executed by processor 902 as it is received, or may be stored in storage device 908 or other non-volatile storage for later execution, or both. In this manner, computer system 900 may obtain application program code in the form of a signal on a carrier wave.

[0098] Various forms of computer readable media may be involved in carrying one or more sequence of instructions or data or both to processor 902 for execution. For example, instructions and data may initially be carried on a magnetic disk of a remote computer such as host 982. The remote computer loads the instructions and data into its dynamic memory and sends the instructions and data over a telephone line using a modem. A modem local to the computer system 900 receives the instructions and data on a telephone line and uses an infra-red transmitter to convert the instructions and data to a signal on an infra-red a carrier wave serving as the network link 978. An infrared detector serving as communications interface 970 receives the instructions and data carried in the infrared signal and places information representing theinstructions and data onto bus 910. Bus 910 carries the information to memory 904 from which processor 902 retrieves and executes the instructions using some of the data sent with the instructions. The instructions and data received in memory 904 may optionally be stored on storage device 908, either before or after execution by the processor 902.

[0099] FIG. 10 illustrates a chip set 1000 upon which an embodiment of the invention may be implemented. Chip set 1000 is programmed to perform one or more steps of a method described herein and includes, for instance, the processor and memory components described with respect to FIG. 9 incorporated in one or more physical packages (e.g., chips). By way of example, a physical package includes an arrangement of one or more materials, components, and / or wires on a structural assembly (e.g., a baseboard) to provide one or more characteristics such as physical strength, conservation of size, and / or limitation of electrical interaction. It is contemplated that in certain embodiments the chip set can be implemented in a single chip. Chip set 1000, or a portion thereof, constitutes a means for performing one or more steps of a method described herein.

[0100] In one embodiment, the chip set 1000 includes a communication mechanism such as a bus 1001 for passing information among the components of the chip set 1000. A processor 1003 has connectivity to the bus 1001 to execute instructions and process information stored in, for example, a memory 1005. The processor 1003 may include one or more processing cores with each core configured to perform independently. A multi-core processor enables multiprocessing within a single physical package. Examples of a multi-core processor include two, four, eight, or greater numbers of processing cores. Alternatively or in addition, the processor 1003 may include one or more microprocessors configured in tandem via the bus 1001 to enable independent execution of instructions, pipelining, and multithreading. The processor 1003 may also be accompanied with one or more specialized components to perform certain processing functions and tasks such as one or more digital signal processors (DSP) 1007, or one or more application-specific integrated circuits (ASIC) 1009. A DSP 1007 typically is configured to process real-world signals (e.g., sound) in real time independently of the processor 1003. Similarly, an ASIC 1009 can be configured to performed specialized functions not easily performed by a general purposed processor. Other specialized components to aid in performing the inventive functions described herein include one or more field programmable gate arrays(FPGA) (not shown), one or more controllers (not shown), or one or more other special-purpose computer chips.

[0101] The processor 1003 and accompanying components have connectivity to the memory 1005 via the bus 1001. The memory 1005 includes both dynamic memory (e.g., RAM, magnetic disk, writable optical disk, etc.) and static memory (e.g., ROM, CD-ROM, etc.) for storing executable instructions that when executed perform one or more steps of a method described herein. The memory 1005 also stores the data associated with or generated by the execution of one or more steps of the methods described herein. 5. Alternatives, Deviations and modifications

[0102] In the foregoing specification, the invention has been described with reference to specific embodiments thereof. It will, however, be evident that various modifications and changes may be made thereto without departing from the broader spirit and scope of the invention. The specification and drawings are, accordingly, to be regarded in an illustrative rather than a restrictive sense. Throughout this specification and the claims, unless the context requires otherwise, the word “comprise” and its variations, such as “comprises” and “comprising,” will be understood to imply the inclusion of a stated item, element or step or group of items, elements or steps but not the exclusion of any other item, element or step or group of items, elements or steps. Furthermore, the indefinite article “a” or “an” is meant to indicate one or more of the item, element or step modified by the article.

[0103] Notwithstanding that the numerical ranges and parameters setting forth the broad scope are approximations, the numerical values set forth in specific non-limiting examples are reported as precisely as possible. Any numerical value, however, inherently contains certain errors necessarily resulting from the standard deviation found in their respective testing measurements at the time of this writing. Furthermore, unless otherwise clear from the context, a numerical value presented herein has an implied precision given by the least significant digit. Thus, a value 1.1 implies a value from 1.05 to 1.15. The term “about” is used to indicate a broader range centered on the given value, and unless otherwise clear from the context implies a broader range around the least significant digit, such as “about 1.1” implies a range from 1.0 to 1.2. If the least significant digit is unclear, then the term “about” implies a factor of two, e.g., “about X” implies a value in the range from 0.5X to 2X, for example, about 100 implies a value in a range from 50 to 200. Moreover, all ranges disclosed herein are to be understood to encompass any and all sub-ranges subsumed therein. For example, a range of "less than 10" for a positive only parameter can include any and all sub-ranges between (and including) the minimum value of zero and the maximum value of 10, that is, any and all sub-ranges having a minimum value of equal to or greater than zero and a maximum value of equal to or less than 10, e.g., 1 to 4. 6. References

[0104] All the references listed here are hereby incorporated by reference as if fully set forth herein except for terminology inconsistent with that used herein. [1] J. A. Detre, J. S. Leigh, D. S. Williams, and A. P. Koretsky, "Perfusion imaging," Magnetic Resonance in Medicine, vol. 23, no. 1, pp. 37-45, 1992. [2] D. S. Williams, J. A. Detre, J. S. Leigh, and A. P. Koretsky, "Magnetic resonance imaging of perfusion using spin inversion of arterial water," Proceedings of the National Academy of Sciences, vol. 89, no. 1, pp. 212-216, 1992. [3] J. A. Detre, J. Wang, Z. Wang, and H. Rao, "Arterial spin-labeled perfusion MRI in basic and clinical neuroscience," Current Opinion in Neurology, vol. 22, no. 4, pp. 348-355, 2009. [4] M. A. Fernández‐Seara et al., "Imaging mesial temporal lobe activation during scene encoding: comparison of fMRI using BOLD and arterial spin labeling," Human Brain Mapping, vol. 28, no. 12, pp. 1391-1400, 2007. [5] Y. V. Chang, M. Vidorreta, Z. Wang, and J. A. Detre, "3D‐accelerated, stack‐of‐spirals acquisitions and reconstruction of arterial spin labeling MRI," Magnetic Resonance in Medicine, vol. 78, no. 4, pp. 1405-1419, 2017. [6] Z. Wang et al., "Assessment of functional development in normal infant brain using arterial spin labeled perfusion MRI," Neuroimage, vol. 39, no. 3, pp. 973-978, 2008. [7] D. C. Alsop et al., "Recommended implementation of arterial spin‐labeled perfusion MRI for clinical applications: a consensus of the ISMRM perfusion study group and the European consortium for ASL in dementia," Magnetic Resonance in Medicine, vol. 73, no. 1, pp. 102-116, 2015. [8] J. A. Detre, H. Rao, D. J. Wang, Y. F. Chen, and Z. Wang, "Applications of arterial spin labeled MRI in the brain," Journal of Magnetic Resonance Imaging, vol. 35, no. 5, pp. 1026- 1037, 2012.[9] Y. Chen et al., "Quantification of cerebral blood flow as biomarker of drug effect: arterial spin labeling phMRI after a single dose of oral citalopram," Clinical Pharmacology & Therapeutics, vol. 89, no. 2, pp. 251-258, 2011.

[0010] G. Zheng et al., "Cerebral blood flow measured by arterial-spin labeling MRI: a useful biomarker for characterization of minimal hepatic encephalopathy in patients with cirrhosis," European Journal of Radiology, vol. 82, no. 11, pp. 1981-1988, 2013.

[0011] G. Zheng et al., "Changes in cerebral blood flow after transjugular intrahepatic portosystemic shunt can help predict the development of hepatic encephalopathy: an arterial spin labeling MR study," European Journal of Radiology, vol. 81, no. 12, pp. 3851-3856, 2012.

[0012] T. Le, N. Fischbein, J. André, C. Wijman, J. Rosenberg, and G. Zaharchuk, "Identification of venous signal on arterial spin labeling improves diagnosis of dural arteriovenous fistulas and small arteriovenous malformations," American Journal of Neuroradiology, vol. 33, no. 1, pp.61-68, 2012.

[0013] T. Blauwblomme et al., "Arterial spin labeling magnetic resonance imaging: toward noninvasive diagnosis and follow-up of pediatric brain arteriovenous malformations," Journal of Neurosurgery: Pediatrics, vol. 15, no. 4, pp. 451-458, 2015.

[0014] G. Zaharchuk, "Arterial spin–labeled perfusion imaging in acute ischemic stroke," Stroke, vol. 45, no. 4, pp. 1202-1207, 2014.

[0015] J. A. Chalela, D. C. Alsop, J. B. Gonzalez-Atavales, J. A. Maldjian, S. E. Kasner, and J. A. Detre, "Magnetic resonance perfusion imaging in acute ischemic stroke using continuous arterial spin labeling," Stroke, vol. 31, no. 3, pp. 680-687, 2000.

[0016] W. Hu, Z. Wang, V.-Y. Lee, J. Trojanowski, J. Detre, and M. Grossman, "Distinct cerebral perfusion patterns in FTLD and AD," Neurology, vol. 75, no. 10, pp. 881-888, 2010.

[0017] S. Dolui et al., "Comparison of PASL, PCASL, and background‐suppressed 3D PCASL in mild cognitive impairment," Human Brain Mapping, vol. 38, no. 10, pp. 5260-5273, 2017.

[0018] Z. Wang, "Characterizing early Alzheimer's disease and disease progression using hippocampal volume and arterial spin labeling perfusion MRI," Journal of Alzheimer's Disease, vol. 42, no. s4, pp. S495-S502, 2014.

[0019] Q. Zhang, R. B. Stafford, Z. Wang, S. E. Arnold, D. A. Wolk, and J. A. Detre, "Microvascular perfusion based on arterial spin labeled perfusion MRI as a measure of vascular risk in Alzheimer's disease," Journal of Alzheimer's Disease, vol. 32, no. 3, pp. 677-687, 2012.

[0020] Z. Wang et al., "Arterial spin labeled MRI in prodromal Alzheimer's disease: a multi-site study," Neuroimage: Clinical, vol. 2, pp. 630-636, 2013.

[0021] A. Camargo, Z. Wang, and A. s. D. N. Initiative, "Longitudinal cerebral blood flow changes in normal aging and the Alzheimer’s disease continuum identified by arterial spin labeling MRI," Journal of Alzheimer's Disease, vol. 81, no. 4, pp. 1727-1735, 2021.

[0022] L. Zhang et al., "Improving Sensitivity of Arterial Spin Labeling Perfusion MRI in Alzheimer's Disease Using Transfer Learning of Deep Learning‐Based ASL Denoising," Journal of Magnetic Resonance Imaging, vol. 55, no. 6, pp. 1710-1722, 2022.

[0023] L. Hernandez‐Garcia et al., "Recent technical developments in ASL: a review of the state of the art," Magnetic Resonance in Medicine, vol. 88, no. 5, pp. 2021-2042, 2022.

[0024] H. Zhu, G. He, and Z. Wang, "Patch-based local learning method for cerebral blood flow quantification with arterial spin-labeling MRI," Medical & Biological Engineering & Computing, vol. 56, pp. 951-956, 2018.

[0025] Z. Wang et al., "Empirical optimization of ASL data analysis using an ASL data processing toolbox: ASLtbx," Magnetic Resonance Imaging, vol. 26, no. 2, pp. 261-269, 2008.

[0026] Z. Li et al., "A two-stage multi-loss super-resolution network for arterial spin labeling magnetic resonance imaging," in International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI), 2019, pp. 12-20.

[0027] Z. Wang, "Arterial spin labeling perfusion MRI signal processing through traditional methods and machine learning," Investigative Magnetic Resonance Imaging, vol. 26, no. 4, p. 220, 2022.

[0028] J. Wang, Z. Wang, G. K. Aguirre, and J. A. Detre, "To smooth or not to smooth? ROC analysis of perfusion fMRI data," Magnetic Resonance Imaging, vol. 23, no. 1, pp. 75-81, 2005.

[0029] Y. Behzadi, K. Restom, J. Liau, and T. T. Liu, "A component based noise correction method (CompCor) for BOLD and perfusion based fMRI," Neuroimage, vol. 37, no. 1, pp. 90- 101, 2007.

[0030] Z. Wang, "Improving cerebral blood flow quantification for arterial spin labeled perfusion MRI by removing residual motion artifacts and global signal fluctuations," Magnetic Resonance Imaging, vol. 30, no. 10, pp. 1409-1415, 2012.

[0031] A. Bibic, L. Knutsson, F. Ståhlberg, and R. Wirestam, "Denoising of arterial spin labeling data: wavelet-domain filtering compared with Gaussian smoothing," Magnetic Resonance Materials in Physics, Biology and Medicine, vol. 23, pp. 125-137, 2010.

[0032] J. A. Wells, D. L. Thomas, M. D. King, A. Connelly, M. F. Lythgoe, and F. Calamante, "Reduction of errors in ASL cerebral perfusion and arterial transit time maps using image de‐ noising," Magnetic Resonance in Medicine, vol. 64, no. 3, pp. 715-724, 2010.

[0033] Z. Wang, "Support vector machine learning-based cerebral blood flow quantification for arterial spin labeling MRI," Human Brain Mapping, vol. 35, no. 7, pp. 2869-2875, 2014.

[0034] S. M. Spann, K. S. Kazimierski, C. S. Aigner, M. Kraiger, K. Bredies, and R. Stollberger, "Spatio-temporal TGV denoising for ASL perfusion imaging," Neuroimage, vol. 157, pp. 81-96, 2017.

[0035] H. Zhu, J. Zhang, and Z. Wang, "Arterial spin labeling perfusion MRI signal denoising using robust principal component analysis," Journal of Neuroscience Methods, vol. 295, pp. 10- 19, 2018.

[0036] A. Buades, B. Coll, and J.-M. Morel, "A non-local algorithm for image denoising," in IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2005, vol. 2, pp. 60-65.

[0037] X. Liang, A. Connelly, and F. Calamante, "Voxel-wise functional connectomics using arterial spin labeling functional magnetic resonance imaging: the role of denoising," Brain Connectivity, vol. 5, no. 9, pp. 543-553, 2015.

[0038] D. Xie et al., "Denoising arterial spin labeling perfusion MRI with deep machine learning," Magnetic Resonance Imaging, vol. 68, pp. 95-105, 2020.

[0039] H. Liu, J. Zhang, and R. Xiong, "CAS: Correlation adaptive sparse modeling for image denoising," IEEE Transactions on Computational Imaging, vol. 7, pp. 638-647, 2021.

[0040] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Image denoising by sparse 3-D transform-domain collaborative filtering," IEEE Transactions on Image Processing, vol. 16, no. 8, pp. 2080-2095, 2007.

[0041] H. Liu, R. Xiong, J. Zhang, and W. Gao, "Image denoising via adaptive soft- thresholding based on non-local samples," in IEEE Conference on Computer Vision and Pattern Recognition, 2015, pp. 484-492.

[0042] M. Tao and X. Yuan, "Recovering low-rank and sparse components of matrices from incomplete and noisy observations," SIAM Journal on Optimization, vol. 21, no. 1, pp. 57-81, 2011.

[0043] H. Liu, R. Xiong, D. Liu, F. Wu, and W. Gao, "Low rank regularization exploiting intra and inter patch correlation for image denoising," in IEEE Visual Communications and Image Processing (VCIP), 2017, pp. 1-4.

[0044] Y. Zhang et al., "Image denoising via structure-constrained low-rank approximation," Neural Computing and Applications, vol. 32, pp. 12575-12590, 2020.

[0045] J. Wright, A. Ganesh, S. Rao, Y. Peng, and Y. Ma, "Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization," Advances in Neural Information Processing Systems, vol. 22, 2009.

[0046] V. Koltchinskii, K. Lounici, and A. B. Tsybakov, "Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion," The Annals of Statistics, vol. 39, no. 5, pp. 2302-2329, 28, 2011.

[0047] K. Mohan and M. Fazel, "Iterative reweighted algorithms for matrix rank minimization," The Journal of Machine Learning Research, vol. 13, no. 1, pp. 3441-3473, 2012.

[0048] Z. Zha et al., "Non-convex weighted ℓp nuclear norm based ADMM framework for image restoration," Neurocomputing, vol. 311, pp. 209-224, 2018.

[0049] H. Liu, X. Zhang, and R. Xiong, "Content-adaptive low rank regularization for image denoising," in IEEE International Conference on Image Processing (ICIP), 2016, pp. 3091-3095.

[0050] D. P. Wipf, B. D. Rao, and S. Nagarajan, "Latent variable Bayesian models for promoting sparsity," IEEE Transactions on Information Theory, vol. 57, no. 9, pp. 6236-6255, 2011.

[0051] H. Liu, R. Xiong, X. Zhang, Y. Zhang, S. Ma, and W. Gao, "Nonlocal gradient sparsity regularization for image restoration," IEEE Transactions on Circuits and Systems for Video Technology, vol. 27, no. 9, pp. 1909-1921, 2016.

[0052] M. Vidorreta, Z. Wang, I. Rodríguez, M. A. Pastor, J. A. Detre, and M. A. Fernández- Seara, "Comparison of 2D and 3D single-shot ASL perfusion fMRI sequences," Neuroimage, vol. 66, pp. 662-671, 2013.

[0053] M. Vidorreta, Z. Wang, Y. V. Chang, D. A. Wolk, M. A. Fernández-Seara, and J. A. Detre, "Whole-brain background-suppressed pCASL MRI with 1D-accelerated 3D RARE Stack- Of-Spirals readout," PLOS ONE, vol. 12, no. 8, p. e0183762, 2017.

[0054] J. Veraart, E. Fieremans, and D. S. Novikov, "Diffusion MRI noise mapping using random matrix theory," Magnetic Resonance in Medicine, vol. 76, no. 5, pp. 1582-1593, 2016.

[0055] J. L. Olesen, A. Ianus, L. Østergaard, N. Shemesh, and S. N. Jespersen, "Tensor denoising of multidimensional MRI data," Magnetic Resonance in Medicine, vol. 89, no. 3, pp. 1160-1172, 2023.

[0056] P. Coupé, J. V. Manjón, E. Gedamu, D. Arnold, M. Robles, and D. L. Collins, "An object-based method for Rician noise estimation in MR images," in International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI), 2009: Springer, pp. 601-608.

[0057] A. Mittal, A. K. Moorthy, and A. C. Bovik, "No-reference image quality assessment in the spatial domain," IEEE Transactions on Image Processing, vol. 21, no. 12, pp. 4695-4708, 2012. ________________________________________________

Claims

CLAIMS What is claimed is:

1. A method for de-noising label-differenced medical images, comprising: operating a medical scanner to collect a control image of a region of interest in a subject; operating the medical scanner to collect a labeled image of the region of interest in the subject; selecting, automatically on a processor, a plurality of patches within the control image, each patch having a similarity metric value better than a threshold when compared to a patch in the labeled image at a same location in the region of interest, wherein each patch of the plurality of patches contains fewer than half a number of voxels in the region of interest in the control image; deriving, automatically on the processor, a control low-rank variability model for the plurality of patches within the control image and producing a de-noised control image in the region of interest based on the control low-rank variability model; deriving, automatically on the processor, a labeled low-rank variability model for the plurality of patches within the labeled image and producing a de-noised labeled image in the region of interest based on the labeled low-rank variability model; and outputting, automatically from the processor to a display device, a difference image in the region of interest based on a difference between the de-noised labeled image in the region of interest and the de-noised control image in the region of interest.

2. The method as recited in claim 1, wherein each patch has a first dimension size in a range from 4 voxels to 32 voxels and a perpendicular second dimension size in a range from 4 voxels to 32 voxels.

3. The method as recited in claim 1, wherein each patch has a square shape.

4. The method as recited in claim 1, wherein the similarity metric is a low vector distance between vectors representing each patch.

5. The method as recited in claim 4, wherein the similarity metric is a low Euclidean distance and the threshold is selected so that a number of patches is at least 8.

6. The method as recited in claim 1, wherein the low-rank variability model for each of the control image and the label image is based on a subset of eigenvectors having the highest eigenvalues for a low-rank matrix having a low-rank first dimension comprising the voxels in each patch and a low-rank second dimension comprising the plurality of patches.

7. The method as recited in claim 6, wherein the eigenvalues are based on a non-convex regularization of the low-rank matrix.

8. The method as recited in claim 7, wherein the non-convex regularization of the low-rank matrix is based on a log-determinant of a square covariance matrix for the low-rank matrix.

9. The method as recited in claim 1, wherein: the medical scanner is a Magnetic Resonance Imager (MRI) scanner; the labeled image is an arterial spin labeling (ASL) image; and the difference image is a perfusion image.

10. The method as recited in claim 9, wherein the region of interest transects a brain of the subject and the perfusion image indicates cerebral blood flow (CBF).

11. A non-transitory computer-readable medium carrying one or more sequences of instructions, wherein execution of the one or more sequences of instructions by one or more processors causes the one or more processors to perform the steps of: operating a medical scanner to collect a control image of a region of interest in a subject; operating the medical scanner to collect a labeled image of the region of interest in the subject; selecting, automatically on a processor, a plurality of patches within the control image, each patch having a similarity metric value better than a threshold when compared to a patch in the labeled image at a same location in the region of interest, wherein each patch of the plurality of patches contains fewer than half a number of voxels in the region of interest in the control image; deriving, automatically on the processor, a control low-rank variability model for the plurality of patches within the control image and producing a de-noised control image in the region of interest based on the control low-rank variability model; andderiving, automatically on the processor, a labeled low-rank variability model for the plurality of patches within the labeled image and producing a de-noised labeled image in the region of interest based on the labeled low-rank variability model; outputting, automatically from the processor to a display device, a difference image in the region of interest based on a difference between the de-noised labeled image in the region of interest and the de-noised control image in the region of interest.

12. A system comprising: a medical scanner; at least one processor; and at least one memory including one or more sequences of instructions, the at least one memory and the one or more sequences of instructions configured to, with the at least one processor, cause the apparatus to perform at least the following, operating the medical scanner to collect a control image of a region of interest in a subject; operating the medical scanner to collect a labeled image of the region of interest in the subject; selecting, automatically on a processor, a plurality of patches within the control image, each patch having a similarity metric value better than a threshold when compared to a patch in the labeled image at a same location in the region of interest, wherein each patch of the plurality of patches contains fewer than half a number of voxels in the region of interest in the control image; deriving, automatically on the processor, a control low-rank variability model for the plurality of patches within the control image and producing a de-noised control image in the region of interest based on the control low-rank variability model; and deriving, automatically on the processor, a labeled low-rank variability model for the plurality of patches within the labeled image and producing a de-noised labeled image in the region of interest based on the labeled low-rank variability model; outputting, automatically from the processor to a display device, a difference image in the region of interest based on a difference between the de-noised labeled image in the region of interest and the de-noised control image in the region of interest.

Images