CLIJ optimisation on GLCM

Hello everyone,

I am developping a plugin to perform 3D texture analysis using GLCM. The concept is to look at pixel values in the neighborhood of the central pixel and count the cooccurence of pairs of values. I did a first version using imglib2 in groovy, and it works pretty smoothly. But my users didn’t have the same performances and I am convinced that it can be even better using CLIJ for such low-level operations !
As discussed with Robert @haesleinhuepf, I investigated the two functions he kindly suggested me,
generateTouchMatrix and generateTouchCountMatrix as ways to perform the GLCM creation.
Unfortunately, the first one create a table with boolean answer (1 if the pair exist on the image, 0 if not) and the second looks promising, but works only on the diamond neigborhood, not a cubic one, and is not designed to deal with cases where the pair of pixels share the same pixel values…
So, here I am, asking for new ideas :slight_smile: !
My plan-B would be to look at OpenCL language and use custom code, but I dont know how hard it could be (if you have some good OpenCL Tutorial to share…).
My plan-C would be to wait for a CLIJ function closer to what I would like to do and advise my user to not perform image analysis on a toaster !
Thank you for your help

Julien

1 Like

Hey @jdumont,

great to hear from you. Generating two-dimensional histograms is tricky on GPUs because of the parallel nature of GPU based processing. That’s why such an algorithm I would preferrably implement on CPUs. For a quick test, I just copied the generateTouchCountMatrix command you mentioed (which is implemented on the CPU as well) and updated it for you so that it also counts neighboring pixels of the same intensity and in a half-box neighborhood. (for every pixel neighbors towards positive direction in each dimension. That’s 3 neighbors in 2D and 7 in 3D). The command comes via CLIJx to your Fiji if you update it. I gave it the short and simple name generateIntegerGreyValueCooccurrenceCountMatrixHalfBox :slight_smile:

If this command is reasonably fast even though it runs on the CPU, I would argue that your groovy imglib2 implementation might be the better way to go, because it doen’t push/pull pixels to the GPU in between It might just need some tweeking. Feel free to share your code and the imglib2 experts here in the forum may guide you.

Anyway, I’d love to hear if the method above does what you expect it to do.

Furthermore, the clij website has a couple of links to OpenCL tutorials for beginners and advanced developers, just scroll down to “Further reading”. If you find other great resources, I’m hapy to put them in this list.

Cheers,
Robert

2 Likes

I didn’t expect a function by itself so quickly ! Thank you very much @haesleinhuepf, I will look into that.

Here is my groovy implementation of GLCM :

import net.imagej.ops.image.cooccurrenceMatrix.MatrixOrientation3D
import net.imglib2.img.array.ArrayImgs
import net.imglib2.img.display.imagej.ImageJFunctions
import net.imglib2.type.numeric.real.DoubleType
import net.imglib2.Cursor
import ij.plugin.ImageCalculator
import ij.IJ
import loci.plugins.BF
import loci.formats.ImageReader
import loci.formats.MetadataTools
import loci.plugins.in.ImporterOptions
import net.imglib2.type.numeric.integer.LongType

def process(src, dst) {
   //Calculus of the co-occurence matrix in the 13 directions. Return stack of GLCM and number of grey-level
   img = imps[0]
   orientation = [MatrixOrientation3D.HORIZONTAL, MatrixOrientation3D.VERTICAL, MatrixOrientation3D.DIAGONAL,
   				MatrixOrientation3D.ANTIDIAGONAL, MatrixOrientation3D.HORIZONTAL_VERTICAL, MatrixOrientation3D.HORIZONTAL_DIAGONAL,
   				MatrixOrientation3D.VERTICAL_VERTICAL, MatrixOrientation3D.VERTICAL_DIAGONAL, MatrixOrientation3D.DIAGONAL_VERTICAL,
   				MatrixOrientation3D.DIAGONAL_DIAGONAL, MatrixOrientation3D.ANTIDIAGONAL_VERTICAL, MatrixOrientation3D.ANTIDIAGONAL_DIAGONAL,
   				MatrixOrientation3D.DEPTH]
   long[] dimensions = [nGL, nGL, orientation.size()]
   GLCM = ArrayImgs.doubles(dimensions)
   cursor = GLCM.randomAccess()

   for (z = 0; z < orientation.size(); z++){
   	result = ops.run("image.cooccurrenceMatrix", img, nGL, distance, orientation[z])
   	for(y = 0; y < nGL; y++){
   		for(x = 0; x < nGL; x++){
   			cursor.setPosition(x, 0)
   			cursor.setPosition(y, 1)
   			cursor.setPosition(z, 2)
   			value = result[x][y]
   			cursor.get().set(value)

   		}
   	}
   }
   img.close()
   return [GLCM, nGL]
}

