Skip to content

Commit

Permalink
Decoupling of matrix computation from image-registration API (#16)
Browse files Browse the repository at this point in the history
* Created an API to hook-in custom linear algebra solver implementations.
* Created a SolverProvider as service hook.
* Provided 4 equivalent solver implementations.
* Removed the Jama and La4J specific implementations.
* Introduced the SolverProvider as linear algebra lookup service.
* Updated all calculation, transform and equation classes accordingly.
* Fixed typo in POM.
  • Loading branch information
Oliver-Loeffler committed Aug 1, 2020
1 parent 01ca901 commit a2c3dfa
Show file tree
Hide file tree
Showing 23 changed files with 722 additions and 524 deletions.
8 changes: 7 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,13 @@

Image registration is the process of finding the transform to match a given image with a desired reference. This library supports rigid body transforms (translation and rotation) as well as affine transforms (anisotropic scaling, anisotropic rotation/shear aka. non-orthogonality). Currently non-linear transforms are not supported.
Versions up to and including 0.0.4 run with Java-8, all later version starting with
0.0.5 require Java-11.
0.0.6 require Java-11.

Version 0.0.5 will introduce a method to select different linear algebra frameworks to
be used in `image-registration`. Here it will be possible to choose either `gov.nist.math.jama`,
`org.la4j`, `org.ejml` or `org.apache.commons.math3`. There will be a core or api JAR
for `image-registration` and there will be a group of JARs providing a solver implementation.
The solver will be selected using the Java SPI (Service Provider Interface) mechanism.

These methods are used e.g. in photomask manufacturing, medical imaging or geospatial applications.
Control point or feature based methods have only limited scope of use in medical imaging, there intensity based or voxel based methods are preferred due to the natural structure of medical image data. These additional intensity and voxel based methods are not supported by this library.
Expand Down
19 changes: 17 additions & 2 deletions pom.xml
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
<groupId>net.raumzeitfalle.registration</groupId>
<artifactId>image-registration</artifactId>
<version>0.0.5-SNAPSHOT</version>
<description>A colletion of functions to perform rigid and affine transforms for image registration.</description>
<description>A collection of functions to perform rigid and affine transforms for image registration.</description>
<inceptionYear>2019</inceptionYear>
<organization>
<name>Oliver Loeffler, Raumzeitfalle.net</name>
Expand All @@ -22,7 +22,9 @@

<org.la4j.version>0.6.0</org.la4j.version>
<gov.nist.math.jama.version>1.1.1</gov.nist.math.jama.version>

<org.ejml.version>0.39</org.ejml.version>
<org.apache.math.commons.version>3.6.1</org.apache.math.commons.version>

<junit.jupiter.version>5.6.0</junit.jupiter.version>

<maven.enforcer.plugin.version>3.0.0-M3</maven.enforcer.plugin.version>
Expand Down Expand Up @@ -255,6 +257,19 @@
<version>${gov.nist.math.jama.version}</version>
<scope>compile</scope>
</dependency>

<dependency>
<groupId>org.ejml</groupId>
<artifactId>ejml-simple</artifactId>
<version>${org.ejml.version}</version>
<scope>compile</scope>
</dependency>

<dependency>
<groupId>org.apache.commons</groupId>
<artifactId>commons-math3</artifactId>
<version>${org.apache.math.commons.version}</version>
</dependency>

<dependency>
<groupId>org.junit.jupiter</groupId>
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
package net.raumzeitfalle.registration;

public interface ModelEquation extends Orientable {
public double getDeltaValue();
}
Original file line number Diff line number Diff line change
@@ -1,136 +1,122 @@
/*-
* #%L
* Image-Registration
* %%
* Copyright (C) 2019 Oliver Loeffler, Raumzeitfalle.net
* %%
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
* #L%
*/
package net.raumzeitfalle.registration.alignment;

import java.util.Collection;
import java.util.Iterator;

import Jama.Matrix;
import Jama.QRDecomposition;
import net.raumzeitfalle.registration.Dimension;
import net.raumzeitfalle.registration.Orientation;
import net.raumzeitfalle.registration.Orientable;

/**
*
* Calculates translation and rotation for the given set of equations. A rigid
* transform preserves all distances between points.
*
* <ul>
* <li>Translation</li>
* <li>Rotation</li>
* </ul>
*
* @author oliver
*
*/
final class RigidBodyJamaModel implements RigidBodyModel {

private void prepare(Collection<RigidModelEquation> equations, Matrix references, Matrix deltas, Orientation direction) {
int row = 0;
Iterator<RigidModelEquation> it = equations.iterator();
while (it.hasNext()) {
row = addEquation(references, deltas, row, it.next(), direction);
}
}

private int addEquation(Matrix references, Matrix deltas, int row, RigidModelEquation eq, Orientation direction) {

if (Orientation.X.equals(direction)) {
references.set(row, 0, eq.getXf());
references.set(row, 1, eq.getDesignValue());
}

if (Orientation.Y.equals(direction)) {
references.set(row, 0, eq.getYf());
references.set(row, 1, eq.getDesignValue());
}

if (Orientation.XY.equals(direction)) {
references.set(row, 0, eq.getXf());
references.set(row, 1, eq.getYf());
references.set(row, 2, eq.getDesignValue());
}

deltas.set(row, 0, eq.getDeltaValue());
row++;
return row;
}

private RigidTransform solve(Matrix references, Matrix deltas, Orientation direction) {

QRDecomposition qr = new QRDecomposition(references);

Matrix rInverse = qr.getR().inverse();
Matrix qTransposed = qr.getQ().transpose();
Matrix solved = rInverse.times(qTransposed).times(deltas);

return createTransform(solved,direction);
}

private RigidTransform createTransform(Matrix solved, Orientation direction) {

if (Orientation.X.equals(direction)) {
return SimpleRigidTransform
.with(solved.get(0, 0), 0.0, solved.get(1, 0));
}

if (Orientation.Y.equals(direction)) {
return SimpleRigidTransform
.with(0.0, solved.get(0, 0), solved.get(1, 0));
}

return SimpleRigidTransform
.with(solved.get(0, 0), solved.get(1, 0), solved.get(2, 0));
}

@Override
public <T extends Orientable> RigidTransform solve(Collection<RigidModelEquation> equations, Dimension<T> dimension) {

int rows = equations.size();
int cols = dimension.getDimensions()+1;

Orientation direction = dimension.getDirection();

Matrix deltas = new Matrix(rows, 1);
Matrix references = new Matrix(rows, cols);
prepare(equations, references, deltas, direction);

// escape here before singular matrix exception can be thrown
if (Orientation.XY.equals(direction) && equations.size() == 2) {
double tx = deltas.get(0, 0);
double ty = deltas.get(1, 0);
return SimpleRigidTransform.translation(tx, ty);
}

if (Orientation.X.equals(direction) && equations.size() == 1) {
double tx = deltas.get(0, 0);
return SimpleRigidTransform.shiftX(tx);
}

if (Orientation.Y.equals(direction) && equations.size() == 1) {
double ty = deltas.get(0, 0);
return SimpleRigidTransform.shiftY(ty);
}

return solve(references, deltas, direction);
}

}
/*-
* #%L
* Image-Registration
* %%
* Copyright (C) 2019 Oliver Loeffler, Raumzeitfalle.net
* %%
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
* #L%
*/
package net.raumzeitfalle.registration.alignment;

import java.util.Collection;
import java.util.Iterator;

import net.raumzeitfalle.registration.Dimension;
import net.raumzeitfalle.registration.Orientable;
import net.raumzeitfalle.registration.Orientation;
import net.raumzeitfalle.registration.solver.Deltas;
import net.raumzeitfalle.registration.solver.References;
import net.raumzeitfalle.registration.solver.Solution;
import net.raumzeitfalle.registration.solver.Solver;
import net.raumzeitfalle.registration.solver.SolverProvider;

/**
*
* Calculates translation and rotation for the given set of equations. A rigid
* transform preserves all distances between points.
*
* <ul>
* <li>Translation</li>
* <li>Rotation</li>
* </ul>
*
* @author oliver
*
*/
final class DefaultRigidBodyModel implements RigidBodyModel {

private void prepare(Collection<RigidModelEquation> equations, References references, Deltas deltas, Orientation direction) {
int row = 0;
Iterator<RigidModelEquation> it = equations.iterator();
while (it.hasNext()) {
row = addEquation(references, deltas, row, it.next(), direction);
}
}

private int addEquation(References references, Deltas deltas, int row, RigidModelEquation eq, Orientation direction) {

references.set(row, eq, direction);
deltas.set(row, eq);

row++;
return row;
}

private RigidTransform solve(References references, Deltas deltas, Orientation direction) {

Solver solver = new SolverProvider().getSolver();
Solution solution = solver.apply(references, deltas);

return createTransform(solution,direction);
}

private RigidTransform createTransform(Solution solved, Orientation direction) {

if (Orientation.X.equals(direction)) {
return SimpleRigidTransform
.with(solved.get(0), 0.0, solved.get(1));
}

if (Orientation.Y.equals(direction)) {
return SimpleRigidTransform
.with(0.0, solved.get(0), solved.get(1));
}

return SimpleRigidTransform
.with(solved.get(0), solved.get(1), solved.get(2));
}

@Override
public <T extends Orientable> RigidTransform solve(Collection<RigidModelEquation> equations, Dimension<T> dimension) {

int rows = equations.size();
int cols = dimension.getDimensions()+1;

Orientation direction = dimension.getDirection();

Deltas deltas = new Deltas(rows);
References references = new References(rows, cols);
prepare(equations, references, deltas, direction);

// escape here before singular matrix exception can be thrown
if (Orientation.XY.equals(direction) && equations.size() == 2) {
double tx = deltas.get(0);
double ty = deltas.get(1);
return SimpleRigidTransform.translation(tx, ty);
}

if (Orientation.X.equals(direction) && equations.size() == 1) {
double tx = deltas.get(0);
return SimpleRigidTransform.shiftX(tx);
}

if (Orientation.Y.equals(direction) && equations.size() == 1) {
double ty = deltas.get(0);
return SimpleRigidTransform.shiftY(ty);
}

return solve(references, deltas, direction);
}

}
Loading

0 comments on commit a2c3dfa

Please sign in to comment.