Convert 2D into a 3D object: radiotherapy treatment planning system

henrikschachner · September 11, 2016, 8:40pm

Dear all,

my data consist of a list of lists of points in 3D. After running the code

ClearAll["Global`*"]
SetDirectory[NotebookDirectory[]];
sliceData = << "SliceData.txt";
Graphics3D[Line /@ sliceData, Boxed -> True, Axes -> True]

one gets:

Clearly the data describe a 3D object (as “wire frame”).

QUESTION: How can those data be converted into a single 3D Mathematica object (graphics, mesh, …)? Is there an already implemented way (a routine) I am missing?

NOTE on DATA:

The data stem from our radiotherapy treatment planning system. There we are working with contours which are drawn on any single CT slice. Organs at risk and target volumes are defined by those contours. One special contour is the “BODY contour”; this is what is shown in my example data. Dose will be calculated only inside this BODY contour, therefore e.g. the volume around the ears and nose appears to be exaggerated (for being on the safe side). Those treatment plans can be exported as DICOM files and nicely imported in Mathematica.

When high energy radiation is applied to a body it turns out that the dose next to the skin is highly diminished; this is due to the buildup effect of the dose. When full dose at the skin is wanted, one has to anticipate this buildup effect, and this can be realized by putting a “flab” onto the skin: A flab is a layer made of some tissue equivalent material.

Optimal flabs can have quite irregular shapes. I recently learned at a conference that there is the possibility to 3D print those individual flabs - if one can provide the data … So - probably to everybody’s disappointment - I do not want to print any BODY contour, but I imagined a BODY contour might serve here as kind of a “honey pot”.

Best regards and many thanks! – Henrik

mthiel · September 11, 2016, 9:33pm

Hi Henrik,

there is this:

ListSurfacePlot3D[Flatten[sliceData, 1], MaxPlotPoints -> 75]

but there are obvious artefacts. I have another approach which I will implement next.

Cheers,

Marco

mthiel · September 11, 2016, 11:00pm

The second idea was to transform everything into slices and then use Image3D.

