Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update to ejml-0.41 #1

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Update to ejml-0.41 #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9dba05f
Add minimal test for Ellipse
imagejan Oct 18, 2023
72d165d
[wip] Update to ejml-0.41
imagejan Oct 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 19 additions & 5 deletions pom.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,13 +55,22 @@

<!-- NB: Deploy releases to the SciJava Maven repository. -->
<releaseProfiles>sign,deploy-to-scijava</releaseProfiles>

<ejml.version>0.41</ejml.version>
</properties>

<dependencies>
<dependency>
<groupId>net.imglib2</groupId>
<artifactId>imglib2-ij</artifactId>
</dependency>
<exclusions>
<exclusion>
<!-- declare the exclusion here -->
<groupId>jitk</groupId>
<artifactId>jitk-tps</artifactId>
</exclusion>
</exclusions>
</dependency>
<dependency>
<groupId>mpicbg</groupId>
<artifactId>mpicbg</artifactId>
Expand All @@ -76,10 +85,15 @@
<version>0.1</version>
</dependency>
<dependency>
<groupId>
com.googlecode.efficient-java-matrix-library
</groupId>
<artifactId>ejml</artifactId>
<groupId>org.ejml</groupId>
<artifactId>ejml-all</artifactId>
<version>${ejml.version}</version>
</dependency>

<dependency>
<groupId>junit</groupId>
<artifactId>junit</artifactId>
<scope>test</scope>
</dependency>
</dependencies>

Expand Down
53 changes: 27 additions & 26 deletions src/main/java/fit/circular/Ellipse.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,12 +36,13 @@
import java.util.ArrayList;
import java.util.Collection;

import org.ejml.data.DenseMatrix64F;
import org.ejml.factory.DecompositionFactory;
import org.ejml.factory.LinearSolverFactory;
import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.factory.DecompositionFactory_DDRM;
import org.ejml.dense.row.factory.LinearSolverFactory_DDRM;
import org.ejml.interfaces.decomposition.EigenDecomposition;
import org.ejml.interfaces.linsol.LinearSolver;
import org.ejml.ops.CommonOps;
import org.ejml.interfaces.linsol.LinearSolverDense;
import org.ejml.dense.row.CommonOps_DDRM;

import fit.util.TransformUtil;
import ij.ImageJ;
Expand Down Expand Up @@ -213,26 +214,26 @@ public void fitFunction( final Collection< Point > points )
// @author Peter Abeles (released under http://www.apache.org/licenses/LICENSE-2.0)

// qudratic part of design matrix
final DenseMatrix64F D1 = new DenseMatrix64F( 3, 1 );
final DMatrixRMaj D1 = new DMatrixRMaj( 3, 1 );
// linear part of design matrix
final DenseMatrix64F D2 = new DenseMatrix64F( 3, 1 );
final DMatrixRMaj D2 = new DMatrixRMaj( 3, 1 );

// quadratic part of scatter matrix
final DenseMatrix64F S1 = new DenseMatrix64F( 3, 3 );
final DMatrixRMaj S1 = new DMatrixRMaj( 3, 3 );
// combined part of scatter matrix
final DenseMatrix64F S2 = new DenseMatrix64F( 3, 3 );
final DMatrixRMaj S2 = new DMatrixRMaj( 3, 3 );
// linear part of scatter matrix
final DenseMatrix64F S3 = new DenseMatrix64F( 3, 3 );
final DMatrixRMaj S3 = new DMatrixRMaj( 3, 3 );
// Reduced scatter matrix
final DenseMatrix64F M = new DenseMatrix64F( 3, 3 );
final DMatrixRMaj M = new DMatrixRMaj( 3, 3 );

// storage for intermediate steps
final DenseMatrix64F T = new DenseMatrix64F( 3, 3 );
final DenseMatrix64F Ta1 = new DenseMatrix64F( 3, 1 );
final DenseMatrix64F S2_tran = new DenseMatrix64F( 3, 3 );
final DMatrixRMaj T = new DMatrixRMaj( 3, 3 );
final DMatrixRMaj Ta1 = new DMatrixRMaj( 3, 1 );
final DMatrixRMaj S2_tran = new DMatrixRMaj( 3, 3 );

final LinearSolver< DenseMatrix64F > solver = LinearSolverFactory.linear( 3 );
final EigenDecomposition< DenseMatrix64F > eigen = DecompositionFactory.eig( 3, true, false );
final LinearSolverDense<DMatrixRMaj> solver = LinearSolverFactory_DDRM.linear( 3 );
final EigenDecomposition< DMatrixRMaj > eigen = DecompositionFactory_DDRM.eig( 3, true, false );

final int N = points.size();

Expand All @@ -256,23 +257,23 @@ public void fitFunction( final Collection< Point > points )
}

