DGL • Thema anzeigen - Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Moderator: DGL-Team

Seite 1 von 1

[ 7 Beiträge ]

Vorheriges Thema | Nächstes Thema

Autor

Nachricht

Wölfchen

Betreff des Beitrags: Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Verfasst: So Jan 17, 2010 16:28

DGL Member

Registriert: Do Mär 05, 2009 20:17
Beiträge: 284
Wohnort: Kaiserslautern

Huhu mal wieder,

heute möchte ich mich einem neuen Thema widmen, der Volumenberechnung eines beliebigen (geschlossenen) Mesh bestehend aus Dreiecksflächen.

Mehrstündiges Googeln hat mir nur wenig Treffer geliefert. Das hier sieht mir derzeit am vielversprechendsten aus:
http://www.geometrictools.com/Documentation/PolyhedralMassProperties.pdf

Dort heisst es unter 5:

Zitat:

5 Pseudocode
The pseudocode for computing the integrals is quite simple. The polyhedron vertices are passes as the array
p[]. The number of triangles is tmax. The array index[] has tmax triples of integers that are indices into
he vertex array. The return values are the mass, the center of mass, and the inertia tensor relative to the
center of mass. The code assumes that the rigid body has constant density 1. If your rigid body has constant
density D, then you need to multiply the output mass by D and the output inertia tensor by D.

Bis hierher hört sich das ziemlich brauchbar an, denn ein array of TVertex und ein Indexarray aus dem sich integer triples ziehen lassen habe ich ja sowieso zum Rendern.
Leider überfordert mich der Pseudocode dann doch:

Zitat:

MACRO Subexpressions(w0,w1,w2,f1,f2,f3,g0,g1,g2)
{
temp0 = w0+w1; f1 = temp0+w2; temp1 = w0*w0; temp2 = temp1+w1*temp0;
f2 = temp2+w2*f1; f3 = w0*temp1+w1*temp2+w2*f2;
g0 = f2+w0*(f1+w0); g1 = f2+w1*(f1+w1); g2 = f2+w2*(f1+w2);
}
void Compute (Point p[], int tmax, int index[], Real& mass, Point& cm, Matrix& inertia)
{
constant Real mult[10] = {1/6,1/24,1/24,1/24,1/60,1/60,1/60,1/120,1/120,1/120};
Real intg[10] = {0,0,0,0,0,0,0,0,0,0}; // order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx
for (t = 0; t < tmax; t++)
{
// get vertices of triangle t
i0 = index[3*t]; i1 = index[3*t+1]; i2 = index[3*t+2];
x0 = p[i0].x; y0 = p[i0].y; z0 = p[i0].z;
x1 = p[i1].x; y1 = p[i1].y; z1 = p[i1].z;
x2 = p[i2].x; y2 = p[i2].y; z2 = p[i2].z;
// get edges and cross product of edges
a1 = x1-x0; b1 = y1-y0; c1 = z1-z0; a2 = x2-x0; b2 = y2-y0; c2 = z2-z0;
d0 = b1*c2-b2*c1; d1 = a2*c1-a1*c2; d2 = a1*b2-a2*b1;
// compute integral terms
Subexpressions(x0,x1,x2,f1x,f2x,f3x,g0x,g1x,g2x);
Subexpressions(y0,y1,y2,f1y,f2y,f3y,g0y,g1y,g2y);
Subexpressions(z0,z1,z2,f1z,f2z,f3z,g0z,g1z,g2z);
// update integrals
intg[0] += d0*f1x;
intg[1] += d0*f2x; intg[2] += d1*f2y; intg[3] += d2*f2z;
intg[4] += d0*f3x; intg[5] += d1*f3y; intg[6] += d2*f3z;
intg[7] += d0*(y0*g0x+y1*g1x+y2*g2x);
intg[8] += d1*(z0*g0y+z1*g1y+z2*g2y);
intg[9] += d2*(x0*g0z+x1*g1z+x2*g2z);
}
for (i = 0; i < 10; i++)
intg[i] *= mult[i];
mass = intg[0];
// center of mass
cm.x = intg[1]/mass;
cm.y = intg[2]/mass;
cm.z = intg[3]/mass;
// inertia tensor relative to center of mass
inertia.xx = intg[5]+intg[6]-mass*(cm.y*cm.y+cm.z*cm.z);
inertia.yy = intg[4]+intg[6]-mass*(cm.z*cm.z+cm.x*cm.x);
inertia.zz = intg[4]+intg[5]-mass*(cm.x*cm.x+cm.y*cm.y);
inertia.xy = -(intg[7]-mass*cm.x*cm.y);
inertia.yz = -(intg[8]-mass*cm.y*cm.z);
inertia.xz = -(intg[9]-mass*cm.z*cm.x);
}

