Is there an option to make Image3DProjection respect the underlying data type of the image array?

Question

Even though Image3D supports a variety of data types (Bit, Byte, Bit16, Real32, Real64), the Image3DProjection function discards the quantitative information and rescales everything to the range from 0 to 1. This is problematic when dealing with calibrated image data where the numerical values are meaningful. X-ray computed tomography data is an example of an image type where the signal intensity values are calibrated and meaningful. The Hounsfield unit scale is defined with ‑1000 HU being air and 0 being water. Various biological tissues have known ranges in the Hounsfield unit scale that are invariant between scans, patients, different vendors' equipment, etc.

For the sake of demonstration, let's define a simulated low resolution CT image of a phantom object:

airbackground = Table[-1000., {10}, {10}, {10}];
softtissue = RandomReal[NormalDistribution[125, 100], {7, 5, 5}];
hardobject = RandomReal[NormalDistribution[2000, 500], {3, 3, 3}];

data = airbackground;
data[[2 ;; 8, 4 ;; 8, 3 ;; 7]] = softtissue;
data[[4 ;; 6, 7 ;; 9, 1 ;; 3]] = hardobject;

img = Image3D[data, "Real64", ColorSpace -> "Grayscale", 
    ColorFunction -> "XRay", Background -> Black];

and display it with a window setting from ‑1000 HU (air) to 3000 HU (cortical bone)

ImageAdjust[img, 0.,{-1000., 3000.}]

We can check and see that Mathematica has retained the signal intensity information for this 3D image in the original scale that it was specified in (Hounsfield units, in this case).

stillrealdata = ImageData[img];
stillrealdata[[3 ;; 7, 3 ;; 7, 5]] // MatrixForm

The output is something like this (actual numbers will vary since we defined the phantom with RandomReal values:

\begin{array}{ccccc} -1000. & -122.221 & 94.8954 & 199.379 & 202.886 \\ -1000. & 362.572 & 7.79696 & 156.627 & 251.862 \\ -1000. & 59.939 & 8.30324 & 263.878 & 247.806 \\ -1000. & 285.41 & 132.313 & 183.714 & 6.53021 \\ -1000. & 261.366 & 57.2342 & 102.978 & 254.443 \\ \end{array}

Now, let's compute a mean intensity projection of the 3D image along a specified direction:

theta1 = 60 Degree;
theta2 = 40 Degree;

projimg = Image3DProjection[img,
    {
        Cos[theta1],
        Sin[theta1]*Cos[theta2],
        Sin[theta1]*Sin[theta2]
    }, "MeanIntensity"];

and display it.

ImageResize[projimg, 300]

Unfortunately, Image3DProjection has rescaled everything to the range 0 to 1. Thus, the resulting projection image data, which would be meaningful if it was in Hounsfield units, is not.

ImageData[projimg] // MatrixForm

\begin{array}{ccccccccccccccc} 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0.0555556 & 0.0555556 & 0.0555556 & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0.111111 & 0.333333 & 0.333333 & 0.333333 & 0.166667 & 0.111111 & 0.0555556 & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0.222222 & 0.444444 & 0.444444 & 0.388889 & 0.277778 & 0.166667 & 0.0555556 & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0.166667 & 0.388889 & 0.444444 & 0.277778 & 0.388889 & 0.166667 & 0.111111 & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0.111111 & 0.222222 & 0.277778 & 0.444444 & 0.277778 & 0.222222 & 0.0555556 & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0.111111 & 0.222222 & 0.333333 & 0.444444 & 0.333333 & 0.222222 & 0.111111 & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0.111111 & 0.277778 & 0.333333 & 0.277778 & 0.333333 & 0.111111 & 0.111111 & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0.111111 & 0.166667 & 0.222222 & 0.222222 & 0.166667 & 0.111111 & 0.0555556 & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0.0555556 & 0.111111 & 0.0555556 & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. & 0. \\ \end{array}

---

Is there a way to turn off this counterproductive automatic rescaling behavior?

Answer 1

I think what you are looking for is to use "Mean" instead of "MeanIntensity":

img1 = Image3DProjection[
  img, 
  {Cos[theta1], Sin[theta1]*Cos[theta2], Sin[theta1]*Sin[theta2]}, 
  "Mean"
  ];

Now one has

ImageData[img1] // MinMax

{-555.556, 372.955}

which seems to preserve the units.