Estimate Euclidean transformation with python

I want to do something similar to what in image analysis would be a standard 'image registration' using features.

I want to find the best transformation that transforms a set of 2D coordinates A in another one B. But I want to add an extra constraint being that the transformation is 'rigid/Euclidean transformation' Meaning that there is no scaling but only translation and rotation. Normally allowing scaling I would do:

 from skimage import io, transform
 destination = array([[1.0,2.0],[1.0,4.0],[3.0,3.0],[3.0,7.0]])
 source = array([[1.2,1.7],[1.1,3.8],[3.1,3.4],[2.6,7.0]])
 T = transform.estimate_transform('similarity',source,destination)

I believeestimate_transform under the hood just solves a least squares problem. But I want to add the constraint of no scaling.

Are there any function in skimage or other packages that solve this? Probably I need to write my own optimization problem with scipy, CVXOPT or cvxpy. Any help to phrase/implement this optimization problem?

EDIT: My implementation thanks to Stefan van der Walt Answer

from matplotlib.pylab import *
from scipy.optimize import *

def obj_fun(pars,x,src):
    theta, tx, ty = pars
    H = array([[cos(theta), -sin(theta), tx],\
         [sin(theta), cos(theta), ty],
         [0,0,1]])
    src1 = c_[src,ones(src.shape[0])]
    return sum( (x - src1.dot(H.T)[:,:2])**2 )

def apply_transform(pars, src):
    theta, tx, ty = pars
    H = array([[cos(theta), -sin(theta), tx],\
         [sin(theta), cos(theta), ty],
         [0,0,1]])
    src1 = c_[src,ones(src.shape[0])]
    return src1.dot(H.T)[:,:2]

res = minimize(obj_fun,[0,0,0],args=(dst,src), method='Nelder-Mead')

With that extra constraint you are no longer solving a linear least squares problem, so you'll have to use one of SciPy's minimization functions. The inner part of your minimization would set up a matrix H:

H = np.array([[np.cos(theta), -np.sin(theta), tx],
              [np.sin(theta), np.cos(theta), ty],
              [0, 0, 1]])

Then, you would compute the distance

|x_target - H.dot(x_source)|

for all data-points and sum the errors. Now, you have a cost function that you can send to the minimization function. You probably will also want to make use of RANSAC, which is available as skimage.measure.ransac, to reject outliers.