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

GROUPS:

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

Attachments:
1 year ago
19 Replies
 Marco Thiel 8 Votes 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
1 year ago
 Marco Thiel 10 Votes 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
1 year ago
 Chip Hurst 4 Votes 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"] 
1 year ago
 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
1 year ago
 Vitaliy Kaurov 9 Votes 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] 
1 year ago
 Chip Hurst 4 Votes 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"] 
1 year ago
 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.
1 year ago
 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?
1 year ago
 Chip Hurst 1 Vote It's towards the bottom of the Details and Options section in the documentation for MeshRegion and BoundaryMeshRegion.
1 year ago
 Oooops .... thanks again! Henrik
1 year ago
 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.
1 year ago
 Henrik Schachner 1 Vote 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
1 year ago
 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!
1 year ago
 Valeriu Ungureanu 2 Votes Congratulations @Henrik Schachner! Congratulations all the contributors! Impressive common result!
1 year ago
 Henrik Schachner 1 Vote 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
1 year ago
 There are impressive both the result and the Community contribution. So, I think your new posts may have similar resonance.
 Henrik Schachner 2 Votes 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):
 Henrik Schachner 1 Vote 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 Attachments: