Re: Fläche sich überlagernder Objekte

[Udo und Susanne Krause <su.krause@XXXXXXX.ch>]

Hallo Dominik,

meine Idee mit der Triangulation ist etwas mühsam, die Idee von Rolf, 
Boole[] zu nutzen, ist besser.
Da ein Wegweiser nicht in die falsche Richtung zeigen sollte, habe ich 
diese Hausaufgabe kurz implementiert:

In[1] :=  <<DiscreteMath`ComputationalGeometry`
In[3]:= Clear[edgeCrossQ];
(* do these two edges cross, i.e. do they have a common point not being 
an end point *)
edgeCrossQ[e1_List, e2_List] :=
Module[{x1, x2, x3, x4, y1, y2, y3, y4, m},
 {{x1, y1}, {x2, y2}} = e1;
 {{x3, y3}, {x4, y4}} = e2;
 m = {{x1 - x2, x4 - x3}, {y1 - y2, y4 - y3}};
 If[(x1 == x2 && y1 == y3) || (x3 == x4 && y3 == y4) || Chop[Det[m]] == 0,
   Return[False],
   Length[Select[LinearSolve[m, {x1 - x3, y1 - y3}], (0 < # && # < 
1)&]] == 2
 ]
] /; Dimensions[e1] == {2, 2} && Dimensions[e2] == {2, 2} && (And @@ 
(NumericQ[#]& /@ Flatten[Join[e1, e2]]))

In[7]:= Clear[selfCrossingQ];
(* contains inefficiency, connected edges must not be tested in a 
non-toy implementation *)
selfCrossingQ[p_Polygon] := Or @@ Apply[edgeCrossQ, Subsets[Partition[#, 
2, 1, 1], {2}]& @@ p, {1}]

In[10]:= Clear[hezelConvexQ];
(* this allows lower dimensional artefacts and colinear edges *)
hezelConvexQ[p_Polygon] := (Length[ConvexHull @@ p] == Length[Union @@ 
p]) && Not[selfCrossingQ[p]]

In[13]:= Clear[hezelEdgeInequality]
hezelEdgeInequality[e_List, s_List] :=
Module[{x1, x2, y1, y2},
 {{x1, y1}, {x2, y2}} = e;
 {sx, sy} = s;
 If[Chop[N[x2 - x1]] == 0,
   If[sx >= x1, GreaterEqual[x, x1], LessEqual[x, x1]],
   If[sy >= ((y2 - y1)/(x2 - x1)) (sx - x1) + y1,
     GreaterEqual[y, ((y2 - y1)/(x2 - x1)) (x - x1) + y1],
     LessEqual[y, ((y2 - y1)/(x2 - x1)) (x - x1) + y1]
   ]
 ]
]

In[15]:= Clear[hezelPolygonalInequalities]
(* uses ugly implicit x and y as veriables *)
hezelPolygonalInequalities[p_Polygon] :=
Module[{edges, xS, yS},
 (* the center of mass gives orientation *)
 {xS, yS} = (Plus @@ (Identity @@ p))/(Length @@ p);
 (* one inequality per edge *)
 And @@ (hezelEdgeInequality[#, {xS, yS}]& /@ (Partition[#, 2, 1, 1]& 
@@ p))
]

In[20]:= Clear[hezelCoveredArea];
(* takes a list of polygons, stops if a non-convex polygon is there *)
hezelCoveredArea[l0_List] :=
Module[{escapeIt = False, o},
 (* input check, possibly a module condition *)
 For[i = 1, i <= Length[l0], i++,
   If[Not[hezelConvexQ[l0[[i]]]], Print["Non-convex: ", l0[[i]]]; 
escapeIt = True]
 ];
 If[escapeIt, Return[$Failed]];
 (* to get integration domain *)
 o = Transpose[Flatten[l0 /. Polygon -> Identity, 1]];
 (* generate inequalities and get result*)
 Integrate[Boole[Or @@ (hezelPolygonalInequalities[#]& /@ l0)],
   {x, Min[First[o]], Max[First[o]]}, {y, Min[Last[o]], Max[Last[o]]}]
] /; Length[Cases[l0, _Polygon]] == Length[l0]

Das war's schon. Die Beschränkung auf konvexe Polygone ist drin, weil 
bei einem
nicht-konvexen Polygon die Ungleichungen, die man von den Seiten 
(Kanten) bezieht,
nicht logisch mit AND verknüpft werden dürfen; verschiedene konvexe 
Bereiche (deren
gibt es mindestens 2 in einer nicht konvexen ebenen Figur) müssen 
logisch mit OR verknüpft
werden. Das würde hezelPolygonalInequalities[] schwieriger machen, da 
man eben in konvexe
Teilfiguren zerlegen muss. Just do it yourself.

Ein paar Beispiele,  nicht-konvexe Polygone
In[24]:=
hezelCoveredArea[{Polygon[{{0, 0}, {0, 1}, {1, 0}, {1, 1}}],
 Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}],
 Polygon[{{0, 0}, {1, 0}, {1/4, 1/4}, {0, 1}}]}]

From In[24]:=
Non-convex: Polygon[{{0,0},{0,1},{1,0},{1,1}}]
Non-convex: Polygon[{{0, 0}, {1, 0}, {1/4, 1/4}, {0, 1}}]

Out[24]=
$Failed

diese beiden sind gar nicht zusammenhängend:
In[23]:= hezelCoveredArea[{Polygon[{{0, 0}, {2, 0}, {1, 1}}], 
Polygon[{{0, -1}, {1, -2}, {2, -1}}]}]

Out[23]=
2

die Überlagerung von Objekten geht:
In[30]:=
hezelCoveredArea[{Polygon[{{0, 0}, {2, 0}, {1, 1}}], Polygon[{{0, 1}, 
{1, 0}, {2, 1}}]}]

Out[30]=
3/2

und ein leeres Quadrat in der Mitte wird auch nicht mitgerechnet:
In[29]:=
hezelCoveredArea[{Polygon[{{0, 0}, {2, 0}, {1, 1}}],
                 Polygon[{{2, 0}, {3, -1}, {2, -2}}],
                 Polygon[{{0, -2}, {1, -3}, {2, -2}}],
                 Polygon[{{0, 0}, {0, -2}, {-1, -1}}]}]

Out[29]=
4

Versuchen Sie, diese Implementation zu crashen, verbessern und erweitern 
Sie sie und lassen Sie allenfalls die Gruppe von Fortschritten hören.

Mit den besten Grüssen
Udo.

Udo und Susanne Krause wrote:

> Hallo Dominik,
>
> das würde ich in 3 Schritten machen (lassen)
> (1) Aussenkontour feststellen, sofern die überdeckte Fläche einfach 
> zusammenhängend ist
> (2) Triangulieren
> (3) die Flächeninhalte der Dreiecke der Triangulation summieren
>
> Falls der Triangulierer nur konvexe Gebiete korrekt bearbeitet, dann
> (2.1) in konvexe Gebiete zerlegen (dazu müssen Teilmengenränder für 
> ConvexHull[] ein  Fixpunkt sein, sozusagen, das kann sehr kompliziert 
> werden, es wirklich auszuführen. d.h. Teilmengen zu Testen und die 
> ausgewählten Ränder zu ändern)
> (2.2) jedes Gebiet einzeln triangulieren
> (3) die Flächeninhalte der Dreiecke der Triangulationen summieren.
>
> Mathematica hält für 2D in DiscreteMath`ComputationalGeometry` 
> wesentliche Algorithmen bereit, DelaunayTriangulation[], ConvexHull[]. 
> Man kann auch im Internet nach Triangulierern suchen, selbstverständlich.
>
> Mit den besten Grüssen
> Udo.
>
> Hans.Dolhaine@XXXXXXX.com wrote:
>
> > Im Jahr 2000 hat es bereits eine Diskussion über Boolsche Operatoren auf
> > Polygonen gegeben. Ich habe damals eine einfache Version zur 
> > Realisierung
> > solcher Dinge an dmug geschickt.
> >
> > Mit freundlichen Grüßen
> >
> > Hans Dolhaine
> > _________________________________
> >
> > VTR-TS
> > Phone:      +49-211-797-4809
> > Fax:        +49-211-798-1853
> > Mobile:     0171 97 17 049
> > E-Mail:     Hans.Dolhaine@XXXXXXX.com
> >
> >
> >
> >
> >
> >