Add Hessian eigenvalues #29

hanslovsky · 2016-08-08T21:19:55Z

This pull request features:

Computation of Hessian matrices (second derivatives)
Computation of eigenvalues for square and symmetric second order tensors in n-dimensional spaces.

For this purpose, these classes and methods were added:

Package net.imglib2.algorithm.corner:
- RealComposite*Matrix: Implementation of apache RealMatrix as a View on imglib RealComposite to avoid unnecessary allocation and copy
- HessianMatrix: Static methods to compute Hessian matrix from n-dimensional input.
- TensorEigenValues: Static methods to compute tensor from n+1-dimensional input, where the last dimension stores linear representation of the tensor.
Package net.imglib2.algorithm.gradient:
- PartialDerivative: Added static method gradientCentralDifferenceParallel to speed up gradient computation.

Before merging this, please comment on the package structure:

Should RealComposite*Matrix get their own package?
Should TensorEigenValues be in a package of its own, e.g. net.imglib2.algorithm.tensor?
Should the corner package be renamed?

Once the pull request is ready to be merged, I will rebase in my local repository and create a new pull request.

- Unified RealCompositeMatrix classes - eigenvalues for non-symmetric matrices possible now

That's what the interface requires: http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/linear/RealMatrix.html#copy()

…uesTest e7b31ff added separate test for Hessian matrix making the existing consistency test obsolete.

@tpietzsch

Use [email protected] as contact. Change if requested by @tpietzsch.

hanslovsky · 2017-02-08T22:33:50Z

@tinevez's commit tinevez/CircleSkinner@84e4569 reminded me of this outstanding pull request. I should have some time next week to clean-up and simplify the interface. Specifically, I will change the interface such that all inputs must be supplied by the caller. That might be as extreme as having an interface that accepts only the gradient image. That way the number of overloads would be greatly reduced. Also, I would like to get rid of one of

func( args..., int nThreads, ExecutorService es)

and

func( args..., int nThreads )

Currently, I lean towards keeping the first option and removing the second one. Another option would be to get rid of the parallel computing in its entirety as this can be done by the caller by calculating the Hessian on Chunks of data. That would greatly reduce the number of overloads.

Thinking about the ImageJ picture, this should probably go into ops @ctrueden @dietzc . If you agree, I will write an op once this is merged in imglib2-algorithm. I will approach you with questions once I get there.

dietzc · 2017-02-08T22:41:23Z

Hi @hanslovsky,

I'd keep the variant with ExecutorService which is a pattern used in other imglib2-algorithms, too (for example DoG). Having a wrapper in Ops would of course be great and I'm happy help to help you with that!

tinevez · 2017-02-09T10:02:03Z

Hi @hanslovsky ,
Great news! I am sure that at the very least Hessian and partial derivatives will be a of great interest for people developing algorithms.

Here is a question regarding the Hessian calculation.
In the HessianMatrix class, you offer the possibility to specify a non-isotropic smoothing. This is great because then we accommodate e.g. 3D images for which the Z pixel size is different ant the X and Y pixel size (95% of confocal images).

But for the Hessian to be calculated correctly, the partial derivatives should be scaled accordingly too right?

hanslovsky · 2017-02-16T17:06:20Z

@dietzc What I don't like about the ExecutorService pattern is that we have to somehow arbitrarily pick the number of tasks, e.g. as in DoG:

// FIXME find better heuristic?
final int numThreads = Runtime.getRuntime().availableProcessors();
final int numTasks = numThreads <= 1 ? 1 : numThreads * 20;

I will keep the ExecutorService pattern but add a nTasks parameter so the caller has to decide how many tasks he'd like to run.

@tinevez I will have to think about that!

hanslovsky · 2017-02-22T13:42:48Z

@tinevez I think the appropriate way to get partial derivatives right is to re-sample the (smoothed) input image for isotropy. I would leave it to the caller to do that. What I add instead is a View that takes a Hessian matrix as input and scales the second derivatives according to a sigma parameter (the same you would use for smoothing), something like this:

