Linear Least Squares using Householder transform

Hi everyone,
I am trying to find the scale, rotation and translation difference between two images. The coordinate transformation is defined by the equations:
x’ = Ax + By + C
y’ = Dx + Ey + F
Where (x’,y’) are points from the reference image, and (x,y) are corresponding points in the target image that needs to be scaled, rotated and translated to match the reference.

I believe this problem is usually solved with Linear Least Squares, using HouseHolder transformations. Is this possible within Hipparchus?

I tried using the Hipparchus Optimization LeastSquares, (Least squares | Hipparchus::Optim), but it was unclear to me how to handle the C or F term.

Thanks for your help, John

I would suggest you look at some of the validation tests in the org.hipparchus.optim.nonlinear.vector.leastsquare package. The first step would be to look at one of the implementations of test problems (say CircleProblem for example) and how it is used in LevenbergMarquardtOptimizerTest to create a LeastSquaresBuilder.
The two main things to implement are the model function and its Jacobian whith respect to the parameters:

f = \left\{ \begin{gather} A x_1 + B y_1 + C\\ D x_1 + E y_1 + F\\ \ldots\\ A x_n + B y_n + C\\ D x_n + E y_n + F\\ \end{gather} \right.
J = \begin{bmatrix} \frac{\partial x'_1}{\partial A} & \frac{\partial x'_1}{\partial B} & \frac{\partial x'_1}{\partial C} & \frac{\partial x'_1}{\partial D} & \frac{\partial x'_1}{\partial E} & \frac{\partial x'_1}{\partial F}\\ \frac{\partial y'_1}{\partial A} & \frac{\partial y'_1}{\partial B} & \frac{\partial y'_1}{\partial C} & \frac{\partial y'_1}{\partial D} & \frac{\partial y'_1}{\partial E} & \frac{\partial y'_1}{\partial F}\\ \ldots&\ldots&\ldots&\ldots&\ldots&\ldots\\ \frac{\partial x'_n}{\partial A} & \frac{\partial x'_n}{\partial B} & \frac{\partial x'_n}{\partial C} & \frac{\partial x'_n}{\partial D} & \frac{\partial x'_n}{\partial E} & \frac{\partial x'_n}{\partial F}\\ \frac{\partial y'_n}{\partial A} & \frac{\partial y'_n}{\partial B} & \frac{\partial y'_n}{\partial C} & \frac{\partial y'_n}{\partial D} & \frac{\partial y'_n}{\partial E} & \frac{\partial y'_n}{\partial F}\\ \end{bmatrix} = \begin{bmatrix} x_1 & y_1 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & x_1 & y_1 & 1\\ \ldots&\ldots&\ldots&\ldots&\ldots&\ldots\\ x_n & y_n & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & x_n & y_n & 1 \end{bmatrix}

When you have implemented these two mathematical functions (for example in a class that provides two methods getModelFunction() and getModelFunctionJacobian()), you can put them in a LeastSquaresProblem thanks to a LeastSquaresBuilder, using something along these lines:

LeastSquaresProblem myProblem =
  new LeastSquaresBuilder()
      .checkerPair(new SimpleVectorValueChecker(1e-12, 1e-12))
      .maxEvaluations(100)
      .maxIterations(100)
      .model(imageTransform.getModelFunction(),
             imageTransform.getModelFunctionJacobian())
      .target(targetArrayContainingXprimesYprimes)
      .start(startArrayContainingInitialABCDEFWhichMayBeFullOfZeros)
      .build();

The last step is to solve this least squares problem. This is done by calling:

Optimum optimum =
    new LevenbergMarquardtOptimizer().optimize(myProblem);

The optimum instance will contain a lot of things related to the optimum solution of your problem. This includes the number of iterations, the value of the cost function, covariances, residuals for all the test points, RMS and of course using getPoint() the vector containing the parameters (in your case A, B, C, D, E, and F).

There is no need to bother about Householder transformations. Everything is taken care of in the Levenberg-Marquardt algorithm. The Householder transformations are really low level parts and there is a lot of stuff on top of that to build a least squares solver. Levenberg-Marquardt is a good solver. It is a kind of hybrid method between a gradient method and a Gauss-Newton, its most renown property is that it is able to start even far from the solution (which a Gauss-Newton alone can’t do). This is why I suggest to just use an array full of zeros for the start point. If your two images are only slightly different, you could also start with a point where only A and E are set to 1.0, so your initial model corresponds to the identity transform.

1 Like

Hi Luc,
Thank you for your very clear explanation. I have implemented it, and it works extremely well.

I was at the start of a steep learning curve. I had studied the Least squares | Hipparchus::Optim web page, but I now realize that I had failed to grasp some key concepts. This caused my first attempts to fail.

After reading your clear explanation, I now have a much better understanding, so not only is my current problem solved, I think I should be able to solve future ones as well.

I could not have done it without your help. Thank you!