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

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

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]

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.

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

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.

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!

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

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

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:

Bolus3D_complex.txt

in.txt

out.txt