public static < T extends RealType< T > > AbstractConvertedRandomAccessible< T, T > scaleHessian( final RandomAccessible< T > hessian, final T t, final double... sigma )
{
	final double minSigma = Arrays.stream( sigma ).reduce( Double.MAX_VALUE, ( d1, d2 ) -> Math.min( d1, d2 ) );
	final double minSigmaSq = minSigma * minSigma;
	final double[] scales = new double[ sigma.length * ( sigma.length + 1 ) / 2 ];
	for ( int i1 = 0, k = 0; i1 < sigma.length; ++i1 )
		for ( int i2 = i1; i2 < sigma.length; ++i2, ++k )
			scales[ k ] = sigma[ i1 ] * sigma[ i2 ] / minSigmaSq;
		final AbstractConvertedRandomAccessible< T, T > wrapped = new AbstractConvertedRandomAccessible< T, T >( hessian )
	{
		@Override
		public AbstractConvertedRandomAccess< T, T > randomAccess()
		{
			return new AbstractConvertedRandomAccess< T, T >( hessian.randomAccess() )
			{
				private final T tt = t;
				@Override
				public T get()
				{
					tt.set( source.get() );
					tt.mul( scales[ source.getIntPosition( sigma.length ) ] );
					return tt;
				}

				@Override
				public AbstractConvertedRandomAccess< T, T > copy()
				{
					// TODO Auto-generated method stub
					return null;
				}
			};
		}

		@Override
		public AbstractConvertedRandomAccess< T, T > randomAccess( final Interval interval )
		{
			return randomAccess();
		}
	};
	return wrapped;
}

Note that the scaling factors are normalized by the smallest sigma squared. The scaling has to be quadratic because the hessian matrix contains second derivatives. This should have minimal overhead and you could even write the contents of the view back into the hessianMatrix in a single pass without allocation of extra memory. Let me know what you think.

Non-symmetric matrices may have complex eigenvalues. This interface change is thus required for compatibility with use cases other than Hessian matrix (symmetric) eigenvalues.

…ompute Now, matrix based EigenValue implementations hold a RealCompositeMatrix object with exchangable data Composite instead of creating an instance at each call of compute. As a consequence EigenValues are not thread safe in general and a copy method was added to the interface.

Replace the checks by asserts and if the caller has input that might be out of bounds the check should be done outside getEntry.

…rom RealComposite to Composite. We do not care about the kind of composite that is passed. We only need to call Composite.get.

hanslovsky · 2017-02-24T21:59:35Z

I added a convenience method to scale the hessian matrix based on a sigma[] parameter passed by the caller (cf the last commit @tinevez ). I also added the missing license headers (how did I miss them?)

Is this ready to be merged now?

tinevez · 2017-04-03T13:35:08Z

I could use this PR, once merged, in imagej-ops to port the tubenessop.

…-ojalgo

ojAlgo v43 was just released so we can migrate to ojAlgo.

hanslovsky · 2017-04-24T16:28:26Z

ojAlgo v43 was just released: optimatika/ojAlgo@5487fa8
I will update the tensor eigenvalues and matrices to use ojAlgo as backend instead of commons-math and then we should be able to merge this (after further review).

Update: Once this is confirmed for a merge I will make a rebase into a single commit.

Hessian eigenvalues ojalgo

Eclipse (neon.2) seems to have issues indenting parameters following generic wildcards.

axtimwalde · 2017-05-15T16:23:52Z

@/all Is this ready for a merge?

hanslovsky · 2017-05-15T19:15:11Z

It is from my side (after a rebase, of course).

tpietzsch · 2017-05-16T10:22:28Z

I'm looking at it now

tpietzsch

I left a few inline suggestions/comments, but from my side this is in principle ready for merging.

Cool that you are using ojAlgo.

In general, I would be interested in the trade-offs, with respect to performance and expressiveness, of Vector/Matrix NativeTypes vs. the Composite approach. Unfortunately I have no time to look into it now...

tpietzsch · 2017-05-16T11:43:57Z

src/main/java/net/imglib2/algorithm/gradient/PartialDerivative.java

+		final int nDim = source.numDimensions();
+		if ( nDim < 2 )
+		{
+			gradientCentralDifference2( source, gradient, dimension );


Why are you calling gradientCentralDifference2, not gradientCentralDifference in this case?
Is it faster? Did you test?

No particular reason: I will change that to gradientCentralDifference to be consistent with what happens in HessianMatrix and use the faster version (according to the comments in PartialDerivative.

tpietzsch · 2017-05-16T12:24:00Z

src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java

+		assert sigma.length * ( sigma.length + 1 ) / 2 == hessian.dimension( sigma.length );
+		final int maxD = sigma.length;
+
+		final double minSigma = Arrays.stream( sigma ).reduce( Double.MAX_VALUE, ( d1, d2 ) -> Math.min( d1, d2 ) );


final double minSigma = Util.min( sigma ) would be more readable. (net.imglib2.util.Util)

tpietzsch · 2017-05-16T12:25:52Z

src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java

+		return hessianMatrix;
+	}
+
+	public static < T extends RealType< T > > IntervalView< T > scaleHessianMatrix( final RandomAccessibleInterval< T > hessian, final double[] sigma )


Can you add javadoc? It's not obvious what this method does

tpietzsch · 2017-05-16T13:19:54Z

src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeMatrix.java

+		this( data, nRows, nCols, nRows * nCols );
+	}
+
+	public RealCompositeMatrix( final Composite< T > data, final int nRows, final int nCols, final int length )


