diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py
new file mode 100644
index 00000000..6d3ea72f
--- /dev/null
+++ b/doc/examples/plot_matching.py
@@ -0,0 +1,140 @@
+"""
+============================
+Robust matching using RANSAC
+============================
+
+In this simplified example we first generate two synthetic images as if they
+were taken from different view points.
+
+In the next step we find interest points in both images and find
+correspondencies based on a weighted sum of squared differences of a small
+neighbourhood around them. Note, that this measure is only robust towards
+linear radiometric and not geometric distortions and is thus only usable with
+slight view point changes.
+
+After finding the correspondencies we end up having a set of source and
+destination coordinates which can be used to estimate the geometric
+transformation between both images. However, many of the correspondencies are
+faulty and simply estimating the parameter set with all coordinates is not
+sufficient. Therefore, the RANSAC algorithm is used on top of the normal model
+to robustly estimate the parameter set by detecting outliers.
+
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+
+from skimage import data
+from skimage.feature import corner_harris, corner_subpix, corner_peaks
+from skimage.transform import warp, AffineTransform
+from skimage.exposure import rescale_intensity
+from skimage.color import rgb2gray
+from skimage.measure import ransac
+
+
+# generate synthetic checkerboard image and add gradient for the later matching
+checkerboard = data.checkerboard()
+img_orig = np.zeros(list(checkerboard.shape) + [3])
+img_orig[..., 0] = checkerboard
+gradient_r, gradient_c = np.mgrid[0:img_orig.shape[0], 0:img_orig.shape[1]]
+img_orig[..., 1] = gradient_r
+img_orig[..., 2] = gradient_c
+img_orig = rescale_intensity(img_orig)
+img_orig_gray = rgb2gray(img_orig)
+
+# warp synthetic image
+tform = AffineTransform(scale=(0.9, 0.9), rotation=0.2, translation=(20, -10))
+img_warped = warp(img_orig, tform.inverse, output_shape=(200, 200))
+img_warped_gray = rgb2gray(img_warped)
+
+# extract corners using Harris' corner measure
+coords_orig = corner_peaks(corner_harris(img_orig_gray), threshold_rel=0.001,
+                           min_distance=5)
+coords_warped = corner_peaks(corner_harris(img_warped_gray),
+                             threshold_rel=0.001, min_distance=5)
+
+# determine subpixel corner position
+coords_orig_subpix = corner_subpix(img_orig_gray, coords_orig, window_size=10)
+coords_warped_subpix = corner_subpix(img_warped_gray, coords_warped,
+                                     window_size=10)
+
+
+def gaussian_weights(window_ext, sigma=1):
+    y, x = np.mgrid[-window_ext:window_ext+1, -window_ext:window_ext+1]
+    g = np.zeros(y.shape, dtype=np.double)
+    g[:] = np.exp(-0.5 * (x**2 / sigma**2 + y**2 / sigma**2))
+    g /= 2 * np.pi * sigma * sigma
+    return g
+
+
+def match_corner(coord, window_ext=5):
+    r, c =  np.round(coord)
+    window_orig = img_orig[r-window_ext:r+window_ext+1,
+                           c-window_ext:c+window_ext+1, :]
+
+    # weight pixels depending on distance to center pixel
+    weights = gaussian_weights(window_ext, 3)
+    weights = np.dstack((weights, weights, weights))
+
+    # compute sum of squared differences to all corners in warped image
+    SSDs = []
+    for cr, cc in coords_warped:
+        window_warped = img_warped[cr-window_ext:cr+window_ext+1,
+                                   cc-window_ext:cc+window_ext+1, :]
+        SSD = np.sum(weights * (window_orig - window_warped)**2)
+        SSDs.append(SSD)
+
+    # use corner with minimum SSD as correspondency
+    min_idx = np.argmin(SSDs)
+    return coords_warped_subpix[min_idx]
+
+
+# find correspondencies using simple weighted sum of squared differences
+src = []
+dst = []
+for coord in coords_orig_subpix:
+    src.append(coord)
+    dst.append(match_corner(coord))
+src = np.array(src)
+dst = np.array(dst)
+
+
+# estimate affine transform model using all coordinates
+model = AffineTransform()
+model.estimate(src, dst)
+
+# robustly estimate affine transform model with RANSAC
+model_robust, inliers = ransac((src, dst), AffineTransform, min_samples=3,
+                               residual_threshold=2, max_trials=100)
+outliers = inliers == False
+
+
+# compare "true" and estimated transform parameters
+print tform.scale, tform.translation, tform.rotation
+print model.scale, model.translation, model.rotation
+print model_robust.scale, model_robust.translation, model_robust.rotation
+
+
+# visualize correspondencies
+img_combined = np.concatenate((img_orig_gray, img_warped_gray), axis=1)
+
+fig, ax = plt.subplots(nrows=2, ncols=1)
+plt.gray()
+
+ax[0].imshow(img_combined, interpolation='nearest')
+ax[0].axis('off')
+ax[0].axis((0, 400, 200, 0))
+ax[0].set_title('Correct correspondencies')
+ax[1].imshow(img_combined, interpolation='nearest')
+ax[1].axis('off')
+ax[1].axis((0, 400, 200, 0))
+ax[1].set_title('Faulty correspondencies')
+
+
+for ax_idx, (m, color) in enumerate(((inliers, 'g'), (outliers, 'r'))):
+    ax[ax_idx].plot((src[m, 1], dst[m, 1] + 200), (src[m, 0], dst[m, 0]), '-',
+                    color=color)
+    ax[ax_idx].plot(src[m, 1], src[m, 0], '.', markersize=10, color=color)
+    ax[ax_idx].plot(dst[m, 1] + 200, dst[m, 0], '.', markersize=10,
+                    color=color)
+
+plt.show()
diff --git a/doc/examples/plot_ransac.py b/doc/examples/plot_ransac.py
new file mode 100644
index 00000000..26fd78c0
--- /dev/null
+++ b/doc/examples/plot_ransac.py
@@ -0,0 +1,55 @@
+"""
+=========================================
+Robust line model estimation using RANSAC
+=========================================
+
+In this example we see how to robustly fit a line model to faulty data using
+the RANSAC algorithm.
+
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+
+from skimage.measure import LineModel, ransac
+
+
+np.random.seed(seed=1)
+
+# generate coordinates of line
+x = np.arange(-200, 200)
+y = 0.2 * x + 20
+data = np.column_stack([x, y])
+
+# add faulty data
+faulty = np.array(30 * [(180, -100)])
+faulty += 5 * np.random.normal(size=faulty.shape)
+data[:faulty.shape[0]] = faulty
+
+# add gaussian noise to coordinates
+noise = np.random.normal(size=data.shape)
+data += 0.5 * noise
+data[::2] += 5 * noise[::2]
+data[::4] += 20 * noise[::4]
+
+# fit line using all data
+model = LineModel()
+model.estimate(data)
+
+# robustly fit line only using inlier data with RANSAC algorithm
+model_robust, inliers = ransac(data, LineModel, min_samples=2,
+                               residual_threshold=1, max_trials=1000)
+outliers = inliers == False
+
+# generate coordinates of estimated models
+line_x = np.arange(-250, 250)
+line_y = model.predict_y(line_x)
+line_y_robust = model_robust.predict_y(line_x)
+
+plt.plot(data[inliers, 0], data[inliers, 1], '.b', alpha=0.6,
+         label='Inlier data')
+plt.plot(data[outliers, 0], data[outliers, 1], '.r', alpha=0.6,
+         label='Outlier data')
+plt.plot(line_x, line_y, '-k', label='Line model from all data')
+plt.plot(line_x, line_y_robust, '-b', label='Robust line model')
+plt.legend(loc='lower left')
+plt.show()
diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py
index 226c4235..89f707c1 100755
--- a/skimage/measure/__init__.py
+++ b/skimage/measure/__init__.py
@@ -2,10 +2,16 @@ from .find_contours import find_contours
 from ._regionprops import regionprops, perimeter
 from ._structural_similarity import structural_similarity
 from ._polygon import approximate_polygon, subdivide_polygon
