# A simple way to discretize a linear region

I have a few hundred linear regions of dimensions 0, 1, 2 and 3 embedded in $R^3$ that I need to discretise. DiscretizeRegion does a good job with most of them, but there are some (of dimensions 2 and 3) for which it fails. As far as I can tell, the failures come down to the region being very thin in one dimension, which makes it tricky for DiscretizeRegion to figure out cell sizes. (But that's really just an uninformed guess.) Messing with MaxCellMeasure seems to fix the problem in 2D, but in 3D it just gives me a choice between failure and freezing the notebook. It seems like these regions should be easy to discretise because they're all linear -- I can't imagine any of them would need more than ten simplexes (which would be fine for my purpose of sampling random points).

What I'm looking for is either: A way to get DiscretizeRegion working on these regions; or an alternative approach to building a simple (perhaps containing the fewest possible simplexes) MeshRegion from the defining constraints or the corresponding ImplicitRegion.

Here's an example:

In:= cons = -(3/5) + a < c < -(8/5) + a + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) && 1/5 + a + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) < 0 && 4/5 + b + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) < a;
reg = ImplicitRegion[cons, {{a, -5, 5}, {b, -10, 10}, {c, -5, 5}}];
DiscretizeRegion[reg]

During evaluation of In:= DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region ImplicitRegion[<<2>>]. >>

Out= DiscretizeRegion[ImplicitRegion[-(3/5) + a <  c < -(8/5) + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) + a && 1/5 + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) + a < 0 && 4/5 + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) + b < a && -5 <= a <= 5 && -10 <= b <= 10 && -5 <= c <= 5, {a, b, c}]]


It appears to be a sensible region:

In:= RegionQ[reg]
RegionDimension[reg]
N@RegionMeasure[reg]

Out= True

Out= 3

Out= 0.331639


but I've tried everything I can think of using the built-in Region functionality without success. I would have thought that there would be some straightforward way to get the bounding vertices of a linear region, and then maybe use DelaunayMesh, or MeshRegion to discretise, but I haven't had any luck down that route, either.

Turns out it's not as hard as I thought to get a basic function working:

linearDiscretizeRegion[fullconstraints_, vars_] :=
Block[{closedregion, equations, points},
closedregion = ImplicitRegion[fullconstraints /. {Less -> LessEqual, Greater -> GreaterEqual}, Evaluate@vars];
equations = List@@LogicalExpand[fullconstraints] /. {Less -> Equal, Greater -> Equal, LessEqual -> Equal, GreaterEqual -> Equal};
points = vars /. (Solve[#, vars] & /@ Subsets[equations, {Length[vars]}] /. {} -> Nothing)[[;;, 1]];
DelaunayMesh[Pick[points, RegionMember[closedregion] /@ points]]
]


Which yields

In:= cons = -(3/5) + a < c < -(8/5) + a + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) && 1/5 + a + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) < 0 && 4/5 + b + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) < a;
fullcons = cons && -5 <= a <= 5 && -10 <= b <= 10 && -5 <= c <= 5;
mesh = linearDiscretizeRegion[fullcons, {a, b, c}];
RegionMeasure[mesh]
RegionDimension[mesh]

Out= 0.331639

Out= 3


But I'm sure there are better/faster implementations, or entirely different approaches.

Update: This is what I ended up using. It's seems to work well as far as it goes, but it's a bit messy and long-winded for my liking. The function colFast gives the interior points in a set of colinear points (which can then be romoved), and credit for that goes to @yarchik and @ciao from this post. The function linearDiscretizeRegion takes as its first argument a simply connected ImplicitRegion of dimensions 0, 1, 2 or 3, with embedding dimension 3, and for which the defining equations are linear and do not involve any strict inequalities (ie, the region is closed). Because of these limitations it may not be widely applicable, but in cases where it is applicable it can actually give a more accurate result than DiscretizeRegion (as shown in an example below).

The function:

colFast[pts_] := Module[{a = MapThread[Append, {pts, ConstantArray[1, Length[pts]]}]},
If[MatrixRank[a] == 2, Sort[pts][], {}]
];