I would make this constructor (and all constructors of derived classes having a length argument) protected. The constructor only makes sense for being able to derive classes. In all other scenarios the length element is redundant.

expectedLength() can be removed then.

tpietzsch · 2017-05-16T13:26:14Z

src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java

+ * @author Philipp Hanslovsky
+ *
+ */
+public class HessianMatrix


I wonder whether it is worth adding a version that takes gaussian and gradient arguments that have appropriate size so that an OutOfBounds factory is not required. This would be useful if you want to do the Hessian computation blockwise to avoid artifacts at block boundaries. Of course, this can also be added later when it is needed.

I have not thought about that but I like the idea. I did not like the extra OufOfBounds so we should keep that in our minds. I have only one concern which actually might be a non-issue: If you already have the gradients pre-computed, you can just call the overload takes only the gradients as input. If you do not have them pre-computed, you will probably want to re-use them as features elsewhere and you might want to have them at the same size as the input.

Yes, I fully agree

tpietzsch · 2017-05-16T13:26:43Z

src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java

+	 *            n+1-dimensional {@link RandomAccessibleInterval} for storing
+	 *            the gradients along all axes of the smoothed source (size of
+	 *            last dimension is n)
+	 * @param hessianMatrix


Can you add a javadoc description of the upper triangular layout of the hessianMatrix (can be adapted from TensorEigenValues javadoc)?

tpietzsch · 2017-05-16T16:54:50Z

src/main/java/net/imglib2/algorithm/linalg/eigen/TensorEigenValues.java

+		{
+			// TODO Auto-generated catch block
+			e.printStackTrace();
+		}


e.printStackTrace() is probably not what should happen here.
This is a problem also in other algorithms (e.g., DifferenceOfGaussian, SeparableSymmetricConvolution, ...).
So I would merge this anyway and open an issue to fix it in all places.

For InterruptedException the right thing to do would be to restore the interruption status of the current thread and return.

For ExecutionException I don't know ...

@tpietzsch

Quoting @tpietzsch "The constructor only makes sense for being able to derive classes. In all other scenarios the length element is redundant. expectedLength() can be removed then."

@tpietzsch

@tpietzsch rightfully mentioned that it was not clear what scaleHessianMatrix does. Added some JavaDoc for clarification.

…sian-eigenvalues

hanslovsky · 2017-05-16T19:36:42Z

#41 is the same pull request but from a rebased branch to have a single commit.

hanslovsky added 18 commits August 3, 2016 14:06

added first version of hessian eigenvalues

e1678fb

implemented copy method

032c525

Changes to Matrix classes and eigenvalue methods

fda376d

- Unified RealCompositeMatrix classes - eigenvalues for non-symmetric matrices possible now

copy() returning deep copy now

ee213c3

That's what the interface requires: http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/linear/RealMatrix.html#copy()

added methods for passing result RandomAccesisbleIntervals

17dded2

Added parallel version for gradient.

9aecda0

Use parallel gradient, pass number of threads

fbe5282

removed time logs

9b76979

added parrallel option to eigenvalues

b96cd39

added author tags (PartialDerivative unclear)

00fc4a9

Merge branch 'master' into hessian-eigenvalues

a515974

added documentation/comments

dc8c079

Add HessianMatrixEigenValuesTest

3074bee

Add test for HessianMatrix (2D and 3D)

e7b31ff

Remove consistency test for Hessian matrix from HessianMatrixEigenVal…

3b7b075

…uesTest e7b31ff added separate test for Hessian matrix making the existing consistency test obsolete.

Add @tpietzsch as author of PartialDerivative.java

7fa9b69

Use [email protected] as contact. Change if requested by @tpietzsch.

Merge branch 'master' into hessian-eigenvalues

d5871fd

removed e-mails from author tags

b72af5f

dietzc mentioned this pull request Aug 24, 2016

Sobel and hessian filter imagej/imagej-ops#349

Merged

stelfrich mentioned this pull request Sep 8, 2016

MatrixTypes imglib/imglib2#137

Open

hanslovsky added 6 commits February 16, 2017 14:02

Merge remote-tracking branch 'origin/master' into hessian-eigenvalues

80a2b9e

Simplify interfaces

a7dd4d8

Rename packages (more appropriately)

5273d3c

Use nTasks instead of nThreads

c925680

Add convenience methods for gaussian smoothing to Hessian matrix

58f24bd

Remove number of threads from nTasks JavaDoc

a4e427c

hanslovsky added 6 commits February 22, 2017 11:49

