teem	/	nrrd	/	Visible Female	Inter-slice Brightness Problem

The real problem with the Visible Human Female color data lies in how the brightness changes suddenly from slice to slice. This problem was previously described in a paper by Máquez and Schmitt presented at the 1996 Visible Human Project Conference. That work was focused on the lungs of the Visible Human Male color data, but the problem is the same as described here. Those authors also described a method for correcting the problem. My pages here serve more to precisely visualize the problem as it occurs in the head of the Visible Human Female.

The easiest way to see the brightness problem is with simple summed projections, using the half-size head images created in an earlier step:

unu join -i *.headsm.ppm -a 3 \
  | unu project -a 2 -m sum | unu swap -a 1 2 -o xsum.nrrd
unu join -i *.headsm.ppm -a 3 \
  | unu project -a 1 -m sum | unu swap -a 1 2 -o ysum.nrrd
unu quantize -i xsum.nrrd -b 8 | unu save -f pnm | topng doc/ysum.png
unu quantize -i ysum.nrrd -b 8 | unu save -f pnm | topng doc/xsum.png

• unu project ... The axis ordering of the unu join output is (RGB,X,Y,Z). To create an image showing something about X versus Z, we sum along axis 2 (Y), to produce ordering (RGB,X,Z). To show Y versus Z, w sum along axis 1 (X), to produce ordering (RGB,Y,Z).

xsum.png: Summation along X: Y versus Z

Graph of R,G,B sums (see below)

The problem is the vertical bands in the RGB projection images. The fact that these bands appear in both images (regardless of projection direction), and that they appear uniformly in the images (regardless of the non-projected position) suggest that they are not a property of the tissue or the background, but of the image intensity itself. We can also graph the sum of component values over the whole slice. Nrrd doesn't (yet) have any facility for graphing, so we'll start with nrrd ...

unu project -i xsum.nrrd -a 2 -m sum -o colsum.txt

• unu project ... The axis ordering in xsum.nrrd is (RGB,Z,X), so summing along axis 2 gives us a 2-D nrrd with ordering (RGB,Z).
• ... -o colsum.txt: All the unu commands (except unu make and unu save ) guess output format based on filename extension. ".txt" signifies a plain headerless text file.

>> rgb = dlmread('colsum.txt');
>> d = 0:1:854;
>> plot(d,rgb(:,1),'r',  d,rgb(:,2),'g',  d,rgb(:,3),'b')
>> axis tight                                          

Because the variations seem to apply to all colors, we can get another kind of representation which starts with summing the three color components together, and then histogramming them on a per-slice basis:

unu join -i *.headsm.ppm -a 3 | unu project -a 0 -m sum \
  | unu reshape -s 92338 855 \
  | unu histax -a 0 -b 324 -t float -o histax.nrrd
cat > colormap.txt
  0.3 0.3 0
  0 0 1
  1 0 0
  1 1 0
  1 1 1
(Ctrl-D)
unu heq -i histax.nrrd -b 1024 -s 2 -a 0.7 \
  | unu swap -a 0 1 | unu flip -a 1 \
  | unu rmap -m colormap.txt | unu quantize -b 8 -o histax.ppm
convert histax.ppm doc/histax.png

• unu project ... This sums along the RGB axis, reducing the (RGB,X,Y,Z) volume to (X,Y,Z).
• unu reshape ... Just like Matlab's reshape: the data isn't changed, but the logical axis sizes are. This turns the (X , Y , Z) volume into a 2-dimensional (XY , Z) image, in which every (gray-scale) slice has been stretched into one long scanline.
• unu histax ... This replaces every scanline with a histogram of the scanline. The dimension of the nrrd stays the same, but the option for number of bins ("-b") determines the new size. Why 324 bins? A previous run of the command found the min and max values along the histogram axis to be 33 and 682; this number of bins produces the least banding by putting two values per bin. The output axis ordering is (H,Z).
• unu heq ... The goal is the put the axis histogram on screen in an intelligible way, and one way of improving the contrast of an image (so as to better utilitize the 8-bits of screen display) is to apply some amount of histogram equalization to the image. (Never mind the fact that the image in question is already a histogram!)
• unu swap ... Changes the axis ordering from (H,Z) to (Z,H).
• unu flip ... The axis ordering stays the same, but the samples along the histogram axis (the second axis, axis 1) are flipped end-to-end, so that the histogram has the same sense as the RGB graphs above.
• unu rmap ... This applies the colormap stored in colormap.txt to the image generated by the previous steps. The colormap makes it much easier to discern positional discontinuities in the histogram. The nrrd reader can read plain text files as 2-D arrays of floats.

