The posts about colour deconvolution today have reminded me of a question I had some time ago…

Color deconvolution, as described by Ruifrok and Johnston, involves generating a 3x3 stain matrix using three stain vectors.

I understand that if two stains are available, then the remaining elements can be created by generating a third (pseudo)stain that is orthogonal to the first two.

As far as I can tell, this third stain is generated using the cross product in several places:

I wasn’t able to tell what approach CellProfiler takes (I figure this code is relevant, but I got lost ).

However, I understand that this is not exactly how it is implemented in the @gabriel’s ImageJ/Fiji plugin. From a quick look at the code, this may be because of negative values being clipped to 0:

In any case, the following Groovy script written for Fiji shows that the third stain vector is *not* orthogonal using the Fiji plugin, i.e. the dot product is not zero:

```
import sc.fiji.colourDeconvolution.*
def mat = new StainMatrix()
mat.init("H&E", 0.644211000,0.716556000,0.266844000,0.09278900,0.95411100,0.28311100,0.00000000,0.00000000,0.0000000)
mat.compute(true, false, new ij.ImagePlus("Anything", new ij.process.ColorProcessor(10, 10)))
double[] stain1 = [mat.cosx[0], mat.cosy[0], mat.cosz[0]]
double[] stain2 = [mat.cosx[1], mat.cosy[1], mat.cosz[1]]
double[] stain3 = [mat.cosx[2], mat.cosy[2], mat.cosz[2]]
println 'Stain 1: ' + stain1
println 'Stain 2: ' + stain2
println 'Stain 3: ' + stain3
println 'Dot product stain 1 x stain 2: ' + dot(stain1, stain2)
println 'Dot product stain 1 x stain 3: ' + dot(stain1, stain3)
println 'Dot product stain 2 x stain 3: ' + dot(stain2, stain3)
double dot(double[] v1, double[] v2) {
double s = 0
for (int i = 0; i < v1.length; i++)
s += v1[i] * v2[i]
return s
}
```

It’s not completely clear to me that negative values must be avoided in the stain matrix, and that this is more important than orthogonality.

It’s also not clear to me if/how much this matters.

I think probably all of us benefited from @gabriel’s implementation – I know I did, and I’ve seen it referred to a lot of time in other people’s code. But even though it seems to be pretty much the standard reference for many (in the absence of the original macro), I’m not sure it’s widely recognized that other implementations seem to have deviated a bit in this detail.

In any case, I’d be really interested to understand if there is a ‘right’ way to do it.

I’m also very interested in whether @phaub might have any more best practice suggestions from this