Change EigenValue interface to accept ComplexType

0ff0e19

Non-symmetric matrices may have complex eigenvalues. This interface change is thus required for compatibility with use cases other than Hessian matrix (symmetric) eigenvalues.

Remove bounds check for RealCompositeMatrix.getEntry

1cc9f93

Replace the checks by asserts and if the caller has input that might be out of bounds the check should be done outside getEntry.

Relax interface for RealCompositeMatrix.data and EigenValue.compute f…

7f928eb

…rom RealComposite to Composite. We do not care about the kind of composite that is passed. We only need to call Composite.get.

Add license header.

6489ec5

Add convenience method for scaling Hessian according to sigma.

de07cf4

Add eigenvalue decomposition using ojAlgo

fca300c

hanslovsky mentioned this pull request Mar 8, 2017

Re-use eigenvalue/-vector storage arrays in RawEigenvalue.doDecomposeGeneral when applicable optimatika/ojAlgo#33

Closed

hanslovsky added 2 commits April 24, 2017 11:47

Merge remote-tracking branch 'origin/master' into hessian-eigenvalues…

6d91ca2

…-ojalgo

Use ojAlgo instead of commons-math

0b6b80f

ojAlgo v43 was just released so we can migrate to ojAlgo.

hanslovsky and others added 2 commits April 24, 2017 13:44

Merge pull request #1 from hanslovsky/hessian-eigenvalues-ojalgo

ba28c27

Hessian eigenvalues ojalgo

Fix formatting issues

f709b25

Eclipse (neon.2) seems to have issues indenting parameters following generic wildcards.

tpietzsch reviewed May 16, 2017

View reviewed changes

hanslovsky added 7 commits May 16, 2017 13:54

Use gradientCentralDifference instead of gradientCentralDifference2

cbfc6b5

Use more concise and readable Util.min over streams.

c6a1ad0

Make constructors with length parameters protected

1e85a8c

Quoting @tpietzsch "The constructor only makes sense for being able to derive classes. In all other scenarios the length element is redundant. expectedLength() can be removed then."

Remove unused import java.util.Arrays

0d6e352

Specify linear representation of upper triangular Hessian in JavaDoc

1b63be8

Add JavaDoc to HessianMatrix.scaleHessianMatrix

7366220

@tpietzsch rightfully mentioned that it was not clear what scaleHessianMatrix does. Added some JavaDoc for clarification.

Merge branch 'master' of github.com:imglib/imglib2-algorithm into hes…

c7608c4

…sian-eigenvalues

hanslovsky mentioned this pull request May 16, 2017

Add HessianMatrix and TensorEigenValues #41

Merged

hanslovsky closed this May 16, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add Hessian eigenvalues #29

Add Hessian eigenvalues #29

hanslovsky commented Aug 8, 2016 •

edited

Loading

hanslovsky commented Feb 8, 2017

dietzc commented Feb 8, 2017

tinevez commented Feb 9, 2017

hanslovsky commented Feb 16, 2017

hanslovsky commented Feb 22, 2017 •

edited

Loading

hanslovsky commented Feb 24, 2017

tinevez commented Apr 3, 2017

hanslovsky commented Apr 24, 2017 •

edited

Loading

axtimwalde commented May 15, 2017

hanslovsky commented May 15, 2017

tpietzsch commented May 16, 2017

tpietzsch left a comment

tpietzsch May 16, 2017

hanslovsky May 16, 2017

tpietzsch May 16, 2017

tpietzsch May 16, 2017

tpietzsch May 16, 2017

tpietzsch May 16, 2017

hanslovsky May 16, 2017

tpietzsch May 17, 2017

tpietzsch May 16, 2017

tpietzsch May 16, 2017

tpietzsch May 16, 2017

hanslovsky commented May 16, 2017

Add Hessian eigenvalues #29

Add Hessian eigenvalues #29

Conversation

hanslovsky commented Aug 8, 2016 • edited Loading

hanslovsky commented Feb 8, 2017

dietzc commented Feb 8, 2017

tinevez commented Feb 9, 2017

hanslovsky commented Feb 16, 2017

hanslovsky commented Feb 22, 2017 • edited Loading

hanslovsky commented Feb 24, 2017

tinevez commented Apr 3, 2017

hanslovsky commented Apr 24, 2017 • edited Loading

axtimwalde commented May 15, 2017

hanslovsky commented May 15, 2017

tpietzsch commented May 16, 2017

tpietzsch left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hanslovsky commented May 16, 2017

hanslovsky commented Aug 8, 2016 •

edited

Loading

hanslovsky commented Feb 22, 2017 •

edited

Loading

hanslovsky commented Apr 24, 2017 •

edited

Loading