Frangi Vesselness Filter Feedback

frangi
blood-vessel
imagej-ops
Tags: #<Tag:0x00007fd69c142650> #<Tag:0x00007fd69c142510> #<Tag:0x00007fd69c142380>

#1

Hello! I have worked to implement an ImageJ2 op implementation of the Frangi Vesselness Filter, inspired by Mark Longair’s ImageJ1 plugin and adapted from the original paper. This filter is not yet released in ImageJ2 but we would like to release it soon; for now, however, it is only a prototype. Before the filter is finished we would love to get any feedback and/or suggestions from those for whom this op might be relevant (@tferr?).

The filter is currently able to take in an image composed of lighter vessels upon a darker background (left) and is able to return an image of only the vessel-like structures (highlighted in green, superimposed on original image in red, right) (image retrieved from https://en.wikipedia.org/wiki/Magnetic_resonance_angiography#/media/File:Mra-mip.jpg):

As of now those who would like to try out the filter can do so as follows, using a terminal or Git Bash:

  1. Clone branch issue490 from imagej/imagej-ops
  2. Checkout the branch (issue490)
  3. Enter into the terminal: mvn -Dimagej.app.directory=PATH_TO_Fiji.app/
  4. Launch Fiji and use the Groovy script posted below. Make sure to use a grayscale input image and to open the image in ImageJ before running the Groovy script.

Here is the script to run the filter:

#@OpService op
#@UIService ui
#@Img input
#@String(label="scales: e.g. \"1, 2, 3\"",value="3, 5") scale
#@String(label="spacing: e.g. \"1, 1, 1\"",value="1, 1, 1") spacing
//scales refers the the pixel distance at which the filter operates. Larger scales measure larger structures.
//spacing refers to the physical distance between pixels. The default setting is {1, 1, 1...}
import java.util.Arrays
import net.imglib2.img.Img
import net.imglib2.img.array.ArrayImgs
import net.imglib2.type.numeric.real.DoubleType
import net.imglib2.type.numeric.integer.IntType
import net.imagej.ops.filter.vesselness.DefaultFrangi
//parse the strings

double[] scales = Arrays.stream(scale.split(",")).map{
	t -> Double.parseDouble(t.trim())
}.collect() as double[]
double[] spacings = Arrays.stream(spacing.split(",")).map{
	t -> Double.parseDouble(t.trim())
}.collect() as double[]
long[] dims = new long[input.numDimensions()]
input.dimensions(dims)
actualOutput = ArrayImgs.doubles(dims)
op.run(DefaultFrangi.class, actualOutput, input, spacings, scales)
ui.show("result", actualOutput)

We plan to merge this op by mid-September. In the meantime, we welcome all feedback and suggestions!


Issue on plugin: Multi-Scale Oriented Flux Tubularity Measure - Fiji crashes
#2

@gselzer this is really really nice.
Didn’t have much time to look into it, but tried on a demo image and looks great. Will try some more tests and will let you know.

One question: If I recall correctly Mark was skeptical of the accuracy of his implementation as it did not yield the same result of some other MATLAB code. Did you had a chance to look into it?

I had to use mvn -Denforcer.skip=true -DskipTests -Dimagej.app.directory=PATH_TO_Fiji.app/ otherwise it wouldn’t compile. (BTW: Compilation took forever: almost 30min!)


#3

It’s great that you have worked on this @gselzer! :thumbsup:

The results look promising for the 2D case, but somehow weird in 3D (although that might be on my side and due to anisotropic data). So I am trying to understand what the issue is with the following simple test data generated from an IJ1 macro:

newImage("Untitled", "8-bit black", 400, 100, 7);
setForegroundColor(255, 255, 255);

run("Specify...", "width=300 height=1 x=50 y=52 slice=2");
run("Fill", "slice");

run("Specify...", "width=300 height=3 x=50 y=51 slice=3");
run("Fill", "slice");

run("Specify...", "width=300 height=5 x=50 y=50 slice=4");
run("Fill", "slice");

run("Specify...", "width=300 height=3 x=50 y=51 slice=5");
run("Fill", "slice");

run("Specify...", "width=300 height=1 x=50 y=52 slice=6");
run("Fill", "slice");

run("Select None");

run("Gaussian Blur...", "sigma=1 stack");

If you apply the op to slice 4 with scale “5” and spacing “1,1” everything looks as expected. If the op is applied to the isotropic 3D data with scale “5” and spacing “1,1,1”, I get a stack with empty images (which I did not expect…).

Any idea?

Best,
Stefan


#4

Unfortunately, this is a known issue with the javac compiler acting upon the imagej-ops codebase. Compiling it in the Eclipse IDE does not have this problem. The most recent java 1.8.0_144 seems a bit faster—but still very slow—compared to older versions of 1.8.0. The problem has to do with interface multiple inheritance and use of default methods. The problem may become moot with the next iteration of ImageJ Ops, which will simplify the interface inheritance hierarchy for core ops interfaces. In the meantime, compiling from the CLI is indeed painful.


#5

I did get a chance to look through Mark’s code. For the most part it was correct but there were a few flaws in his implementation, which I worked to understand and fix with @awalter17.


#6

This is interesting as I did not take much of a look at testing 3D images. I wonder if it has something to do with the larger scale as, with such a small image, that scale will only allow the largest vessel structures through (not to mention 3 dimensional analysis should be pickier than 2 dimensional). Have you tried varying the scale at all to see if that is the issue? If that does not fix it I will take a look and see if I can determine the issue.

Thanks,

Gabe


#7

I did: still an all empty output for the test image. Even if I scale the image up to 70 z-slices I still get an empty image for scales between 1 and 10…

newImage("Untitled", "8-bit black", 400, 100, 7);
setForegroundColor(255, 255, 255);

run("Specify...", "width=300 height=1 x=50 y=52 slice=2");
run("Fill", "slice");

run("Specify...", "width=300 height=3 x=50 y=51 slice=3");
run("Fill", "slice");

run("Specify...", "width=300 height=5 x=50 y=50 slice=4");
run("Fill", "slice");

run("Specify...", "width=300 height=3 x=50 y=51 slice=5");
run("Fill", "slice");

run("Specify...", "width=300 height=1 x=50 y=52 slice=6");
run("Fill", "slice");

run("Select None");

run("Gaussian Blur...", "sigma=1 stack");

run("Scale...", "x=2 y=2 z=10 width=800 height=200 depth=70 interpolation=Bilinear average process create");

Thanks for looking into it, @gselzer!


#8

Sorry for the 2+ week delay in my answer, but I had to put the filter on hold while starting college!

After playing with various parameters I finally got results by ignoring the third eigenvalue. This particular value should be negative upon all vessels in the image but for some reason, the value returns positive for the whole vessel. I plan to run additional 3D images through the filter to see if I get the same result. The question, however, is whether or not ignoring this parameter is the okay thing to do, even though it gives results.

I pushed my code change to the issue490 branch if you would like to test again. I found that 1 was the only scale required for the unscaled version of the macro image . Although the output does not seem to return the entire vessel boundary, as you can see in most of the slices, but this is definitely a start of a solution.


#9

I looked at the results as well as your implementation today, @gselzer:

  1. It don’t see the reasoning for DefaultFrangi.java#L222, since the publication uses a value of S^2 not S^4. Maybe I am missing something, so just correct me if I am wrong
  2. The second part of the conditional in DefaultFrangi.java#L246 is not mentioned in the publication
  3. I guess the conditionals should be connected by || in DefaultFrangi.java#L269; at least according to the definition in the publication
  4. The results look different even in 2D, when I compare the example image that you have posted; the MATLAB implementation as well as the scikit-image implementation emphasize more regions when using the same scales

Gotta run now, but will post an update with images tomorrow!


#10

Alright, you can find the sources for my tests at https://github.com/stelfrich/frangi-vesselness.

Ops:

Scikit-Image:
result-skimage

Matlab:


#11
  1. Whoops. That was supposed to be deleted. When I was first ironing out the kinks in 2D I messed around with all of the different parameters to see their effect and get a better idea of how the filter works. I’ll remove that.

  2. While you are right in that the publication does not use this calculation, on the third page of the paper, there are three (or two, if you are working in 2 dimensions) characteristics of an ideal tubular structure identified. I took a snapshot of them for reference. The second characteristic is the one that you saw in the conditional. It is helpful to have in that conditional as it removes noise, as blob-like structures (and noisy pixels) will have similar values for l1 and l2, and we don’t want to let those pass through.

IdealTubeness

  1. Yes, according to the publication, my implementation is wrong; both l2 and l3 should technically be negative for a vessel (according to the publication, if either l2 or l3 are greater than 0, the vesselness should be 0). However, when I had the conditional as thus, we get the terrible black slices, and only by removing one of the two l-conditionals was I able to get some filtered data. I kept l3 over l2 because it made for a better image in my opinion. However, I am not yet happy with the 3D filter version, so consider it a WIP.

  2. Wow. Those look different. It looks as if MATLAB, which is the one I am currently comparing to the ops implementation, has a Gaussian of some sort applied to it. Mark’s filter also applied a Gaussian in the ImageJ1 plugin. @ctrueden and I decided against doing it in this implementation because not everyone might want it. If they do want a Gaussian on top of this filter they can always do so afterward. I’ll play around a little with the 2D parameters if you think they could use some fine-tuning; just let me know what you think it needs.

EDIT: I just pushed a branch fixing 1). The 3D results, when viewed in 3D Viewer, look pretty good. I edited your macro a little and got good results when viewing the filter on its own. Note that for some reason ImageJ thinks that the TIFF has time slices instead of a third dimension; if this is not changed in the image properties 3D Viewer does not work right.

