Method and System for Unsupervised Cross-Modal Medical Image Synthesis

[0001] This application claims the benefit of U.S. Provisional Application No. 62/077,523, filed Nov. 10, 2014, the disclosure of which is herein incorporated by reference.

[0002] The present invention relates to synthesis of medical image data, and more particularly, to synthesizing subject-specific medical image data across image domains or image modalities.

[0003] A multitude of imaging modalities, such as X-ray, computed tomography (CT), diffuser tensor imaging (DTI), T1-weighted magnetic resonance imaging (MRI), T2-weighted MRI, etc., can be used for medical image analysis of a of a patient. Each of these imaging modalities captures different characteristics of the underlying anatomy and the relationship between any two modalities is highly nonlinear.

[0004] In many practical medical image analysis problems, a situation is often encountered in which medical image data available for training, for example for machine learning based anatomical object detection, has a different distribution or representation than the medical image data given during testing due to modality heterogeneity or domain variation. Due to variations in the image characteristics across modalities, medical image analyses algorithms trained with data from one modality may not work well when applied to medical image data from a different modality. A straightforward way to address this issue is to collect large amounts of training data from each imaging modality. However, this solution is impractical since collecting medical images is often time consuming and expensive.

[0005] Cross-modal synthesis generates medical images in a desired target modality from given source modality images. The ability to synthesize medical images without actual acquisition has many potential applications, such as atlas construction, multi-modal registration, super-resolution, and building virtual models. Various approaches for cross-modal synthesis have been proposed, but such approaches are either tailored to specific applications or work under a supervised setting in which training data from the same set of subjects in both the source and target modalities is required. Availability of such paired data is often limited and collecting such paired data is not desirable because each subject must be scanned multiple times. Accordingly, an unsupervised cross-modal medical image synthesis approach that generates target modality images without the need for paired training data is desirable.

[0006] The present invention provides a method and system for unsupervised cross-modal synthesis of medical images. Embodiments of the present invention can be used with any pair of imaging modalities and do not require paired training data from the source and target modalities.

[0007] In one embodiment of the present invention, a source modality medical image is received. Multiple candidate target modality intensity values are generated for each of a plurality of voxels of a target modality medical image based on corresponding voxels in the source modality medical image. A synthesized target modality medical image is generated by selecting, jointly for all of the plurality of voxels in the target modality medical image, intensity values from the multiple candidate target modality intensity values generated for each of the plurality of voxels. The synthesized target modality medical image can be refined using coupled sparse representation.

[0008] These and other advantages of the invention will be apparent to those of ordinary skill in the art by reference to the following detailed description and the accompanying drawings.

[0014] The present invention relates to a method and system for unsupervised cross-modal synthesis of medical images. Embodiments of the present invention are described herein to give a visual understanding of the medical image synthesis method. A digital image is often composed of digital representations of one or more objects (or shapes). The digital representation of an object is often described herein in terms of identifying and manipulating the objects. Such manipulations are virtual manipulations accomplished in the memory or other circuitry/hardware of a computer system. Accordingly, is to be understood that embodiments of the present invention may be performed within a computer system using data stored within the computer system.

[0015] Embodiments of the present invention provide a generalized and robust framework for cross-modal synthesis of medical images. Embodiments of the present invention can be used to synthesize medical images in a target modality from available images in a source modality without having to perform image acquisition in the target modality. Embodiments of the present invention may be used to synthesize target modality medical images in order to create large training set of target modality medical images for training machine learning based classifiers for anatomical object detection, segmentation, tracking, and classification, without having to perform additional image acquisition on a large number of subjects. In addition, embodiments of the present invention may be used to synthesize target modality medical images for other applications, such as to create visualization tools for virtual domains, to perform cross-modality registration, or to up-sample the resolution of image data. As used herein, cross-model synthesis refers to synthesis of medical images across medical imaging modalities, such as synthesizing a CT image from an MR image, as well as synthesis of images across an image domain, such MR images with different protocols (e.g., T1 and T2), contrast CT images and non-contrast CT images, CT image captured with low kV and CT images captured with high kV, or any type of low resolution medical image to a corresponding high resolution medical image. That is, different modalities may refer to different domains or protocols within the same overall imaging modalities, and the “source modality” and “target modality” may be completely different medical imaging modalities or different image domains or protocols within the same overall imaging modality. FIG. 1 illustrates images of the human brain in three different imaging modalities. It is to be understood that the present invention is not limited to the above examples and may be performed to synthesize any type of target modality image from any type of source modality image.