ysum.png (repeated for comparison)

These histograms are a little hard to understand at first. As arranged here, every vertical scanling corresponds to one Z slice. The bottom of the histogram represents low values, the top represents high values. Bright regions indicate where (both in terms of slice and in terms grayscale value) there were the most pixels. When the slices suddenly get dark near the eyes, there is a jog downward in the upper band of the histogram. The isolated dark slices are visible as small spikes along the bottom edge of the histogram. The value in looking at this histogram, as opposed to simply looking at the R,G,B sum plots (above), is that it provides a sense of how the whole distribution of values is changing.

You can make a color version of this histogram too, by not initially summing the color components, so that the red, green, and blue channels of the axis histogram portray positional variations in the corresponding color component. The result may have more aesthetic value than informational value:

unu join -i *.headsm.ppm -a 3 | unu reshape -s 3 92338 855 \
  | unu histax -a 0 -b 256 -t float -o histaxcol.nrrd
unu heq -i histaxcol.nrrd -b 1024 -s 2 -a 0.7 \
  | unu swap -a 1 2 | unu flip -a 2 \
  | unu quantize -b 8 | topng doc/histaxcol.png

• unu join ... Creates a (RGB , X , Y , Z) volume
• unu reshape ... Converts it to an (RGB , XY , Z) "image"
• unu histax ... Generates an (RGB , H , Z) image.
• unu swap, flip ... Ends up an (RGB , Z , H) image.

A better way to get a component-specific visualization of how the brightness variations is to simply do axis histogramming for the red, green, and blue channels in isolation:

unu slice -i histaxcol.nrrd -a 0 -p 0 | unu heq -b 1024 -s 2 -a 0.7 \
  | unu swap -a 0 1 | unu flip -a 1 | unu rmap -m colormap.txt \
  | unu quantize -b 8 -o histaxR.ppm
unu slice -i histaxcol.nrrd -a 0 -p 1 | unu heq -b 1024 -s 2 -a 0.7 \
  | unu swap -a 0 1 | unu flip -a 1 | unu rmap -m colormap.txt \
  | unu quantize -b 8 -o histaxG.ppm
unu slice -i histaxcol.nrrd -a 0 -p 2 | unu heq -b 1024 -s 2 -a 0.7 \
  | unu swap -a 0 1 | unu flip -a 1 | unu rmap -m colormap.txt \
  | unu quantize -b 8 -o histaxB.ppm
convert histaxR.ppm doc/histaxR.png
convert histaxG.ppm doc/histaxG.png
convert histaxB.ppm doc/histaxB.png

• unu slice ... The axis histogram histaxcol.nrrd had axes (RGB , H , Z), and we want to see the (H,Z) image for red, green, and blue seperately, so we slice along the first axis (axis 0).

histaxG.png

histaxB.png

These component-specific axis histograms tell a great deal:

• The location (in Z) of the head is evident from from the location of the upper histogram bands in red and green. The blue histogram has only one major band.
• The largest step discontinuities appear in all components.
• There is more high-frequency variation in the green component than the other colors.
• On the right side of the histograms, single slices which are darker than their neighbors are more visible in the green and blue histograms than the red. These slices are indicated by spikes pointing downwards from the lower histogram band.
• Two features associated with blue voxels are visible in the blue histogram: the sinuses, and the missing chin.