linearDiscretizeRegion[closedregion_, vars_] :=
Block[{d = RegionDimension[closedregion], equations, points, inpoints},
equations = DeleteDuplicates[Simplify[Flatten[
BooleanMinimize[LogicalExpand[closedregion[]], "DNF"] /. {Or -> List, And -> List, LessEqual -> Equal, GreaterEqual -> Equal}
]]];

(* Solve defining equations *)
points = DeleteDuplicates[Pick[#, RegionMember[closedregion] /@ #], Norm[#1 - #2] < 10^-12 &] &@
(vars /. Select[
Simplify[Quiet[Solve[#, vars]] & /@ Subsets[equations, {3}] /. {} -> Nothing][[;; , 1]],
Length[#] == 3 &]);

inpoints = colFast /@ Subsets[points, {3}] /. {} -> Nothing;
points = Select[points, ! MemberQ[inpoints, #] &];

Which[
d == 0,
MeshRegion[points, {Point[{1}]}],

d == 1,
MeshRegion[points, {Line[{1, 2}]}],

d == 2, (* Mainly concerned with ordering the points by rotation around the centroid *)
Block[{numpts = Length[points], normal, centre, sortpoints},
normal = #/Norm[#] &@ FirstCase[
Cross @@@ Subsets[Subtract @@@ Subsets[points, {2}], {2}],
Except[{0, 0, 0}]];
centre = N@RegionCentroid[closedregion];
sortpoints = With[{mp = Mean[points[[{1, 2}]]]},
SortBy[points, ArcTan[Cross[# - centre, mp - centre].normal, (mp - centre).(# - centre)] &]];
MeshRegion[sortpoints, Triangle /@ Partition[Riffle[Range[Ceiling[numpts/2]], Range[numpts, Ceiling[numpts/2] + 1, -1]], 3, 1]]
],

d == 3, (* Mainly concerned with removing unnecessary points and invalid cells from a DelaunayMesh *)
Block[{validcells, newpoints},
validcells = With[{dm = DelaunayMesh[points]},
Pick[MeshPrimitives[dm, 3], # > 10^-12 & /@
PropertyValue[{dm, 3}, MeshCellMeasure]]];
newpoints = DeleteDuplicates[Flatten[validcells /. Tetrahedron -> List, 2]];
MeshRegion[newpoints, With[{cellcoords = #[]},
Tetrahedron[First@FirstPosition[newpoints, #] & /@ cellcoords]] & /@  validcells]
]
]
];


Applying it to the test case in the question gives:

closedcons = -(3/5) + a <= c <= -(8/5) + a + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) && 5 Sqrt[-2 + Sqrt] + 2 (1 + 5 a) ArcCsc[1 + 2 GoldenRatio] <= 0 && 4/5 + b + Sqrt[-2 + Sqrt]/(2 ArcCsc[1 + 2 GoldenRatio]) <= a && b >= -10 && c >= -5;
reg = ImplicitRegion[closedcons, {a, b, c}];
discreg = DiscretizeRegion[reg];
mydiscreg = linearDiscretizeRegion[reg, {a, b, c}];
HighlightMesh[mydiscreg, Style[1, Black, Thick]]

(*During evaluation of In:= DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region ImplicitRegion[<<2>>]. >>*) In the following example DiscretizeRegion does work but misses out points that should lie in the region due to cutting of sharp edges. (This could be mitigated to some extent by playing around with the options for DiscretizeRegion.)

closedcons = -10. <= a <= -1.2193719504789793 && -20. <= b <= 0.20387439009579586 (-8.923984761519574 + 4.904980951899468 a) && 0.20387439009579586 (6.961992380759787 + 4.904980951899468 a) <= c <= 10.;
impreg = ImplicitRegion[closedcons, {a, b, c}];
discreg = DiscretizeRegion[impreg];
mydiscreg = linearDiscretizeRegion[impreg, {a, b, c}];
GraphicsRow[{discreg, HighlightMesh[mydiscreg, Style[1, Black, Thick]]}] So both functions work, but with very different results. Moreover, discreg is actually missing some points near the edges which are included in mydiscreg. Taking a test point {a0, b0, c0}, satisfying the constraints closedcons, it turns out that the point is not contained in discreg but is contained in mydiscreg.

In:= testpoint = {a0, b0, c0} = {-1.2193719504789793, -3.0387439009579587, 10};
closedcons /. {a -> a0, b -> b0, c -> c0}
RegionMember[discreg, testpoint]
RegionMember[mydiscreg, testpoint]

Out= True

Out= False

Out= True


So while linearDiscretizeRegion may have a limited range of applicability, it appears to have its uses for a certain class of simple regions.