Deshalb dachte ich, ihr habt vielleicht ein ähnliches Problem schon gelöst oder schonmal gesehen...
Für sachdienliche Hinweise bin ich wie immer dankbar!

*** EDIT ***
habe eben noch eine Lösung für das Programm Blender gefunden... sieht garnicht so schwierig aus....

Code:

selection = Blender.Object.GetSelected()
l = len(selection)
 
if l == 0:
    print 'No Objects Selected'
elif l > 1:
    print 'More than one Object Selected'
else:
    obj = selection[0]
    obtype = obj.getType()
    if obtype == 'Mesh':
        obname = obj.getName()
        mymesh = NMesh.GetRawFromObject(obname)
        vtot = 0
        for f in mymesh.faces:
            fzn = f.normal[2]
            x1 = f.v[0].co[0]
            y1 = f.v[0].co[1]
            z1 = f.v[0].co[2]
            x2 = f.v[1].co[0]
            y2 = f.v[1].co[1]
            z2 = f.v[1].co[2]
            x3 = f.v[2].co[0]
            y3 = f.v[2].co[1]
            z3 = f.v[2].co[2]
            pa = 0.5*abs((x1*(y3-y2))+(x2*(y1-y3))+(x3*(y2-y1)))
            volume = ((z1+z2+z3)/3.0)*pa
            if fzn < 0:
                fzn = -1
            elif fzn > 0:
                fzn = 1
            else:
                fzn = 0
            vtot = vtot + (fzn * volume)
        print 'Volume = ', vtot, ' BU'
        print '-----------------------------------------------'
    else:
        print 'Object must be a Mesh'
 

Gruß
Wölfchen

Zuletzt geändert von Wölfchen am So Jan 17, 2010 17:03, insgesamt 1-mal geändert.

Nach oben

Coolcat

Betreff des Beitrags: Re: Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Verfasst: So Jan 17, 2010 17:00

DGL Member

Registriert: Do Dez 29, 2005 12:28
Beiträge: 2249
Wohnort: Düsseldorf
Programmiersprache: C++, C#, Java

Öhm, der Pseudocode ist doch recht einleuchtend, also der sollte sich nach Delphi (oder was auch immer) übersetzen lassen? Die Syntax ist ziemlich ähnlich zu C/C++. Gut, was der da rechnet verstehe ich auch nicht, aber das muss man ja auch nicht unbedingt wissen. Falls doch wird das aber doch sicher im Paper erklärt?

Das einzige was vielleicht als Delphianer nicht ganz klar ist, ist dieses Marco Subexpressions. Ein Marco funktioniert so, dass der Inhalt des Marcos von einem Präprozessor vor dem compilieren einfach an der entsprechenden Stelle eingefügt wird. Die Parameter des Marcos werden dabei einfach ersetzt. Also:

Code:

Subexpressions(x0,x1,x2,f1x,f2x,f3x,g0x,g1x,g2x);

wird vom Präprozessor in das folgende umgewandelt:

Code:

temp0 = x0+x1;
f1x = temp0+x2;
temp1 = x0*x0;
temp2 = temp1+x1*temp0;
f2x = temp2+x2*f1x;
f3x = x0*temp1+x1*temp2+x2*f2x;
g0 = f2x+x0*(f1x+x0);
g1x = f2x+x1*(f1x+x1);
g2x = f2x+x2*(f1x+x2);

_________________
Yeah! :mrgreen:

Nach oben

Wölfchen

Betreff des Beitrags: Re: Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Verfasst: So Jan 17, 2010 17:16

DGL Member

Registriert: Do Mär 05, 2009 20:17
Beiträge: 284
Wohnort: Kaiserslautern

Coolcat hat geschrieben:

Code:

Subexpressions(x0,x1,x2,f1x,f2x,f3x,g0x,g1x,g2x);

wird vom Präprozessor in das folgende umgewandelt:

Code:

temp0 = x0+x1;
f1x = temp0+x2;
temp1 = x0*x0;
temp2 = temp1+x1*temp0;
f2x = temp2+x2*f1x;
f3x = x0*temp1+x1*temp2+x2*f2x;
g0 = f2x+x0*(f1x+x0);
g1x = f2x+x1*(f1x+x1);
g2x = f2x+x2*(f1x+x2);

die ganzen Marcos sind eigentlich Macros, oder?

ist intg denn integer oder integral? und was bedeutet der ausdruck "+=" ist das soviel wie ":="?

Ich werde mal versuchen den blender scriptcode, den ich im ursprungsposting noch hinzueditiert hab, irgendwie in delphi zu bringen der erscheint mir irgendwie einfacher... (ich mach das ja alles (Also Programmieren) erst wieder seit 6 wochen und für mich sind so sachen echt harte Nüsse :mrgreen:

Deswegen bin ich ja so froh um jedes Einsteigerforum das ich finde )

Nach oben

Coolcat

Betreff des Beitrags: Re: Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Verfasst: So Jan 17, 2010 17:24

DGL Member

Registriert: Do Dez 29, 2005 12:28
Beiträge: 2249
Wohnort: Düsseldorf
Programmiersprache: C++, C#, Java

Zitat:

die ganzen Marcos sind eigentlich Macros, oder?

Äh, ja, natürlich...

Zitat:

ist intg denn integer oder integral?

Da steht "Real intg[10] = {0,0,0,0,0,0,0,0,0,0};"
Es handelt sich um ein Array des Typs "Real" mit der Länge 10 (von 0 bis 9). Der Name das Arrays ist "intg", was wahrscheinlich für "integral" stehen soll.

Zitat:

und was bedeutet der ausdruck "+=" ist das soviel wie ":="?

"a += b;" wäre in Delphi "a := a + (b);". Beachte die Klammern, also was immer "b" ist, es wird zuerst ausgewertet. Das gleiche gilt für "-=", "*=" usw...

Ein einfaches "=" ist eine Zuweisung, also in Delphi ":=".

_________________
Yeah! :mrgreen:

Nach oben

Wölfchen

Betreff des Beitrags: Re: Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Verfasst: So Jan 17, 2010 19:40

DGL Member

Registriert: Do Mär 05, 2009 20:17
Beiträge: 284
Wohnort: Kaiserslautern

so ich hab jetzt mal versucht den blender script dingens in delphi zu übertragen...

schaut so aus:

Code:

procedure TCW3D.Volumenberechnung(Sender:TObject);
var i,k:integer;
var fzn:double;
var x1,x2,x3:double;
var y1,y2,y3:double;
var z1,z2,z3:double;
var pa:double;
var volume:double;
var vtot:double;
 begin
 k:=0;
  for i:=0 to  length(TCW3d(instanzliste[0]).varray)-1 do begin
             fzn:=TCW3d(instanzliste[0]).narray[k].vz;
             if fzn<0 then fzn:=-1.0;
             if fzn>0 then fzn:=1.0;
             if fzn=0 then fzn:=0.0;
             x1 := TCW3d(instanzliste[0]).varray[k].vx;
             x2 := TCW3d(instanzliste[0]).varray[k+1].vx;
             x3 := TCW3d(instanzliste[0]).varray[k+2].vx;
             y1 := TCW3d(instanzliste[0]).varray[k].vy;
             y2 := TCW3d(instanzliste[0]).varray[k+1].vy;
             y3 := TCW3d(instanzliste[0]).varray[k+2].vy;
             z1 := TCW3d(instanzliste[0]).varray[k].vz;
             z2 := TCW3d(instanzliste[0]).varray[k+1].vz;
             z3 := TCW3d(instanzliste[0]).varray[k+2].vz;
             pa := 0.5*abs((x1*(y3-y2))+(x2*(y1-y3))+(x3*(y2-y1)));
             volume := ((z1+z2+z3)/3.0)*pa;
             vtot := vtot + (fzn * volume);
             k:=k+3;
 end;
  form1.panel5.caption:=floattostr(vtot);
 
 end;