+from .fit import LineModel, CircleModel, EllipseModel, ransac
+
 
 __all__ = ['find_contours',
            'regionprops',
            'perimeter',
            'structural_similarity',
            'approximate_polygon',
-           'subdivide_polygon']
+           'subdivide_polygon',
+           'LineModel',
+           'CircleModel',
+           'EllipseModel',
+           'ransac']
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
new file mode 100644
index 00000000..87febf70
--- /dev/null
+++ b/skimage/measure/fit.py
@@ -0,0 +1,653 @@
+import math
+import numpy as np
+from scipy import optimize
+
+
+def _check_data_dim(data, dim):
+    if data.ndim != 2 or data.shape[1] != dim:
+        raise ValueError('Input data must have shape (N, %d).' % dim)
+
+
+class BaseModel(object):
+
+    def __init__(self):
+        self._params = None
+
+
+class LineModel(BaseModel):
+
+    """Total least squares estimator for 2D lines.
+
+    Lines are parameterized using polar coordinates as functional model::
+
+        dist = x * cos(theta) + y * sin(theta)
+
+    This parameterization is able to model vertical lines in contrast to the
+    standard line model ``y = a*x + b``.
+
+    This estimator minimizes the squared distances from all points to the
+    line::
+
+        min{ sum((dist - x_i * cos(theta) + y_i * sin(theta))**2) }
+
+    The ``_params`` attribute contains the parameters in the following order::
+
+        dist, theta
+
+    A minimum number of 2 points is required to solve for the parameters.
+
+    """
+
+    def estimate(self, data):
+        """Estimate line model from data using total least squares.
+
+        Parameters
+        ----------
+        data : (N, 2) array
+            N points with ``(x, y)`` coordinates, respectively.
+
+        """
+
+        _check_data_dim(data, dim=2)
+
+        X0 = data.mean(axis=0)
+
+        if data.shape[0] == 2:  # well determined
+            theta = np.arctan2(data[1, 1] - data[0, 1],
+                               data[1, 0] - data[0, 0])
+        elif data.shape[0] > 2:  # over-determined
+            data = data - X0
+            # first principal component
+            _, _, v = np.linalg.svd(data)
+            theta = np.arctan2(v[0, 1], v[0, 0])
+        else:  # under-determined
+            raise ValueError('At least 2 input points needed.')
+
+        # angle perpendicular to line angle
+        theta = (theta + np.pi / 2) % np.pi
+        # line always passes through mean
+        dist = X0[0] * math.cos(theta) + X0[1] * math.sin(theta)
+
+        self._params = (dist, theta)
+
+    def residuals(self, data):
+        """Determine residuals of data to model.
+
+        For each point the shortest distance to the line is returned.
+
+        Parameters
+        ----------
+        data : (N, 2) array
+            N points with ``(x, y)`` coordinates, respectively.
+
+        Returns
+        -------
+        residuals : (N, ) array
+            Residual for each data point.
+
+        """
+
+        _check_data_dim(data, dim=2)
+
+        dist, theta = self._params
+
+        x = data[:, 0]
+        y = data[:, 1]
+
+        return dist - (x * math.cos(theta) + y * math.sin(theta))
+
+    def predict_x(self, y, params=None):
+        """Predict x-coordinates using the estimated model.
+
+        Parameters
+        ----------
+        y : array
+            y-coordinates.
+        params : (2, ) array, optional
+            Optional custom parameter set.
+
+        Returns
+        -------
+        x : array
+            Predicted x-coordinates.
+
+        """
+
+        if params is None:
+            params = self._params
+        dist, theta = params
+        return (dist - y * math.sin(theta)) / math.cos(theta)
+
+    def predict_y(self, x, params=None):
+        """Predict y-coordinates using the estimated model.
+
+        Parameters
+        ----------
+        x : array
+            x-coordinates.
+        params : (2, ) array, optional
+            Optional custom parameter set.
+
+        Returns
+        -------
+        y : array
+            Predicted y-coordinates.
+
+        """
+
+        if params is None:
+            params = self._params
+        dist, theta = params
+        return (dist - x * math.cos(theta)) / math.sin(theta)
+
+
+class CircleModel(BaseModel):
+
+    """Total least squares estimator for 2D circles.
+
+    The functional model of the circle is::
+
+        r**2 = (x - xc)**2 + (y - yc)**2
+
+    This estimator minimizes the squared distances from all points to the
+    circle::
+
+        min{ sum((r - sqrt((x_i - xc)**2 + (y_i - yc)**2))**2) }
+
+    The ``_params`` attribute contains the parameters in the following order::
+
+        xc, yc, r
+
+    A minimum number of 3 points is required to solve for the parameters.
+
+    """
+
+    def estimate(self, data):
+        """Estimate circle model from data using total least squares.
+
+        Parameters
+        ----------
+        data : (N, 2) array
+            N points with ``(x, y)`` coordinates, respectively.
+
+        """
+
+        _check_data_dim(data, dim=2)
+
+        x = data[:, 0]
+        y = data[:, 1]
+        # pre-allocate jacobian for all iterations
+        A = np.zeros((3, data.shape[0]), dtype=np.double)
+        # same for all iterations: r
+        A[2, :] = -1
+
+        def dist(xc, yc):
+            return np.sqrt((x - xc)**2 + (y - yc)**2)
+
+        def fun(params):
+            xc, yc, r = params
+            return dist(xc, yc) - r
+
+        def Dfun(params):
+            xc, yc, r = params
+            d = dist(xc, yc)
+            A[0, :] = -(x - xc) / d
+            A[1, :] = -(y - yc) / d
+            # same for all iterations, so not changed in each iteration
+            #A[2, :] = -1
+            return A
+
+        xc0 = x.mean()
+        yc0 = y.mean()
+        r0 = dist(xc0, yc0).mean()
+        params0 = (xc0, yc0, r0)
+        params, _ = optimize.leastsq(fun, params0, Dfun=Dfun, col_deriv=True)
+
+        self._params = params
+
+    def residuals(self, data):
+        """Determine residuals of data to model.
+
+        For each point the shortest distance to the circle is returned.
+
+        Parameters
+        ----------
+        data : (N, 2) array
+            N points with ``(x, y)`` coordinates, respectively.
+
+        Returns
+        -------
+        residuals : (N, ) array
+            Residual for each data point.
+
+        """
+
+        _check_data_dim(data, dim=2)
+
+        xc, yc, r = self._params
+
+        x = data[:, 0]
+        y = data[:, 1]
+
+        return r - np.sqrt((x - xc)**2 + (y - yc)**2)
+
+    def predict_xy(self, t, params=None):
+        """Predict x- and y-coordinates using the estimated model.
+
+        Parameters
+        ----------
+        t : array
+            Angles in circle in radians. Angles start to count from positive
+            x-axis to positive y-axis in a right-handed system.
+        params : (3, ) array, optional
+            Optional custom parameter set.
+
+        Returns
+        -------
+        xy : (..., 2) array
+            Predicted x- and y-coordinates.
+
+        """
+        if params is None:
+            params = self._params
+        xc, yc, r = params
+
+        x = xc + r * np.cos(t)
+        y = yc + r * np.sin(t)
+
+        return np.concatenate((x[..., None], y[..., None]), axis=t.ndim)
+
+
+class EllipseModel(BaseModel):
+
+    """Total least squares estimator for 2D ellipses.
+
+    The functional model of the ellipse is::
+
+        xt = xc + a*cos(theta)*cos(t) - b*sin(theta)*sin(t)
+        yt = yc + a*sin(theta)*cos(t) + b*cos(theta)*sin(t)
+        d = sqrt((x - xt)**2 + (y - yt)**2)
+
+    where ``(xt, yt)`` is the closest point on the ellipse to ``(x, y)``. Thus
+    d is the shortest distance from the point to the ellipse.
+
+    This estimator minimizes the squared distances from all points to the
+    ellipse::
+
+        min{ sum(d_i**2) } = min{ sum((x_i - xt)**2 + (y_i - yt)**2) }
+
+    Thus you have ``2 * N`` equations (x_i, y_i) for ``N + 5`` unknowns (t_i,
+    xc, yc, a, b, theta), which gives you an effective redundancy of ``N - 5``.
+
+    The ``_params`` attribute contains the parameters in the following order::
+
+        xc, yc, a, b, theta
+
+    A minimum number of 5 points is required to solve for the parameters.
+
+    """
+
+    def estimate(self, data):
+        """Estimate circle model from data using total least squares.
+
+        Parameters
+        ----------
+        data : (N, 2) array
+            N points with ``(x, y)`` coordinates, respectively.
+
+        """
+
+        _check_data_dim(data, dim=2)
+
+        x = data[:, 0]
+        y = data[:, 1]
+
+        N = data.shape[0]
+
+        # pre-allocate jacobian for all iterations
+        A = np.zeros((N + 5, 2 * N), dtype=np.double)
+        # same for all iterations: xc, yc
+        A[0, :N] = -1
+        A[1, N:] = -1
+
+        diag_idxs = np.diag_indices(N)
+
+        def fun(params):
+            xyt = self.predict_xy(params[5:], params[:5])
+            fx = x - xyt[:, 0]
+            fy = y - xyt[:, 1]
+            return np.append(fx, fy)
+
+        def Dfun(params):
+            xc, yc, a, b, theta = params[:5]
+            t = params[5:]
+
+            ct = np.cos(t)
+            st = np.sin(t)
+            ctheta = math.cos(theta)
+            stheta = math.sin(theta)
+
+            # derivatives for fx, fy in the following order:
+            #       xc, yc, a, b, theta, t_i
+
+            # fx
+            A[2, :N] = - ctheta * ct
+            A[3, :N] = stheta * st
+            A[4, :N] = a * stheta * ct + b * ctheta * st
+            A[5:, :N][diag_idxs] = a * ctheta * st + b * stheta * ct
+            # fy
+            A[2, N:] = - stheta * ct
+            A[3, N:] = - ctheta * st
+            A[4, N:] = - a * ctheta * ct + b * stheta * st
+            A[5:, N:][diag_idxs] = a * stheta * st - b * ctheta * ct
+
+            return A
+
+        # initial guess of parameters using a circle model
+        params0 = np.empty((N + 5, ), dtype=np.double)
+        xc0 = x.mean()
+        yc0 = y.mean()
+        r0 = np.sqrt((x - xc0)**2 + (y - yc0)**2).mean()
+        params0[:5] = (xc0, yc0, r0, 0, 0)
+        params0[5:] = np.arctan2(y - yc0, x - xc0)
+
+        params, _ = optimize.leastsq(fun, params0, Dfun=Dfun, col_deriv=True)
+
+        self._params = params[:5]
+
+    def residuals(self, data):
+        """Determine residuals of data to model.
+
+        For each point the shortest distance to the ellipse is returned.
+
+        Parameters
+        ----------
+        data : (N, 2) array
+            N points with ``(x, y)`` coordinates, respectively.
+
+        Returns
+        -------
+        residuals : (N, ) array
+            Residual for each data point.
+
+        """
+
+        _check_data_dim(data, dim=2)
+
+        xc, yc, a, b, theta = self._params
+
+        ctheta = math.cos(theta)
+        stheta = math.sin(theta)
+
+        x = data[:, 0]
+        y = data[:, 1]
+
+        N = data.shape[0]
+
+        def fun(t, xi, yi):
+            ct = math.cos(t)
+            st = math.sin(t)
+            xt = xc + a * ctheta * ct - b * stheta * st
+            yt = yc + a * stheta * ct + b * ctheta * st
+            return (xi - xt)**2 + (yi - yt)**2
+
+        # def Dfun(t, xi, yi):
+        #     ct = math.cos(t)
+        #     st = math.sin(t)
+        #     xt = xc + a * ctheta * ct - b * stheta * st
+        #     yt = yc + a * stheta * ct + b * ctheta * st
+        #     dfx_t = - 2 * (xi - xt) * (- a * ctheta * st
+        #                                - b * stheta * ct)
+        #     dfy_t = - 2 * (yi - yt) * (- a * stheta * st
+        #                                + b * ctheta * ct)
+        #     return [dfx_t + dfy_t]
+
+        residuals = np.empty((N, ), dtype=np.double)
+
+        # initial guess for parameter t of closest point on ellipse
+        t0 = np.arctan2(y - yc, x - xc) - theta
+
+        # determine shortest distance to ellipse for each point
+        for i in range(N):
+            xi = x[i]
+            yi = y[i]
+            # faster without Dfun, because of the python overhead
+            t, _ = optimize.leastsq(fun, t0[i], args=(xi, yi))
+            residuals[i] = np.sqrt(fun(t, xi, yi))
+
+        return residuals
+
+    def predict_xy(self, t, params=None):
+        """Predict x- and y-coordinates using the estimated model.
+
+        Parameters
+        ----------
+        t : array
+            Angles in circle in radians. Angles start to count from positive
+            x-axis to positive y-axis in a right-handed system.
+        params : (5, ) array, optional
+            Optional custom parameter set.
+
+        Returns
+        -------
+        xy : (..., 2) array
+            Predicted x- and y-coordinates.
+
+        """
+
+        if params is None:
+            params = self._params
+        xc, yc, a, b, theta = params
+
+        ct = np.cos(t)
+        st = np.sin(t)
+        ctheta = math.cos(theta)
+        stheta = math.sin(theta)
+
+        x = xc + a * ctheta * ct - b * stheta * st
+        y = yc + a * stheta * ct + b * ctheta * st
+
+        return np.concatenate((x[..., None], y[..., None]), axis=t.ndim)
+
+
+def ransac(data, model_class, min_samples, residual_threshold,
+           is_data_valid=None, is_model_valid=None,
+           max_trials=100, stop_sample_num=np.inf, stop_residuals_sum=0):
+    """Fit a model to data with the RANSAC (random sample consensus) algorithm.
+
+    RANSAC is an iterative algorithm for the robust estimation of parameters
+    from a subset of inliers from the complete data set. Each iteration
+    performs the following tasks:
+
+    1. Select `min_samples` random samples from the original data and check
+       whether the set of data is valid (see `is_data_valid`).
+    2. Estimate a model to the random subset
+       (`model_cls.estimate(*data[random_subset]`) and check whether the
+       estimated model is valid (see `is_model_valid`).
+    3. Classify all data as inliers or outliers by calculating the residuals
+       to the estimated model (`model_cls.residuals(*data)`) - all data samples
+       with residuals smaller than the `residual_threshold` are considered as
+       inliers.
+    4. Save estimated model as best model if number of inlier samples is
+       maximal. In case the current estimated model has the same number of
+       inliers, it is only considered as the best model if it has less sum of
+       residuals.
+
+    These steps are performed either a maximum number of times or until one of
+    the special stop criteria are met. The final model is estimated using all
+    inlier samples of the previously determined best model.
+
+    Parameters
+    ----------
+    data : [list, tuple of] (N, D) array
+        Data set to which the model is fitted, where N is the number of data
+        points and D the dimensionality of the data.
+        If the model class requires multiple input data arrays (e.g. source and
+        destination coordinates of  ``skimage.transform.AffineTransform``),
+        they can be optionally passed as tuple or list. Note, that in this case
+        the functions ``estimate(*data)``, ``residuals(*data)``,
+        ``is_model_valid(model, *random_data)`` and
+        ``is_data_valid(*random_data)`` must all take each data array as
+        separate arguments.
+    model_class : object
+        Object with the following object methods:
+
+         * ``estimate(*data)``
+         * ``residuals(*data)``
+
+    min_samples : int
+        The minimum number of data points to fit a model to.
+    residual_threshold : float
+        Maximum distance for a data point to be classified as an inlier.
+    is_data_valid : function, optional
+        This function is called with the randomly selected data before the
+        model is fitted to it: `is_data_valid(*random_data)`.
+    is_model_valid : function, optional
+        This function is called with the estimated model and the randomly
+        selected data: `is_model_valid(model, *random_data)`, .
+    max_trials : int, optional
+        Maximum number of iterations for random sample selection.
+    stop_sample_num : int, optional
+        Stop iteration if at least this number of inliers are found.
+    stop_residuals_sum : float, optional
+        Stop iteration if sum of residuals is less equal than this threshold.
+
+    Returns
+    -------
+    model : object
+        Best model with largest consensus set.
+    inliers : (N, ) array
+        Boolean mask of inliers classified as ``True``.
+
+    References
+    ----------
+    .. [1] "RANSAC", Wikipedia, http://en.wikipedia.org/wiki/RANSAC
+
+    Examples
+    --------
+
+    Generate ellipse data without tilt and add noise:
+
+    >>> t = np.linspace(0, 2 * np.pi, 50)
+    >>> a = 5
+    >>> b = 10
+    >>> xc = 20
+    >>> yc = 30
+    >>> x = xc + a * np.cos(t)
+    >>> y = yc + b * np.sin(t)
+    >>> data = np.column_stack([x, y])
+    >>> np.random.seed(seed=1234)
+    >>> data += np.random.normal(size=data.shape)
+
+    Add some faulty data:
+
+    >>> data[0] = (100, 100)
+    >>> data[1] = (110, 120)
+    >>> data[2] = (120, 130)
+    >>> data[3] = (140, 130)
+
+    Estimate ellipse model using all available data:
+
+    >>> model = EllipseModel()
+    >>> model.estimate(data)
+    >>> model._params
+    array([  4.85808595e+02,   4.51492793e+02,   1.15018491e+03,
+             5.52428289e+00,   7.32420126e-01])
+
+    Estimate ellipse model using RANSAC:
+
+    >>> ransac_model, inliers = ransac(data, EllipseModel, 5, 3, max_trials=50)
+    >>> # ransac_model._params, inliers
+
+    Should give the correct result estimated without the faulty data:
+
+        [ 20.12762373, 29.73563061, 4.81499637, 10.4743584, 0.05217117]
+        [ 4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
+         21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,
+         38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]
+
+    Robustly estimate geometric transformation:
+
+    >>> from skimage.transform import SimilarityTransform
+    >>> src = 100 * np.random.random((50, 2))
+    >>> model0 = SimilarityTransform(scale=0.5, rotation=1,
+    ...                              translation=(10, 20))
+    >>> dst = model0(src)
+    >>> dst[0] = (10000, 10000)
+    >>> dst[1] = (-100, 100)
+    >>> dst[2] = (50, 50)
+    >>> model, inliers = ransac((src, dst), SimilarityTransform, 2, 10)
+    >>> inliers
+    array([ 3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
+       20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
+       37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49])
+
+
+    """
+
+    best_model = None
+    best_inlier_num = 0
+    best_inlier_residuals_sum = np.inf
+    best_inliers = None
+
+    if not isinstance(data, list) and not isinstance(data, tuple):
+        data = [data]
+
+    # make sure data is list and not tuple, so it can be modified below
+    data = list(data)
+    # number of samples
+    N = data[0].shape[0]
+
+    for _ in range(max_trials):
+
+        # choose random sample set
+        samples = []
+        random_idxs = np.random.randint(0, N, min_samples)
+        for d in data:
+            samples.append(d[random_idxs])
+
+        # check if random sample set is valid
+        if is_data_valid is not None and not is_data_valid(*samples):
+            continue
+
+        # estimate model for current random sample set
+        sample_model = model_class()
+        sample_model.estimate(*samples)
+
+        # check if estimated model is valid
+        if is_model_valid is not None and not is_model_valid(sample_model,
+                                                             *samples):
+            continue
+
+        sample_model_residuals = np.abs(sample_model.residuals(*data))
+        # consensus set / inliers
+        sample_model_inliers = sample_model_residuals < residual_threshold
+        sample_model_residuals_sum = np.sum(sample_model_residuals**2)
+
+        # choose as new best model if number of inliers is maximal
+        sample_inlier_num = np.sum(sample_model_inliers)
+        if (
+            # more inliers
+            sample_inlier_num > best_inlier_num
+            # same number of inliers but less "error" in terms of residuals
+            or (sample_inlier_num == best_inlier_num
+                and sample_model_residuals_sum < best_inlier_residuals_sum)
+        ):
+            best_model = sample_model
+            best_inlier_num = sample_inlier_num
+            best_inlier_residuals_sum = sample_model_residuals_sum
+            best_inliers = sample_model_inliers
+            if (
+                best_inlier_num >= stop_sample_num
+                or best_inlier_residuals_sum <= stop_residuals_sum
+            ):
+                break
+
+    # estimate final model using all inliers
+    if best_inliers is not None:
+        # select inliers for each data array
+        for i in range(len(data)):
+            data[i] = data[i][best_inliers]
+        best_model.estimate(*data)
+
+    return best_model, best_inliers
diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py
new file mode 100644
index 00000000..bde1fcc7
--- /dev/null
+++ b/skimage/measure/tests/test_fit.py
@@ -0,0 +1,208 @@
+import numpy as np
+from numpy.testing import assert_equal, assert_raises, assert_almost_equal
+from skimage.measure import LineModel, CircleModel, EllipseModel, ransac
+from skimage.transform import AffineTransform
+
+
+def test_line_model_invalid_input():
+    assert_raises(ValueError, LineModel().estimate, np.empty((5, 3)))
+
+
+def test_line_model_predict():
+    model = LineModel()
+    model._params = (10, 1)
+    x = np.arange(-10, 10)
+    y = model.predict_y(x)
+    assert_almost_equal(x, model.predict_x(y))
+
+
+def test_line_model_estimate():
+    # generate original data without noise
+    model0 = LineModel()
+    model0._params = (10, 1)
+    x0 = np.arange(-100, 100)
+    y0 = model0.predict_y(x0)
+    data0 = np.column_stack([x0, y0])
+
+    # add gaussian noise to data
+    np.random.seed(1234)
+    data = data0 + np.random.normal(size=data0.shape)
+
+    # estimate parameters of noisy data
+    model_est = LineModel()
+    model_est.estimate(data)
+
+    # test whether estimated parameters almost equal original parameters
+    assert_almost_equal(model0._params, model_est._params, 1)
+
+
+def test_line_model_residuals():
+    model = LineModel()
+    model._params = (0, 0)
+    assert_equal(abs(model.residuals(np.array([[0, 0]]))), 0)
+    assert_equal(abs(model.residuals(np.array([[0, 10]]))), 0)
+    assert_equal(abs(model.residuals(np.array([[10, 0]]))), 10)
+    model._params = (5, np.pi / 4)
+    assert_equal(abs(model.residuals(np.array([[0, 0]]))), 5)
+    assert_equal(abs(model.residuals(np.array([[np.sqrt(50), 0]]))), 5)
+
+
+def test_line_model_under_determined():
+    data = np.empty((1, 2))
+    assert_raises(ValueError, LineModel().estimate, data)
+
+
+def test_circle_model_invalid_input():
+    assert_raises(ValueError, CircleModel().estimate, np.empty((5, 3)))
+
+
+def test_circle_model_predict():
+    model = CircleModel()
+    r = 5
+    model._params = (0, 0, r)
+    t = np.arange(0, 2 * np.pi, np.pi / 2)
+
+    xy = np.array(((5, 0), (0, 5), (-5, 0), (0, -5)))
+    assert_almost_equal(xy, model.predict_xy(t))
+
+
+def test_circle_model_estimate():
+    # generate original data without noise
+    model0 = CircleModel()
+    model0._params = (10, 12, 3)
+    t = np.linspace(0, 2 * np.pi, 1000)
+    data0 = model0.predict_xy(t)
+
+    # add gaussian noise to data
+    np.random.seed(1234)
+    data = data0 + np.random.normal(size=data0.shape)
+
+    # estimate parameters of noisy data
+    model_est = CircleModel()
+    model_est.estimate(data)
+
+    # test whether estimated parameters almost equal original parameters
+    assert_almost_equal(model0._params, model_est._params, 1)
+
+
+def test_circle_model_residuals():
+    model = CircleModel()
+    model._params = (0, 0, 5)
+    assert_almost_equal(abs(model.residuals(np.array([[5, 0]]))), 0)
+    assert_almost_equal(abs(model.residuals(np.array([[6, 6]]))),
+                        np.sqrt(2 * 6**2) - 5)
+    assert_almost_equal(abs(model.residuals(np.array([[10, 0]]))), 5)
+
+
+def test_ellipse_model_invalid_input():
+    assert_raises(ValueError, EllipseModel().estimate, np.empty((5, 3)))
+
+
+def test_ellipse_model_predict():
+    model = EllipseModel()
+    r = 5
+    model._params = (0, 0, 5, 10, 0)
+    t = np.arange(0, 2 * np.pi, np.pi / 2)
+
+    xy = np.array(((5, 0), (0, 10), (-5, 0), (0, -10)))
+    assert_almost_equal(xy, model.predict_xy(t))
+
+
+def test_ellipse_model_estimate():
+    # generate original data without noise
+    model0 = EllipseModel()
+    model0._params = (10, 20, 15, 25, 0)
+    t = np.linspace(0, 2 * np.pi, 100)
+    data0 = model0.predict_xy(t)
+
+    # add gaussian noise to data
+    np.random.seed(1234)
+    data = data0 + np.random.normal(size=data0.shape)
+
+    # estimate parameters of noisy data
+    model_est = EllipseModel()
+    model_est.estimate(data)
+
+    # test whether estimated parameters almost equal original parameters
+    assert_almost_equal(model0._params, model_est._params, 0)
+
+
+def test_line_model_residuals():
+    model = EllipseModel()
+    # vertical line through origin
+    model._params = (0, 0, 10, 5, 0)
+    assert_almost_equal(abs(model.residuals(np.array([[10, 0]]))), 0)
+    assert_almost_equal(abs(model.residuals(np.array([[0, 5]]))), 0)
+    assert_almost_equal(abs(model.residuals(np.array([[0, 10]]))), 5)
+
+
+def test_ransac_shape():
+    np.random.seed(1)
+
+    # generate original data without noise
+    model0 = CircleModel()
+    model0._params = (10, 12, 3)
+    t = np.linspace(0, 2 * np.pi, 1000)
+    data0 = model0.predict_xy(t)
+
+    # add some faulty data
+    outliers = (10, 30, 200)
+    data0[outliers[0], :] = (1000, 1000)
+    data0[outliers[1], :] = (-50, 50)
+    data0[outliers[2], :] = (-100, -10)
+
+    # estimate parameters of corrupted data
+    model_est, inliers = ransac(data0, CircleModel, 3, 5)
+
+    # test whether estimated parameters equal original parameters
+    assert_equal(model0._params, model_est._params)
+    for outlier in outliers:
+        assert outlier not in inliers
+
+
+def test_ransac_geometric():
+    np.random.seed(1)
+
+    # generate original data without noise
+    src = 100 * np.random.random((50, 2))
+    model0 = AffineTransform(scale=(0.5, 0.3), rotation=1,
+                             translation=(10, 20))
+    dst = model0(src)
+
+    # add some faulty data
+    outliers = (0, 5, 20)
+    dst[outliers[0]] = (10000, 10000)
+    dst[outliers[1]] = (-100, 100)
+    dst[outliers[2]] = (50, 50)
+
+    # estimate parameters of corrupted data
+    model_est, inliers = ransac((src, dst), AffineTransform, 2, 20)
+
+    # test whether estimated parameters equal original parameters
+    assert_almost_equal(model0._matrix, model_est._matrix)
+    assert np.all(np.nonzero(inliers == False)[0] == outliers)
+
+
+def test_ransac_is_data_valid():
+    np.random.seed(1)
+
+    is_data_valid = lambda data: data.shape[0] > 2
+    model, inliers = ransac(np.empty((10, 2)), LineModel, 2, np.inf,
+                            is_data_valid=is_data_valid)
+    assert_equal(model, None)
+    assert_equal(inliers, None)
+
+
+def test_ransac_is_model_valid():
+    np.random.seed(1)
+
+    def is_model_valid(model, data):
+        return False
+    model, inliers = ransac(np.empty((10, 2)), LineModel, 2, np.inf,
+                            is_model_valid=is_model_valid)
+    assert_equal(model, None)
+    assert_equal(inliers, None)
+
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 448829de..1ff8e9d5 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -41,6 +41,28 @@ class GeometricTransform(object):
         """
         raise NotImplementedError()
 
+    def residuals(self, src, dst):
+        """Determine residuals of transformed destination coordinates.
+
+        For each transformed source coordinate the euclidean distance to the
+        respective destination coordinate is determined.
+
+        Parameters
+        ----------
+        src : (N, 2) array
+            Source coordinates.
+        dst : (N, 2) array
+            Destination coordinates.
+
+        Returns
+        -------
+        residuals : (N, ) array
+            Residual for coordinate.
+
+        """
+
+        return np.sqrt(np.sum((self(src) - dst)**2, axis=1))
+
     def __add__(self, other):
         """Combine this transformation with another.