Teem: nrrd: Visible Female Color Data Processing: Cropping


teem	/	nrrd	/	Visible Female	Cropping

Supposing that we'd like to crop the slices so to only contain the head, we need to know where to crop. We can project through a small subset of the slices in order to see where the "footprint" of the head is, within the frame of the color images. One kind of projection which should work is simply maximum intensity projection, since any tissue will be brighter in red and green than the blue surrounding material. The subset of slices we'll use is *{0,5}a.raw.Z (every 15th slice), but only up to avf1215a, since that is where the shoulders begin:

ls -1 data/avf11?{0,5}a.raw.Z data/avf12{0,1}{0,5}a.raw.Z | wc -l
            24

Knowing that the number of slices to project through is 24, we can now create the color volume and project it:

zcat data/avf11?{0,5}a.raw.Z data/avf12{0,1}{0,5}a.raw.Z \
  | unu make -i - -t uchar -s 2048 1216 3 24 -e raw \
  | unu project -a 3 -m max \
  | unu permute -p 2 0 1 -o headmax.ppm

unu make ... Running zcat on a sequence of files (each of which is an image for a different Z coordinate) will generate raw data in which the Z coordinate changes the slowest. Since the axis ordering for one image is (X,Y,RGB), the axis ordering for the zcat output should be (X,Y,RGB,Z), so the axis sizes given to "-s" are in that order.
unu project ... Projecting along the Z axis means that the Z axis is eliminated from the volume. The dimension is decreased by one, and the axis ordering of the output is (X,Y,RGB). Each scaline along Z is replaced by a single scalar which measures some aspect of the scanline. We're using maximum intensity projection: "-m max".
unu permute ... The output of projection has to be axis permuted just like an individual slice.

The resulting image (headmax.ppm) is as big as one of the slices; here's a downsampled version:

unu resample -i headmax.ppm -s = x0.3 x0.3 | topng doc/headmaxsm.png

By inspecting this image, a bound on the head was found with the help of xv (middle button gives pixel coordinates):

unu crop -i headmax.ppm -min 0 771 271 -max M 1319 945 \
  | unu resample -k cubic:0,0.5 -s = x0.5 x0.5 \
  | topng doc/headmaxcrop.png

unu crop ... This selects a sub-array from the input, bounded by the coordinates given to "-min" and "-max". "M" means "highest coordinate along this axis".

Besides cropping down to the head, we would also like to see if there are any problems with spatial displacements between the images, which would create sharp edges in the data where none existed in the tissue. This is facilitated with the fiducial markers which appear in every slice. Inspection of individual slices (with the commands shown earlier) determined 80x80 pixel windows which bound the circle marker, and the I-beam marker on the opposite side of the slice. Here they are for slice avf1130a:

zcat data/avf1130a.raw.Z \
  | unu make -i - -t uchar -s 2048 1216 3 -e raw \
  | unu permute -p 2 0 1 -o tmp.ppm
unu crop -i tmp.ppm -min 0 365 945 -max M m+79 m+79 \
  | topng doc/circle.png
unu crop -i tmp.ppm -min 0 1660 950 -max M m+79 m+79 \
  | topng doc/ibeam.png


`circle.png`	`ibeam.png`

Having identified the region for the head and markers in each image, we go through all 855 slices and save out these regions as PPMs. We also save a half-size version of the head images, as these will be easier to handle in a later processing step.

foreach SLC ( `(cd data; echo avf*.raw.Z)` )
  echo $SLC
  zcat data/$SLC \
     | unu make -i - -t uchar -s 2048 1216 3 -l X Y RGB -e raw \
     | unu permute -p 2 0 1 -o tmp.ppm
  unu crop -i tmp.ppm -min 0 365 945 -max M m+79 m+79 -o $SLC.circle.ppm
  unu crop -i tmp.ppm -min 0 1660 950 -max M m+79 m+79 -o $SLC.ibeam.ppm
  unu crop -i tmp.ppm -min 0 771 271 -max M 1319 945 -o $SLC.head.ppm
  unu resample -i $SLC.head.ppm -k tent -s = x0.5 x0.5 -o $SLC.headsm.ppm
end

unu make ... The only new thing here is the use of the "-l" option to give a "label" for each axis. Axis labels are short descriptive strings which follow the axis through various nrrd transformations. These will be used below to keep track of which axes are which.

Go grab lunch will this runs, it can take over an hour, a significant part of which is simply running zcat over and over.

One useful aspect of the NRRD writer is that when it writes PPM images, it stores information for the NRRD headers in PPM comments. Since the NRRD file format can contain many different fields, its handy to have them stored in non-NRRD files. It is easy to read these headers, too. Instead of using head to look at the PPM header (and risk clobbering your terminal settings with bizarro binary data), use unu head:

unu save -i avf1130a.raw.Z.head.ppm -f nrrd | unu head -

NRRD0001
content: crop(permute(???,2,0,1),[0,2]x[771,1319]x[271,945])
type: unsigned char
dimension: 3
sizes: 3 549 675
labels: "RGB" "X" "Y"
encoding: raw

The content field shows that this image is the result of permuting and then cropping some original array, and there is enough information here to recreate the result. The labels field shows what became of the information originally given to unu make's "-l" option.