Suche Matrix Libraries

[Degush]

Hey,

ich suche fuer ein aktuelles Projekt eine gute Bibliothek zum Rechnen mit Vektoren im dreidimensionalen Raum. Wichtig sind mir z.B. Methoden, die ein lineares Gleichungssystem loesen.
Schoen waeren auch Methoden zur Abstandsbestimmung von Gerade/Gerade, Punkt/Punkt, Gerade/Ebene usw.
Also quasi eine Bibliothek fuer lineare Algebra.
Ich koennte so erwas natuerlich selbst schreiben, aber das ist viel Arbeit und ich bin mir sicher, dass es schon solche Bibliotheken gibt, die mir diese Arbeiten ersparen.

MfG,
Schmalte

 

[Flown]

Vielleicht suchst du genau sowas?

JAMA

 

[kama]

Hi,

schon mal JScience angeschaut ?

JScience

Gruß
Karl-Heinz Marbaise

 

[Guybrush Threepwood]

matrix-toolkits-java - a comprehensive collection of matrix data structures, linear solvers, least squares methods, eigenvalue and singular value decompositions. - Google Project Hosting

 

[Marco13]

Es gibt mehrere Bibliotheken, die Teilbereiche davon (mehr oder weniger gut) Abbilden. Eine Lib ganz konkret für solche Abstandsberechnungen kenne ich nicht (ist schwieriger als man im ersten Moment meinen könnte...)

 

[Degush]

Danke fuer eure Hilfe.
Habe jetzt la4j verwendet. Die anderen Bibliotheken sind teilweise verdammt unuebersichtlich - mit 300 Klassen, von denen jede irgendeine komische Abkuerzung als Namen hat.

Die Abstandsbestimmung zweier windschiefer Geraden musste ich selbst implementieren.
Das habe ich mithilfe des Vektorprodukts und der Hesseschen Normalform gemacht:

    public static double getDistance(Vector pointA, Vector directionA, Vector pointB, Vector directionB)
    {
        double[] normal = new double[3];
        normal[0] = directionA.get(1)*directionB.get(2)-directionA.get(2)*directionB.get(1);
        normal[1] = directionA.get(2)*directionB.get(0)-directionA.get(0)*directionB.get(2);
        normal[2] = directionA.get(0)*directionB.get(1)-directionA.get(1)*directionB.get(0);
        return ( normal[0]*(pointA.get(0)-pointB.get(0)) + 
                normal[1]*(pointA.get(1)-pointB.get(1)) + 
                normal[2]*(pointA.get(2)-pointB.get(2)) ) / 
                (Math.sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]));
    }

 

[Degush]

Habe zu diesem Thema auch gleich noch eine Frage, die sich auf Lineare Algebra bezieht.
Ich will die Lotfusspunkte der beiden Geraden bestimmen.
Theoretisch sollte ich den Abstand der Geraden doch auch bekommen, wenn ich ein LGS loese mit folgenden Bedingungen:

Geraden sind g und h. In g wird der RIchtungsvektor mit u multipliziert, in h mit k.
Abstandsvektor d = g-h
d*g = 0
d*h = 0

Also eingesetzt
(g-h)*g = 0
(g-h)*h = 0

Und umgeformt

g-h liefert fuer beliebige Parameter eine Gerade, die durch g und h geht. An der Stelle, an der der Abstand am kleinsten ist, sollten g und h orthogonal sein (oder irre ich mich?), das Skalarprodukt also 0.
Wenn ich also das LGs loese, sollten nach meiner Ueberlegung die Lotfusspunkte dadurch bestimmbar sein, dass ich die Werte h und k, die im LGS bestimmt worden sind, wieder in die Parameterglg eingesetzt werden.