// Compute scatter matrix
CommonOps.multTransA( D1, D1, S1 ); // S1 = D1'*D1
CommonOps.multTransA( D1, D2, S2 ); // S2 = D1'*D2
CommonOps.multTransA( D2, D2, S3 ); // S3 = D2'*D2
CommonOps_DDRM.multTransA( D1, D1, S1 ); // S1 = D1'*D1
CommonOps_DDRM.multTransA( D1, D2, S2 ); // S2 = D1'*D2
CommonOps_DDRM.multTransA( D2, D2, S3 ); // S3 = D2'*D2

// for getting a2 from a1
// T = -inv(S3)*S2'
if ( !solver.setA( S3 ) )
throw new IllDefinedDataPointsException( "Could not fit ellipse, failed at T = -inv(S3)*S2'" );

CommonOps.transpose( S2, S2_tran );
CommonOps.changeSign( S2_tran );
CommonOps_DDRM.transpose( S2, S2_tran );
CommonOps_DDRM.changeSign( S2_tran );
solver.solve( S2_tran, T );

// Compute reduced scatter matrix
// M = S1 + S2*T
CommonOps.mult( S2, T, M );
CommonOps.add( M, S1, M );
CommonOps_DDRM.mult( S2, T, M );
CommonOps_DDRM.add( M, S1, M );

// Premultiply by inv(C1). inverse of constraint matrix
for (int col = 0; col < 3; col++)
Expand All @@ -289,12 +290,12 @@ public void fitFunction( final Collection< Point > points )
if ( !eigen.decompose( M ) )
throw new IllDefinedDataPointsException( "Could not fit ellipse, failed at eigen.decompose( M )" );

final DenseMatrix64F a1 = selectBestEigenVector( eigen );
DMatrixRMaj a1 = selectBestEigenVector( eigen );
if ( a1 == null )
throw new IllDefinedDataPointsException( "Could not fit ellipse, could not find best eigenvector." );

// ellipse coefficients
CommonOps.mult( T, a1, Ta1 );
CommonOps_DDRM.mult( T, a1, Ta1 );

this.a = a1.data[ 0 ];
this.b = a1.data[ 1 ] / 2;
Expand Down Expand Up @@ -381,14 +382,14 @@ public void intersectsAt( final double[] p, final double[] i )
}
}

protected DenseMatrix64F selectBestEigenVector( final EigenDecomposition< DenseMatrix64F > eigen )
protected DMatrixRMaj selectBestEigenVector( final EigenDecomposition< DMatrixRMaj > eigen )
{
int bestIndex = -1;
double bestCond = Double.MAX_VALUE;

for (int i = 0; i < eigen.getNumberOfEigenvalues(); i++)
{
final DenseMatrix64F v = eigen.getEigenVector( i );
final DMatrixRMaj v = eigen.getEigenVector( i );

if ( v == null ) // TODO WTF?!?!
continue;
Expand Down
52 changes: 52 additions & 0 deletions src/test/java/fit/circular/EllipseTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
package fit.circular;

import org.junit.Test;

import mpicbg.models.IllDefinedDataPointsException;
import mpicbg.models.NotEnoughDataPointsException;

import static org.junit.Assert.assertTrue;


public class EllipseTest {

@Test
public void test() throws NotEnoughDataPointsException, IllDefinedDataPointsException {
final double[][] points = new double[ 8 ][ 2 ];

points[ 0 ][ 0 ] = 320;
points[ 0 ][ 1 ] = 443;

points[ 0 ][ 0 ] = 0; // non-axis-aligned ellipse
points[ 0 ][ 1 ] = 17;

points[ 1 ][ 0 ] = 377;
points[ 1 ][ 1 ] = 377;

points[ 2 ][ 0 ] = 507;
points[ 2 ][ 1 ] = 350;

points[ 3 ][ 0 ] = 640;
points[ 3 ][ 1 ] = 378;

points[ 4 ][ 0 ] = 694;
points[ 4 ][ 1 ] = 444;

points[ 5 ][ 0 ] = 639;
points[ 5 ][ 1 ] = 511;

points[ 6 ][ 0 ] = 508;
points[ 6 ][ 1 ] = 538;

points[ 7 ][ 0 ] = 376;
points[ 7 ][ 1 ] = 511;

//final ShapePointDistanceFactory< Ellipse, ?, ? > factory = new BruteForceShapePointDistanceFactory< Ellipse >( 0.001 );
final ShapePointDistanceFactory< Ellipse, ?, ? > factory = new EllipsePointDistanceFactory( 10 );

final Ellipse ellipse = new Ellipse( factory );
ellipse.fitFunction( Ellipse.toPoints( points ) );

assertTrue(ellipse.isEllipse());
}
}