setForegroundColor(255, 255, 255);

run("Specify...", "width=300 height=1 x=50 y=52 slice=2");
run("Fill", "slice");
run("Specify...", "width=1 height=50 x=50 y=52 slice=2");
run("Fill", "slice");

run("Specify...", "width=300 height=3 x=50 y=51 slice=3");
run("Fill", "slice");

run("Specify...", "width=3 height=50 x=50 y=51 slice=3");
run("Fill", "slice");

run("Specify...", "width=300 height=5 x=50 y=50 slice=4");
run("Fill", "slice");

run("Specify...", "width=5 height=50 x=50 y=50 slice=4");
run("Fill", "slice");

run("Specify...", "width=300 height=3 x=50 y=51 slice=5");
run("Fill", "slice");

run("Specify...", "width=3 height=50 x=50 y=51 slice=5");
run("Fill", "slice");

run("Specify...", "width=300 height=1 x=50 y=52 slice=6");
run("Fill", "slice");

run("Specify...", "width=1 height=50 x=50 y=52 slice=6");
run("Fill", "slice");

run("Select None");

run("Gaussian Blur...", "sigma=1 stack");

run("Scale...", "x=2 y=2 z=2 width=800 height=200 depth=14 interpolation=Bilinear average process create");selectWindow("Untitled");

