Symbolic Computing with Geometric Regions

Adam Strzeboński, Wolfram Research

Paweł Strzeboński

Join the Conversation #WolframTechConf

Region types

A region is a subset R⊆(ℝ^n) for some n∈(SubPlus[ℤ]).

Basic (special) regions

In Mathematica 10.3 there are 36 special region types supported.

These 10 special region types are new in Mathematica 10.3:

• RegularPolygon
• DiskSegment
• Annulus
• StadiumShape
• SphericalShell
• CapsuleShape
• AffineSpace
• AffineHalfSpace
• Hyperplane
• HalfSpace

The remaining 26 special region types:

• Point
• Line
• Interval
• Triangle
• Rectangle
• Parallelogram
• Polygon
• RegularPolygon
• Disk
• Circle
• DiskSegment
• Annulus
• StadiumShape
• Tetrahedron
• Pyramid
• Prism
• Hexahedron
• Cylinder
• Cone
• SphericalShell
• CapsuleShape
• Simplex
• Cuboid
• Parallelepiped
• Ball
• Sphere
• Ellipsoid
• InfiniteLine
• HalfLine
• InfinitePlane
• HalfPlane
• AffineSpace
• AffineHalfSpace
• Hyperplane
• HalfSpace
• ConicHullRegion

Formula regions

• ImplicitRegion
• ParametricRegion

Derived regions

• RegionUnion
• RegionIntersection
• RegionDifference
• RegionSymmetricDifference
• BooleanRegion
• RegionProduct
• TransformedRegion
• InverseTransformedRegion
• RegionBoundary

Mesh regions

• MeshRegion
• BoundaryMeshRegion

Symbolic computation with regions -- general vs. specialized methods

Region properties

• Element
• RegionMember
• RegionEmbeddingDimension
• RegionDimension
• RegionMeasure
• RegionCentroid
• ArcLength
• Area
• Volume
• RegionDistance
• SignedRegionDistance
• RegionNearest
• BoundedRegionQ
• RegionBounds

Regions as input to symbolic computation functions

• Reduce
• Solve
• FindInstance
• Resolve
• Minimize
• Maximize
• MinValue
• MaxValue
• ArgMin
• ArgMax
• Integrate
• assumptions

General methods

The general methods for symbolic computation with regions work by finding systems of equations and inequalities describing the region and then using general methods for solving systems of equations and inequalities, quantifier elimination, or constrained optimization. In general, derived regions do not autosimplify, however we can construct systems of equations and inequalities describing them (if the region involves ParametricRegion or TransformedRegion the corresponding system may involve quantifiers).

In principle, all problems involving semialgebraic regions can be solved using the cylindrical algebraic decomposition (CAD) algorithm. In practice, in some cases the algorithm may be very slow. Also, in computations involving integration, algebraic function bounds produced by the CAD algorithm may not be covenient for Integrate.

Specialized methods

For basic regions we often have either closed-form formulas for region properties or there are special case algorithms that can compute the property much faster than the general algorithm. Symbolic computation functions, notably Integrate, can also use faster special case algorithms for basic regions. Derived regions with basic region components often can autosimplify to basic regions, hence making further computation with the region easier.

This computes the volume of the unit ball using the general algorithm:

R=ImplicitRegion[x^2+y^2+z^2≤1,{x,y,z}];
Volume[R]//Timing

This uses a closed-form formula to find the result:

Volume[Ball[3]]//Timing

Implementation of specialized methods

Mathematica 10.3 contains a collection of formulas and special case algorithms computing properties of basic regions, operations on basic regions, and relations between basic regions.

Region properties

There are 20 properties of basic regions and 12 properties of basic region boundaries currently supported. Some of these directly correspond to region property functions, and some give properties that are helpful in symbolic computation algorithms.

Properties of the region:

“Assumptions”
“SpaceDimension”
“RegionDimension”
“ImplicitDescription”
“ParametricDescription”
“InjectiveParametricDescription”
“BoundedQ”
“ConvexQ”
“BoundingBox”
“Instance”
“Measure”
“Centroid”
“Distance”
“Nearest”
“SignedDistance”
“Type”
“Primitive”
“Boundary”
“SimpleBoundary”
“SimplicialDecomposition”

