Using CLIJ to remove outliers / noise

Dear CLIJ users,

I was searching for a function to remove the outliers (such as dead/bright pixels) in the GPU, similar to the one that is available on FIJI.
Is there anyway of replicating it using CLIJ commands?

Many thanks to all!

1 Like

Dear @anpl,

this particular functions is not implemented in CLIJ yet. But for noise removal, a median filter with a small radius (e.g. 1) might do the job. Would you mind trying?

Furthermore, I added your idea to the wishist I’ll keep you posted if I find a volunteer implementing it :wink:

Cheers,
Robert

1 Like

Hello @haesleinhuepf,

Thank you for the reply.
Yes, I saw it and thought of implementing it. I will probably use it if there is no other solution.

I am a happy volunteer. :smiley:
I did some gpgpu coding in the past (cuda), which is pretty similar (coding structure) to opencl.
Tell me what can I do to start testing.

Cheers!

Hey @anpl,

do I read you right, you want to help programming it? That is so awesome! I’m happy to give you a quick tour. First of all, CLIJ2 plugins are also just Fiji plugins. You can learn in this video how to build and maintain them:

In the video, Carina shows how to do this with Eclipse. I personally use IntelliJ. But it may not matter. Furthermore, you see in the video how an ImageJ plugin is built. Of course, you want to build a CLIJ2 plugin. Thus, instead of downloading the repository mentioned in the video, download this one:

The rest of the procedure should be the same! I would say after you managed building this Template Plugin and tested it successfully, we chat again on how to proceed.

You may in the meantime though have a look at the RemoveOutliers code. As you have experience with CUDA already, you know that when translating this to OpenCL, we basically need to get rid of some of the for-loops to parallelize it. But again, we can have a look at details at the right time. Furthermore, I just CC @axtimwalde who wrote the original RemoveOutliers code. Maybe he wants to point us to relevant resources. Maybe there is an imglib2 version somewhere?

Let me know how it goes! I’m happy to assist!

Cheers,
Robert

1 Like

The only additional resource that could be helpful in addition to the open source code is may be this page

summarizing all the integral image methods used to calculate STD and MEAN in very large boxes around each individual pixel quickly. This is to identify the outlier pixels. The code to fill them in is iterative diffusion which is interesting for larger regions but not so much for individual pixels.

2 Likes

Oh, keep in mind that this is not ImageJ’s remove outliers filter (which is very slow), but the Fiji plugin Integral Image Filters > Remove Outliers (which is very fast).

2 Likes

Hi @haesleinhuepf!

So, I’m also interested in a faster version of ImageJ’s Remove Outliers. I have tried the Integral Image Filters version, but it does not particularly do what I want it to.

Essentially, I am trying to remove deep valleys along a line profile while keeping small, high frequency peaks intact. I essentially use the largest radius in the ImageJ Remove Outliers, 0.0 threshold, and choose the Dark option. This works fantastic but takes forever!

Any suggestions in doing this in a roundabout way in CLIJ2?

Could you share an example before/after image? :upside_down_face:

1 Like

Yep! I’ll work up an example and try to post it sooner rather than later.

Calcite-1_OriginalImage_16bits-U_895x898x1516_45.26um0000.tif (1.5 MB)

@haesleinhuepf
So, here’s an example of a 2D reconstructed X-ray computed tomography slice of some calcite. Not my real data, but it can provide a quick example of what I’m trying to do.

From ImageJ Macro Recorder…

makeLine(4, 471, 878, 471);
run(“Plot Profile”);
selectWindow(“Calcite-1_OriginalImage_16bits-U_895x898x1516_45.26um0000.tif”);
run(“Remove Outliers…”, “radius=100 threshold=0 which=Dark”);

The line profile is there just to see the effect using the Live and Preview options.


^The original data


^The original data’s line profile


^After performing Remove Outlier


^The line profile, same location as the original, of the Remove Outlier image

This is not quite like my actual data, but is close enough. I have some deep valleys (negative Gaussian peaks) that are quite wide at the base (say 100-200 pixels wide). With close to 2000 slice images, I’m looking at over 2 hour processing time using the Remove Outlier function. Also, my “baseline” is not quite straight like this data, there is a long arc of negative curvature (corresponding to X-ray beam hardening cupping effect, which is what I’m actually trying to correct). So, once the negative peaks are removed, there is a very nice interpolation of the curved baseline.

Anyways, here is the example :smiley:

1 Like

Hey @Nik,

thanks for this nice challenge :slight_smile:

Alright, I read a bit about the algorithm. “Remove outliers” is a special median filter. If I understood it correctly, threshold 0 leads to this filter being a pure median filter. But I’m not sure if I interpret that part correctly. Your image was made with threshold 0 apparently and doesn’t look like a pure median filtered image.

Unfortunately, the median filter is not very fast on the GPU. The master student fixing this hasn’t been found yet. But why don’t we try it with a mean filter instead :wink:

close("*");

// Init GPU
run("CLIJ2 Macro Extensions", "cl_device=");
Ext.CLIJ2_clear();

radius = 30;
threshold = 2000;

// Load image from disc 
open("C:/Users/rober/Downloads/Calcite-1_OriginalImage_16bits-U_895x898x1516_45.26um0000.tif");
input = getTitle();
Ext.CLIJ2_push(input);

makeLine(4, 471, 878, 471);
run("Plot Profile");