Thanks for all your help in this!


#12

That’s true. I was just wondering why we actually need to explicitly check that those properties hold. From my understanding, those assumptions are implicitly encoded in the vesselness value (as rb and ra) and don’t need to be checked individually. All of the implementations (scikit-image, ITK, MATLAB) I had a look at are operating on smoothed versions of the input image and seem to be doing fine without those kind of checks…

I might not have thought that through properly, but I guess it’s not equivalent to what’s implemented in e.g. scikit-image. In that implementation (same for MATLAB and ITK) the images are Gauss filtered at each scale and the Hessian is computed on that filtered image instead of the input.


EDIT: That’s the ITK implementation I was talking about: https://itk.org/Doxygen/html/classitk_1_1FrangiTubularnessImageFilter.html


#13

You’re right. They shouldn’t. I just ran a few comparison tests and it did not change the output visually when I omitted that check from the conditional. I don’t know why I found that useful earlier but it does not seem useful now, so I will remove it from the conditional.

I am very confused by what you mean here. Are you saying that the whole image should be run through a gaussian before the filter runs? If so, I misunderstood the paper, and I can work to implement it. However, looking at the ITK and SciKit implementations, it doesn’t seem like they are doing anything drastically different from this implementation (in the SciKit they use sigma instead of scale, but they square their sigma, which is strange. I will check that out, though I cannot find anything in the paper to back that up. In the ITK implementation I cannot find any reference to a scale or sigma at all, though I am not familiar with C++, so maybe I missed it). If you are seeing anything that I might be missing let me know.

