0
$\begingroup$

This question is in reference to this paper. In Fig. 1b, the joint histogram of data value and gradient magnitude is plotted. What I understand is that the values in the plot represent the frequencies of voxels having the corresponding data value and gradient. However, I am confused how to implement this in Mathematica. Below is what I have tried so far:

img = ImageRotate[Import["ChapelHillCThead.tif", "Image3D"], {Pi, {1, 0, 0}}];
data = ImageData[img];
imggrad = GradientFilter[img, 1];
grad = ImageData[imggrad];
dataf = Flatten@data;
gradf = Flatten@grad;
nd = BinCounts[dataf, {0, 1, 1/255}];
ng = BinCounts[gradf, {0, 1, 1/255}];
(* This is to calculate the joint histogram *)
jh = Table[Min[nd[[i]], ng[[j]]], {i, 1, 255}, {j, 1, 255}];

I know this cannot give me the desired result. However, I don't find any alternative way apart from checking the data value and gradient for each voxel and count them.

How can I get the desired result in a fast and efficient way?

[The Chapel Hill CT dataset can be found here]

$\endgroup$
7
  • $\begingroup$ I find .pvn files in the link you mention, how did you convert it to .tif? Can you upload a directly importable file? $\endgroup$
    – Ruud3.1415
    Nov 10, 2017 at 17:11
  • $\begingroup$ Also have a look at the Histogram function $\endgroup$
    – Ruud3.1415
    Nov 10, 2017 at 17:17
  • $\begingroup$ @Ruud3.1415 Please check this link[drive.google.com/file/d/1sUjrMR7V76fOK3uMv5_Cc3bkNCZPmEra/… for the file. $\endgroup$
    – user36426
    Nov 10, 2017 at 17:22
  • $\begingroup$ Yeah, thanks. How does ImageHistogram[Import["ChapelHillCThead.tif", "Image3D"]] work for you? $\endgroup$
    – Ruud3.1415
    Nov 10, 2017 at 17:26
  • $\begingroup$ @Ruud3.1415 It works fine. Do you want the screenshot? $\endgroup$
    – user36426
    Nov 10, 2017 at 17:29

1 Answer 1

2
$\begingroup$

This is more of an extended comment than an answer.

Because there are around 7.4 million observations, taking a large random sample can make things simpler.

(* Take a random sample of the 7.4 million observations *)
d = Transpose[{dataf, gradf}];
r = RandomSample[d, 100000];

(* Plot the values of the random sample *)
ListPlot[r]

Plot of random sample of data

To obtain a smooth version of a histogram you can use the SmoothKernelDistribution function:

skd = SmoothKernelDistribution[r];
ContourPlot[PDF[skd, {x, y}], {x, 0, 1}, {y, 0, 0.12},
 PlotRange -> All, Contours -> 20, PlotPoints -> 100]

SmoothKernelDistribution contour plot

There's also the Histogram3D function:

Histogram3D[r, Automatic, "PDF"]

3D histogram

But in all of these figures, I don't see any relationship for the two variables similar to what is displayed in the referenced article.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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