|
![]() |
|
Registriert seit: 14. Apr 2009 673 Beiträge |
#1
Also noch mal zur Triangulation.
So einfach wie Dejan Vu das meint, war es nicht. Man muß auch irgendwie den Rand/Umriss festlegen. Ebenso die Umrisse der Löcher. Das kann ganz schön verwirrend sein. War es jedenfalls bei mir. Ich konnte zum Schluß aber Flächenmomente von beliebigen Querschnitten ausrechnen. Dies Konzept "Fortscheitende Front" von jfheins würde mich auch interessieren. Hatte ich noch nicht gehört. -zumal es sich wie eine rus Militärtaktik anhört- Wofür brauchst Du das jetzt genau? Finite Elemente?
Achtung: Bin kein Informatiker sondern komme vom Bau.
|
![]() |
Registriert seit: 28. Feb 2011 Ort: Mannheim 1.384 Beiträge Delphi 10.4 Sydney |
#2
Ja, FE. Hab zur Zeit nur eine Lösung für Rechtecke. Wie hats du denn die gemeinsamen Lagerpunkte hingekriegt?
|
![]() |
Registriert seit: 14. Apr 2009 673 Beiträge |
#3
![]() Wie hats du denn die gemeinsamen Lagerpunkte hingekriegt?
![]() Das mit den Flächenmomenten habe ich einfach nur so aus Spaß dazuprogrmmiert, als ich mit dem Programmteil fertig war. War auch kaum noch Arbeit. Ich glaube aber, dass sich Flächenmomente anders einfacher rechnen lassen als mit Triangulationen.
Achtung: Bin kein Informatiker sondern komme vom Bau.
|
![]() |
Registriert seit: 28. Feb 2011 Ort: Mannheim 1.384 Beiträge Delphi 10.4 Sydney |
#4
Jens, kannst du mir sagen was der Autor hier macht?
Delphi-Quellcode:
Function TDelaunay.Triangulate(nvert: integer): integer;
//Takes as input NVERT vertices in arrays Vertex() //Returned is a list of NTRI triangular faces in the array //Triangle(). These triangles are arranged in clockwise order. var Complete: PComplete; Edges: PEdges; Nedge: integer; //For Super Triangle xmin: double; xmax: double; ymin: double; ymax: double; xmid: double; ymid: double; dx: double; dy: double; dmax: double; //General Variables i: integer; j: integer; k: integer; ntri: integer; xc: double; yc: double; r: double; inc: Boolean; begin //Allocate memory GetMem(Complete, sizeof(Complete^)); GetMem(Edges, sizeof(Edges^)); //Find the maximum and minimum vertex bounds. //This is to allow calculation of the bounding triangle xmin := Vertex^[1].X; ymin := Vertex^[1].Y; xmax := xmin; ymax := ymin; For i := 2 To nvert do begin if Vertex^[i].X < xmin then xmin := Vertex^[i].X; if Vertex^[i].X > xmax then xmax := Vertex^[i].X; if Vertex^[i].Y < ymin then ymin := Vertex^[i].Y; if Vertex^[i].Y > ymax then ymax := Vertex^[i].Y; end; dx := xmax - xmin; dy := ymax - ymin; if dx > dy then dmax := dx Else dmax := dy; xmid := Trunc((xmax + xmin) / 2); ymid := Trunc((ymax + ymin) / 2); //Set up the supertriangle //This is a triangle which encompasses all the sample points. //The supertriangle coordinates are added to the end of the //vertex list. The supertriangle is the first triangle in //the triangle list. Vertex^[nvert + 1].X := (xmid - 2 * dmax); Vertex^[nvert + 1].Y := (ymid - dmax); Vertex^[nvert + 2].X := xmid; Vertex^[nvert + 2].Y := (ymid + 2 * dmax); Vertex^[nvert + 3].X := (xmid + 2 * dmax); Vertex^[nvert + 3].Y := (ymid - dmax); Triangle^[1].vv0 := nvert + 1; Triangle^[1].vv1 := nvert + 2; Triangle^[1].vv2 := nvert + 3; Triangle^[1].Precalc := 0; Complete[1] := False; ntri := 1; //Include each point one at a time into the existing mesh For i := 1 To nvert do begin Nedge := 0; //Set up the edge buffer. //if the point (Vertex(i).X,Vertex(i).Y) lies inside the circumcircle then the //three edges of that triangle are added to the edge buffer. j := 0; repeat j := j + 1; if Complete^[j] <> True then begin inc := InCircle(Vertex^[i].X, Vertex^[i].Y, Vertex^[Triangle^[j].vv0].X, Vertex^[Triangle^[j].vv0].Y, Vertex^[Triangle^[j].vv1].X, Vertex^[Triangle^[j].vv1].Y, Vertex^[Triangle^[j].vv2].X, Vertex^[Triangle^[j].vv2].Y, xc, yc, r, j); //Include this if points are sorted by X if (xc + r) < Vertex[i].X then complete[j] := True Else if inc then begin Edges^[1, Nedge + 1] := Triangle^[j].vv0; Edges^[2, Nedge + 1] := Triangle^[j].vv1; Edges^[1, Nedge + 2] := Triangle^[j].vv1; Edges^[2, Nedge + 2] := Triangle^[j].vv2; Edges^[1, Nedge + 3] := Triangle^[j].vv2; Edges^[2, Nedge + 3] := Triangle^[j].vv0; Nedge := Nedge + 3; Triangle^[j].vv0 := Triangle^[ntri].vv0; Triangle^[j].vv1 := Triangle^[ntri].vv1; Triangle^[j].vv2 := Triangle^[ntri].vv2; Triangle^[j].PreCalc := Triangle^[ntri].PreCalc; Triangle^[j].xc := Triangle^[ntri].xc; Triangle^[j].yc := Triangle^[ntri].yc; Triangle^[j].r := Triangle^[ntri].r; Triangle^[ntri].PreCalc := 0; Complete^[j] := Complete^[ntri]; j := j - 1; ntri := ntri - 1; end; end; until j >= ntri; // Tag multiple edges // Note: if all triangles are specified anticlockwise then all // interior edges are opposite pointing in direction. For j := 1 To Nedge - 1 do begin if Not (Edges^[1, j] = 0) And Not (Edges^[2, j] = 0) then begin For k := j + 1 To Nedge do begin if Not (Edges^[1, k] = 0) And Not (Edges^[2, k] = 0) then begin if Edges^[1, j] = Edges^[2, k] then begin if Edges^[2, j] = Edges^[1, k] then begin Edges^[1, j] := 0; Edges^[2, j] := 0; Edges^[1, k] := 0; Edges^[2, k] := 0; end; end; end; end; end; end; // Form new triangles for the current point // Skipping over any tagged edges. // All edges are arranged in clockwise order. For j := 1 To Nedge do begin if Not (Edges^[1, j] = 0) And Not (Edges^[2, j] = 0) then begin ntri := ntri + 1; Triangle^[ntri].vv0 := Edges^[1, j]; Triangle^[ntri].vv1 := Edges^[2, j]; Triangle^[ntri].vv2 := i; Triangle^[ntri].PreCalc := 0; Complete^[ntri] := False; end; end; end; //Remove triangles with supertriangle vertices //These are triangles which have a vertex number greater than NVERT i := 0; repeat i := i + 1; if (Triangle^[i].vv0 > nvert) Or (Triangle^[i].vv1 > nvert) Or (Triangle^[i].vv2 > nvert) then begin Triangle^[i].vv0 := Triangle^[ntri].vv0; Triangle^[i].vv1 := Triangle^[ntri].vv1; Triangle^[i].vv2 := Triangle^[ntri].vv2; i := i - 1; ntri := ntri - 1; end; until i >= ntri; Triangulate := ntri; //Free memory FreeMem(Complete, sizeof(Complete^)); FreeMem(Edges, sizeof(Edges^)); end; |
![]() |
Registriert seit: 14. Apr 2009 673 Beiträge |
#5
Ich glaube, er macht erst so ein Supertriangle. Das soll wahrscheinlich den äußeren Rand beschreiben, was aber ein Problem ist, wenn Du eine bestimmte Außenkontur haben willst und Punkte für ein Loch hast. Ansonsten wird er da irgendwie den Algo durchlaufen. Das sehe ich auch nicht so auf den ersten Blick.
Achtung: Bin kein Informatiker sondern komme vom Bau.
|
![]() |
Registriert seit: 28. Feb 2011 Ort: Mannheim 1.384 Beiträge Delphi 10.4 Sydney |
#6
Ich habe es zwar schon gemacht, aber so komplett ohne Ansatz und Ausgangslage ist sinnvolle Hilfe nicht möglich. Fang an, und bei konkreten Problemen sind wir sicher gerne dabei - ich zumindest, weil ich finde solche Themen spannend
![]() Ich glaube, er macht erst so ein Supertriangle. Das soll wahrscheinlich den äußeren Rand beschreiben, was aber ein Problem ist, wenn Du eine bestimmte Außenkontur haben willst und Punkte für ein Loch hast.[..]
Delphi-Quellcode:
unit uDelaunay;
interface uses SysUtils, Dialogs, Graphics, Types, Math, uUtils, uClasses; type TDelaunayTriangle = record public Indices: array [0..2] of integer; // Indices in der Knotenliste FPoints; procedure Clear; end; TDelaunayTriangles = class private procedure SetCount(const Value: integer); procedure Flip(const Index1, Index2, Edge1, Edge2: integer); public Item: array of TDelaunayTriangle; function Center(const Index: integer; var R: double; Nodes: TFloatPoints): TFLoatPoint; function IndexOf(const A, B, C: integer): integer; function IndexOfPtInCircumcircle(const NodeIndex: integer; Nodes: TFloatPoints): integer; function IndexOfPtInCircumcircleEx(const NodeIndex, StarIndex: integer; Nodes: TFloatPoints): integer; function IndexOfEdge(const NodeIndex, Index: integer; const Nodes: TFloatPoints): integer; procedure Add(const A, B, C: integer); procedure AddItem(const Value: TDelaunayTriangle); procedure Delete(const Index: integer); procedure Assign(Value: TDelaunayTriangles); procedure Clear; function Count: integer; destructor Destroy; override; end; TDelaunayTriangulation = class private FBorder, FPoints: TFloatPoints; FTriangles: TDelaunayTriangles; procedure CheckCircumcircles; procedure Triangulate; procedure TriangulatePoint(const Index: integer); public procedure Clear; procedure Init; procedure AddPoint(const Value: TFloatPoint); procedure AddBorderPoint(const Value: TFloatPoint); procedure AddPoints(Value: TFloatPoints); procedure AddBorderPoints(Value: TFloatPoints); procedure Draw(Canvas: TCanvas; const ppMM: double); property Points: TFloatPoints read FPoints; // mm; property Border: TFloatPoints read FBorder; // mm; property Triangles: TDelaunayTriangles read FTriangles; constructor Create; destructor Destroy; override; end; implementation function CircumcircleCenter(const A, B, C: TFloatPoint; var R: double): TFloatPoint; var M1, M2, T1, T2: TFloatPoint; begin M1 := FloatPoint((A.X + B.X) / 2, (A.Y + B.Y) / 2); // Mitte AB; M2 := FloatPoint((A.X + C.X) / 2, (A.Y + C.Y) / 2); // Mitte AC; T1 := Util_RotateFloatPoint(A, M1, 0, 1); // A um M1 90° drehen; 0 = Cos(90); 1 = Sin(90); T2 := Util_RotateFloatPoint(A, M2, 0, 1); // A um M2 90° drehen; .. if Util_IntersectLinesEx(M1, T1, M2, T2, Result) then // Schnittpunkt M1T1, M2T2; R := Util_FloatPointDistance(A, Result) else R := 0; end; function PtInCircle(const P, Center: TFloatPoint; const R: double): boolean; begin Result := Util_PtInEllipse(P, Center, R, R); end; { TDelaunayTriangle } procedure TDelaunayTriangle.Clear; begin Indices[0] := -1; Indices[1] := -1; Indices[2] := -1; end; { TDelaunayTriangles } destructor TDelaunayTriangles.Destroy; begin Clear; inherited Destroy; end; procedure TDelaunayTriangles.SetCount(const Value: integer); begin SetLength(Item, Value); end; procedure TDelaunayTriangles.Assign(Value: TDelaunayTriangles); var I: integer; begin Clear; for I := 0 to Value.Count - 1 do AddItem(Value.Item[I]); end; procedure TDelaunayTriangles.Clear; begin SetCount(0); end; function TDelaunayTriangles.Count: integer; begin Result := Length(Item); end; procedure TDelaunayTriangles.Add(const A, B, C: integer); begin SetCount(Count + 1); Item[Count - 1].Indices[0] := A; Item[Count - 1].Indices[1] := B; Item[Count - 1].Indices[2] := C; end; procedure TDelaunayTriangles.AddItem(const Value: TDelaunayTriangle); begin SetCount(Count + 1); Item[Count - 1] := Value; end; procedure TDelaunayTriangles.Delete(const Index: integer); var I: integer; begin for I := Index to Count - 2 do Item[I] := Item[I + 1]; SetCount(Count - 1); end; function TDelaunayTriangles.IndexOf(const A, B, C: integer): integer; var I: integer; B1, B2, B3, B4, B5, B6: boolean; begin Result := -1; for I := 0 to Count - 1 do begin B1 := (Item[I].Indices[0] = A) and (Item[I].Indices[1] = B) and (Item[I].Indices[2] = C); B2 := (Item[I].Indices[0] = A) and (Item[I].Indices[1] = C) and (Item[I].Indices[2] = B); B3 := (Item[I].Indices[0] = B) and (Item[I].Indices[1] = A) and (Item[I].Indices[2] = C); B4 := (Item[I].Indices[0] = B) and (Item[I].Indices[1] = C) and (Item[I].Indices[2] = A); B5 := (Item[I].Indices[0] = C) and (Item[I].Indices[1] = A) and (Item[I].Indices[2] = B); B6 := (Item[I].Indices[0] = C) and (Item[I].Indices[1] = B) and (Item[I].Indices[2] = A); if B1 or B2 or B3 or B4 or B5 or B6 then begin Result := I; Break; end; end; end; function TDelaunayTriangles.IndexOfEdge(const NodeIndex, Index: integer; const Nodes: TFloatPoints): integer; var IndexA, IndexB, IndexC: integer; begin Result := -1; try IndexA := Item[Index].Indices[0]; IndexB := Item[Index].Indices[1]; IndexC := Item[Index].Indices[2]; if Util_SameFloatPoint(Nodes[NodeIndex], Nodes[IndexA]) then Result := 0 else if Util_SameFloatPoint(Nodes[NodeIndex], Nodes[IndexB]) then Result := 1 else if Util_SameFloatPoint(Nodes[NodeIndex], Nodes[IndexC]) then Result := 2; except on E: Exception do ShowMessage('Fehler vom Typ: ' + E.ClassName + ', Meldung: ' + E.Message); end; end; procedure TDelaunayTriangles.Flip(const Index1, Index2, Edge1, Edge2: integer); begin try except on E: Exception do ShowMessage('Fehler vom Typ: ' + E.ClassName + ', Meldung: ' + E.Message); end; end; function TDelaunayTriangles.IndexOfPtInCircumcircle(const NodeIndex: integer; Nodes: TFloatPoints): integer; begin Result := IndexOfPtInCircumcircleEx(NodeIndex, 0, Nodes); end; function TDelaunayTriangles.Center(const Index: integer; var R: double; Nodes: TFloatPoints): TFLoatPoint; var IndexA, IndexB, IndexC: integer; A, B, C: TFloatPoint; begin R := 0; Result.Clear; try IndexA := Item[Index].Indices[0]; IndexB := Item[Index].Indices[1]; IndexC := Item[Index].Indices[2]; A := Nodes[IndexA]; B := Nodes[IndexB]; C := Nodes[IndexC]; Result := CircumcircleCenter(A, B, C, R); except on E: Exception do ShowMessage('Fehler vom Typ: ' + E.ClassName + ', Meldung: ' + E.Message); end; end; function TDelaunayTriangles.IndexOfPtInCircumcircleEx(const NodeIndex, StarIndex: integer; Nodes: TFloatPoints): integer; var I: integer; C: TFloatPoint; R: double; begin Result := -1; try for I := StarIndex to Count - 1 do begin C := Center(I, R, Nodes); if PtInCircle(Nodes[NodeIndex], C, R) then Result := I; // No Break; end; except on E: Exception do ShowMessage('Fehler vom Typ: ' + E.ClassName + ', Meldung: ' + E.Message); end; end; { TDelaunayTriangulation } constructor TDelaunayTriangulation.Create; begin inherited; FPoints := TFloatPoints.Create; FBorder := TFloatPoints.Create; FTriangles := TDelaunayTriangles.Create; end; destructor TDelaunayTriangulation.Destroy; begin FPoints.Free; FBorder.Free; FTriangles.Free; inherited; end; procedure TDelaunayTriangulation.Init; // need Border; var xMin, xMax, yMin, yMax: double; begin // 0....1....2; // . . . . // . . . . // . . . . // 3.........4; xMin := FBorder.Left; yMin := FBorder.Top; xMax := FBorder.Right; yMax := FBorder.Bottom; FPoints.Clear; FPoints.AddXY(xMin, yMin); // 0; FPoints.AddXY(xMax / 2, yMin); // 1; FPoints.AddXY(xMax, yMin); // 2; FPoints.AddXY(xMin, yMax); // 3; FPoints.AddXY(xMax, yMax); // 4; FTriangles.Clear; FTriangles.Add(0, 1, 3); FTriangles.Add(1, 4, 3); FTriangles.Add(1, 2, 4); end; procedure TDelaunayTriangulation.Clear; begin FBorder.Clear; FPoints.Clear; FTriangles.Clear; end; procedure TDelaunayTriangulation.AddPoint(const Value: TFloatPoint); // Need Border; begin if (FPoints.IndexOf(Value) < 0) and (FBorder.PtInPolygon(Value)) then begin FPoints.Add(Value); TriangulatePoint(FPoints.Count - 1); end; end; procedure TDelaunayTriangulation.AddBorderPoint(const Value: TFloatPoint); begin if FBorder.IndexOf(Value) < 0 then FBorder.Add(Value); end; procedure TDelaunayTriangulation.AddPoints(Value: TFloatPoints); var I: integer; begin for I := 0 to Value.Count - 1 do AddPoint(Value[I]); end; procedure TDelaunayTriangulation.AddBorderPoints(Value: TFloatPoints); var I: integer; begin for I := 0 to Value.Count - 1 do if FBorder.IndexOf(Value[I]) < 0 then FBorder.Add(Value[I]); end; procedure TDelaunayTriangulation.CheckCircumcircles; // ??? var I, J, Index1, Index2, Edge1, Edge2: integer; begin EXIT; I := 0; J := 0; while I < FPoints.Count do begin while J < FTriangles.Count - 1 do begin Index1 := FTriangles.IndexOfPtInCircumcircleEx(I, J, FPoints); if Index1 > -1 then begin Index2 := FTriangles.IndexOfPtInCircumcircleEx(I, J + 1, FPoints); if Index2 > -1 then begin Edge1 := FTriangles.IndexOfEdge(I, Index1, FPoints); Edge2 := FTriangles.IndexOfEdge(I, Index2, FPoints); FTriangles.Flip(Index1, Index2, Edge1, Edge2); I := 0; J := -1; end; end; Inc(J); end; Inc(I); end; end; procedure TDelaunayTriangulation.Triangulate; var TriangleIndex, I, IndexA, IndexB, IndexC: integer; begin for I := 0 to FPoints.Count - 1 do begin TriangleIndex := FTriangles.IndexOfPtInCircumcircle(I, FPoints); if TriangleIndex > -1 then begin IndexA := FTriangles.Item[TriangleIndex].Indices[0]; IndexB := FTriangles.Item[TriangleIndex].Indices[1]; IndexC := FTriangles.Item[TriangleIndex].Indices[2]; FTriangles.Add(IndexA, IndexB, I); FTriangles.Add(IndexB, IndexC, I); FTriangles.Add(IndexC, IndexA, I); end; end; CheckCircumcircles; end; procedure TDelaunayTriangulation.TriangulatePoint(const Index: integer); var TriangleIndex, IndexA, IndexB, IndexC: integer; begin TriangleIndex := FTriangles.IndexOfPtInCircumcircle(Index, FPoints); if TriangleIndex > -1 then begin IndexA := FTriangles.Item[TriangleIndex].Indices[0]; IndexB := FTriangles.Item[TriangleIndex].Indices[1]; IndexC := FTriangles.Item[TriangleIndex].Indices[2]; FTriangles.Add(IndexA, IndexB, Index); FTriangles.Add(IndexB, IndexC, Index); FTriangles.Add(IndexC, IndexA, Index); end; CheckCircumcircles; end; procedure TDelaunayTriangulation.Draw(Canvas: TCanvas; const ppMM: double); var R: double; I, IndexA, IndexB, IndexC: integer; A, B, C, D: TPoint; Center: TFloatPoint; begin Canvas.Brush.Color := clWhite; Canvas.Brush.Style := bsSolid; Canvas.FillRect(Canvas.ClipRect); For I := 0 To FTriangles.Count - 1 do begin IndexA := FTriangles.Item[I].Indices[0]; IndexB := FTriangles.Item[I].Indices[1]; IndexC := FTriangles.Item[I].Indices[2]; A := Util_CanvasPoint(FPoints[IndexA], ppMM); B := Util_CanvasPoint(FPoints[IndexB], ppMM); C := Util_CanvasPoint(FPoints[IndexC], ppMM); Center := FTriangles.Center(I, R, FPoints); D := Util_CanvasPoint(Center, ppMM); // ShowMessage(Format('%d %d %d', [IndexA, IndexB, IndexC])); Canvas.Brush.Style := bsClear; Canvas.Polygon([A, B, C]); Canvas.Brush.Color := clWhite; Canvas.Brush.Style := bsSolid; Canvas.TextOut(A.X, A.Y, IntToStr(IndexA)); Canvas.TextOut(B.X, B.Y, IntToStr(IndexB)); Canvas.TextOut(C.X, C.Y, IntToStr(IndexC)); // Canvas.TextOut(D.X, D.Y, IntToStr(I)); A := Util_CanvasPoint(FloatPoint(Center.X - R, Center.Y - R), ppMM); B := Util_CanvasPoint(FloatPoint(Center.X + R, Center.Y + R), ppMM); Canvas.Brush.Style := bsClear; // Canvas.Ellipse(A.X, A.Y, B.X, B.Y); end; end; end. |
![]() |
Registriert seit: 10. Jun 2004 Ort: Garching (TUM) 4.579 Beiträge |
#7
Ist es eine Option, einen vorhandenen Netzgenerator herzunehmen?
Da habe ich zum Beispiel ![]() Ansonsten: Advancing-Frint geht wie folgt vor: Du unterteilst deine Außenkontur zunächst in Knoten und Kanten. Da kannst du ganz stumpf machen, dass du eine Wunschlänge definierst, und dann an jeder Polygonkante schaust, wie viele Unterknoten da denn eingezogen werden müssen. Advancing-Front geht jetzt schrittweise die Knoten durch, die auf der aktuellen Front liegen und verkleinert die Front (bzw. genauer: das eingeschlossene Gebiet). Wenn die Front leer ist, hast du das Gebiet vollständig vernetzt. Es gibt dabei drei Möglichkeiten:
![]() Das Netz, dass daraus hervoirgeht, kannst du dann natürlich auch noch auf Delaunay überprüfen und Kanten ggf. flippen. Allgemein würde ich dir aber empfehlen, einen fertigen Netzgenerator zu verwenden. Edit: Habe gerade noch eine gute Arbeit gefunden: ![]() Ab Seite 29. Geändert von jfheins ( 7. Sep 2014 um 16:46 Uhr) |
![]() |
Ansicht |
![]() |
![]() |
![]() |
ForumregelnEs ist dir nicht erlaubt, neue Themen zu verfassen.
Es ist dir nicht erlaubt, auf Beiträge zu antworten.
Es ist dir nicht erlaubt, Anhänge hochzuladen.
Es ist dir nicht erlaubt, deine Beiträge zu bearbeiten.
BB-Code ist an.
Smileys sind an.
[IMG] Code ist an.
HTML-Code ist aus. Trackbacks are an
Pingbacks are an
Refbacks are aus
|
|
Nützliche Links |
Heutige Beiträge |
Sitemap |
Suchen |
Code-Library |
Wer ist online |
Alle Foren als gelesen markieren |
Gehe zu... |
LinkBack |
![]() |
![]() |