// apply mean instead of median
Ext.CLIJ2_mean2DBox(input, filtered, radius, radius);

// determine abs difference between original and filtered image
Ext.CLIJ2_absoluteDifference(input, filtered, abs_difference);

// differentiate where the abs difference is above/below a given threshold
Ext.CLIJ2_greaterConstant(abs_difference, above, threshold);
Ext.CLIJ2_binaryNot(above, below);

// mask images above/below
Ext.CLIJ2_multiplyImages(input, below, input_below);
Ext.CLIJ2_multiplyImages(filtered, above, filtered_above);

// assemble final image
Ext.CLIJ2_addImages(input_below, filtered_above, result);

Ext.CLIJ2_pull(result);
makeLine(4, 471, 878, 471);
run("Plot Profile");

And it will look like this:

Would you mind testing it on your data and adapting the parameters accordingly? I’m not sure about the threshold…

If you want to try it with the median instead, I’d recommend downsampling the image, before computing the median filtered image and then upsampling the result again. This should make the median fast enough.

Let me know what you think!

Cheers,
Robert

1 Like

Hm, I finally understood the threshold thingy. The original code doesn’t use the absolute difference, but the difference instead.

So here is the updated code with threshold=0:

// To make this script run in Fiji, please activate 
// the clij and clij2 update sites in your Fiji 
// installation. Read more: https://clij.github.io
close("*");


// Init GPU
run("CLIJ2 Macro Extensions", "cl_device=");
Ext.CLIJ2_clear();

radius = 30;
threshold = 0;

// Load image from disc 
open("C:/Users/rober/Downloads/Calcite-1_OriginalImage_16bits-U_895x898x1516_45.26um0000.tif");
input = getTitle();
Ext.CLIJ2_push(input);

makeLine(4, 471, 878, 471);
run("Plot Profile");

// apply mean instead of median
Ext.CLIJ2_mean2DBox(input, filtered, radius, radius);

// determine difference between original and filtered image
Ext.CLIJ2_subtractImages(filtered, input, difference);

// differentiate where the difference is above/below a given threshold
Ext.CLIJ2_greaterConstant(difference, above, threshold);
Ext.CLIJ2_binaryNot(above, below);

// mask images above/below
Ext.CLIJ2_multiplyImages(input, below, input_below);
Ext.CLIJ2_multiplyImages(filtered, above, filtered_above);

// assemble final image
Ext.CLIJ2_addImages(input_below, filtered_above, result);

Ext.CLIJ2_pull(result);
makeLine(4, 471, 878, 471);
run("Plot Profile");

1 Like

Wow!

Yes, this works very, very well! I tried with my data (1000 images, 16 bit grayscale, 786 px x 788 px, 1.2 GB total). Needed to split the dataset into 2 1000 image stacks, as the almost 2.5 GB is larger than what my cards can handle (2 P5000 SLI-bridged). Tried three times, with results for 1000 images in 12.262 s, 9.407 s, and 9.469 s.

I will go through and see what the threshold does, I did vary the radius a bit just to see. Other than the smoothing of the image edges (which can be taken care of with a threshold/mask), this looks very much like the Remove Outlier results to me. I may try with a median, just to see if I can do it (will have to be slice-wise, I can get a median2DSphere, radius = 4, on the whole data set in about 38 s but I’m sure a kernel of 100 is too large)

Thank you for looking into this and coming up with a very nice solution!

1 Like

How long would plain ImageJ need for this? Asking for a friend :smirk:

1 Like

:grin:

I tried this morning on my data (so I lied, the image stack is 1987 images, not 2000, but that’s close enough).

It took 2 hours and 4 minutes for the whole stack. So, that’s about 62 minutes for 1000 images.

1 Like

Awesome! :star_struck:
After you made sure clij is doing the right thing, feel free to tweet about this incredible speedup factor :wink:

1 Like

So, I find this very interesting and spent the weekend trying to understand why a median filter takes so much longer. It’s rather obvious when one thinks about the math, though! But, what if an approximation of a median filter is good enough? This seems to be a common thing…

1 Like

I think it really depends on the use-case. In many cases a median filter with a huge radius can be replaced by some other filter combination. Furthermore, median filters are edge-preserving. If you don’t care about the edges, you definitely can use another approach. In your images, you remove these black holes. May I ask what they are? Air holes?

1 Like

Well, for these images in the example I don’t exactly know. I wish I could share the images I’m actually working on, but I can’t! But sure, let’s assume they are air holes that I want to get rid of with out doing a grayscale-based thresholding (because on my images, I can’t). Let’s also say that I want to preserve the edges as best as I can, and that the filter is resistant to extreme outliers (mean really isn’t). I could do a gaussian blur, but that doesn’t preserve the edges. I’ve been playing with the bilateral filter (very nice, btw!) and that seems to do rather well.

I use Avizo quite a bit for my X-ray tomography images, and there are some nice diffusion-based edge-preserving filters (did something like this ever make the wishlist?).

Do you have any other suggestions?

When I saw your images for the first time, I was thinking of a maximumImages operation of the original image and a mean-blurred version with a huge radius. That should remove dark areas as well while keeping the image in the remaining regions mostly untouched. I btw tried and it didn’t look so well :smiley: In these considerations, we should also take the subsequent analysis into account. Intensity and contrast measurements obviously suffer from filters being applied in advance…

1 Like