Group Abstract Group Abstract

Message Boards Message Boards

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

Attachments:
POSTED BY: Henrik Schachner
19 Replies
Attachments:
POSTED BY: Henrik Schachner
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
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 9 years ago
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 9 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 9 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