Abstract:
The dataset describing an entity in a first modality and of a first, high resolution is used to enhance the resolution of a dataset describing the same entity in a second modality of a lower resolution. The two data sets of different modalities are spatially registered to each other. From this information, a joint histogram of the values in the two datasets is computed to provide a raw analysis of how the intensities in the first dataset correspond to intensities in the second dataset. This is converted into a joint probability of possible intensities for the missing pixels in the low resolution dataset as a function of the intensities of the corresponding pixels in the high-resolution dataset to provide a very rough estimate of the intensities of the missing pixels in the low resolution dataset. Then, an objective function is defined over the set of possible new values that gives preference to datasets consistent with (1) the joint probability distributions, (2) the existing values in the low resolution dataset, and (3) smoothness throughout the data set. Finally, an annealing or similar iterative method is used to minimize the objective function and find an optimal solution over the entire dataset.

Description:
FIELD OF THE INVENTION  
       [0001]     The invention pertains to computer science and, particularly, digital imaging. More particularly, the invention pertains to propagating information between data sets of two different modalities.  
       BACKGROUND OF THE INVENTION  
       [0002]     Often, two or more different data sets of information about an entity are available in two or more different modalities. A classic example of this occurs frequently in the field of medical imaging. For instance, an MRI (magnetic resonance imaging) scan of an anatomical feature of a patient as well as a CT (computer tomography) scan of the same feature may be available. One of these images, most likely the MRI scan in this example, may have a higher resolution than the other image. Nevertheless, the CT image may show certain features with much better clarity than the MRI image (e.g., bone), whereas the MRI image may show a different feature (e.g., soft tissue) with much better clarity than the CT image. Often, this is exactly why two images of the same object are obtained in different modalities.  
         [0003]     The image(s) in one of the modalities may have a higher resolution that image(s) in the other modality.  
         [0004]     Another example of a common situation in which two images (or datasets) are available describing the same entity occurs in geographic mapping. Often a relatively low resolution terrain map of a geographic area is available as well as a higher resolution aerial photograph of the same geographic area. For example, ten-meter resolution terrain elevation data and corresponding one-meter resolution orthophotos are available online for most of the United States from the USGS&#39;s Seamless Data Distribution System.  
         [0005]     As yet another example, visible and thermal IR images of the same scene may have different effective resolutions because the sensors have different resolutions or because one camera is focused on a narrow part of the scene.  
         [0006]     Several techniques are known for computer enhancing the resolution of a given image in a given modality in order to improve the clarity of the image. One of the simpler techniques is interpolation or extrapolation. In such a technique, typically some form of averaging of the values (e.g., intensities) of adjacent pixels is done in order to fill in (or add) pixels between those adjacent pixels. In image processing, a standard interpolation technique is bilinear interpolation. This technique uses a weighted average of the four existing pixel values to the upper left, lower left, upper right, and lower right of the location that is to be filled in. One way to think of interpolation is to fill in values by assuming smoothness.  
         [0007]     Another technique for enhancing resolution of an image is superresolution. Usually in superresolution, multiple images of the same object in the same modality are obtained. The images might be obtained from slightly or vastly different perspectives. The images are registered with respect to each other, i.e., the translation and location parameters that properly align overlapping views of an object are determined so as to reconstruct from these partial views an integrated overall representation of the object. Then, the sub-pixel motions between the images estimated by the registration are used to construct a single image that has a higher resolution than any of the single images.  
         [0008]     This process has been used with multiple images from the Hubble telescope to increase the telescope&#39;s effective resolution. A superresolution approach that uses multiple images of terrain to construct a three-dimensional surface with higher resolution than would be possible using only stereo vision with two images is described in P. Cheeseman, B. Kanefsky, R. Kraft, J. Stutz, and R. Hanson. Super-resolved surface reconstruction from multiple images. Technical Report FIA-94-12, NASA Ames Research Center, December 1994.  
         [0009]     Another technique is texture synthesis. In texture synthesis, a statistical model of an image texture is created from existing images that demonstrate the texture. Then, that texture model is used in a new image to fill holes with a similar, but unique, instance of that texture. For example, if it is necessary to remove a person from a field of grass in an image, one might model the grass texture from the image and then write over the person in the image with newly synthesized “grass” that is visually consistent with the other grass in the image. See, e.g., S. Zhu, Y. Wu, and D. B. Mumford. FRAME: Filters, random field and maximum entropy: towards a unified theory for texture modeling. International Journal of Computer Vision, 27(2):1-20, March/April 1998.  
         [0010]     Surface hole filling is another technique that can be used in connection with 3-D imaging in which holes in a 3-D surface are filled in with a surface patch that is consistent with the hole border&#39;s geometry. See, e.g., J. Davis, S. E. Marschner, M. Garr, and M. Levoy. Filling holes in complex surfaces using volumetric diffusion. In  Proceedings of the First International Symposium on  3 D Data Processing, Visualization, and Transmission , June 2002.  
         [0011]     It also is known to register images from two different modalities with each other. For instance, this is known in the field of surgical navigation in which real-time infrared images of the surgical field containing a patient are obtained and registered to previously obtained images, such as CT images of the patient. This, for instance, is often used to help guide a surgical instrument in real-time towards an anatomical feature, such as a tumor.  
         [0012]     Furthermore, a technique known as shape-from-shading is a process of recovering 3-D shape information from a single image using a known mapping from image intensities to possible object surface slopes. Shape-from-shading was originally developed by B. K. P. Horn and colleagues and is described in his textbook (B. K. P. Horn, Robot Vision, MIT Press, Cambridge, Mass., 1991). While shape-from-shading does infer a data set of one modality (depths) from a data set of a second modality, it does not exploit values of the second modality that do already exist, such as sparse values in a low-resolution image.  
         [0013]     None of the aforementioned image processing techniques involves enhancing the resolution of an image in one modality based on a data set in another modality.  
         [0014]     Accordingly, is an object of the present invention to provide a new and improved image processing technique.  
         [0015]     It is another object of the present invention to provide a method and apparatus for enhancing the resolution of a data set in a first modality based on information in a data set in a second modality.  
       SUMMARY OF THE INVENTION  
       [0016]     The dataset describing an entity in a first modality and of a first, high resolution is used to enhance the resolution of a dataset describing the same entity in a second modality of a lower resolution.  
         [0017]     In accordance with the concepts of the present invention, the two data sets of different modalities are spatially registered to each other, for example, using maximization of mutual information. From this information, a joint histogram of the values in the two datasets is computed. This provides a raw analysis of how the intensities in the first dataset correspond to intensities in the second dataset. This is converted into a joint probability of possible intensities for the missing pixels in the low resolution dataset (i.e., the pixels that will be inserted between the existing pixels in order to increase the resolution) as a function of the intensities of the corresponding pixels in the high-resolution dataset. This provides a very rough estimate of the intensities of the “missing” pixels in the low resolution dataset. Then, an objective function is defined over the set of possible new values that gives preference to datasets consistent with (1) the joint probabilities, (2) the existing values in the low resolution dataset, and (3) smoothness throughout the data set. Finally, an annealing or similar iterative method is used to minimize the objective function and find an optimal solution over the entire dataset.  
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0018]      FIGS. 1A-1F  are diagrams illustrating the propagation of detail from a data set of one modality to a data set of another modality using the principles of the present invention.  
         [0019]      FIGS. 2A-2H  are a series of diagrams illustrating the recovery of terrain normals from a terrain image and subsampled normals using the principles of the present invention.  
         [0020]      FIG. 3  is a flow diagram illustrating the basic steps of a process in accordance with one particular embodiment of the present invention.  
     
    
     DETAILED DESCRIPTION OF THE INVENTION  
       [0021]     The invention concerns a method and apparatus for propagating higher-resolution detail from a dataset of one modality to a spatially registered, lower-resolution (subsampled) dataset of a second modality. A Gibbs joint posterior distribution for the missing values in the low-resolution dataset is defined, incorporating (1) an MRF (Markov Random Fields) smoothness prior and (2) likelihoods for the observed high-resolution values derived from empirical joint statistics relating the two modalities. This posterior is maximized using Gibbs sampling, producing a maximum a posteriori (MAP) estimate for the missing values.  
         [0022]     The goal of the present invention is to improve the resolution or coverage of an existing dataset of a first modality using a more complete (e.g., higher resolution) dataset of a second modality. The technique exploits both the data that exists in the lower-quality dataset and the statistical information provided by the values in the higher-quality dataset.  
         [0023]     Our algorithm estimates missing values in a lower-resolution dataset, given a spatially registered, higher-resolution dataset of a different modality.  FIGS. 1A-1F  illustrate this process. In this example, a full-resolution T1 MRI dataset shown in  FIG. 1A  and a low-resolution (sampled) T2 dataset shown in  FIG. 1B  are supplied. Using a smoothness prior and likelihoods for the observed T1 values given the missing T2 values, maximum a posteriori (MAP) estimates for missing values are found using Gibbs sampling. Three different estimates generated by the algorithm assuming different kinds and levels of smoothness are shown in  FIGS. 1C, 1D , and  1 E. The ground truth, full-resolution T2 image is shown in  FIG. 1F .  
         [0000]     Datasets  
         [0024]     In this specification, the terms “dataset” and “image” are used interchangeably, although the datasets are not always optical intensity images. In addition, the following definitions will apply.  
         [0025]     The image f is the enhanced resolution image of modality  1  that is to be generated. It is generated by observing an object or scene with modality  1 . The size of f is r rows by c columns, and its pixels take on discrete values. Most of f is unknown, but we observe part of f as the degraded image f′ described below.  
         [0026]     The image g is the existing higher resolution image in modality  2 . It is generated by observing the same object or scene with modality  2 . Like f, g has r rows and c columns and takes on discrete values. Image g is known and spatially registered with f.  
         [0027]     A subsampled image f′ is the existing lower resolution image in modality  1  that will be resolution enhanced in accordance with the principles of the present invention. It gives values of f at a few points. In this example, the samples are at every fifth grid point, (1, 1), (1, 6), . . . (6, 1), (6, 6), . . . , so that we have approximately 4% of the information of f in f′. f′ has size r×c but is only defined as the known sample points.  
         [0028]     s 1 , . . . , s r×c  specify the two-dimensional (row, column) pixel locations, or “sites” or “datapoints”, in these images. These subsets s d1 , s d2 , . . . and s u1 , s u2 , . . . give the sites where f′ is defined and undefined, respectively.  
         [0000]     Gibbs Distribution and MRFs  
         [0029]     In subsequent sections, we will perform maximum a posteriori (MAP) estimation of the unknown values of f, by finding the mode of a posterior distribution incorporating both a smoothness prior for f and our observations g. The posterior and prior distributions will both be Gibbs distributions, or equivalently, Markov random fields (MRFs). We present a brief summary of Gibbs distributions and MRFs that closely follows that in S. Geman and D. Geman. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  IEEE Transactions on Pattern Analysis and Machine Intelligence , PAMI-6(6):721-741, November 1984, incorporated herein by reference, but is specialized to the present problem.  
         [0030]     The Gibbs distributions have the form  
               π   ⁡     (   f   )       =       1   Z     ⁢     ⅇ       -     U   ⁡     (   f   )         /   T                 (   1   )             
 
 where Z is a normalizing constant that brings the volume under π to 1,  
             Z   =       ∑   f     ⁢     ⅇ       -     U   ⁡     (   f   )         /   T                 (   2   )             
 
 and the “temperature” T is a constant with respect to f; T&#39;s effect is explained in the subsection entitled Gibbs sampling below. The energy function U has the form  
               U   ⁡     (   f   )       =       ∑     C   ∈   χ       ⁢       V   C     ⁡     (   f   )                 (   3   )             
 
 where, for our problem, X is a set of sets C including: 
 
 singleton sets C={s u     i   } representing a single site where the value of f is unknown, and 
 
 pair sets C={s i ,s j }, each representing a pair of eight-neighboring sites V C (f) assigns an energy to a singleton set C={s u     i   } or pair set C={s i ,s j } based on the assignments in f . The V C (f) for our prior and posterior are described in the subsections below. 
 
         [0031]     The local characteristics of the distribution are the conditional distributions for the value f (S i ) at a single site given the values at all other sites. For the distribution (1), the local characteristics are:  
               P   ⁡     (         f   ⁡     (     s   i     )       |     f   ⁡     (     s   j     )         ,       u   i     ≠   j       )       =       Z     s   i       -   1       ⁢   exp   ⁢     {       -     1   T       ⁢       ∑       C   ⁢     :     ⁢           ⁢     s   i       ∈   C       ⁢       V   C     ⁡     (   f   )           }               (   4   )             
 
 where Z s     i    is again a normalizing constant that brings the volume under the distribution to one. Taken together, the local characteristics constitute the MRF equivalent to the original Gibbs distribution. 
 
 Prior Distribution 
 
         [0032]     To derive a smoothness prior on f, we assign a joint distribution to pairs of values at eight-neighbors based on the approximate spatial derivatives of f. Specifically, if s i  and s j  are eight-neighbors, then  
               Δ     i   ,   j       =         f   ⁡     (     s   j     )       -     f   ⁡     (     s   i     )                  s   j     -     s   i                      (   5   )                   P     i   ,   j       ⁡     (   f   )       =         N     0   ,     σ   s         ⁡     (     Δ     i   ,   j       )       =       1     2   ⁢     πσ   s   2         ⁢   exp   ⁢     {       -     Δ     i   ,   j     2         2   ⁢     σ   s   2         }                 (   6   )             
 
 where N 0,σ , is the normal distribution with zero mean and standard deviation σ s . We assume that the P i,j  for each pair of neighboring sites are independent, including pairs that include a common site. So,  
               P   ⁡     (   f   )       =       ∏       {       s   i     ,     s   j       }     ∈   χ       ⁢       P     i   ,   j       ⁡     (   f   )                 (   7   )             
 
         [0033]     If we write  
                 V     {       s   i     ,     s   j       }     prior     ⁡     (   f   )       =       -     log   ⁡     (     1     2   ⁢     πσ   s   2         )         +       Δ     i   ,   j     2       2   ⁢     σ   s   2                   (   8   )                   U   prior     ⁡     (   f   )       =       ∑       {       s   i     ,     s   j       }     ∈   χ       ⁢       V     {       s   i     ,     s   j       }     prior     ⁡     (   f   )                 (   9   )             
 
 then substituting (9) for U(f) in (1), (2) with T=1 produces (7). That is, the prior distribution (7) is a Gibbs distribution with the energy function given by (9). 
 
         [0034]     This smoothness prior has one parameter, σ s , which will be called the spatial smoothness parameter hereinafter.  
         [0035]     This smoothness prior is a relatively blunt instrument, favoring values f(s u     i   ) close to f(s j ) for nearby s u     i   , s j , without regard for the corresponding values of g. A second smoothness assumption, which will be called weak smoothness hereinafter, is that f(s u     i   ) is close to f(s j ) if s u     i   is near s j  and g(s u     i   ) is similar to g(s j ). Weak smoothness can be deployed in those situations where the output dataset is not expected to exhibit much pixel-to-pixel spatial smoothness. In the next subsection, a weak smoothness assumption is introduced into the estimates via the observation likelihood function.  
         [0000]     Likelihoods  
         [0036]     To derive a likelihood P(g|f) for the observations g, let us assume that the likelihoods of the individual g values given the corresponding f values are independent, so that:  
               P   ⁡     (     g   |   f     )       =       ∏     s     u   i         ⁢     P   ⁡     (       g   ⁡     (     s     u   i       )       |     f   ⁡     (     s     u   i       )         )                 (   10   )             
 
         [0037]     If we then write  
                 V     {     s     u   i       }     like     ⁡     (   f   )       =       -   log     ⁢           ⁢     P   ⁡     (       g   ⁡     (     s     u   i       )       |     f   ⁡     (     s     u   i       )         )                 (   11   )                   U   like     ⁡     (   f   )       =       ∑       {     s     u   i       }     ∈   χ       ⁢       V     {     s     u   i       }     like     ⁡     (   f   )                 (   12   )             
 
 then substituting (12) for U(f) in (1), (2) produces (10). This relationship will be used below in showing that our posterior distribution is a Gibbs distribution. 
 
         [0038]     The likelihood functions P(g(s u     i   )|f(s u     i   )) for the individual sites s u     i    are constructed by extracting them from joint distributions P(g(s u     i   ),f(s u     i   )). Two ways of choosing the joint probability distributions are explored.  
         [0039]     One possible way is to use a single global joint distribution at every site s u     i   . This global distribution might be generated empirically by taking the joint histogram of known values in f and the corresponding values in g, or by taking the joint histogram of corresponding values in a larger ensemble of registered modality  1  and modality  2  example images. Alternately, the global distribution might be found analytically using a known model of the modality  1  and modality  2  observation processes. For instance, the intensity values of certain tissue types in MRIs and CT scans, respectively, may be well known such that it is not necessary to generate a histogram based on the actual images in question.  
         [0040]     A second way is to use a different joint distribution for every site s u     i   . For example, the joint distribution for s u     i   , can be chosen by augmenting the global distribution with an additional probability mass around the known values of f near s u     i   and the corresponding values of g. This is a type of weak smoothness, as described in the previous subsection.  
         [0041]     In experiments, a global joint distribution was generated empirically, and then generated a local joint distribution for each s u     i   by adding Gaussian-shaped masses to the joint distribution near the values of the known samples in f′ near s u     i   and the corresponding values in g. These masses are specified by a shape covariance matrix C W  and weight m w .  
         [0000]     Posterior Distribution  
         [0042]     Using Bayes&#39; theorem, the posterior distribution P(f|g) that we will maximize can be written as  
               P   ⁡     (     f   |   g     )       =         P   ⁡     (     g   |   f     )       ⁢     P   ⁡     (   f   )           P   ⁡     (   g   )                 (   13   )             
 
 where P(g|f) and P(f) are the likelihood (10) and prior (7), respectively. Combining the energies U prior  and U like  for the prior and likelihood distributions, we have:
 
 U   post   =U   prior   +U   like   (14)
 
 which substituted into (1), (2) produces the distribution (1). That is, the posterior is a Gibbs distribution with energy U post . So, the joint posterior can be maximized using Gibbs sampling as described below. 
 
 Gibbs Sampling 
 
         [0043]     Geman and Geman describe a stochastic relaxation algorithm called the Gibbs Sampler for finding the modes of a Gibbs distribution. The Gibbs Sampler&#39;s operation when applied to the posterior distribution in accordance with the present invention is as follows.  
         [0044]     On each iteration k, the estimate is updated for each unknown pixel value f(s u     i   ), sequentially and in raster scan order. (That is, each iteration k includes a full raster scan). The new estimate for f(s u     i   ) is found by randomly sampling from the local characteristics (4) of the posterior for site s u     i   . This operation is the same as the Gibbs Sampler operation described by Geman and Geman for image restoration, except that they perform an update at every site (i.e., pixel location) s i , whereas, in accordance with the present invention, updates are performed only at the sites s u     i    where the value of f is unknown.  
         [0045]     On each iteration k, the Gibbs distribution temperature T is reduced according to a preset annealing schedule. During early iterations, the high temperature (T&gt;1.0) flattens the distribution, reducing the prominence of the modes. Sampling from this high temperature distribution encourages site updates that lower the probability of the joint estimate f (as measured by the distribution for T=1.0) and allows the algorithm to explore widely separated parts of the solution space. During later iterations, after the temperature has been gradually reduced to below T=1.0, the lower temperatures increase the prominence of the modes that have been discovered. Sampling from these lower temperature distributions increasingly focuses the search in the high probability regions of the T=1.0 distribution.  
         [0046]     Geman and Geman prove that for a sufficiently slow annealing schedule that they describe, this procedure “forces the system into the lowest energy states,” but also point out that this schedule is too slow to be used in practice. Instead, the authors suggested adopting a faster annealing schedule with the same parametric form. In particular,  
               T   ⁡     (   k   )       =     C     log   ⁡     (     1   +   k     )                 (   15   )             
 
 where T(k) is the temperature during all of iteration (i.e., raster scan) k. The authors used C=3 and C=4 in their experiments. 
 
         [0047]     For the present application, a similar form for the annealing schedule may be used that can be easily fit to a specified number of iterations k final , initial temperature T initial , and final temperature T final :  
               T   ⁡     (   k   )       =     C     log   ⁡     (     1   +   ɛ   +     δ   ⁡     (     k   -   1     )         )                 (   16   )             
 
         [0048]     Choosing:  
             ɛ   =     10     -   6               (   17   )               C   =       T   initial     ⁢     log   ⁡     (     1   +   ɛ     )                 (   18   )               δ   =         ⅇ     C   /     T   final         -     (     1.0   +   ɛ     )           k   final     -   1               (   19   )             
 
 fits (16) to T(1)=T initial , T(k final )=T final  and produces an appropriate (e.g., smooth) temperature curve for a wide range of k final , T initial , and T final . 
 
         [0049]     In experiments, the initial temperature has been chosen high enough that the sample image resulting from the first raster scan is a random, salt-and-pepper image. Typically, as σ s  decreases, T initial  must increase to achieve this.  
         [0000]     Estimating Coupled Datasets  
         [0050]     The previous subsections describe a Gibbs sampler for inferring missing values in a low-resolution, subsampled dataset f with the help of a spatially registered, high-resolution dataset g of a different modality.  FIG. 1  shows an example of such an f, g pair, and the performance of the above-described algorithm on this dataset pair is described in the Medical Images section below.  
         [0051]     In an alternate embodiment, the algorithm may be generalized slightly to simultaneously estimate missing values in two coupled datasets f 1 , f 2 , with the help of a third dataset g of a different modality. In the section below entitled Terrain, this generalized algorithm&#39;s performance in simultaneously estimating x and Y components of terrain surface normals from a simple image of the terrain is investigated. In this example, the output x and Y surface normals at each point should not be estimated independently, but should be consistent with each other as well as the observed image intensity at that point.  
         [0052]     The generalized algorithm for estimating coupled datasets is a natural extension of the base algorithm. The necessary extensions are: 
    1. The joint probability function used to compute the observation likelihood is now three-dimensional, describing the co-occurrence statistics of f 1 , f 2 , and g. Similarly, the weak smoothness mass is three-dimensional.     2. During each raster scan of the simulated annealing process, values of both f 1 , and f 2  at each of the missing data sites s u     i    are generated by sampling their respective local characteristics. As an ordering, let us choose to generate a sample value for f 2  at s u     i    immediately after generating the corresponding sample for f 1 (s u     i   ). However, because the likelihood for g(s u     i   ) depends on both f 1 (s u     i   ) and f 2 (s u     i   ), the local characteristics for f 2 (s u     i   ) that we sample from depend on the current sample for f 1 (s u     i   ), and vice-versa. 
 
 Results 
   
 
         [0055]     Medical Images  
         [0056]     This subsection summarizes the algorithm&#39;s performance in estimating full-resolution T2 MRI images from subsampled T2 and full-resolution T1 images. These experiments use simulated datasets from the McConnell Brain Imaging Centre&#39;s online database. These datasets are generated using realistic models of the brain and MRI imaging process, are available with a range of systematic error and random noise levels, and are spatially registered.  
         [0057]     The datasets selected for the experiments have size 181×217×181 voxels, with individual voxels having size 1 mm 3  and taking on greyscale values in the range [0, 255]. In each of the experiments described below, saggital image  108  was arbitrarily selected for processing from the 181 saggital images available. As described below, the datasets selected include different levels of random, pixelwise intensity noise and different levels of intensity non-uniformity (INU), which is a spatially varying error in the intensities.  
         [0058]     As described above, the algorithm&#39;s behavior is determined by several parameters, including the spatial smoothness parameter σ s , the weak smoothness mass shape C w  and weight m w , the choice of joint probability function, the initial annealing temperature T initial , final annealing temperature T final , and number of annealing iterations k final . Together, the dataset and algorithm parameters describe a large space of possible scenarios. Below, two sets of experiments that each explore a part of this space are summarize.  
         [0059]     In the first set of experiments, the algorithm&#39;s performance is investigated as a function of (1) the INU and random noise levels, and (2) the adopted joint probability function. The INU and random noise levels for the four datasets considered are shown in Table 7. In brief, the four datasets cover all four combinations of {INU, no INU} and {random noise, no random noise}, with the noise and INU levels chosen to be the McConnell database defaults. The joint probability functions considered are each generated by applying the Parzen window technique to one of the following joint histograms: 
    1. An “ideal image” joint histogram, generated from the given T1 full-resolution image and the ground truth, full-resolution T2 image. Of course, the ground truth T2 image would not be available in a normal application of the algorithm.     2. A “subsampled image” joint histogram, generated from the given full-resolution T1 image and the given subsampled T2 image. Both of these images would be available in a normal application of the algorithm.     3. A “subsampled volume” joint histogram, generated from the given full-resolution T1 volume and the given subsampled T2 volume. Both of these volumes would be available in a normal application of the algorithm to images extracted from medical volumes.    
 
         [0063]     In each experiment, the spatial smoothness parameter σ s  is chosen to be the actual standard deviation σ s  of the approximate spatial derivatives in the ground truth, full-resolution T2 image. The weak smoothness mass has weight m w =50 in each experiment, and a shape C w  chosen based on the random noise levels in the T1 and T2 images; these values are summarized in Table 1. In each experiment in this set, an annealing schedule with T initial =15, T final =0.075, and k final =250 was used.  
         [0064]     Table 2 gives the average absolute difference in greylevel for the pixels in the estimated, full-resolution T2 image and ground truth T2 image. In this experiment, it appears that (1) random image noise has a larger adverse effect on the algorithm&#39;s accuracy than spatially varying errors, and (2) the subsampled image probability function, which is the probability function expected in practice, provides almost the same accuracy as an ideal probability function generated using the full (and usually unknown) T2 image.  
         [0065]     In the second suite of experiments, the effects of the spatial smoothness parameter σ s  and the weight m w  of the weak smoothness mass are considered. Specifically, we consider σ s  ε{σ s /4,σ s /2,σ s ,2σ s ,4σ s }, and m w  ε{0,25,50}. In each case we have used dataset  4 , which includes both random noise and INU, as described in Table 1. Accordingly, the same weak smoothness mass shape C w  was used as that described for dataset  4  in the first set of experiments.  
         [0066]     As described in the section above entitled Gibbs Sampling, the initial annealing temperature T initial  required to produce a random, salt and pepper initial estimate for the annealing differs depending on the level of smoothness assumed in the MRF prior. More smoothness requires a higher starting temperature. So, the annealing schedule is manually adjusted for different choices of σ s . For σ s =σ s /4, T initial =120.0, T final =0.6; for σ s =σ s /2, T initial =45.0, T final =0.225; and for σ s ≧σ s , T initial =15.0, T final =0.075 are used. In each case, k final =250 iterations.  
         [0067]     The resulting errors are summarized in Table 3. While the errors in Table 3 might suggest a relatively small change in the accuracy as σ s . and m w  vary, the visual differences between these results are significant. In particular, choosing σ s  far smaller than the ideal value σ s  results in images that do not reflect the fine detail in the ground truth image, and tuning m w  can sometimes remove qualitative errors in the output image.  
         [0068]     One would not expect or want the algorithm to estimate random noise in the full-resolution output images. So, for each experiment, σ s  and the average absolute errors are computed using a ground truth image without random noise.  
         [0069]     Terrain  
         [0070]     This subsection describes initial experiments in estimating full-resolution x and Y components of terrain surface normals from a simple, full-resolution image of the surface and subsampled normals. Using the algorithm for coupled datasets described above, once the full-resolution x and Y components of a surface have been recovered, each full normal vector can be reconstructed by appending the z component value that brings the total vector magnitude to 1.0. A full-resolution terrain elevation map could then be estimated by integrating these reconstructed normals.  
         [0071]     An elevation dataset representing a small patch of New Mexico was taken from the ⅓ arc second (10 meter) resolution USGS National Elevation Dataset (NED), available online via the USGS&#39;s Seamless Data Distribution. These elevation values were numerically differentiated to generate ground truth n(s i )=[n x (s i ) n y (s i ) n z (s i )] T  unit surface normals for each site s i , and an image I(s i ) of the surface was generated from these normals using I(s i )=n(s i )·I, where I is the unit lighting direction. In rendering terms, this simple image ignores shadowing, and assumes Lambertian reflection, constant surface albedo, and a distant light source. The resulting n x , n y , and I were then independently linearly transformed with transformations t x , t y , and t I  to fill the range [0,149] and discretized. As with the medical images described in subsection 4.1 above, subsampled normal images n x , and n y  were generated by selecting every fifth pixel from n x  and n y , as described in subsection 3.2. So, the subsampled normal images contain about 4% of the full-resolution normal data.  
         [0072]     The full-resolution normal images n x  and n y  were estimated by applying the method for coupled datasets, described above to I, n x , and n y . For this experiment, the lighting direction is [0 0 1], so that the image I reduces to a map of the surface normal z components. The spatial smoothness parameter was chosen to be σ s =21.56, which is the average of actual smoothness values σ s  calculated from x and Y ground truth normal images. After some experimentation, the shape C w  and weight m w  of the weak smoothness mass was chosen to be 40 2 ·I 3  and 200, and the annealing start temperature, final temperature, and number of iterations to be T initial =30.0, T final =0.2, and k final =250.  
         [0073]     For use in computing the likelihoods, a joint probability function was computed by applying the Parzen window technique to an ideal histogram computed using the ground truth x normal image n x , Y normal image n y , and terrain image I. This histogram is similar to the “ideal image” histogram described for MRI images above in the subsection entitled Medical Images, but has three dimensions rather than two.  
         [0074]     The input images and resulting outputs are shown graphically in  FIGS. 2A-2H . The average absolute errors in estimated x and Y normal components, computed from the [0, 149] greyscale values, are 5.68 and 6.93, respectively. When these values are mapped back to unit normal vectors via t x   −1  and t y   −1 , the average absolute error in angle between the estimated normal vector and the ground truth normal vectors is 6.17 degrees.  
         [0075]     This specification has described a Gibbs sampler for estimating high-resolution datasets from a lower-resolution (subsampled) dataset of a first modality and a spatially registered, higher-resolution dataset of a second modality. Initial experiments indicate that while the resulting estimated datasets are not exact, they do correctly capture most of the high-resolution detail in the corresponding ground truth data.  
         [0000]     Flowchart  
         [0076]     Having thus described the principles of the present in invention,  FIG. 3  illustrates basic steps in accordance with one particular embodiment. A lower resolution image of a first modality  301  and a higher resolution image of a second modality  302  comprise the basic input datasets. In step  305 , the two images are spatially registered with each other.  
         [0077]     Registration can be obtained by using a registration-by-mutual-information technique. Registration-by-mutual-information techniques are known for datasets of the same resolution. Such techniques can be adapted to accommodate datasets of different resolutions. Furthermore, a registration technique might be used in which registration is performed iteratively and/or simultaneously with the filling in of the missing pixels. Alternately, a simple initial estimate for the registration in such an iterative and/or simultaneous process could be found by interpolating between existing values in the low-resolution datasets using bilinear interpolation, and then performing a mutual information registration between the high-resolution dataset and the interpolated dataset.  
         [0078]     Next, in step  307 , a joint histogram is generated.  
         [0079]     In step  309 , the global joint probability is developed using, for example, a Parzen window technique.  
         [0080]     Finally, in step  311 , an enhanced resolution image of the first modality is generated by iteratively maximizing the posterior distribution function P(f|g).  
         [0081]     The invention has been described hereinabove with respect to certain examples in which the datasets are images, e.g., photographs, CT scans, MRIs. However, the concepts of the present invention should be generally applicable to an y  two or more datasets of different modalities describing the same entity.  
         [0082]     Having thus described a few particular embodiments of the invention, various alterations, modifications, and improvements will readily occur to those skilled in the art. Such alterations, modifications and improvements as are made obvious by this disclosure are intended to be part of this description though not expressly stated herein, and are intended to be within the spirit and scope of the invention. Accordingly, the foregoing description is by way of example only, and not limiting. The invention is limited only as defined in the following claims and equivalents thereto.