[0016] FIG. 1 illustrates images of the human brain in three different imaging modalities. As shown in FIG. 1, image 102 is a T1 MR image of the brain, image 104 is a T2 MR image of the brain, and image 106 is a diffusion tensor imaging-fractional anisotropy (DTI-FA) image of the brain. FIG. 2 illustrates an overview of cross-modal synthesis of medical images according to an embodiment of the present invention. In the exemplary embodiment of FIG. 2, cross-modal synthesis is used to synthesize a T2 MR image from an acquired T1 MR image. The left side (200) of FIG. 2 shows the acquisition of T1 and T2 MR images without using cross-modal synthesis. In this case, the patient must be scanned twice to acquire T1 and T2 MR images of the patient. The right side (210) of FIG. 2 shows the use of cross-modal synthesis to generate the T2 MR image. In this case, the patient is scanned once to acquire a T1 MR image (source modality image) and the T2 MR image (target modality image) is synthesized from the acquired T1 MR image.

[0017] Embodiments of the present invention provide a general fully-unsupervised method for cross-modal synthesis of subject-specific medical images. Embodiments of the present invention can be used with any pair of imaging modalities and do not require paired training data from the source and target imaging modalities. Since synthesizing a full medical image is a fairly complex task, embodiments of the present invention break this task into steps of candidate generation and candidate selection. Given a source modality image, multiple target modality candidate values are generated for each voxel using a cross-modal nearest neighbor search. The best candidate values are then selected for all the voxels jointly by solving an optimization problem that simultaneously maximizes a global mutual information cost function and minimizes a local spatial consistency cost function, resulting in a synthesized full target modality image. Coupled sparse representation can then be used to further refine the synthesized target modality image.

[0018] FIG. 3 illustrates a method for unsupervised cross-modal synthesis of medical images according to an embodiment of the present invention. The method of FIG. 3 transforms medical image data of a first modality (source modality) to generate a new medical image of a second modality (target modality). The method of FIG. 3 can be used for any type of cross-modal modality medical image synthesis, including but not limited to, synthesizing an MR image of one type of MR protocol (e.g., T2 MR image) from an MR image of another type of MR protocol (e.g., T1 MR image), synthesizing a contrast CT image from a non-contrast CT image, synthesizing a high-kV CT image from a low-kV CT image, synthesizing a high resolution medical image from a low resolution medical image of the same modality, synthesizing a CT image from an MR image, or for any other combination of imaging modalities (e.g. MR, CT, PET, DynaCT, etc.).

[0019] At step 302, a source modality medical image is received. The source modality medical image will also be referred to herein as the “source image”. The source image can be acquired using any type of imaging modality, such as MR, CT, Ultrasound, X-ray fluoroscopy, DynaCT, positron emission tomography (PET), etc.

[0020] The medical image can be a 2D or 3D medical image. It is to be understood that although the medical image can be 2D or 3D, we use the term “voxel” herein to refer to elements of the medical image, regardless of the dimensionality of the medical image. In one possible implementation, the medical image can be a previously acquired medical image that is stored on a memory or storage of a computer system, or stored remotely on a server or other network device, and the medical image is received by loading the medical image to a computer system performing the method of FIG. 3. In another possible implementation, the medical image can be received directly from an image acquisition device, such as an MR scanner, CT scanner, etc.

[0021] At step 304, multiple candidate values for each voxel in the target modality medical image are generated based on the corresponding voxels in the source modality medical image. The target modality medical image will also be referred to herein as the “target image”. Let Φ^v denote the set consisting of voxel v and its neighbors. In an advantageous implementation, the six voxels which are at a unit distance from v are used as neighbors. More neighbors can be added to the set Φ^v without significant changes to method described herein. We use the notation Φ^v(p,q,r) to represent the elements (voxels) of Φ^v. Here, Φ^v(p,q,r) refers to the voxel (v+(p,q,r)). We represent the l₀ and l₂ norms using ∥ ∥₀ and ∥ ∥₂, respectively. The notation v˜v′ is used to indicate that the voxels v and v′ are neighbors. We use II to denote the indicator function, which is equal to 1 when true and equal to 0 when false.