Properties of boundary of the region :

“RegionDimension”
“ImplicitDescription”
“ParametricDescription”
“InjectiveParametricDescription”
“ConvexQ”
“Instance”
“Measure”
“Centroid”
“Distance”
“Nearest”
“SignedDistance”
“SimplicialDecomposition”

Region operations

There are two binary operations (“Intersection” and “Product”) and two unary operations (“Image” and “Preimage”) on basic regions currently supported. The operations are used in autosimplification of the corresponding derived regions.

Region relations

There are three binary relations (“Subset”, “Disjoint”, “Equal”) between basic regions currently supported. The relations are used in autosimplification of derived regions. Future versions will contain region relation predicates and basic region relation information will be used in autosimplification of these predicates.

Integration rules

Currently we have specialized integration rules for 33 basic and derived regions. These include code for fast integration of polynomials as well as special parametrization formulas that are more convenient for Integrate than algebraic parametrizations that would be computed by the general algorithm.

Implementation statistics

Currenly the implementation of specialized basic region methods consists of more than 17,000 lines of Mathematica code. The code implements 1752 different formulas and special case algorithms.

Autosimplification of Boolean operations

Examples using region intersection rules

R=HalfPlane[{2,3},{1,5},{-1,1}];
S=StadiumShape[{{2,2},{3,2}},1];
RS=RegionIntersection[R,S]

Graphics[{{Opacity[0.3],Yellow,R,Blue,S},{Green,RS}}]

R=Cuboid[];
S=InfinitePlane[{1/2,1/2,1/2},{{3,2,1},{1,1,0}}];
RS=RegionIntersection[R,S]

Graphics3D[{{Opacity[0.3],Yellow,R,Blue,S},{Green,RS}},PlotRange->{{-0.5,1.5},{-0.5,1.5},{-0.5,1.5}}]

Examples using region relation rules

R=Disk[{1/2,1/2}];
S=RegularPolygon[{0,0},2,7];
RS=RegionUnion[R,S]

Graphics[{{Opacity[0.3],Red,R,Yellow,S,Green,RS}}]

R=Ellipsoid[{3,2,1},{1,2,3}];
S=Cuboid[];
RS=RegionDifference[R,S]

Graphics3D[{{Opacity[0.3],Yellow,R,Blue,S},{Green,RS}}]

Autosimplification makes further computation with the region much easier. This uses a specialized algorithm computing distance an ellipsoid:

ans1=RegionDistance[RegionDifference[R,S],{1,2,3}]//Timing

Here RegionDifference does not autosimplify and RegionDistance uses the general method based on constrained optimization:

iS=ImplicitRegion@@{RegionMember[S,{x,y,z}],{x,y,z}};
ans2=RegionDistance[RegionDifference[R,iS],{1,2,3}]//Timing

ans2/ans1//RootReduce

Autosimplification of RegionProduct

R=Disk[{1,2},3];
S=Interval[{0,1}];
RegionProduct[R,S]

Autosimplification makes further computation with the region much easier. This uses a specialized algorithm for integration of polynomials over a cylinder:

ans1=Integrate[x^2y z+x y^2z+x y z^2,Element[{x,y,z},RegionProduct[R,S]]]//Timing

Here RegionProduct does not autosimplify and Integrate uses the general method.

iS=ImplicitRegion[0≤x≤1,x];
ans2=Integrate[x^2y z+x y^2z+x y z^2,Element[{x,y,z},RegionProduct[R,iS]]]//Timing

ans2/ans1

Autosimplification of TransformedRegion

R=Ball[3];
F=AffineTransform[{{{1,2,3},{4,7,5},{9,8,6}},{3,2,1}}];
TransformedRegion[R,F]

Volume[TransformedRegion[R,F]]//Timing

Here TransformedRegion does not autosimplify and the volume computation uses the general method:

iR=ImplicitRegion[(x^2)+(y^2)+(z^2)≤1,{x,y,z}];
Volume[TransformedRegion[iR,F]]//Timing

Distance properties

R=CapsuleShape[{{0,0,0},{1,2,3}},1];

Manipulate[With[{d=N[SignedRegionDistance[R,{x,y,z}]],np=N[RegionNearest[RegionBoundary[R],{x,y,z}]]},Show[Graphics3D[{{Opacity[0.5],Green,R},{Text[Style["signed distance = "<>ToString[PaddedForm[d,{3,3}]],Large],{0,2,5.5}],{Red,PointSize[Large],Point[{{x,y,z},np}]},{Red,Line[{np,{x,y,z}}]}}}],PlotRange->{{-1,4},{-1,4},{-1,4}}]],{{x,1},-1,4},{{y,1},-1,4},{{z,1},-1,4}]

Here we use specialized algorithms for signed distance and the nearest boundary point for CapsuleShape:

SignedRegionDistance[R,{1,1,1}]//Timing

RegionNearest[RegionBoundary[R],{1,1,1}]//Timing

With the general method the signed distance computation takes much longer:

iR=ImplicitRegion@@{RegionMember[R,{x,y,z}],{x,y,z}};
SignedRegionDistance[iR,{1,1,1}]//Timing

Computation of the nearest point on the boundary fails, since in Mathematica 10.3 we cannot compute equation and inequality description of the boundary of an arbitrary 3D semialgebraic region. We will be able to do it in a future version using and algorithm computing CAD cell adjacencies.

RegionNearest[RegionBoundary[iR],{1,1,1}]//Timing

The specialized algorithm can quickly compute a signed distance formula for a symbolic point. The general method does not finish in 12 hours.

SignedRegionDistance[R,{x,y,z}]//Timing

Measure properties

RegionMeasure

This finds the measure of a 5D ellipsoid using a closed-form formula:

R=Ellipsoid[{0,0,0,0,0},{1,2,3,4,5}];
RegionMeasure[R]//Timing

Computation of the measure using the general method takes much longer:

iR=ImplicitRegion@@{RegionMember[R,{x,y,z,u,v}],{x,y,z,u,v}};
RegionMeasure[iR]//Timing

This finds the area of the boundary of a 3D ellipsoid using a closed-form formula. Unfortunately, we have no formula for higher dimensions.

RegionMeasure[RegionBoundary[Ellipsoid[{0,0,0},{1,2,3}]]]//Timing

This finds an approximation of the measure of the boundary of R using the general method and numeric integration:

RegionMeasure[RegionBoundary[N[R]]]//Timing

RegionCentroid

This finds the centroid of a disk segment using a closed-form formula:

R=DiskSegment[{1,2},{3,4},{5,6}];
(c=RegionCentroid[R])//Timing

Graphics[{Green,R,Red,PointSize[Large],Point[c]}]

Computation of the centroid using the general method takes much longer and gives a more complicated result:

iR=ImplicitRegion@@{RegionMember[R,{x,y}],{x,y}};
(c1=RegionCentroid[iR])//Timing

N[c-c1,20]//Quiet

Integration

This uses a specialized algorithm for integration of polynomials over a cone:

R=Cone[{{1,2,3},{4,5,6}},7];
(∫_({x,y,z}∈R))((x^2)+(y^2)+(z^2))//Timing

The general integration method applied to an ImplicitRegion version of R fails because Integrate cannot compute the generated integrals with algebraic function bounds.

SeedRandom[1234];
pts=RandomInteger[{-100,100},{6,5}];
R=Simplex[pts];
(∫_({x,y,z,t,u}∈R))(x+2 (y^2)+3 (z^3)+4 (t^4)+5 (u^5))//Timing

The general integration method didn’t finish overnight.

Integrate a polynomial over a 7D unit ball:

R=Ball[7];
(∫_({x,y,z,t,u,v,w}∈R))((x^2) (y^4)+(z^2) (t^2) (u^2)+(v^4) (w^2))//Timing

The general integration method is much slower.

iR=ImplicitRegion@@{RegionMember[R,{x,y,z,t,u,v,w}],{x,y,z,t,u,v,w}};
(∫_({x,y,z,t,u,v,w}∈iR))((x^2) (y^4)+(z^2) (t^2) (u^2)+(v^4) (w^2))//Timing