I’m quite sure it can be improved, I learned groovy and Fiji plugin during the french lockdown. Any suggestions are welcome. Please note I pasted only the GLCM calculus part, if one need more, feel free to ask.

I will check the Further reading part, and look if i find anything better.

Cheers
Julien

1 Like

Hello again @haesleinhuepf ,

I tested your solution and it works, and quite fast ! In order to be consistent in the comparaison, I reduced my initial script posted above from 13 directions to 7, just like yours, and the start and stop time took into account the push. The CLIJ script (written in .ijm) outperformed the Groovy script by a factor of 7-fold !! 2.165s versus 14.950s exactly
Either CLIJ is outstandingly powerful, either my script is suboptimal, or the two options are both true.

Yet, I would have two remarks :

  • In order to make a GLCM, I would need to investigate also antidiagonal and top-dimension, for a total of 13 directions (a transposition afterwards allows to not have to look at the 26 neighborhood). Would you, at some point have the possibility to add them ?
  • I didn’t succeed in calling clij2.generateIntegerGreyValueCooccurrenceCountMatrixHalfBox(); in groovy, saying there is No signature of method. But i might have failed badly on this one, I didn’t spent too much time on it.

Anyway, thank you already for all your help ! I didn’t expected such support :smiley:

Have a good day

Julien

I’m glad to hear that. I assume that’s because the implementation in clij runs in parallel on the CPU. And the implementation is really not very sophisticated. :slight_smile:
I’m sure that parallelization is also possible with imglib2 in groovy. I’ll just tag @imagejan and @maarzt here, because they may know how a “parallel loop-builder” and/or image-ops call could do the grey-level coocurrence matrix operation you’re looking for. Let’s see what they say.

I’m happy to adapt the code. But I’m not sure if I understand. Which are these 13 neighbors? At the moment, all neighbor relationships in 2D and 3D are counted. Going from 7 to 13 in 3D would mean that some are counted twice. Assume a 1-D image with pixels A, B, C. I take the two pixels values of A and B and increase the value in the corresponding matrix entry. Doing this for B and A would just double the same value in our triangular matrix. So why do you want to go from 7 neighbors to 13 in 3D and not to 26? You need a full matrix, right? Are you interested in directional relationships? Curious!

Yes. That’s because this method is not available in CLIJ2. It lives in its experimental sibling CLIJx. It should show up in the auto-completion, e.g. if you enter “cooc”:

You can import CLIJx and use it analogously to CLIJ2, they are fully compatible:

import net.haesleinhuepf.clijx.CLIJx;

clijx = CLIJx.getInstance();

I collect new interesting methods (this is one for sure) in CLIJx, work a bit on them with collaborators and when we are confident that it does something reasonable, we can consider it for the next stable release.

You’re welcome. Thanks for the suggestion! I think this method can become quite useful.

Cheers,
Robert

2 Likes

Yes, I’m sure it is possible, but, at the time and at the level of my skills, it was not a priority. Again, this performance question came in my mind quite recently thanks to the feedback of the user of my plugin. I would have dig it much more if CLIJ wasn’t an option.

I already use an image-ops to perform the GLCM in my groovy code :

result = ops.run("image.cooccurrenceMatrix", img, nGL, distance, orientation[z])

The parallelization really does the trick.

Your questions are completely right, and i’m afraid to not be the best one to answer it, because I never thought about it ! Actually, the 13-directions came from what I read in multiple publications (1, 2, 3). As we want to investigate the possible directions from the central voxel, we have to look at in 13 ways. We usually dont do the 26 because, I think, you would lose some time computing A=>B and B=>A.The preferred strategy is to transpose the 13-directions matrix to obtain a 26-directions. I am personally not interested in directional relationship, as, at the end, I sum the GLCM of all directions before processing, to make a direction-invariant GLCM. I might have missed many concepts behind that, and I’m sure @dietzc would be more suited to answer, as he designed the corresponding imagej-ops .

It makes sense ! Ok, I didn’t import the right one so.

Cheers
Julien

1 Like

Again, I’m happy to adapt the code. However, I’m having issues understanding what directions should be implemented. Your code example above lists those:

MatrixOrientation3D.HORIZONTAL, MatrixOrientation3D.VERTICAL, MatrixOrientation3D.DIAGONAL,
   				MatrixOrientation3D.ANTIDIAGONAL, MatrixOrientation3D.HORIZONTAL_VERTICAL, MatrixOrientation3D.HORIZONTAL_DIAGONAL,
   				MatrixOrientation3D.VERTICAL_VERTICAL, MatrixOrientation3D.VERTICAL_DIAGONAL, MatrixOrientation3D.DIAGONAL_VERTICAL,
   				MatrixOrientation3D.DIAGONAL_DIAGONAL, MatrixOrientation3D.ANTIDIAGONAL_VERTICAL, MatrixOrientation3D.ANTIDIAGONAL_DIAGONAL,
   				MatrixOrientation3D.DEPTH

