1. A method for registering digitized images using Markov random Fields, said method 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−1 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 .