public static BoundingBox calculateOrientedBoundingBox(Object3D object) 
	{
		// calculate mean (average) position
		float avgX = 0;
		float avgY = 0;
		float avgZ = 0;

		final float[] vertices = object.getVertices();
		final int N = (int) (vertices.length / 4);
		for ( int i = 0 ; i < vertices.length ; ) {
			avgX += vertices[i++];
			avgY += vertices[i++];
			avgZ += vertices[i++];
			i++; // ignore W
		}

		avgX /= N;
		avgY /= N;
		avgZ /= N;

		// construct 3x3 covariance matrix
		float c11 = 0;
		float c12 = 0;
		float c22 = 0;
		float c13 = 0;
		float c33 = 0;
		float c23 = 0;
		for ( int i = 0 ; i < vertices.length ; i+=4 ) 
		{
			final float x = vertices[i];
			final float y = vertices[i+1];
			final float z = vertices[i+2];

			c11 += ( x - avgX ) * ( x - avgX );
			c12 += ( x - avgX ) * ( y - avgY );
			c22 += ( y - avgY ) * ( y - avgY );
			c13 += ( x - avgX ) * ( z - avgZ );
			c33 += ( z - avgZ ) * ( z - avgZ );
			c23 += ( y - avgY ) * ( z - avgZ );
		}
		c11 /= N;
		c12 /= N;
		final float c21 = c12;
		c22 /= N;
		c13 /= N;
		final float c31 = c13;
		c33 /= N;
		c23 /= N;
		final float c32 = c23;

		final Matrix c = new Matrix( new Vector4(c11,c12,c13),
				new Vector4(c21,c22,c23),
				new Vector4(c31,c32,c33),
				new Vector4(0,0,0,0)
				);

		/* The eigenvalues (c11, c22 , c33) are the roots of the characteristic polynomial:
		 *  
		 *                           | c11 - lambda  c21         c31         |                  
		 * det( C - lambda * I ) =   | c12           c22-lambda  c32         |
		 *                           | c13           c23         c32-lambda  |
		 * 
		 */

		// calculate eigenvectors

		final Matrix eigenVectors = new Matrix();
		calculateEigenSystem( c , eigenVectors );

		Vector4 r = eigenVectors.getRow(0); // .multiply( -1 );
		Vector4 s = eigenVectors.getRow(1);
		Vector4 t = eigenVectors.getRow(2); // .multiply( -1 );

		float minPR = 0;
		float maxPR = 0;

		float minPS = 0;
		float maxPS = 0;

		float minPT = 0;
		float maxPT = 0;

		for ( Iterator<Vector4> iterator = object.getVertexIterator(); iterator.hasNext() ; ) {

			final Vector4 v = iterator.next();

			/*
			 * If b is a unit vector, then the dot product a*b gives |a| * cos(theta), i.e., 
			 * the magnitude of the projection of a in the direction of b, with a minus sign if the direction is opposite
			 */
			float dot1 = v.dotProduct( r );
			float dot2 = v.dotProduct( s );
			float dot3 = v.dotProduct( t );

			if ( dot1 < minPR ) {
				minPR = dot1;
			}
			if ( dot2 < minPS ) {
				minPS = dot2;
			}
			if ( dot3 < minPT ) {
				minPT = dot3;
			}		

			if ( dot1 > maxPR ) {
				maxPR = dot1;
			}	
			if ( dot2 > maxPS ) {
				maxPS = dot2;
			}	
			if ( dot3 > maxPT ) {
				maxPT = dot3;
			}				
		}

		// calculate center of bounding box
		final float extendR = (minPR + maxPR) / 2;
		final float extendS = (minPS + maxPS) / 2;
		final float extendT = (minPT + maxPT) / 2;

		final Vector4 center = r.multiply( extendR ).plus( s.multiply( extendS ) ).plus( t.multiply( extendT ) );

		System.out.println("Center: "+center);

		System.out.println("R = "+r);
		System.out.println("S = "+s);
		System.out.println("T = "+t);

		final float width = maxPR - minPR;
		final float height = maxPS - minPS;
		final float depth = maxPT - minPT;

		System.out.println("Extent R: "+minPR+" , "+maxPR+" (width "+width+")");
		System.out.println("Extent S: "+minPS+" , "+maxPS+" (height "+height+")");
		System.out.println("Extent T: "+minPT+" , "+maxPT+" (depth "+depth+")");

		return new BoundingBox( center , r , s , t , width , height , depth , false );
	}