imagej / imagej-ops

ImageJ Ops: "Write once, run anywhere" image processing
https://imagej.net/libs/imagej-ops
BSD 2-Clause "Simplified" License
90 stars 42 forks source link

Ridge Detection #535

Closed gselzer closed 6 years ago

gselzer commented 6 years ago

This PR adds support for a Ridge Detection algorithm, inspired by Thorsten Wagner's original implementation and adapted from Carsten Steger's paper. The PR adds multiple ops intended to be used in a single workflow however separated to accommodate user preferences. The package includes a Ridge Detection algorithm, used to find bright "Curvilinear Structures" against a dark background, and a Junction Detection algorithm used to find junctions or near junctions (as determined by the threshold parameter) between a list of polylines. As of now the algorithm package is not yet complete, however this PR is being made for the forum discussion about the scientific need for the current and potentially future ops. Any input on the ops would be greatly appreciated there. Before this PR is merged the snapshot dependency on the Imglib2 ROIs will need to be sorted out as well.

gselzer commented 6 years ago

@imagejan Yes, you are right. I was trying to hammer out the main points before I left yesterday. I will make sure to make it separate from the UnaryOp.

@dietzc maybe I can elaborate on why @ctrueden and I decided on this method. Sometimes, as in Ridge Detection or in Frangi Vesselness it is necessary to compile the information of a bunch of different partial derivative images, for example, when you need to create a hessian matrix to find the eigenvectors and eigenvalues at that point. When this is the case you have the same input image, but you need to vary the int[] array that indicates the order of the partial derivative to be taken in each dimension. In the case of Ridge Detection the metadata generation requires five calls to this op, each with a different int array.

Now, we could just create a public void setPartials(int[] partials){} setter method, but this is not really done anywhere else in Ops, and this int array isn't really just used in calculations like other parameters are - this array is instead defining what the op does and how the op functions, more like the input image does than the sigma used for the gauss. This is why we chose to make it a BinaryFunctionOp instead of adding a setter method. Hope this helps.

gselzer commented 6 years ago

@imagejan after thinking about it for a while today and discussing with @ctrueden, I don't actually think it is necessary to have both a UnaryComputerOp and a BinaryComputerOp. I really could not think of an example where it would be helpful to have that UnaryComputerOp instead of the binary version. Furthermore, @ctrueden suggested that all BinaryComptuerOps can run as UnaryComputerOps as is. I will lean on his experience with Ops for this part of the explanation if you would like clarification, but between the two of us we decided that a singular BinaryComputerOp is the way to go.

gselzer commented 6 years ago

@hadim @ctrueden I compressed the Git history down to make it simpler and easier to read. All commits build and all tests pass on each commit.

hadim commented 6 years ago

@gselzer did you have a chance to look at a way to improve the performance of the implementation? Or do you think that would be too much work at the moment?

gselzer commented 6 years ago

@hadim when I refactored the Derivative Gaussian Op into a BinaryComputerOp I made an instance of that op in Ridge Detection's Metadata Generation phase, which brought down the number of times that the op went through Op matching down by 4. This should save some time, however I am not sure how much is being saved. With that being said I have not put any more work into improving the speed of this op. I might have a little more time to investigate this later.

ctrueden commented 6 years ago

I added a commit that improves the time performance by more than 2x on my system. There are still inefficiencies which would not be super hard to address. For example, the use of JAMA eigenvalue decomposition invokes new Matrix(...) internally—that code could be extracted and made allocation-free. There are likely also algorithmic aspects that could be improved, and use of sub-ops to optimize case by case when the image is e.g. an ArrayImg or PlanarImg. But the latter would be more effort. Already what is here is hopefully useful, so: merging it! 😁