Teem: nrrd: Visible Female Color Data Processing: Fiducial Inspection


teem	/	nrrd	/	Visible Female	Fiducial Inspection

Having cropped out images of the fiducials, we need a way of differentiating the pixels inside and outside the fiducial markers, in order to analyze how they move from slice to slice. Since the markers are gray, the the background is blue, the saturation of the pixel color should be a good indicator. "unu" doesn't know anything about saturation or color, but it can compute the standard deviation of the three color components. More saturated colors will have higher standard deviation. Here are images of the standard deviations of the the circle and I-beam images shown above:

unu project -i avf1130a.raw.Z.circle.ppm -a 0 -m sd \
  | unu quantize -b 8 | topng doc/circle-sd.png
unu project -i avf1130a.raw.Z.ibeam.ppm -a 0 -m sd \
  | unu quantize -b 8 | topng doc/ibeam-sd.png

unu quantize ... Takes integral and floating point data down to a smaller number of bits. By not explicitly specifying a min and max, the min and max present in the input are mapped to 0 and 255, respectively.


`circle-sd.png`	`ibeam-sd.png`

This looks very promising, and it certainly beats trying to do some sort of segmentation! It also demonstrates a benefit of nrrd's philosophy that "everything is a scalar". That is, nrrd has no notion of vector-valued sample values. There is no such thing as a 3D array of colors, there is only a 4D array of color components. While it might seem strange to make a seperate axis for color information, since its not a spatial axis like X, Y, and Z, in this case it enabled us to use one of nrrd's many projection measures to collapse the RGB color into a useful scalar.

Now we can produce what are essentially graphs of the fiducials' X and Y positions, going through the 855 slices of the head:

unu join -i *circle.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 1 -m min \
  | unu quantize -b 8 | unu resample -k box -s x2 = \
  | unu swap -a 0 1 | topng doc/circle-X.png
unu join -i *ibeam.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 1 -m min \
  | unu quantize -b 8 | unu resample -k box -s x2 = \
  | unu swap -a 0 1 | topng doc/ibeam-X.png
unu join -i *circle.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 0 -m min \
  | unu quantize -b 8 | unu resample -k box -s x2 = \
  | unu swap -a 0 1 | topng doc/circle-Y.png
unu join -i *ibeam.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 0 -m min \
  | unu quantize -b 8 | unu resample -k box -s x2 = \
  | unu swap -a 0 1 | topng doc/ibeam-Y.png

There is some non-trivial unu cleverness going on here; each step of the first command (ending with circle-X.png) is:

unu join ... This collects all the PPM images into a volume, a 4-D array of color values. The "-a 3" flag says that the arrays are being joined along axis 3, which doesn't exist in the input arrays, so unu join knows that the output has one higher dimension. The axis ordering of the output of this step is (RGB,X,Y,Z).
1st unu project ... This projection along the color axis gives us an indicator of being inside or outside the fiducial. Projecting along an axis removes it from the axis ordering, going from (RGB,X,Y,Z) to (X,Y,Z).
2nd unu project ... Now we want to project along Y in order to get a picture of how the X position of the marker changes as a function of Z. Because Y is now the second axis, we project along axis 1. We use the "-m min" measure because the marker is darker than the background. Again, projection removes an axis from the ordering, so now the output ordering is (X,Z). To verify this, we can look at the labels of the output of this stage:
```
unu join -i *circle.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 1 -m min \
  | unu head - | grep labels

labels: "X" ""
```
The label for the second axis (Z) is empty because we've never specified one (and currently unu join doesn't have a way of doing that).
unu quantize ... The output of the first projection (standard dev of color) was a floating point volume, and all subsequent steps have also therefore been on floating point volumes. At some point, however, we need to get back to the world of fixed point.
unu resample ... This uses nearest-neighbor interpolation ("-k box") to double the number of samples ("x2") along the first axis (X), while leaving the second axis entirely untouched ("="). This will have the effect of exaggerating any displacements between slices.
unu swap ... The input axis ordering is still (X,Z), which if put on screen would have Z along the vertical, but we want to have these drawn on screen as a more of a horizontal graph, with Z along the horizontal. The easiest way to do this is to just switch the axis ordering to (Z,X). Keep in mind, however, that when the data is shown in usual raster order on a screen, the X values increase as you go downward, not upward. (If this was a problem, we could follow unu swap with "unu flip -a 1".)

circle-X.png: X versus Z

ibeam-X.png: X versus Z

circle-Y.png: Y versus Z

ibeam-Y.png: Y versus Z

What do these pictures tell us? The motivation in making these pictures was to see if there were sharp discontinuities. There are some, such as the shifts downward (toward higher X values) between Z locations 88 and 89, 405 and 406, and 712 and 713, which are visible to some extent in both the circle and I-beam X-versus-Z graphs. But these shifts are small, only about one pixel. For now, we won't try to correct for these; we'll make a volume dataset and then see if any artifacts we notice can be plausibly associated with the displacements measured by the graphs above.

The other thing to notice is that both markers moved in X and Y, and more in Y than in X. We can also see this in differences between images near the top and bottom. To improve image quality, we'll average three images at each location.

unu join -i avf100{7,8,9}a.raw.Z.circle.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 2 -m mean -o top.nrrd
unu join -i avf128{3,4,5}a.raw.Z.circle.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 2 -m mean -o bot.nrrd
unu 2op - top.nrrd bot.nrrd \
  | unu quantize -b 8 | topng doc/diff-circle.png
unu join -i avf100{7,8,9}a.raw.Z.ibeam.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 2 -m mean -o top.nrrd
unu join -i avf128{3,4,5}a.raw.Z.ibeam.ppm -a 3 \
  | unu project -a 0 -m sd | unu project -a 2 -m mean -o bot.nrrd
unu 2op - top.nrrd bot.nrrd \
  | unu quantize -b 8 | topng doc/diff-ibeam.png

2nd unu project ... The axis ordering in the output of the first projection is (X,Y,Z), so then projecting along axis 2 will eliminate Z, to give us (X,Y). We measure the average with "-m mean". The order of projections here matters: taking the standard deviation of the averages wouldn't be as useful. We save the output of the projection as a floating point nrrd.
unu 2op ... A number of unary, binary, and ternery arithmetic operations are available on nrrds; this just does a component-wise difference between the /tmp/top.nrrd and /tmp/bot.nrrd.


`diff-circle.png`	`diff-ibeam.png`

Looking carefully at the circle and I-beam Y-versus-Z graphs, it seems that I-beam moved downward more than the circle did, which implies there was some small amount of rotatation, in addition to the translation. This seems to be a very small effect, so we won't worry about it, but at least we had the means of detecting it.