[0022] Given the source image I_s, multiple target modality candidate intensity values are generated for the respective set Φ^v representing the neighborhood at each voxel independently. To generate the target intensity values for Φ^v, a d₁×d₁×d₁ patch centered on v is extracted from the received source modality medical image. If paired source-target images were available during training, it would be possible to learn a predictor/regressor to predict normal target modality candidate voxels for Φ^v from the source modality patch at voxel v. However, since such paired training data may not be available, embodiments of the present invention do not use such a trained predictor/regressor. In an advantageous embodiment of the present invention, a cross-modal nearest neighbor search is used to target modality candidate intensity values. For each d₁×d₁×d₁ image patch for the source image, K nearest d₁×d₁×d₁ target patches are obtained by searching across a set of target modality training images. In a possible implementation, K can be set equal 10, but the present invention is not limited to any particular value of K. The intensity values of the center voxel and the neighboring voxels from these K nearest image patches from the target modality training images provide the target modality candidate intensity values the set Φ^v.

[0023] The cross-modal nearest neighbor search compares each patch of the source image with patches of the target modality training images using a similarity measure. The similarity measure should be robust to changes in modality. In an advantageous embodiment of the present invention, voxel-intensity based mutual information is used as the cross-modality similarity measure. Given two image patches A and B, their mutual information is given by:

MI(A,B)=H(X_a)+H(X_b)−H(X_a,X_b),   (1)

where H denotes the Shannon entropy function, and X_a and X_b are random variables representing the voxel intensities in patches A and B, respectively. The mutual information similarity measure in Equation (1) measures the consistency between the intensity distributions of the image patches in the source and target domains. The Shannon entropy function is a well-known measure of entropy or uncertainty. H(X_a) represents and uncertainty associated with the intensity value X_a occurring in image patch A, H(X_b) represents an uncertainty associated with the intensity value X_b occurring in image patch B, and H(X_a,X_b) represents an uncertainty associated with X_a and X_b occurring as intensity values of corresponding voxels in image patches A and B.

[0024] In order to generate the candidate target modality intensity values for the target image, a plurality of image patches are extracted from the source image. In a possible implementation, a respective patch can be centered at each voxel in the source image can be extracted. In another possible implementation, a predetermined number of voxels can be skipped in each direction between voxels at which the image patches are centered. In another possible implementation, the source image can be divided into a plurality of non-overlapping image patches that cover all of the voxels of the source image. Each image patch extracted from the source image is then compared with a large number of image patches of the target modality training images by calculating the mutual information similarity measure between the image patch of the source image and each image patch of the target modality training images, and K image patches of target modality training images having the highest mutual information similarity measure values are selected for each image patch of the source image. The intensity values for the center voxel. The intensity values of the center voxel and its neighboring voxels in each of the K image patches of the target modality training images selected for a particular source image are assigned to be candidate target modality intensity values for the voxels in the set Φ^v in the target image, where the voxel v in the target image is the voxel in the target image that corresponds to the center voxel of the image patch extracted from the source image (i.e., voxel v is located at the same location in the target image as the center voxel of the image patch is located in the source image). This results in a plurality of candidate target modality intensity values for each voxel in the target image.

[0025] Although the embodiment of the present invention described herein uses mutual information as a similarity measure, the present invention is not limited thereto and other cross-modality similarity measures can be used instead.

[0026] Returning to FIG. 3, at step 306, the full target modality medical image is synthesized by jointly selecting the best candidate values for all of the voxels in the target modality medical image. Given the K candidate target modality intensity values for the set Φ^v centered at each voxel in the target image, the full target modality medical image I_t is synthesized by selecting one of the K candidates at each voxel. The value of Φ^v(0,0,0) (i.e., the intensity value of the center voxel) from the selected candidate is used to synthesize voxel v in I_t.