das liefert für einige Teile die ich habe ganz glaubwürdige Resultate, aber für eine kugel mit r0,5
die laut allgemeiner Formelsammlung 4/3*pi*r^3 = 0,523598... haben sollte kommt da raus:

Zitat:

4.6265025410383E46

was er nach einem trunc(vtot) aber zu 0 macht :shock:

wenns ja jetzt 0,46 gewesen wär hätt ichs ja eventuell noch geschluckt wegen der ungenauigkeiten mit den Dreiecksflächen.
Mehr Tests dazu kann ich leider erst morgen machen, wenn ich Teile zum Testen erstellen und Ergebnisse vergleichen kann.

Aber was ich jetzt auch nicht hinbekommen habe ist die

Zitat:

4.6265025410383E46

einfach mal als Zahl mit 3 Nachkommastellen auszugeben... wie macht man das in delphi?
**** EDIT Sonntag Abend ****
huiiii wieder ne Stunde gegoogelt und alles Mögliche gelesen wie das angeblich gehen soll...
Das hier hat dann doch tatsächlich aus

Zitat:

4.6265025410383E46

ein nettes kleines

Zitat:

0,52

gezaubert (wieso auch immer)...

Code:

form1.panel5.caption:=floattostrF(vtot,ffFixed,8,2);

Damit steigt die Hoffnung das diese Volumenberechnung irgendwie richtig zu sein scheint.
Dann muß ich jetzt nur noch rausfinden wieso das ganze bei einigen Teilen zu einer Access Violation führt... es bleibt spannend! :mrgreen:

**** EDIT Montag Morgen ****
Hmm da hab ich mich wohl zu früh gefreut, irgendwie war das mit den 0,52 wohl zufall... mein PC rechnet jedesmal was andres... mal steht da NAN mal -INF mal

Zitat:

4.6265025410383E46

oder ähnlich sinnloser kram....

Nach oben

zoiX

Betreff des Beitrags: Re: Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Verfasst: Mo Jan 18, 2010 09:25

DGL Member

Registriert: Mo Aug 31, 2009 13:19
Beiträge: 151

Ich rate mal ins Blaue: Du hast vtot nicht initialisiert, richtig? Dann passiert sowas

Füg mal gleich hinter dem begin ein vtot := 0 ein.
Sollte Delphi dir auch als Hinweis angeben.

Nach oben

Wölfchen

Betreff des Beitrags: Re: Volumenberechnung eines beliebigen Mesh aus Dreiecksflächen.

Verfasst: Mo Jan 18, 2010 18:49

DGL Member

Registriert: Do Mär 05, 2009 20:17
Beiträge: 284
Wohnort: Kaiserslautern

zoiX hat geschrieben:

Ich rate mal ins Blaue: Du hast vtot nicht initialisiert, richtig? Dann passiert sowas

Füg mal gleich hinter dem begin ein vtot := 0 ein.
Sollte Delphi dir auch als Hinweis angeben.

Schade,
hatte echt gehofft das wärs - aber hat nix gebracht :cry:

ich muß wohl weiter suchen

**EDIT**

ich glaub ich habs! :mrgreen:

ich bin zu oft durch die Schleife gelaufen, und zwar basiert der code aus dem Blender Script ja darauf die Schleife nur für jede Fläche einmal zu durchlaufen, ich bin aber 3 mal so oft (nämlich für jeden Vertex) durchgehuscht...

das hier hats jetzt korrigiert und alles was ich bis jetzt getestet hab, schaut sehr gut aus:

Code:

anzflaechen:=round(length(TCW3d(instanzliste[0]).varray)/3);
for i:=0 to  anzflaechen-1 do begin

Nach oben

Seite 1 von 1

[ 7 Beiträge ]

Foren-Übersicht » Programmierung » Einsteiger-Fragen

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 3 Gäste

Du darfst keine neuen Themen in diesem Forum erstellen.
Du darfst keine Antworten zu Themen in diesem Forum erstellen.
Du darfst deine Beiträge in diesem Forum nicht ändern.
Du darfst deine Beiträge in diesem Forum nicht löschen.
Du darfst keine Dateianhänge in diesem Forum erstellen.