The documentation is also not very helpful unfortunately:
https://javadoc.scijava.org/ImageJ/net/imagej/ops/image/cooccurrenceMatrix/MatrixOrientation3D.html

I mean, also to simplify: what neighbors do these three stand for?

MatrixOrientation3D.HORIZONTAL, 
MatrixOrientation3D.VERTICAL, 
MatrixOrientation3D.HORIZONTAL_VERTICAL, 

Any hint is appreciated :slightly_smiling_face:

Cheers,
Robert

1 Like

It might be more obvious from the source code, where you have x, y, z coordinates:


In a 26-neighborhood in 3D, I’d also think you would need to look at 13 directions to count each pair exactly once (and in 2D you’d look at 4 directions for an 8-neighborhood, respectively).

Counting only 7 neighbors in 3D would be the equivalent to counting only 3 neighbors in 2D, i.e. you would miss the “backwards diagonal” directions in both ways.

image

3 Likes

Yes! You’re right! Also a very good visualisation. I’m still not fully understanding the terminology (horizontal_vertical corresponds to 1, 0, 1) but I think I know how to implement it. Thanks for your support @imagejan !

1 Like

Hey @imagejan and @jdumont,

just to make sure that I’m thinking right. Looking from the red pixel, we need to take the blue pixels in 3D into account, right? These are the 13…

Thanks!

Cheers,
Robert

1 Like

Hey @haesleinhuepf,
To me, it should more looks like that :

At least, this is what I understood from the source code @imagejan shared.

Cheers
Julien

2 Likes

Hehe, if you tilt that by 90 degrees along the Y-axis and the Z-axis, it’s the same as I posted… correct?

1 Like

You’re right ! But, according to what you posted, would’nt that mean that if the red pixel is on the last image of the stack it would be biased because the majority of measure would be done on non-existent pixel ? And if the red pixel was on the first slice of the stack would have stronger weight, no ?

Well, pixels outside the image are not taken into account. Thus, on the last slices, the red pixel only takes the neighbors in its 2D slice into account.

It’s really a bit tricky to get this right. Good to chat about it.

Btw. do you think a triangual matrix is right or the fully filled matrix?

That’s right, it would make like 4 neighbors. But for the first image, the 13 would be measured…Sounds odd to me, it is kind of unbalanced !

I would say you need the fully filled matrix : nothing prevent you to have a pair (15,14) somewhere and (14,15) somewhere else. How the triangular matrix would deal with that ?

Not really, because from the second plane no pixels are taken into account toward the first, while from the second last 9 pixels are taken into account towards the last. That’s 13 in both cases.

In a triangular matrix neighboring pixels with (14,15) and (15,14) would both be counted in the same matrix element: 14,15.

I’m now working on this example image which has three planes:

0 0 0 
0 1 0
0 0 0

2 2 2
2 2 2
2 2 2

0 0 0 
0 0 0 
0 0 0   

The corresponding cooccurrence matrix is

32 4 31
 4 0 9
40 0 20

To me it sounds correct. Would you mind cross-checking?

If you turn it by any direction, the “unbalancedness” translates to the right/left/top/bottom border of each plane, but essentially stays the same here. But we’re looking at the occurrence of pairs of pixels here, so for counting, all that matters is that you count each pair of neighbors, which will be true regardless of the direction of your cube.

That’s correct. The order of the pairs plays no role for counting the coocurrence. A full matrix would/should be diagonally symmetric.

That was exactly my thought! I would also expect it to be symetric.

Thus, also the code in ImageJ-Ops should be changed from

matrix[pixels[z][y][x]][pixels[sz][sy][sx]]++;

to

matrix[pixels[z][y][x]][pixels[sz][sy][sx]]++;
matrix[pixels[sz][sy][sx]][pixels[z][y][x]]++;

Correct?

Thanks guys! I’m enjoying this discussion.

I am afraid I can’t do it so easily, inside the ops there is a normalisation per the number of pairs in the image, and i can’t retrieve it.
Still, here is the normalised result, performed with the same code than posted above :

2.111   0.2778  4.0833
0.2778  0       0.4167
4.0278  0.4722  1.3333

Could you normalize yours the same way to compare both ?

And I learned from it as much as you like it !

1 Like

I would expect the normalised result has only numbers between 0 and 1 though…

1 Like