sliceData = ToExpression[Import["/Users/thiel/Desktop/SliceData.txt"]];
range = {MinMax[#[[All, 1]]], MinMax[#[[All, 2]]]} &@Flatten[sliceData[[All, All, {1, 2}]], 1];
imgs = Table[Image@Graphics[Polygon[sliceData[[k, All, {1, 2}]]], PlotRange -> range], {k, 1, Length[sliceData]}];

This gives slices like these:

Partition[imgs[[1 ;; ;; 10]], 3] // TableForm

Then something like this should work:

Image3D[imgs, "Bit", BoxRatios -> {1, 1, 3}, Background -> White, 
 ColorFunction -> "Grayscale"]

But it doesn’t. So I tried something more like this:

graphics = Table[Graphics3D[Polygon[sliceData[[k, All, {1, 2, 3}]]]], {k, 1, Length[sliceData]}];
Show[graphics]

The thing is that this is not a “real” 3D object but just a bunch of 2D planes in 3D space. This here would “fill the gaps”.

Show[Table[
  Graphics3D[
   Polygon[Riffle[
     Append[#, sliceData[[k, 1, -1]] + 1.5] & /@ 
      sliceData[[k, All, {1, 2}]], 
     Append[#, sliceData[[k, 1, -1]] - 1.5] & /@ 
      sliceData[[k, All, {1, 2}]]]]], {k, 1, Length[sliceData]}]]

Cheers,

Marco

henrikschachner · September 12, 2016, 10:43am

Dear Marco,

thank you very much for taking the time looking into this problem!

Actually ListSurfacePlot3D in this case would be the perfect function - which I missed completely (shame on me!). It is a pity that it shows those artifacts. I imagine that someone from WR knows some undocumented Method options …

Otherwise your second solution seems to show a way to go. Probably I have to modify it a bit, because it can happen that there are several disjunct areas in one slice.
My final goal here is to end up with something which can be exported for 3D printing.

Again many thanks and best regards from Germany! – Henrik

Vitaliyk · September 12, 2016, 11:08am

In addition to @Marco Thiel nice suggestions, I’d like to mention a possibility of resampling and building a spline surface. You have nice dense data. But number of points in each loop is different. A high resolution resampled 3D array of 3D points where number of points in each loop is the same could be your final output object. This could help 3D printing too. We start from creating a spline curve in 3D for every loop:

loopsBSF=BSplineFunction/@sliceData;
ParametricPlot3D[#[t]&/@loopsBSF,{t,0,1},
	PlotTheme->"Detailed",PlotLegends->None,PlotStyle->Thickness[0]]

where the loops are already parametric spline functions and not data points. Now we build a higher resolution array with 1000 points per every loop and a spline surface:

higHreSdatA=Table[loopsBSF[[k]][t],{k,Length[sliceData]},{t,0,1,.001}];
surfBSF=BSplineFunction[higHreSdatA,SplineClosed->{False,True}]

ParametricPlot3D[surfBSF[u,v],{u,0,1},{v,0,1},
	PlotRange->All,Mesh->All,PlotStyle->{Orange,Specularity[White,10]},
	SphericalRegion->True,PlotTheme->"Detailed",PlotLegends->None]

You can increase resampling generally for a higher resolution. You can also ask for surface thickness and make .STL file for 3D priniting:

ParametricPlot3D[surfBSF[u,v],{u,0,1},{v,0,1},
	PlotRange->All,Mesh->All,PlotStyle->{Orange,Specularity[White,10]},
	SphericalRegion->True,PlotTheme->"ThickSurface",PlotLegends->None]
Import[Export["test.stl",%]]

And this is object build now from only vertical lines - the new information derived from your data:

Graphics3D[{Opacity[.2], Line[Transpose[higHreSdatA]]}, Boxed -> False, SphericalRegion -> True]

Vitaliyk · September 12, 2016, 11:47am

BTW @Henrik Schachner, where are this data from? What project are you working on and what is your final goal? It is always interesting to know the story behind the data.

ghurst · September 12, 2016, 3:51pm

Nice idea to use Image3D. As of version 11, one can smooth things out a bit with ImageMesh:

(* imgs is defined from the above post *)
im = ImageRotate[ColorNegate[Image3D[imgs]], {?, {0, 1, 0}}];

(* reconstruct data with marching cubes *)
mesh = ImageMesh[im, Method -> "DualMarchingCubes"];

(* visualize with nice box ratios *)
BoundaryMeshRegion[mesh, BoxRatios -> {9, 6, 7}]

For visualization purposes only, we can smooth things out a bit more with PlotTheme -> "SmoothShading":

BoundaryMeshRegion[mesh, BoxRatios -> {9, 6, 7}, PlotTheme -> "SmoothShading"]

ghurst · September 12, 2016, 6:31pm

Very neat idea!

If OP wants to print the solid head, he can fill the holes using RepairMesh.

plot = ParametricPlot3D[surfBSF[u, v], {u, 0, 1}, {v, 0, 1}, PlotRange -> All, PlotPoints -> 100];
RepairMesh[DiscretizeGraphics[plot], PlotTheme -> "SmoothShading"]

Vitaliyk · September 12, 2016, 6:35pm

Awesome tip, @Chip Hurst! RepairMesh is indeed a nice function to take advantage of. Repairing defects in the mesh regions is a crucial step for 3D printing.

henrikschachner · September 12, 2016, 9:10pm

Dear @Vitaliy Kaurov, dear @Chip Hurst,

thank you very much for your valuable contributions - I learned a lot! And I think BSplineFunction does the trick. (Image3D is less suitable because of distortions.)

The data stem from our radiotherapy treatment planning system. There we are working with contours which are drawn on any single CT slice. Organs at risk and target volumes are defined by those contours. One special contour is the “BODY contour”; this is what is shown in my example data. Dose will be calculated only inside this BODY contour, therefore e.g. the volume around the ears and nose appears to be exaggerated (for being on the safe side). Those treatment plans can be exported as DICOM files and nicely imported in Mathematica.

When high energy radiation is applied to a body it turns out that the dose next to the skin is highly diminished; this is due to the buildup effect of the dose. When full dose at the skin is wanted, one has to anticipate this buildup effect, and this can be realized by putting a “flab” onto the skin: A flab is a layer made of some tissue equivalent material.

Sorry for this lengthy explanation! My point is: Optimal flabs can have quite irregular shapes. I recently learned at a conference that there is the possibility to 3D print those individual flabs - if one can provide the data …
So - probably to everybody’s disappointment - I do not want to print any BODY contour, but I imagined a BODY contour might serve here as kind of a “honey pot”.

Best regards and many thanks! – Henrik

Vitaliyk · September 12, 2016, 9:33pm

Happy to help, @Henrik Schachner ! And there is no “lengthy explanations!”, just the detailed ones - and that’s a good thing. We added the data description in the top post, - this is what makes this case interesting and potentially useful to other professionals. Thank you for sharing!

henrikschachner · September 13, 2016, 5:06pm

RepairMesh[DiscretizeGraphics[plot], PlotTheme → “SmoothShading”]

@Chip Hurst, thank you again for this nice hint! I think the option PlotTheme -> "SmoothShading" is important, but I cannot find it anywhere in the documentation. It would be a pity if it were just undocumented - or am I again missing something?

ghurst · September 13, 2016, 5:22pm

It’s towards the bottom of the Details and Options section in the documentation for MeshRegion and BoundaryMeshRegion.

henrikschachner · September 13, 2016, 5:30pm

Oooops … thanks again! Henrik

vaungureanu · September 14, 2016, 11:53am

Congratulations @Henrik Schachner! Congratulations all the contributors! Impressive common result!

henrikschachner · September 14, 2016, 8:36pm

Thank you for the congratulations @Valeriu Ungureanu! I did not expect my little question would get such a resonance. In fact radiotherapy is full of very nice applications for Mathematica - but most of them are much too special to be presented in this community, unfortunately. E.g. our linear accelerator produces very detailed logfiles. These can be read in into Mathematica (excellent for performing checks!) and furthermore can be visualized with Mathematica. If you want get an impression, see my Youtube channel (explanations still only in German, sorry). But nearly nobody else has those logfiles, so it makes no sense showing here any code …

Regards – Henrik

vaungureanu · September 15, 2016, 4:39am

There are impressive both the result and the Community contribution. So, I think your new posts may have similar resonance.

henrikschachner · September 18, 2016, 3:36pm

I would like to make a follow-up to this discussion:

The data I supplied above have the problem that the “wires” in each slice are shifted relative to each other. The following image should make clear what I am talking about:

This of course is a problem when building a regular mesh. I am trying a correction by shifting the wires in a way that point distances within neighboring slices become minimal like so:

Our complete code now reads:

ClearAll["Global`*"]
SetDirectory[NotebookDirectory[]];
sliceData = << "SliceData.txt";

bsp = BSplineFunction[#, SplineClosed -> False] & /@ sliceData;

highResData = Table[bsp[[h]][t], {h, 1, Length[bsp]}, {t, 0, 1, .003}];

(* definition of difference between to sets of "wire data" - to be minimized: *)    
totalPointDist[data1_List, data2_List, n_?NumericQ] := 
 Total[EuclideanDistance[#1, #2] & @@@ Transpose[{data1, RotateLeft[data2, n]}]]

(* minimization and shift: *)
Monitor[
 Do[highResData[[n]] = 
   RotateLeft[highResData[[n]], \[FormalN] /. Last@Minimize[totalPointDist[highResData[[n - 1]], 
        highResData[[n]], \[FormalN]], \[FormalN], Integers]], {n, 2, Length[highResData]}], n]

So with the friendly help of @Marco Thiel, @Vitaliy Kaurov, @Chip Hurst one ends up with:

I.e.now there appears to be no more mesh defects!
The effect can also be seen by showing the vertical lines only (according to Vitaliy):

Vitaliyk · September 18, 2016, 4:08pm

Great idea, @Henrik Schachner and a good choice for a function to minimize.

henrikschachner · November 13, 2016, 3:10pm

I still have to make an addition:

When data of more realistic structures are used, e.g.:

the method above simply does not work! One reason is the possible big difference in the length of neighboring “wires”. This graphic shows as an example an imprint (i.e. negative form) of an human ear - with some typical “complications” added like show by these slices:

The problem here is that those slices contain more than just a single line/wire - and they have different meaning. What I want to end up with is of course:

Hence a completely different approach was needed:

Basically I create from a stack of those images a 3D array;
from this an interpolated function is defined;
using RegionPlot3D a graphic is generated;
finally this graphic is converted into a mesh.

The function mkImage does the relevant job of creating the correct images. Here I had to use high level functions like FillingTransform or SelectComponents. My full (but short) code now reads:

ClearAll["Global`*"]

mkImage[wire2D_, xyRange_] := Module[{imgGraphics, img0, img1},
  imgGraphics = 
   Binarize[
    ColorNegate[(ColorConvert[#, "Grayscale"] &) /@ 
      Image[Graphics[{White, EdgeForm[Black], Polygon[wire2D]}, 
        PlotRange -> xyRange]]]];
  img0 = FillingTransform[imgGraphics];
  img1 = FillingTransform@
    SelectComponents[imgGraphics, #EnclosingComponentCount == 1 &];
  ImageSubtract[img0, img1]
  ]

SetDirectory[NotebookDirectory[]]

sliceData = << "Bolus3D_complex.txt";

planarAssoc = 
  GroupBy[sliceData, #[[1, 3]] &] /. {x_, y_, z_} :> {x, y};

border = 2;
range3D = {{-border, border}, {-border, border}, {0, 0}} + 
  CoordinateBounds@sliceData

imgs = mkImage[#, Most[range3D]] & /@ Values[planarAssoc];

data3D = Transpose[ImageData /@ imgs, {3, 2, 1}];

func3D = ListInterpolation[data3D, {#1, Reverse[#2], #3} & @@ range3D,
   InterpolationOrder -> 3]

gr = RegionPlot3D[func3D[\[FormalX], \[FormalY], \[FormalZ]] > .5, 
  Evaluate[
   Sequence @@ 
    MapThread[
     Prepend, {range3D, {\[FormalX], \[FormalY], \[FormalZ]}}]],
  BoxRatios -> Automatic, PlotPoints -> {100, 100, 100}, 
  ImageSize -> 700]

mesh = DiscretizeGraphics[gr, ImageSize -> 700]

Export["bolus.stl", mesh, "STL"]

And this now seems to work! I end up with:

Many thanks again to all who have helped!

Regards – Henrik

Topic		Replies	Views
Converting data/graph into a 3D printable file (.stl or similar) Questions engineering , import-and-export , modeling	19	37866	January 18, 2017
Making images into 3D cylindrical shape Staff Picks image-processing , geometry , graphics-and-visualization , 3d-printing , staff-pick , wolfram-repositories	1	4506	October 8, 2024
Can I Plot Multiple 2D Graph Slices as a Single Integrated 3D Graph? Questions mathematica , graphics-and-visualization	11	19589	July 6, 2014
3D Thinning algorithm in Wolfram Language Staff Picks geometry , graphics-and-visualization , graphs-and-networks , optimization , 3d-printing , staff-pick , wolfram-repositories	1	16140	November 30, 2018
Plot several 2D plots on one 3D graph? Questions graphics-and-visualization	13	8408	March 8, 2019

Convert 2D into a 3D object: radiotherapy treatment planning system

NOTE on DATA:

Related topics