[0027] Let X_s and X_t be two random variables with support Ψ=(l₁, l₂, . . . , l_L), representing the voxel intensity values of the source and target images I_s and I_t, respectively, where l₁, l₂, . . . , l_L are intensity values sampled from the intensity distribution of the source and target images I_s and I_t. Let I_s(v) and I_t(v) denote the intensity values of voxel v in images I_s and I_t, respectively. Let V represent the set of all voxels with cardinality |V|=N. Let {φ^v1, φ^v2, . . . , φ^vK} denote the K candidate target modality intensity values for the set Φ^v at voxel v. Let W_vk=II[Candidate φ^vk is selected at voxel v]. That is, w_vk equals 1 when the candidate φ^vk is selected at voxel v and equals 0 when the candidate φ^vk is not selected at voxel v. In an advantageous embodiment of the present invention, since the candidates have been obtained for each voxel independently using the nearest neighbor search, a the candidate intensity values are selected jointly for all of the voxels of the target image by solving a selection (optimization) problem based on the following two criteria: (i) Mutual information maximization, which is a global criterion; and (ii) Spatial consistency maximization, which is a local criterion.

[0028] It can be assumed that regions of similar tissue (and hence similar intensity values) in one image would correspond to regions in the other image that also have similar intensity values (though probably different values from the intensity values in the first image). Based on this assumption, mutual information is used as a cost function for cross-modal medical image synthesis. Since we are interested in generating synthesized subject specific scans, the synthesized target image I_t should have high mutual information with the given source image I_s. That is the amount of information I_s and I_t contain about each other should be maximal. This global criterion helps in transferring the image level structure across modalities. The mutual information between images I_s and I_t is given by MI(I_s, I_t)=H(X_s)+H(X_t)−M_s, X_t). Since the entropy H(X_s) is constant for a given source image, maximizing mutual information is equivalent to maximizing H(X_t)−H_s, X_t), where:

$\begin{array}{cc}H\left(X_t\right)=-\sum _{b=1}^LP\left(X_t=l_b\right)\log \left[P\left(X_t=l_b\right)\right],P\left(X_t=l_b\right)=\frac{1}{N}\sum _{v\in V}\sum _{k=1}^Kw_{\mathrm{vk}}\mathrm{II}\left[{\phi }^{\mathrm{vk}}\left(0,0,0\right)=l_b\right],H\left(X_s,X_t\right)=-\sum _{a=1}^L\sum _{b=1}^LP_{\mathrm{ab}}\log \left[P_{\mathrm{ab}}\right],P_{\mathrm{ab}}=P\left(X_s=l_a,X_t=l_b\right)=\frac{1}{N}\sum _{v\in V}\sum _{k=1}^Kw_{\mathrm{vk}}\mathrm{II}\left[I_s\left(v\right)=l_a,{\phi }^{\mathrm{vk}}\left(0,0,0\right)=l_b\right].& \left(2\right)\end{array}$

[0029] Regarding local spatial consistency maximization, let v, v′ ∈ V be two neighboring voxels of the target image. Note that if a candidate φ^vi is selected at voxel v, along with assigning the value φ^vi(0,0,0) to voxel v, the candidate φ^vi can also assign the value φ^vi(v′−v) to the neighboring voxel v′. Similarly, if a candidate φ^v′j is selected at voxel along with assigning the value φ^v′j(0,0,0) to voxel v′, the candidate φ^v′j can also assign the value φ^v′j(v−v′) to the neighboring voxel v. In this case, we would ideally like to have:

φ^vi(0,0,0)=φ^v′j(v−v′), φ^v′j(0,0,0)=φ^vi(v′−v).   (3)

Hence, in an advantageous implementation, to promote spatial consistency among the selected candidate target modality intensity values, the following cost function can be minimized:

$\begin{array}{cc}\mathrm{SC}\left(W\right)=\sum _{\underset{v\sim v^{\prime }}{v,v^{\prime }\in V}}\left[w_{v1}\dots w_{\mathrm{vk}}\right]\left[\begin{array}{ccc}{C}_{11}^{{\mathrm{vv}}^{\prime }}& \dots & C_{1K}^{{\mathrm{vv}}^{\prime }}\\ ⋮& \ddots & ⋮\\ C_{K1}^{{\mathrm{vv}}^{\prime }}& \dots & C_{\mathrm{KK}}^{{\mathrm{vv}}^{\prime }}\end{array}\right]\left[\begin{array}{c}{W}_{v^{\prime }1}\\ ⋮\\ W_{v^{\prime }K}\end{array}\right],\mathrm{where}C_{\mathrm{ij}}^{{\mathrm{vv}}^{\prime }}=\sqrt{{\left({\phi }^{\mathrm{vi}}\left(0,0,0\right)-{\phi }^{v^{\prime }j}\left(v-v^{\prime }\right)\right)}^2+{\left({\phi }^{v^{\prime }j}\left(0,0,0\right)-{\phi }^{\mathrm{vi}}\left(v^{\prime }-v\right)\right)}^2}.& \left(4\right)\end{array}$

[0030] According to an advantageous embodiment of the present invention, the global mutual information cost function and the spatial consistency cost function are combined into a single optimization problem. In particular, the selection of the candidate target modality intensity values for all of the voxels in the target image can be formulated as the following optimization problem:

$\begin{array}{cc}{\mathrm{maximize}}_{\left\{w_{\mathrm{vk}}\right\}}H\left(X_t\right)-H\left(X_s,X_t\right)-\lambda SC\left(W\right)\mathrm{subject}\mathrm{to}\sum _{k=1}^Kw_{\mathrm{vk}}=1,\forall v\in V,w_{\mathrm{vk}}\in \left\{0,1\right\},\mathrm{for}k=1,2,\dots ,K,\forall v\in V,& \left(5\right)\end{array}$

where λ is a trade-off parameter.

[0031] The optimization problem in Equation (5) is combinatorial in nature due to binary integer constraints on w_vk and is difficult to solve. In an advantageous implementation, the binary integer constraints are relaxed to positivity constraints to obtain the following relaxed optimization problem:

$\begin{array}{cc}{\mathrm{maximize}}_{\left\{w_{\mathrm{vk}}\right\}}H\left(X_t\right)-H\left(X_s,X_t\right)-\lambda SC\left(W\right)\mathrm{subject}\mathrm{to}\sum _{k=1}^Kw_{\mathrm{vk}}=1,\forall v\in V,w_{\mathrm{vk}}\ge 0,\mathrm{for}k=1,2,\dots ,K,\forall v\in V,& \left(6\right)\end{array}$

The cost function H(X_t)−H(X_s, X_t)−λSC(W) is differentiable and its derivative with respect to w_vk can be calculated using:

