Message Boards Message Boards

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

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:

enter image description here

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:
POSTED BY: Henrik Schachner
19 Replies

I still have to make an addition:

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

enter image description here

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:

enter image description here

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:

enter image description here

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:

enter image description here

Many thanks again to all who have helped!

Regards -- Henrik

Attachments:
POSTED BY: Henrik Schachner

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:

enter image description here

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:

enter image description here

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:

enter image description here

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):

enter image description here

POSTED BY: Henrik Schachner

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

POSTED BY: Vitaliy Kaurov

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

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

POSTED BY: Henrik Schachner

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

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.

POSTED BY: Vitaliy Kaurov

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

POSTED BY: Henrik Schachner

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!

POSTED BY: Vitaliy Kaurov

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

enter image description here

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]

enter image description here

enter image description here

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",%]]

enter image description here

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]

enter image description here

POSTED BY: Vitaliy Kaurov
Posted 8 years ago

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"]

enter image description here

POSTED BY: Greg Hurst

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.

POSTED BY: Vitaliy Kaurov
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?

POSTED BY: Henrik Schachner
Posted 8 years ago

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

enter image description here

POSTED BY: Greg Hurst

Oooops .... thanks again! Henrik

POSTED BY: Henrik Schachner

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

POSTED BY: Henrik Schachner

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

enter image description here

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]

enter image description here

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

enter image description here

Cheers,

Marco

POSTED BY: Marco Thiel
Posted 8 years ago

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

enter image description here

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

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

enter image description here

POSTED BY: Greg Hurst

Hi Henrik,

there is this:

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

enter image description here

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

Cheers,

Marco

POSTED BY: Marco Thiel
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract