# Counting and extracting hit circles/triangles for randomly chosen points

Below, you see a triangle in which the incircle was inscribed. As we know, the center of this circle is the intersection of the angle bisectors. At this point, one could go further and find the points, where each of the angle bisectors intersect with the circle and create a tangent of the circle which will divide three new triangles from the big one. Each of the new triangles has again an incircle which can be created and cutting (from now on) only one triangle from each new one and repeating this process brings us to the picture:

First question: Is it possible create this figure in Mathematica? Unfortunately, I have no clue and cannot provide some code myself, but I have asked a similar question at TeX.SE and got some pretty nice answers.

What I would like to do in the end is not only to create the figure; I want to shoot random points at the figure and want to find out where the points hit the figure in a circle or not.

I would be delighted to hear some ideas or if possible to collect sample code for my problem.

-
ad-infinitum tasks aren't well fitted for computers –  belisarius Feb 22 at 1:04
@belisarius They are perfectly fitted for computers! The problem is you, because either you get bored waiting for the answer or you die of exhaustion when running back and forth to the computer store to buy more memory :-) –  halirutan Feb 22 at 1:09
@belisarius what about instead of ad-infinitum just a practical amount of recursion? Like until it cannot be readily seen by the naked eye. –  David Caliri Feb 22 at 1:10
@belisarius. I failed it, too, but I don't know where I fit into the ordering of those who have. –  m_goldberg Feb 22 at 1:28
It's much more clear to state the problem as calculating the ratio of area of the triangle and sum of areas of circles, which obviously form three convergent geometric series of circle areas around a incircle. Frankly, this sounds more like an elementary geometry than Mathematica question. Technically, you can formulate an equation that would perform test even for infinite amount of circles, but this sounds... quite awkward. –  kirma Feb 22 at 9:01

Since this is rather long, some might prefer a teaser of what is coming:

### Introduction

First of all, I don't really know why you make your figure inconsistent. I mean, from the first big triangle you separate three smaller triangles. Why don't you just repeat this process and inscribe a circle in each of the new triangles and again separate three new triangles from their corners?

I will show you how to do that and additionally, how you can create your figure. Depending on your experience some of the code might be a bit harder to read. The solution is that you take my code only as reference and re-implement everything in your own style. Nevertheless, I will use two kinds of data structure: triangle[{a,b,c}] is a triangle having the corner points a,b, and c. Please be aware that this is not a function! It will never do something, it just stays as it is and indicates that we work with a triangle. The same is true for circle[cent,r].

When we divide a triangle into its circle and its sub-triangles, we will change structure from

triangle[{a,b,c}]


to

triangle[{a,b,c}, circle[cent,r], {triangle[{a,..}],triangle[{b,..}],triangle[{c,...}}]


With this approach we keep all our triangles in a tree-like structure and it is really easy to work our way from the biggest triangle deeper and deeper to the smallest parts.

### Building the triangle/circle structure to a certain depth

First, we need to do the hard part, which is thinking about how you can compute the inscribed circle of a triangle when you have given the three points of it. This is pure geometry and the necessary calculations can be either looked up or you can make a small sketch and do it by hand. The following function takes a triangle and returns {{cx,cy}, r} where cx, cy is the center of the circle and r its radius:

inscribedCircle[triangle[{pa_, pb_, pc_}]] := Module[
{a = Norm[pc - pb], b = Norm[pa - pc], c = Norm[pa - pb], pp = 0.0, s = 0.0},
pp = a + b + c;
s = pp/2;
{(a pa + b pb + c pc)/pp, Sqrt[(s - a) (s - b) (s - c)/s]}
];


This function will be used in the next routine where we divide a given triangle into 3 new triangles and its incircle. Again, the calculation of the points where the new triangles are split from the bigger one is pure geometry, depends on the center and radius of the incircle and you probably should do the calculation for this by yourself. The key to the triangle split is given in the Function[{a,b,c},...] part. As result of this function you will get a triangle-circle-children mix as I have showed it in my second code block:

divide[t : triangle[coords : {pa_, pb_, pc_}]] := Block[{r, cc},
{cc, r} = inscribedCircle[t];
triangle[coords, circle[cc, r], #] &[
triangle[Prepend[#2, #1]] & @@@
Transpose[{{pa, pb, pc}, Function[{c, a, b},
a + Normalize[# - a]*(Norm[cc - a] - r)/
Cos[ArcCos[Dot[c - a, b - a]/(Norm[c - a]*Norm[b - a])]/
2] & /@ {c, b}] @@@
NestList[RotateLeft, {pc, pa, pb}, 2]}]]]


The last function buildTree is only for convenience so that we can sub-divide over and over again as often as we want. The level argument will help me track how often I have subdivided the figure

buildTree[t_triangle, level_Integer] := With[{subs = divide[t]},
subs /.
triangle[coords_, circ_, childs_List] :>
triangle[coords, circ, buildTree[#, level - 1] & /@ childs]
];
buildTree[t_triangle, 0] := t;


Let's test the stuff we have so far. I will take an example triangle and divide it exactly one time. Be aware that when you use exact integer coordinates, than Mathematica will calculate exact points which will pretty soon blow your machine because of the large roots. I will use machine numbers here:

buildTree[triangle[{{-1., -1.}, {3., 1.}, {2., 4.}}], 1]
(* triangle[{{-1., -1.}, {3., 1.}, {2., 4.}},
circle[{1.7285, 1.52667}, 1.0397],
{triangle[{{-1., -1.}, {0.435586, 1.39264}, {1.4957, 0.24785}}],
triangle[{{3., 1.}, {2.54056, 0.770278}, {2.83756, 1.48731}}],
triangle[{{2., 4.}, {2.50417, 2.48748}, {1.17972, 2.63286}}]}] *)


It seems to do something. Let's create some graphical output to see whether this is useful.

### Creating graphical output from the triangle structure

Let's start with the easy part which is to decide how exactly a triangle and a circle should look like. For this I make two functions drawTriangle and drawCircle. Since I want to color the the parts differently when they get smaller, these functions will have depth parameter which tells them how deep in the triangle structure we are. Therefore, the next functions draw exactly one triangle/circle and they have the knowledge how small (or deep) they are. Note that I only provide the points to of the triangles/circles to these functions:

With[{style = "BrightBands"},
drawTriangle[pts : {_, _, _}, depth_, maxDepth_] :=
With[{col = ColorData[style, depth/maxDepth]},
{FaceForm[Opacity[0.3, col]], EdgeForm[{Black}], Polygon[pts]}];
drawCircle[center_, r_, depth_, maxDepth_] :=
With[{col = ColorData[style, 1 - depth/maxDepth]},
{FaceForm[Opacity[0.3, col]], EdgeForm[{Darker[Gray]}],
Disk[center, r]}];
]


The next part is the fun part because showTree will draw (actually: create graphics primitives) a triangle structure completely by itself by walking through it and drawing all triangles/circles it finds. How is that possible? Very easy (ignore the first definition line, it only calculates the maxDepth at the beginning): When you want to draw a structure like we have

triangle[{a,b,c}, circle[cent,r], triangle[{a,..}],triangle[{b,..}],triangle[{c,...}}]


we draw a triangle {a,b,c}, then we draw a circle {cent,r} and then we start over with the exact same procedure for all child triangles, because they themselves will again have triangle coordinates, an inscircle and child-triangles. Thats it. When we see a triangle which has no incircle, we know this is the smallest one (because we haven't sub-divided it) and there we stop.

showTree[t_triangle] := Graphics[showTree[t, Depth[t]]];
showTree[arg : triangle[pts_, circle[cent_, r_], childs_],
maxDepth_] := With[{d = Depth[arg]},
{drawTriangle[pts, d, maxDepth], drawCircle[cent, r, d, maxDepth],
showTree[#, maxDepth] & /@ childs}
];
showTree[triangle[pts_], maxDepth_] := drawTriangle[pts, 0, maxDepth]


Let's test it and sub-divide a triangle 6 times and then draw it:

showTree[buildTree[triangle[{{-1., 1.}, {3., 1.}, {2., 4.}}], 6]]


If you want the figure you have presented, the only thing you need to do is stop the complete sub-division of the triangles after the first step and only subdivide the smaller triangles into one smaller triangle. We can implement this by adjusting the divide and buildTree function into

divideRetarded[t : triangle[coords : {pa_, pb_, pc_}]] :=
Block[{r, cc},
{cc, r} = inscribedCircle[t];
triangle[coords, circle[cc, r], {triangle[Prepend[#, pa]]}] &[
Function[{c, a, b},
a + Normalize[# - a]*(Norm[cc - a] - r)/
Cos[ArcCos[Dot[c - a, b - a]/(Norm[c - a]*Norm[b - a])]/
2] & /@ {c, b}][pc, pa, pb]]];

buildTreeRetarded[t_triangle, level_Integer] :=
With[{subs = divide[t]},
Block[{divide = divideRetarded},
subs /.
triangle[coords_, circ_, childs_List] :>
triangle[coords, circ, buildTree[#, level - 1] & /@ childs]
]
];


and then you get:

showTree[buildTreeRetarded[triangle[{{-1., 1.}, {3., 1.}, {2., 4.}}], 6]]


### Shooting points at a figure

Your last question was whether it is possible to test certain random points whether or not they hit a circle. If you have read and understand the showTree function, in particular how it walks by itself through the triangle-structure to draw it, you will understand the next easily, because it is exactly the same approach.

Let's say you want to test a random point p, then you start with the biggest (outer) triangle. p has to lie inside this triangle, otherwise p is outside the whole figure. So you test whether p is in the triangle and you test whether or not p is in the incircle of the biggest triangle. If this is not the case, then you have only a two possibilities:

• p is in one of the smaller sub-triangles
• p is in none of the smaller sub-triangles

In the first case, you know that now you have to test the sub-triangle where p is in with the exact same approach. In the second case, you instantly know that p must lie anywhere in the rest area of the bigger triangle, because it is not in the circle and not in the sub-triangles.

There is one last possibility you have to prepare for: When your triangle looks like triangle[{a,b,c}] then you know you have reached the smallest level of resolution and (if you use a fixed triangle tree) you cannot tell whether it is in the incircle of this triangle because you simply haven't sub-divided it. Then, your outcome is unknown. Let me point out, that it is really easy to implement a dynamically growing triangle structure which sub-divides further when you reach the highest resolution without a definitive result.

Since this article is already far too long, I will not show how to do this here. Nevertheless, if you understood that the following hitTest function walks your structure exactly the same way as your showTree function, the only thing which is left is to implement two routines which help you to test whether or not a point is inside a triangle or circle.

pointInTriangleQ[{{p1x_, p1y_}, {p2x_, p2y_}, {p3x_, p3y_}}, {px_,
py_}] :=
Function[#1 > 0 && #2 > 0 && #1 + #2 < 1][(p1y*p3x - p1x*p3y -
p1y*px + p3y*px + p1x*py - p3x*py)/
((-p1y)*p2x + p1x*p2y + p1y*p3x - p2y*p3x - p1x*p3y +
p2x*p3y),
((-p1y)*p2x + p1x*p2y + p1y*px - p2y*px - p1x*py + p2x*py)/
((-p1y)*p2x + p1x*p2y + p1y*p3x - p2y*p3x - p1x*p3y + p2x*p3y)]

pointInCircleQ[cent_, r_, pt_] := Sqrt[#.#] &[pt - cent] <= r


The following hitTest gives you back the triangle or circle the point is in, $Failed if the point is completely outside the structure and $Failed[triangle[..]] if you have reached on of the smallest triangles and you cannot determine whether it is in its incircle:

hitTest[triangle[pts_, circle[cent_, r_], childs_], pt_] :=
With[{inTri = pointInTriangleQ[pts, pt],
inCir = pointInCircleQ[cent, r, pt]},
If[inCir,
circle[cent, r],
If[inTri,
With[{inchilds = pointInTriangleQ[First[#], pt] & /@ childs},
If[Or @@ inchilds,
hitTest[#, pt] & @@ Pick[childs, inchilds],
triangle[pts]
]
],
$Failed] ]]; hitTest[t : triangle[_], _] :=$Failed[t]


Now you can create a triangle structure and test every point you like:

tree =
buildTree[triangle[{{0., 0.}, {1., 0}, {1., 1.}}], 5];
hitTest[tree, {.7, .6}]
(* triangle[{{0., 0.}, {1., 0}, {1., 1.}}] *)


Or you throw constantly random points and count how many of them are in a circle. Or, you create a small dynamic environment to test the function:

With[{col = Green},
draw[circle[cent_, r_]] :=
{FaceForm[Opacity[.3, col]],
EdgeForm[{Black}], Disk[cent, r]};
draw[triangle[pts_]] :=
{FaceForm[Opacity[.3, col]],
EdgeForm[{Black}], Polygon[pts]};
draw[$Failed[triangle[pts_]]] := {FaceForm[Opacity[.3, Red]], EdgeForm[{Black}], Polygon[pts]}; draw[___] := {} ]; With[{bg = showTree[tree]}, DynamicModule[{pt = {.5, .5}}, LocatorPane[Dynamic[pt], Dynamic[Show[bg, Graphics[draw[hitTest[tree, pt]]]]] ]] ] ` - Lovely! This is like a Shifrin, only with images (which makes it even better) :D – Yves Klett Feb 24 at 7:42 Extremely well written answer! – VividD Jun 20 at 8:34 Shooting points at the diagram and determining the fraction inside any of the nested circles may be the beginnings of a Monte Carlo determination of the area fraction of the infinite series of circles. If so (and it's probably a very big if), consider the triangle with side lengths$a$,$b$,$c$. The inradius$r$of the incircle is$r = \sqrt{(a+b-c)(a-b+c)(-a+b+c)/(a+b+c)} / 2$Set$d_1$equal to the length of the line segment from vertex 1, through the incentre, to the far side of the incircle, which it intersects perpendicularly. Similarly for$d_2$and$d_3$. The total area of all three infinite series of circles is CircleArea$ = \frac{\pi r}{4} \left( \frac{d_1^2}{d_1-r} + \frac{d_2^2}{d_2-r} + \frac{d_3^2}{d_3-r} \right) - 2 \pi r^2$The ratio of total circle area to triangle area is AreaFraction$ = \frac{\pi}{2(a+b+c)} \left( \frac{d_1^2}{d_1-r} + \frac{d_2^2}{d_2-r} + \frac{d_3^2}{d_3-r} - 8 r \right)\$

-