1. A computer implemented method for registering digitized images using Markov random Fields, said method executed by a computer comprising the steps of: providing a source image f and a target image g, each image comprising a plurality of intensities associated with an N-dimensional grid of points; defining a deformation grid G of control points, wherein a resolution of grid G is lower than a resolution of the grids for said source and target images; defining a coordinate transformation on G as T ⁡ ( x ) = x + ∑ p ∈ G ⁢ ⁢ η ⁡ (  x - p  ) ⁢ d u p , where x is an image point on said source image, p is a position vector of control point p, d p is a displacement vector for each control point on G, u p is a label for point p associated with displacement d p , and η( ) is a weighting function for said displacement vector; defining a Markov Random Field (MRF) energy functional to be minimized by T as E t = 1  G  ⁢ ∑ p ∈ G ⁢ ⁢ V p t ⁡ ( u p ) + 1  E  ⁢ ∑ p , q ∈ E ⁢ ⁢ V pq ⁡ ( u p , u q ) , wherein V p t ⁡ ( u p ) = ∫ Ω ⁢ η ^ ⁡ (  x - p  ) ⁢ ρ h ⁡ ( g ⁡ ( x ) , f ⁡ ( d u p · T t - 1 ⁡ ( x ) ) ) ⁢ ⁢ ⅆ x , ⁢ V pq ⁡ ( u p , u q ) = λ pq ⁢  d u p - d u q  , wherein f is the source image, g is the target image, |G| is a number of control points on G, |E| is a number of pairs of neighboring control points on a neighborhood system E defined on G, t is an iteration counter, T t−l is a previous approximation to the coordinate transformation, λ pq is a spatially varying weight, ρ h is a dissimilarity measure, and {circumflex over (η)}( ) is normalization factor defined as η ^ ⁡ (  x - p  ) = η ⁡ (  x - p  ) ∫ Ω ⁢ η ⁡ (  y - p  ) ⁢ ⁢ ⅆ y integrated over a domain of f; and associating said MRF with a primary linear program and solving said primary linear program using a fast primal-dual algorithm to yield a coordinate transformation T that minimizes said MRF energy functional by minimizing said dissimilarity measure ρ h .