Mein Code:

    public static void solve(Vector pointA, Vector directionA, Vector pointB, Vector directionB)
    {
        double[][] matrix = new double[2][2];
        matrix[0][0] = directionB.get(0)*directionA.get(0)+directionB.get(1)*directionA.get(1)+directionB.get(2)*directionA.get(2);
        matrix[1][0] = -(directionA.get(0)*directionA.get(0)+directionA.get(1)*directionA.get(1)+directionA.get(2)*directionA.get(2));
        matrix[0][1] = directionB.get(0)*directionB.get(0)+directionB.get(1)*directionB.get(1)+directionB.get(2)*directionB.get(2);
        matrix[1][1] = -matrix[0][0];
        DenseFactory denseFactory = new DenseFactory();
        Matrix a = denseFactory.createMatrix(matrix);
        double result[] = new double[2];
        result[0] = -((pointB.get(0)-pointA.get(0))*directionA.get(0)+(pointB.get(1)-pointA.get(1))*directionA.get(1)+(pointB.get(2)-pointA.get(2))*directionA.get(2));
        result[1] = -((pointB.get(0)-pointA.get(0))*directionB.get(0)+(pointB.get(1)-pointA.get(1))*directionB.get(1)+(pointB.get(2)-pointA.get(2))*directionB.get(2));
        
        Vector b = denseFactory.createVector(result);
        System.out.println(a.toString());
        LinearSystem system = new LinearSystem(a, b);
        try {
            Vector res = system.solve(new GaussianSolver());
            System.out.println(res.toString());
            Vector lotA = pointA.add(directionA.multiply(res.get(0)));
            System.out.println(lotA.toString());
            Vector lotB = pointB.add(directionB.multiply(res.get(1)));
            System.out.println(lotB.toString());
            System.out.println(getMagnitude(lotA.subtract(lotB)));
        } catch (Exception ex) {
            Logger.getLogger(VectorCalculations.class.getName()).log(Level.SEVERE, null, ex);
        }
    }

Liefert allerdings falsche Ergebnisse.
Wenn mein Ansatz falsch ist: Wie soll ich das sonst machen?

MfG,
Degush

edit; Ansatz klappt. Habe aber die Lsg in der falschen Reihenfolge verwendet. Richtig ist

            Vector lotA = pointA.add(directionA.multiply(res.get(1)));
            System.out.println(lotA.toString());
            Vector lotB = pointB.add(directionB.multiply(res.get(0)));

anstatt

            Vector lotA = pointA.add(directionA.multiply(res.get(0)));
            System.out.println(lotA.toString());
            Vector lotB = pointB.add(directionB.multiply(res.get(1)));

 

[Semox]

Hi Degush/Schmalte

Es wäre hilfreich von Dir, wenn Du den Anwendungsfall beschriebst, damit man weiß ob "mal" was gerechnet werden muß, oder ob der Aufruf millionenfach passiert.

Das könnte etwas für Dich sein, wenn Du auch gern mal C++ für Performancegewinn einsetzen möchtest:

JNA - https://github.com/twall/jna --> zum Binden von Java an C++ Libs mit guter Doku und Beispielen
blitz++ --> superhigh-performance C++ Lib für Vektorrechnung und kleine matrizen ist das nice

Für größere Matrizen mit bereits haufenweise fertigen Funktionen:
C++ Lib Boost mit uBlas als Untermenge.

Auf jeden Fall würde ich wegen Speed davon abraten eigene Implementationen für Matrizenberechnungen zu schaffen, es sei denn Du kannst ubercrack-mäßig Algorithmen entwickeln und hast einen Doktor in Mathematik. ^^ Ich habe leider noch keinen.

Ansonsten ist diese Vektorengeschichte auch mit etwas OpenGL bzw. GLSL Hacks echt prima lösbar. Du kannst dann damit zur Matrizenrechnung die blitzflinke Grafikkarte benutzen. Mußt halt nichts an den GL Context schicken und zeichnen lassen, sondern nur die Matrizen retournieren und auslesen. Insgesamt ist das dann aber auch vom Aufwand eher für komplexe numerische Berechungen lohnend.

Hoffe Dir geholfen zu haben.

Viele Grüße,
Semox

 

[Degush]

Um eine bessere Performance werde ich mich kümmern, sobald das Projekt steht und ich mit dem, was ich habe, arbeiten kann.
Die Schnittpunkte der Geraden werden bestimmt, um aus zwei Bildern einer Kamera die Tiefe für jeden Pixel zu bestimmen.

 

[pappawinni]

