# How to draw cylinders inside a cube

I want to create the following graphic (ignore the unit vectors).

What I do is the following (not the most elegent and smart way I guess:-)!):

1) Use the code below to generate randomly distributed but not intersecting circles (I adopt the code from this forum but I don't remember the thread. Actually I learnt from my recent posts other ways to do this.)

distinct[n_, r_] :=
Module[{d, f, p}, d = {Disk[RandomReal[{-1, 1}, 2], r]};
Do[f = RegionDistance[RegionUnion @@ d];
While[p = RandomReal[{-1, 1}, 2]; f[p] < r];
d = Append[d, Disk[p, r]], {n - 1}];
d]


Generate the circles

circles = distinct[75, 0.1];


Generate the respective cylinders which have this circles as cross sections.

cylinders =
Graphics3D[{Cyan, EdgeForm[Thick],
Cylinder[{{#[[1]], #[[2]], -3}, {#[[1]], #[[2]], 3}}, 0.1] & /@
Map[First, fibers]}, PlotRange -> {{-1, 1}, {-1, 1}, Automatic}, Lighting -> "Neutral"];


Generate the parallelepiped.

par= Graphics3D[{Blue, EdgeForm[Thick], Opacity[0.6],
Cuboid[{-1, -1, -3}, {1, 1, 3}]}, Lighting -> "Neutral"];


And adding them altogether

Show[{cylinders, par}, PlotRange -> {{-1, 1}, {-1, 1}, Automatic},
Boxed -> False]


But the result is quite unsatisfactory.

I will appreciate any help. Thanks in advance!

• A bit related: How to plot contours in the faces of a cube Commented Nov 5, 2015 at 20:00
• You can find almost everything here:-)! Commented Nov 5, 2015 at 20:04

Firstly, let us generate some set of random circles with findPoints from this answer

findPoints =
Compile[{{n, _Integer}, {low, _Real}, {high, _Real}, {minD, _Real}},
Block[{data = RandomReal[{low, high}, {1, 2}], k = 1, rv, temp},
While[k < n, rv = RandomReal[{low, high}, 2];
temp = Transpose[Transpose[data] - rv];
If[Min[Sqrt[(#.#)] & /@ temp] > minD, data = Join[data, {rv}];
k++;];];
data]];

npts = 150;
r = 0.03;
minD = 2.2 r;
low = 0;
high = 1;

pts = findPoints[npts, low, high, minD];
g2d = Graphics[{FaceForm@Lighter[Blue, 0.8],
EdgeForm@Directive[Thickness[0.004], Black], Disk[#, r] & /@ pts},
PlotRange -> {{low, high}, {low, high}}, Background -> Lighter@Blue]


## Method 1: Texture

We can simply use this graphics as a texture of the cube

pad = 0.1;
coords = Tuples[{0, 1}, 3];
cube = Polygon[{{1, 3, 7, 5}, {1, 5, 6, 2}, {5, 7, 8, 6}, {7, 3, 4,
8}, {3, 1, 2, 4}, {6, 8, 4, 2}}];
vtc = pad + (1 - 2 pad) coords[[;; , {1, 3}]];

Graphics3D[{Texture[g2d],
GraphicsComplex[coords, cube, VertexTextureCoordinates -> vtc]},
Lighting -> "Neutral", Boxed -> False, ImageSize -> 500]


## Method 2: MeshRegion

I'm appreciate many upvotes so I want to expand my answer and add a more general approach. Mathematica has very powerful (and still very limited) region functions.

Let's try to use some interesting 2D mask:

mask = BoundaryDiscretizeRegion[#, {{0, 1}, {0, 1}},
MaxCellMeasure -> {1 -> .02}] &@
ImplicitRegion[
0.1 < x < 0.9 && 0.1 < y < 0.9 + 0.05 Sin[20 x], {x, y}];
r2d = DiscretizeGraphics[g2d, MaxCellMeasure -> {1 -> .01},
PlotRange -> All];
inside = RegionIntersection[r2d, mask]


Then I find the edge and points on the edge. Unfortunately RegionIntersection doesn't work with lines and points. Here is workaround

edge = DiscretizeRegion@*Line@*Intersection @@
Round[{Sort /@
MeshPrimitives[RegionIntersection[r2d, mask], 1][[;; , 1]],
Sort /@ MeshPrimitives[RegionDifference[r2d, mask], 1][[;; ,
1]]}, .0001];
points = DiscretizeRegion@*Point@*Intersection @@
Round[{MeshPrimitives[RegionDifference[r2d, mask], 0][[;; , 1]],
MeshPrimitives[RegionDifference[mask, r2d], 0][[;; ,
1]]}, .0001];


Then I want to make RegionProduct to create 3D regions from corresponding 2D regions. I also have to use hand-written workaround

regionProduct[reg_, join_: True, y1_: 0, y2_: 1] :=
Module[{n = MeshCellCount[reg, 0]},
MeshRegion[Join @@ (ArrayFlatten@{{#[[;; , ;; 1]], #2, #[[;; , 2 ;;]]}} &[
MeshCoordinates@reg, #] & /@ {y1, y2}), {MeshCells[reg, _],
MeshCells[reg, _] /. p : {__Integer} :> p + n,
If[join, MeshCells[reg, _] /. {(Polygon | Line)[
p_] :> (Polygon@Join[#, Reverse[#, 2] + n, 2] &@
Partition[p, 2, 1, 1]),
Point[p_] :> Line@{p, p + n}}, ## &[]]}]];
inside3d = regionProduct[inside, False];
edge3d = regionProduct@edge;
points3d = regionProduct@points;


The result is impressive

toGC[reg_, dim_] := GraphicsComplex[MeshCoordinates@reg, MeshCells[reg, dim]];

Graphics3D[{FaceForm@Lighter[Blue, 0.7], toGC[inside3d, 2],
EdgeForm[], toGC[edge3d, 2], toGC[points3d, 1], Lighter@Blue,
ScalingTransform[0.999 {1, 1, 1}, RegionCentroid@mask3d]]},
Lighting -> "Neutral", Boxed -> False]


Also with transparency:

Graphics3D[{FaceForm@Lighter[Blue, 0.7],
toGC[regionProduct[RegionBoundary@inside, False], 1], EdgeForm[],
toGC[regionProduct@inside, 2], toGC[edge3d, 2], toGC[points3d, 1],
Blue, Opacity[0.03],
ScalingTransform[0.999 {1, 1, 1} #, RegionCentroid@mask3d] & /@
Range[0, 1, 0.01]]}, Lighting -> "Neutral", Boxed -> False,
BaseStyle -> {RenderingOptions -> {"DepthPeelingLayers" -> 100}}]


I hope future versions will do it more automatically.

• @dimitris, it is a cross section of the circle near its edge. Circles have black edge with nonzero thickness. After texture interpolation it smooths and becomes gray. You can find more accurate random realization of circles. Commented Nov 5, 2015 at 19:29
• Can you explain coords[[;; , {1, 3}]]? I see it yields {{0, 0}, {0, 1}, {0, 0}, {0, 1}, {1, 0}, {1, 1}, {1, 0}, {1, 1}}, but I don't know why. Thanks! Commented Nov 6, 2015 at 1:17
• @JosephO'Rourke It's projection of coordinates to x-z plane Commented Nov 6, 2015 at 1:46
• @dimitris It is number of coords, which makes faces of the cube. You are right, I rotated the graphics manually. Now regionProduct expands in the y direction for the proper orientation. Commented Nov 6, 2015 at 11:01
• Very very nice +1
– eldo
Commented Nov 6, 2015 at 11:15

Although I think that some of the options presented on the other answers are better, this one is probably the one with the shortest code on the 3D aspect:

Starting from any of the 2D generated graphics, from the other answers, hereafter named g2d:

Image3D[{Rasterize[g2d]}, BoxRatios -> 1]


and a lot of tweaking is possible with the Image3D options.

• Thank you very much! I think my collegue who is afraid of large codes will love this solution:-)! Commented Nov 6, 2015 at 14:59
• @dimitris and I just reduced a few more characters... (I recently participated on the one-liner... :-) Commented Nov 6, 2015 at 15:16
• I don't know why but I cannot reproduce your result. I got error message: Image3D::imgarray: The specified argument... Commented Nov 6, 2015 at 16:27
• @dimitris You have to put the image inside a list. Have you done that? {g2d} I'm on 10.3.0. If it doesn't work, try to repeat the image {g2d, g2d} Commented Nov 6, 2015 at 16:35
• Same error message with Image3D[{g2d}, BoxRatios -> 1] & Image3D[{g2d, g2d}, BoxRatios -> 1]. Strangely the Documentation Center says Introduced in 2012 (9.0) | Updated in 2014 (10.0) for this function (reference.wolfram.com/language/ref/Image3D.html) Commented Nov 6, 2015 at 16:39

How about DensityPlot3D or ListDensityPlot3D.

DensityPlot3D[
Sin[\[Pi] z] + Cos[\[Pi] x], {x, -6, 6}, {y, -4, 4}, {z, -4.5, 4.5},
OpacityFunction -> Function[f, If[f > 1, 1, .5]], OpacityFunctionScaling -> False,
ColorFunction -> Function[f, If[f > 1, Gray, Blue]], ColorFunctionScaling -> False]


• If you notice I want the cylinders randomly distributed. Thanks anyway for the answer. Commented Nov 5, 2015 at 22:45
• It is an interesting idea. May be you can use NearestFunction with some random points. Commented Nov 6, 2015 at 3:29
• @ybektukov: You are right. The solution of Edmund is very nice and I agree it can be further extended. Commented Nov 6, 2015 at 8:47

Nothing new here. Just a workaround following the new things I learn from ybeltukov amazing answer (and copy them:-)!).

Use the code below to generate randomly distributed but not intersecting disks

distinct[n_, r_] :=
Module[{d, f, p}, d = {Disk[RandomReal[{0, 1}, 2], r]};
Do[f = RegionDistance[RegionUnion @@ d];
While[p = RandomReal[{0, 1}, 2]; f[p] < r];
d = Append[d, Disk[p, r]], {n - 1}];
d]

disks // Clear
SeedRandom[159]
disks = distinct[75, 0.03];

g2d = Graphics[{FaceForm@Lighter[Blue, 0.8],
EdgeForm@Directive[Thickness[0.004], Black], disks},
PlotRange -> {{0, 1}, {0, 1}}, Background -> Lighter@Blue]


Texture...

pad = 0.1;
coords = Tuples[{0, 1}, 3];
cube = Polygon[{{1, 3, 7, 5}, {1, 5, 6, 2}, {5, 7, 8, 6}, {7, 3, 4,
8}, {3, 1, 2, 4}, {6, 8, 4, 2}}];
vtc = pad + (1 - 2 pad) coords[[;; , {1, 3}]];


and the 3D final drawing...

Graphics3D[{Texture[g2d],
GraphicsComplex[coords, cube, VertexTextureCoordinates -> vtc]},
Lighting -> "Neutral", Boxed -> False, ImageSize -> 500]


And a bit of cylinders now...(do not expect something advanced here:-)!)

mycylinders =
Cylinder[{{#[[1]], #[[2]], -3}, {#[[1]], #[[2]], 3}}, 0.03] & /@
Map[First, disks];

cub = Graphics3D[{FaceForm@Blue, Opacity[0.2],
Cuboid[{0, 0, -3}, {1, 1, 3}]}, Lighting -> "Neutral"];

Show[{Graphics3D[{FaceForm@Lighter[Blue, 0.7], Lighting -> "Neutral",
BaseStyle -> {RenderingOptions -> {"DepthPeelingLayers" -> 100}},
mycylinders}], cub}, Boxed -> False,
PlotRange -> {{0, 1}, {0, 1}, Automatic}]


I guess the final drawing can be made much better...

• Just a side note: BaseStyle -> {RenderingOptions -> {"DepthPeelingLayers" -> 100}} is a specific option to plot complicated semi-transparent graphics (e.g. in quasi-volumetric rendering). Commented Nov 6, 2015 at 12:05
• Thanks for the comment. I will modify it! Commented Nov 6, 2015 at 12:11