I really appreciate your help in all of this.

Gabe


#14

Looking at the scikit-image implementation, the input image is Gaussian filtered for each scale (corner.py#L159-L160) and from this filtered version a Hessian is computed for each pixel (corner.py#L174-L181).

For the ITK implementation, I looked at an example of how to apply a filter like the Frangi one. In that example, an HessianRecursiveGaussianImageFilter is applied to the input, which uses - according to the documentation - a Gaussian-filtered version of the input image to compute Hessians per pixel. Subsequently, the Hessian per pixel can be fed to a FrangiTubularnessImageFilter which computes the vesselness in an identical fashion to your implementation.

Sorry for holding you up for so long!

Best,
Stefan


#15

Aha. You are right. I did not see this when I looked through. It does appear like they do this.

However, @ctrueden and I talked for a while about the options for this filter, since my current implementation would struggle when we add this Gauss. Currently we input various values through scales parameter, and the filter uses these values to filter the images with differing ranges of sensitivity. This is how we were able to get an image that returns both small and large vessels. If we add this Gauss in, however, then the parameter passing would become complicated as users would have to input a 2D array, of (ND)*y size, containing y scales, each with N sigmas to control the Gauss at each stage of filtering. Furthermore @ctrueden and I discussed as to whether or not the flexibility of the op would be hampered by forcing a Gauss in, as maybe some users would find better results without the Gauss.

We decided that a Command might be the better option in this scenario, since then we can let users decide what they want. We can default all of the paper standards (i.e. the Gauss beforehand), but then we can give them the option to run at multiple scales and/or not use the Gauss, if they so choose. It won’t hamper time either since this was all managed through for loops beforehand. Not to mention it will look nice and the GUI might let more people on to its release :smile:

Any thoughts on this?

Gabe


#16

That all sounds good to me! Would this involve refactoring the Op to operate on one scale and do the handling of multiple scales in the Command? This way, people would be able to do something like this in scripts to emulate the behavior of known implementations while using the op itself will give you the most flexibility:

// ... get input ...
for (scale in 1..5) {
    inputGaussFiltered = ops.filter().gauss(input, scale);
    ops.filter().vesselness(inputGaussFiltered, spacing, scale);
}

Also, we might want to rename the Op to FrangiVesselness since there are also other vesselness measures (e.g. in ITK).

Best,
Stefan


#17

This was pretty much our plan.

Right. The class housing the command is FrangiVesselness, and it shows up in the ImageJ UI as Frangi Vesselness. The op itself is called DefaultFrangi and is housed in the package net.imagej.ops.filter.vesselness. Is there anything that I am missing?

Gabe


#18

Awesome!

Sorry, wasn’t clear about that. It’s actually the method name in the FilterNamespace that’s just vesselness. I’d go for frangiVesselness there as well…

Thanks again for working on this @gselzer!


#19

Got it. Knew I was missing something.

The Command is finished, and there are a few things that I have done with this:

  • The Gauss is run at every scale AS LONG AS the user does not uncheck the Gauss Filter box.
  • I removed the sigmas parameter, since we can just use the scale that we are currently at as the sigma (I am 95% sure that this is what they do in equation 2 of the paper).
  • The return is now a stack of images, consisting of the filter returns at each scale value. If the user wants they can then combine them.

Using scales 1, 2, and 3 in the command and then using the auto settings of the B&C tool I was able to generate this image:

The one thing I am still a little confused about in relationship to the paper is Equation 2. I’m not very experienced with scale space theory and as such am pretty confident, but not sure, that this Gauss/scale/normalized derivatives thing is done correctly. I think this only relates to the Gauss application that we have been discussing, but it would be good to make sure that we are not missing something here.

I am going to post something about this onto the Gitter and see if I can get any more feedback, especially on the command and maybe on Equation 2, but unless you have any other concerns I think we are good to go!

Thanks,

Gabe