Oh, für die Berechnung des Abstandes zweier Geraden braucht es doch keine Solver oder irgendwelche Libs.
Das geht doch gut zu Fuß.
Angenommen du hast zwei Geraden in der Punkt-Richtungsform
g: (a,b,c) + lamda * ( e, f,g ) = (X,Y,Z)
h: (h,j,k) + ny * (m,n,o) = (X,Y,Z)
Dann kannst du aus den Richtungsvektoren (e,f,g) und (m, n, o) das Kreuzprodukt bilden und erhältst damit einen zu beiden Richtungsvektoren senkrechten Vektor (q , r, s ), d.h.
q = (f *o)-(g*n)
r = (g*m)-(e*o)
s= (e*n)-(f*m)
Diesen Vektor könntest du dann als Normalenvektor von 2 Ebenen auffassen mit jeweils den Punkten der Geraden.
Dazu normalisierst du den Vektor erst:
z = Wurzel(q^2+r^2+s^2)
q = q / z
r = r / z
s = s / z
und würdest die folgenden Ebenen erhalten:
E1: q * X + r * Y + s * Z = q * a + r * b + s * c
E2: q * X + r * Y + s * Z = q * h + r * j + s * k
Diese Ebenen hätten die Abstände
d1 = q * a + r * b + s * c
d2 = q * h + r * j + s * k
vom Koordinaten-Ursprung
und die Differenz dieser Abstände sagt dir, welchen Abstand sie untereinander haben.

Voraussetzung hier aber, dass die Geraden nicht parallel sind.

 

[pappawinni]

Ach für die Lotpunkte reichte das Alphabet nicht aus, aber das Prinzip wäre:
Du bildest mit dem Normalenvektor der Ebenen und dem Richtungsvektor einer Geraden wiederum das Kreuzprodukt, hast damit einen weiteren Normalenvektor einer Ebene, die jetzt den den Punkt der gleichen Geraden enthalten soll.
Diese Ebene schneidest du mit der anderen Geraden und hast den ersten Lotpunkt.
Das Gleiche mit der anderen Geraden und du hast den zweiten Lotpunkt, fertig.

 

[pappawinni]

Ich hab jetzt auch mal ne Klasse für Geraden gemacht.
Damit könntest du evtl. deine Lotpunkte bestimmen und auch ein paar andere Dinge:

public class GeradenDemo {

	/**
	 * @param args
	 */
	public static void main(String[] args) {		
		LineGO line1 = new LineGO(new double[] {2.0000,-2.000,0.0000},new double[] {0.0000,0.9487,2.8460});
		LineGO line2 = new LineGO(new double[] {2.0000,2.0000,0.0000},new double[] {-0.6030,-0.6030,1.8091});
		System.out.println("Gerade 1: " + line1);
		System.out.println("Gerade 2: " + line2);		
		System.out.println("Winkel zwischen Gerade 1 und 2 im Bogenmaß: " + line1.getAngleRAD(line2));
		System.out.println("Winkel zwischen Gerade 1 und 2 im Dezimalgrad: " + line1.getAngleDEG(line2));
		System.out.println("Abstand zwischen Geraden 1 und 2 : " + line1.getDistance(line2));
		double[] punktA = line1.getNearestPoint(line2);
		double[] punktB = line2.getNearestPoint(punktA);
		System.out.println("Punkt A = Nähester Punkt auf Gerade 1 zur Geraden 2 ("+punktA[0]+","+punktA[1]+","+punktA[2]+")");
		System.out.println("Punkt B = Lotpunkt von Punkt A auf Gerade 2 ("+punktB[0]+","+punktB[1]+","+punktB[2]+")");
		double[] punktC = line2.getNearestPoint(line1);
		System.out.println("Punkt C = Nähester Punkt von Gerade 2 zur Geraden 1 ("+punktC[0]+","+punktC[1]+","+punktC[2]+")");
	}

}

class LineGO{
	public double orientation_X;
	public double orientation_Y;
	public double orientation_Z;
	public double supportpoint_X;
	public double supportpoint_Y;
	public double supportpoint_Z;	

	LineGO(){
		orientation_X=1;
	}
	