$\begin{array}{cc}\begin{array}{c}\frac{H\left(X_t\right)}{w_{\mathrm{vk}}}=-\sum _{b=1}^L\left(1+\log \left[P\left(X_t=l_b\right)\right]\right)\frac{}{w_{\mathrm{vk}}}P\left(X_t=l_b\right)\\ =-\frac{1}{N}\sum _{b=1}^L\left(1+\log \left[P\left(X_t=l_b\right)\right]\right)\mathrm{II}\left[{\phi }^{\mathrm{vk}}\left(0,0,0\right)=l_b\right]\\ =-\frac{1}{N}(1+\log \left[P\left(X_t={\phi }^{\mathrm{vk}}\left(0,0,0\right)\right]\right),\end{array}\begin{array}{c}\frac{H\left(X_{s,}X_t\right)}{w_{\mathrm{vk}}}=-\sum _{a=1}^L\sum _{b=1}^L\left(1+\log \left[P_{\mathrm{ab}}\right]\right)\frac{P_{\mathrm{ab}}}{w_{\mathrm{vk}}}\\ =-\frac{1}{N}\sum _{a=1}^L\sum _{b=1}^L\left(1+\log \left[P_{\mathrm{ab}}\right]\right)\mathrm{II}\left[I_s\left(v\right)=l_a,{\phi }^{\mathrm{vk}}\left(0,0,0\right)=l_b\right]\\ =-\frac{1}{N}(1+\log \left[P\left(X_s=I_s\left(v\right),X_t={\phi }^{\mathrm{vk}}\left(0,0,0\right)\right]\right),\end{array}\frac{\mathrm{SC}\left(W\right)}{{\mathrm{sw}}_{\mathrm{vk}}}=\sum _{v^{\prime }\sim v}\left(\sum _{p=1}^KC_{\mathrm{kp}}^{{\mathrm{vv}}^{\prime }}w_{v^{\prime }p}\right).& \left(7\right)\end{array}$

[0032] The optimization problem in Equation (6) has a differentiable cost function with linear equality and equity constraints. Accordingly, in an exemplary implementation, this optimization problem can be solved using the reduced gradient ascent approach. Solving the optimization problem provides values for w_vk for each candidate target modality intensity value for each voxel. In an advantageous implementation, the candidate φ^vk* at is selected at each voxel v, where k*=argmax_k w_vk. That is, since the binary constraint for w_vk is relaxed in this optimization problem, the candidate with the maximum value for w_vk is selected for each voxel. In an alternative implementation, a number of top candidates for each voxel having the highest values of w_vk can be averaged to obtain the intensity value for each voxel of the target image.

[0033] Since the cost function in Equation (6) is non-convex, it is not guaranteed to find a global optimum. In an advantageous implementation, the local optimum obtained by initializing all of the variables w_vk with a value of

$\frac{1}{K}$

is used. This initialization can be interpreted as giving equal weight to all K candidates at the beginning of the optimization.

[0034] Returning to FIG. 3, at step 308, the synthesized target modality medical image is refined using coupled sparse representation. Sparse representation is a powerful model for image processing. Sparse representation models a signal as a linear combination of a small number of bases, also known as dictionary atoms: x=Dα, where ∥α∥₀≦T₀. Here, x is the signal, D is the dictionary whose columns are dictionary atoms, α is the sparse coefficient, and ∥α∥₀ is the number of non-zero elements which is constrained to be less than or equal to T₀, which is a predetermined threshold value. This representation can be regarded as a union of subspaces. The dictionary is typically not known a priori. The dictionary D can be learned (trained) to optimally represent the signal of interest using dictionary learning algorithms like K-SVD or method of orthogonal direction (MOD). Given a signal corrupted by different types of noise and a learned dictionary in which the signal can be sparsely represented, the original signal can be recovered using sparse coding techniques by finding the optimal sparse coefficient as α*=argmin_α∥x−Dαμ₂, where ∥αa∥₀≦T₀, wherein the recovered signal is reconstructed from the learned dictionary as Dα*. In practice, sparse coding has been shown to perform well and yield superior reconstruction and noise reduction performance. This is due to the fact that random noises and artifacts are rarely sparse in the dictionary, which was learned to optimally represent the signal of interest.

[0035] Coupled sparse representation finds a coupled sparse code that best reconstructs that both of a pair of signals using respective dictionaries for the signals. To refine the synthesized target image using coupled sparse representation, at each voxel v ∈ V, small d₂×d₂×d₂ image patches are extracted from given source modality image I_s and the synthesized target modality image I_t. Let P_v^s and P_v^t denote the patches centered at voxel v extracted from images I_s and I_t, respectively. Using {((P_v^s, P_v^t)|v ∈ V)} as signal pairs from the source and target modalities, coupled sparse representation can be formulated as the following optimization problem:

$\begin{array}{cc}{\mathrm{minimize}}_{D_s,D_t\left\{{\alpha }_v\right\}}\sum _{v\in V}\left({P_v^s-D_s{\alpha }_v}_2^2+{P_v^t-D_t{\alpha }_v}_2^2\right)\mathrm{subject}\mathrm{to}{{\alpha }_v}_0\le T_0\forall v\in V,& \left(8\right)\end{array}$

where D_s and D_t are over-complete dictionaries with M atoms in the source and target modalities, respectively, α_v is the coupled sparse code for signals P_v^s and P_v^t in their respective dictionaries, and T₀ is the sparsity parameter.

[0036] The dictionaries D_s, D_t and the coupled sparse codes α_v for each voxel are learned directly from the source image and the synthesized target image by the solving the optimization problem shown in Equation (8) using the K-SVD algorithm with explicitly re-normalizing the dictionary atoms of D_s and D_t separately to norm 1 after each iteration. Once the dictionaries and the coupled sparse codes for each voxel are learned from the source image and the synthesized target image, the target modality image patch P_v^t is reconstructed at each voxel using the learned coupled sparse code α_v and the learned target modality dictionary D_t as {circumflex over (P)}_v^t=D_tα_v. The synthesized target image is refined by using the intensity value of the center voxel from {circumflex over (P)}_v^t as the new target modality intensity value for voxel v.

[0037] Returning to FIG. 3, at step 310, the synthesized target modality image is output. For example, the synthesized target medical image can be displayed on a display of a computer system and/or stored in a storage or memory of a computer system.

[0038] The method of FIG. 3 can be used to synthesize medical images across medical imaging modalities, as well as synthesis of images across an image domain. For example, the method of FIG. 3 can be used for applications including but not limited to synthesizing a CT image from an MR image, synthesizing an MR image with one protocol (e.g., T1 and T2) from an MR image with a different protocol, synthesizing contrast CT images from non-contrast CT images, synthesizing a high kV CT image from a low kV CT image, or synthesizing any type of high resolution medical image from a corresponding low resolution medical image. In an exemplary implementation, the present inventors have utilized the method of FIG. 3 to synthesize patient-specific T2-weighted MR brain images from T1-weighted MR brain images and vice versa. FIG. 4 illustrates corresponding T1-weighted MR brain images and T2-weighted MR brain images. As shown in FIG. 4, images 402, 404, and 406 are T1-weighted MR brain images of three different subjects and image 412, 414, and 416 are the respective corresponding T2-weighted MR brain images of the three different subjects.

[0039] In the present inventors experiments for synthesizing T2 MR brain images from T1 MR brain images and synthesizing T1 MR brain images from T2 MR brain images, the training images and the input source images were linearly registered, skull stripped, inhomogeneity corrected, histogram matched within each modality, and resampled to 2 mm resolution. Since the present inventors have access a database containing both T1-weighted and T2-weighted MR scans for a number of subjects, the synthesized images can be directly compared to ground truth target modality images for evaluation of the synthesis method. The present inventors utilized normalized cross-correlation as the evaluation metric.

[0040] Since exhaustively searching the target modality training images to find the nearest neighbors is highly computational, the search region in each target modality training image can be restricted to a h×h×h region surrounding the voxel in the target modality training image corresponding to the voxel of interest at which the image patch is centered in the source image. Since MR scans have a high dynamic range, a mutual information measure calculated using the original voxel intensity values may be highly unreliable. Accordingly, for computing the mutual information, the original voxel intensity values are quantized to L levels. The value trade-off λ in the optimization problem expressed in Equation (6) was selected such that the values of the mutual information and spatial consistency costs are of the same order of magnitude. Exemplary parameters values that were used by the present inventors in the MR brain image synthesis experiments are shown in Table 1.

TABLE 1
Parameter values used in the experiments.
	Parameter	Description	Value
	d1	Size of the image patch used for cross-	9
		modal nearest neighbor search
	h	Size of the search region used for cross-	7
		modal nearest neighbor search
	L	Number of quantization levels used for	32
		computing mutual information
	K	Number of nearest neighbors	10
	d2	Size of the image patch used for sparse	3
		representation
	M	Number of dictionary atoms	500
	T0	Sparsity parameter	5
	λ	Trade-off parameter in the optimization	10−8
		problem 6

[0041] Table 2 shows the normalized cross correlation values between the synthesized target modality images and the ground truth target modality images for 19 subjects. Based on results the MR brain image synthesis performed using the method of FIG. 3, it can be observed that target modality images synthesized using the method of FIG. 3 capture most of the structural information of ground truth target modality images. As shown in Table 2, synthesis of T1 images from T2 images produces more accurate results as compared to the synthesis of T2 images from T1 images. In addition, the results in Table 2 also show that coupled sparse representation (CSR) improves the synthesis results, in particular when synthesizing T2 images from T1 images.

TABLE 2
Normalized cross correlation values between synthe
and ground truth target modality images.
	Source: T1-MRI		Source: T2-MRI
	Target: T2-MRI		Target: T1-MRI
Subject	No CSR	CSR	No CSR	CSR
1	0.862	0.878	0.932	0.936
2	0.839	0.859	0.932	0.936
3	0.881	0.894	0.934	0.938
4	0.841	0.853	0.933	0.937
5	0.814	0.832	0.873	0.870
6	0.841	0.863	0.939	0.943
7	0.792	0.811	0.900	0.905
8	0.833	0.846	0.941	0.944
9	0.856	0.875	0.933	0.937
10	0.848	0.863	0.936	0.941
11	0.871	0.886	0.935	0.939
12	0.822	0.837	0.925	0.930
13	0.838	0.853	0.926	0.931
14	0.861	0.872	0.940	0.943
15	0.791	0.810	0.915	0.921
16	0.830	0.845	0.936	0.940
17	0.851	0.868	0.929	0.934
18	0.859	0.874	0.923	0.928
19	0.811	0.826	0.924	0.928
average	0.839	0.855	0.927	0.931

[0042] Embodiments of the present invention provide an unsupervised approach for cross-modal synthesis of subject specific medical images. The method described herein can be used with any pair of imaging modalities, and works without paired training data from the source and target modalities, thereby alleviating the need for scanning each subject multiple times. Given a source modality image, multiple candidate target modality intensity values are generated for each voxel location of a target modality image independently based on image patches extracted from the source modality image using a cross-modal nearest neighbor search. In an exemplary embodiment, voxel-intensity based mutual information is used as a similarity measure for the cross-modal nearest neighbor search, but the present invention is not limited thereto and other cross-modal similarity measures may be used as well. The best candidate values for all the voxels of the target modality image are selected jointly by simultaneously maximizing a global mutual interest cost and minimizing a local spatial consistency cost. Coupled sparse representation is used for further refinement of the synthesized target modality image.

[0043] The above-described method for unsupervised cross-modal synthesis of medical images may be implemented on a computer using well-known computer processors, memory units, storage devices, computer software, and other components. A high-level block diagram of such a computer is illustrated in FIG. 5. Computer 502 contains a processor 504, which controls the overall operation of the computer 502 by executing computer program instructions which define such operation. The computer program instructions may be stored in a storage device 512 (e.g., magnetic disk) and loaded into memory 510 when execution of the computer program instructions is desired. Thus, the steps of the method of FIG. 3 may be defined by the computer program instructions stored in the memory 510 and/or storage 512 and controlled by the processor 504 executing the computer program instructions. An image acquisition device 520, such as an MR scanning device or a CT scanning device, can be connected to the computer 502 to input image data to the computer 502. It is possible to implement the image acquisition device 520 and the computer 502 as one device. It is also possible that the image acquisition device 520 and the computer 502 communicate wirelessly through a network. In a possible embodiment, the computer 502 can be located remotely with respect to the image acquisition device 520 and the method steps described herein can be performed as part of a server or cloud based service. In this case, the method steps may be performed on a single computer or distributed between multiple networked computers. The computer 502 also includes one or more network interfaces 506 for communicating with other devices via a network. The computer 502 also includes other input/output devices 508 that enable user interaction with the computer 502 (e.g., display, keyboard, mouse, speakers, buttons, etc.). Such input/output devices 508 may be used in conjunction with a set of computer programs as an annotation tool to annotate volumes received from the image acquisition device 520. One skilled in the art will recognize that an implementation of an actual computer could contain other components as well, and that FIG. 5 is a high level representation of some of the components of such a computer for illustrative purposes.

[0044] The foregoing Detailed Description is to be understood as being in every respect illustrative and exemplary, but not restrictive, and the scope of the invention disclosed herein is not to be determined from the Detailed Description, but rather from the claims as interpreted according to the full breadth permitted by the patent laws. It is to be understood that the embodiments shown and described herein are only illustrative of the principles of the present invention and that various modifications may be implemented by those skilled in the art without departing from the scope and spirit of the invention. Those skilled in the art could implement various other feature combinations without departing from the scope and spirit of the invention.