Page MenuHomec4science
Files F92320050

TwoDThreeD.pde
No OneTemporary
Actions

File Metadata

Created
Tue, Nov 19, 09:27

TwoDThreeD.pde
View Options

import java.util.List;
import processing.core.PVector;
import org.opencv.core.Mat;
import org.opencv.core.CvType;
import org.opencv.core.Core;
class TwoDThreeD {
// default focal length, well suited for most webcams
float f = 700;
// intrisic camera matrix
float [][] K = {{f, 0, 0},
{0, f, 0},
{0, 0, 1}};
float [][] invK;
PVector invK_r1, invK_r2, invK_r3;
Mat opencv_A, w, u, vt;
double [][] V;
// Real physical coordinates of the Lego board in mm
//float boardSize = 380.f; // large Duplo board
// float boardSize = 255.f; // smaller Lego board
// the 3D coordinates of the physical board corners, clockwise
float [][] physicalCorners = {
{-128, 128, 0, 1},
{128, 128, 0, 1},
{128, -128, 0, 1},
{-128, -128, 0, 1}
};
//Filtering variables: low-pass filter based on arFilterTrans from ARToolKit v5 */
float[] q;
float sampleRate;
float cutOffFreq;
float alpha;
public TwoDThreeD(int width, int height, float sampleRate) {
// set the offset to the center of the webcam image
K[0][2] = 0.5f * width;
K[1][2] = 0.5f * height;
//compute inverse of K
Mat opencv_K= new Mat(3, 3, CvType.CV_32F);
opencv_K.put(0, 0, K[0][0]);
opencv_K.put(0, 1, K[0][1]);
opencv_K.put(0, 2, K[0][2]);
opencv_K.put(1, 0, K[1][0]);
opencv_K.put(1, 1, K[1][1]);
opencv_K.put(1, 2, K[1][2]);
opencv_K.put(2, 0, K[2][0]);
opencv_K.put(2, 1, K[2][1]);
opencv_K.put(2, 2, K[2][2]);
Mat opencv_invK=opencv_K.inv();
invK = new float[][]{
{ (float)opencv_invK.get(0, 0)[0], (float)opencv_invK.get(0, 1)[0], (float)opencv_invK.get(0, 2)[0] },
{ (float)opencv_invK.get(1, 0)[0], (float)opencv_invK.get(1, 1)[0], (float)opencv_invK.get(1, 2)[0] },
{ (float)opencv_invK.get(2, 0)[0], (float)opencv_invK.get(2, 1)[0], (float)opencv_invK.get(2, 2)[0] }};
invK_r1=new PVector(invK[0][0], invK[0][1], invK[0][2]);
invK_r2=new PVector(invK[1][0], invK[1][1], invK[1][2]);
invK_r3=new PVector(invK[2][0], invK[2][1], invK[2][2]);
opencv_A=new Mat(12, 9, CvType.CV_32F);
w=new Mat();
u=new Mat();
vt=new Mat();
V= new double[9][9];
q=new float[4];
q[3]=1;
this.sampleRate=sampleRate;
if (sampleRate>0) {
cutOffFreq=sampleRate/2;
alpha= (1/sampleRate)/(1/sampleRate + 1/cutOffFreq);
}
}
PVector get3DRotations(List<PVector> points2D) {
// 1- Solve the extrinsic matrix from the projected 2D points
double[][] E = solveExtrinsicMatrix(points2D);
// 2 - Re-build a proper 3x3 rotation matrix from the camera's
// extrinsic matrix E
PVector firstColumn=new PVector((float)E[0][0], (float)E[1][0], (float)E[2][0]);
PVector secondColumn=new PVector((float)E[0][1], (float)E[1][1], (float)E[2][1]);
firstColumn.normalize();
secondColumn.normalize();
PVector thirdColumn=firstColumn.cross(secondColumn);
float [][] rotationMatrix={{firstColumn.x, secondColumn.x, thirdColumn.x},
{firstColumn.y, secondColumn.y, thirdColumn.y},
{firstColumn.z, secondColumn.z, thirdColumn.z}};
if (sampleRate>0)
filter(rotationMatrix, false);
// 3 - Computes and returns Euler angles (rx, ry, rz) from this matrix
return rotationFromMatrix(rotationMatrix);
}
double[][] solveExtrinsicMatrix(List<PVector> points2D) {
// p ~= K · [R|t] · P
// with P the (3D) corners of the physical board, p the (2D)
// projected points onto the webcam image, K the intrinsic
// matrix and R and t the rotation and translation we want to
// compute.
//
// => We want to solve: (K^(-1) · p) X ([R|t] · P) = 0
float[][] projectedCorners = new float[4][3];
if(points2D.size() >= 4)
for (int i=0; i<4; i++) {
// TODO:
// store in projectedCorners the result of (K^(-1) · p), for each
// corner p found in the webcam image.
// You can use PVector dot function for computing dot product between K^(-1) lines and p.
//Do not forget to normalize the result
PVector point =points2D.get(i);
projectedCorners[i][0]=point.dot(invK_r1)/point.dot(invK_r3);
projectedCorners[i][1]=point.dot(invK_r2)/point.dot(invK_r3);
projectedCorners[i][2]=1;
}
// 'A' contains the cross-product (K^(-1) · p) X P
float[][] A= new float[12][9];
for (int i=0; i<4; i++) {
A[i*3][0]=0;
A[i*3][1]=0;
A[i*3][2]=0;
// note that we take physicalCorners[0,1,*3*]: we drop the Z
// coordinate and use the 2D homogenous coordinates of the physical
// corners
A[i*3][3]=-projectedCorners[i][2] * physicalCorners[i][0];
A[i*3][4]=-projectedCorners[i][2] * physicalCorners[i][1];
A[i*3][5]=-projectedCorners[i][2] * physicalCorners[i][3];
A[i*3][6]= projectedCorners[i][1] * physicalCorners[i][0];
A[i*3][7]= projectedCorners[i][1] * physicalCorners[i][1];
A[i*3][8]= projectedCorners[i][1] * physicalCorners[i][3];
A[i*3+1][0]= projectedCorners[i][2] * physicalCorners[i][0];
A[i*3+1][1]= projectedCorners[i][2] * physicalCorners[i][1];
A[i*3+1][2]= projectedCorners[i][2] * physicalCorners[i][3];
A[i*3+1][3]=0;
A[i*3+1][4]=0;
A[i*3+1][5]=0;
A[i*3+1][6]=-projectedCorners[i][0] * physicalCorners[i][0];
A[i*3+1][7]=-projectedCorners[i][0] * physicalCorners[i][1];
A[i*3+1][8]=-projectedCorners[i][0] * physicalCorners[i][3];
A[i*3+2][0]=-projectedCorners[i][1] * physicalCorners[i][0];
A[i*3+2][1]=-projectedCorners[i][1] * physicalCorners[i][1];
A[i*3+2][2]=-projectedCorners[i][1] * physicalCorners[i][3];
A[i*3+2][3]= projectedCorners[i][0] * physicalCorners[i][0];
A[i*3+2][4]= projectedCorners[i][0] * physicalCorners[i][1];
A[i*3+2][5]= projectedCorners[i][0] * physicalCorners[i][3];
A[i*3+2][6]=0;
A[i*3+2][7]=0;
A[i*3+2][8]=0;
}
for (int i=0; i<12; i++)
for (int j=0; j<9; j++)
opencv_A.put(i, j, A[i][j]);
Core.SVDecomp(opencv_A, w, u, vt);
for (int i=0; i<9; i++)
for (int j=0; j<9; j++)
V[j][i]=vt.get(i, j)[0];
double[][] E = new double[3][3];
//E is the last column of V
for (int i=0; i<9; i++) {
E[i/3][i%3] = V[i][V.length-1] / V[8][V.length-1];
}
return E;
}
PVector rotationFromMatrix(float[][] mat) {
// Assuming rotation order is around x,y,z
PVector rot = new PVector();
if (mat[1][0] > 0.998) { // singularity at north pole
rot.z = 0;
float delta = (float) Math.atan2(mat[0][1], mat[0][2]);
rot.y = -(float) Math.PI/2;
rot.x = -rot.z + delta;
return rot;
}
if (mat[1][0] < -0.998) { // singularity at south pole
rot.z = 0;
float delta = (float) Math.atan2(mat[0][1], mat[0][2]);
rot.y = (float) Math.PI/2;
rot.x = rot.z + delta;
return rot;
}
rot.y =-(float)Math.asin(mat[2][0]);
rot.x = (float)Math.atan2(mat[2][1]/Math.cos(rot.y), mat[2][2]/Math.cos(rot.y));
rot.z = (float)Math.atan2(mat[1][0]/Math.cos(rot.y), mat[0][0]/Math.cos(rot.y));
return rot;
}
int filter(float m[][], boolean reset) {
float[] q= new float[4];
float alpha, oneminusalpha, omega, cosomega, sinomega, s0, s1;
mat2Quat(m, q);
if (nomalizeQuaternion(q)<0) return -1;
if (reset) {
this.q[0] = q[0];
this.q[1] = q[1];
this.q[2] = q[2];
this.q[3] = q[3];
} else {
alpha = this.alpha;
oneminusalpha = 1.0 - alpha;
// SLERP for orientation.
cosomega = q[0]*this.q[0] + q[1]*this.q[1] + q[2]*this.q[2] + q[3]*this.q[3]; // cos of angle between vectors.
if (cosomega < 0.0) {
cosomega = -cosomega;
q[0] = -q[0];
q[1] = -q[1];
q[2] = -q[2];
q[3] = -q[3];
}
if (cosomega > 0.9995) {
s0 = oneminusalpha;
s1 = alpha;
} else {
omega = acos(cosomega);
sinomega = sin(omega);
s0 = sin(oneminusalpha * omega) / sinomega;
s1 = sin(alpha * omega) / sinomega;
}
this.q[0] = q[0]*s1 + this.q[0]*s0;
this.q[1] = q[1]*s1 + this.q[1]*s0;
this.q[2] = q[2]*s1 + this.q[2]*s0;
this.q[3] = q[3]*s1 + this.q[3]*s0;
nomalizeQuaternion(this.q);
}
if (quat2Mat(this.q, m) < 0) return (-2);
return (0);
}
int nomalizeQuaternion(float[] q) {// Normalise quaternion.
float mag2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3];
if (mag2==0) return (-1);
float mag = sqrt(mag2);
q[0] /= mag;
q[1] /= mag;
q[2] /= mag;
q[3] /= mag;
return (0);
}
int mat2Quat(float m[][], float q[]) {
float t, s;
t = m[0][0] + m[1][1] + m[2][2] + 1.0;
if (t > 0.0001) {
s = sqrt(t) * 2.0;
q[0] = (m[1][2] - m[2][1]) / s;
q[1] = (m[2][0] - m[0][2]) / s;
q[2] = (m[0][1] - m[1][0]) / s;
q[3] = 0.25 * s;
} else {
if (m[0][0] > m[1][1] && m[0][0] > m[2][2]) { // Column 0:
s = sqrt(1.0 + m[0][0] - m[1][1] - m[2][2]) * 2.0;
q[0] = 0.25 * s;
q[1] = (m[0][1] + m[1][0] ) / s;
q[2] = (m[2][0] + m[0][2] ) / s;
q[3] = (m[1][2] - m[2][1] ) / s;
} else if (m[1][1] > m[2][2]) { // Column 1:
s = sqrt(1.0 + m[1][1] - m[0][0] - m[2][2]) * 2.0;
q[0] = (m[0][1] + m[1][0] ) / s;
q[1] = 0.25 * s;
q[2] = (m[1][2] + m[2][1] ) / s;
q[3] = (m[2][0] - m[0][2] ) / s;
} else { // Column 2:
s = sqrt(1.0 + m[2][2] - m[0][0] - m[1][1]) * 2.0;
q[0] = (m[2][0] + m[0][2] ) / s;
q[1] = (m[1][2] + m[2][1] ) / s;
q[2] = 0.25 * s;
q[3] = (m[0][1] - m[1][0] ) / s;
}
}
return 0;
}
int quat2Mat( float q[], float m[][] )
{
float x2, y2, z2;
float xx, xy, xz;
float yy, yz, zz;
float wx, wy, wz;
x2 = q[0] * 2.0;
y2 = q[1] * 2.0;
z2 = q[2] * 2.0;
xx = q[0] * x2;
xy = q[0] * y2;
xz = q[0] * z2;
yy = q[1] * y2;
yz = q[1] * z2;
zz = q[2] * z2;
wx = q[3] * x2;
wy = q[3] * y2;
wz = q[3] * z2;
m[0][0] = 1.0 - (yy + zz);
m[1][1] = 1.0 - (xx + zz);
m[2][2] = 1.0 - (xx + yy);
m[1][0] = xy - wz;
m[0][1] = xy + wz;
m[2][0] = xz + wy;
m[0][2] = xz - wy;
m[2][1] = yz - wx;
m[1][2] = yz + wx;
return 0;
}
}

Event Timeline

c4science · Help