	LineGO(double[] supportpoint,double[] orientation){
		orientation_X=orientation[0];
		orientation_Y=orientation[1];
		orientation_Z=orientation[2];
		supportpoint_X=supportpoint[0];
		supportpoint_Y=supportpoint[1];
		supportpoint_Z=supportpoint[2];
	}
	public double getAngleRAD(LineGO line){
		double d=Math.sqrt(
				(orientation_X*orientation_X
				+orientation_Y*orientation_Y
				+orientation_Z*orientation_Z) *
				(line.orientation_X*line.orientation_X
				+line.orientation_Y*line.orientation_Y
				+line.orientation_Z*line.orientation_Z) );
		double angle = Math.acos(Math.abs(
				(orientation_X*line.orientation_X
				+orientation_Y*line.orientation_Y
			    +orientation_Z*line.orientation_Z)/d 
				));
		return angle;		
	}
	public double getAngleDEG(LineGO line){
		return(Math.toDegrees(getAngleRAD(line)));		
	}
	public double getDistance(double[] point){
		double len = Math.sqrt(orientation_X*orientation_X
				    +orientation_Y*orientation_Y
				    +orientation_Z*orientation_Z);
		double dist = 0;
		if (len > 0) {
			double dx = orientation_Z*(point[1]-supportpoint_Y)-orientation_Y*(point[2]-supportpoint_Z);
			double dy = orientation_X*(point[2]-supportpoint_Z)-orientation_Z*(point[0]-supportpoint_X);
			double dz = orientation_Y*(point[0]-supportpoint_X)-orientation_X*(point[1]-supportpoint_Y);
			dist = Math.sqrt(dx*dx+dy*dy+dz*dz)/len;
		}
		return (dist);
	}
	public double getDistance(LineGO line){
		double[] normal ={
            				orientation_Y*line.orientation_Z-orientation_Z*line.orientation_Y,
            				orientation_Z*line.orientation_X-orientation_X*line.orientation_Z,
            				orientation_X*line.orientation_Y-orientation_Y*line.orientation_X				
		                 };
	    double len = Math.sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
		double dist = 0;
		if (len > 0) {
		   dist=Math.abs(
				(supportpoint_X-line.supportpoint_X)*normal[0]
			   +(supportpoint_Y-line.supportpoint_Y)*normal[1]
			   +(supportpoint_Z-line.supportpoint_Z)*normal[2]
			    )/len;
		}
		else{
		  double[] point ={line.supportpoint_X,line.supportpoint_Y,line.supportpoint_Z};
		  dist = getDistance(point);				
		}
		return (dist);
	}
	public double[] getNearestPoint(double[] point){
		double pnt[] = {point[0],point[1],point[2]};
		double len = Math.sqrt(orientation_X*orientation_X
			    +orientation_Y*orientation_Y
			    +orientation_Z*orientation_Z);
	    if (len > 0) {
    		double dx = orientation_Z * (point[1]-supportpoint_Y) - orientation_Y * (point[2]-supportpoint_Z);
	    	double dy = orientation_X * (point[2]-supportpoint_Z) - orientation_Z * (point[0]-supportpoint_X);
		    double dz = orientation_Y * (point[0]-supportpoint_X) - orientation_X * (point[1]-supportpoint_Y);
		    pnt[0] -= (orientation_Y * dz - orientation_Z * dy) / len / len;
		    pnt[1] -= (orientation_Z * dx - orientation_X * dz) / len / len;
		    pnt[2] -= (orientation_X * dy - orientation_Y * dx) / len / len;				
	    }
	    return (pnt);
	}
	public double[] getNearestPoint(LineGO line){
		double s1 = 
				orientation_X * line.orientation_X +
				orientation_Y * line.orientation_Y +
				orientation_Z * line.orientation_Z;
		double s2 = line.orientation_X * line.orientation_X
				  + line.orientation_Y * line.orientation_Y
				  + line.orientation_Z * line.orientation_Z;
		double[] vec = {s1 * line.orientation_X - s2 * orientation_X,
				        s1 * line.orientation_Y - s2 * orientation_Y,
				        s1 * line.orientation_Z - s2 * orientation_Z
				        };
		double div = orientation_X * vec[0] + orientation_Y * vec[1] + orientation_Z * vec[2];
		double k = vec[0] * (supportpoint_X - line.supportpoint_X)
				  + vec[1] * (supportpoint_Y - line.supportpoint_Y)
				  + vec[2] * (supportpoint_Z - line.supportpoint_Z);
        if (div != 0){
        	k /= div;
        }
        double[] pnt = {supportpoint_X - k * orientation_X,
        		        supportpoint_Y - k * orientation_Y,
        		        supportpoint_Z - k * orientation_Z};
        return (pnt);		
	}	
	public String toString(){
		String str = "(" + supportpoint_X + "," + supportpoint_Y + "," + supportpoint_Z + ") + k * ("
				+ orientation_X + "," + orientation_Y + "," + orientation_Z + ") = (X,Y,Z)";
		return (str);
	}
}