7
$\begingroup$

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.}]

CT Phantom Object

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]

Projection Image

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?

$\endgroup$
  • 1
    $\begingroup$ Erm, have you tried Image3DProjection[img, {Cos[theta1], Sin[theta1]*Cos[theta2], Sin[theta1]*Sin[theta2]}, "Mean"]? Just by using `"MeanIntensity", you just enforced the rescaling yourself... $\endgroup$ – Henrik Schumacher Mar 6 at 19:42
  • 2
    $\begingroup$ It's annoying that the official documentation doesn't say anything about this rescaling behavior and quite misleading that "Mean" mode, which is described as "mean projection for every color channel" is the appropriate option to choose when dealing with a monochrome image that only has one color channel instead of using "MeanIntensity" mode, which is described as "mean projection of the grayscale volume". Nonetheless, you are right @HenrikSchumacher. That is indeed the behavior of Image3DProjection. Your comment deserves to be made into the answer for this question. $\endgroup$ – Matt Mar 6 at 20:09
  • $\begingroup$ Indeed, I also find the documentation a bit vague. It is just a nuance (and I am not a native speaker), but I would have found "mean projection for each color channel" slightly better... I just tried "Mean" because the meaning of "mean projection of the grayscale volume" remains completely opaque to me. $\endgroup$ – Henrik Schumacher Mar 7 at 8:54
3
+50
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for?Browse other questions tagged or ask your own question.