Introduction

Notations

Vector, matrix and scalar

 * vector (lowercase Gothic alphabets): v, p, m = [m1, m2, ..., mk]', ...
 * matrix (uppercase Gothic alphabets): R, A, ...
 * scalar (Italic alphabets: s, t, X, Y, Z, ...

Indices

   v(1)   // the 1st element of the vector v
   R(1,2) // the 1st row and 2nd col of the matrix R
   A(:,1) // the 1st vector _ from A = [v1, v2, ..., vN]

Some reserved symbols

hp: a homogeneous representation of a 2D point, hp = [u, v, 1]'
hx: a homogeneous representation of a 3D point, hx = [X, Y, Z, 1]'
p: a 2D point, p = [u, v]'
x: a 3D point, x = [X, Y, Z]'
H: (3x3) planar homography matrix
P: (3x4) projective projection matrix

Coordinate transformations

o<W>: the origin of the world frame
o<C>: the origin of the camera frame
o<O>: the origin of the object frame
(e.g., x<C> represents a 3D point in the camera frame)
R: (3x3) rotation matrix
t: (3x1) translation vector
T<n, m>: transformation from m to n coordinate system (e.g., o<C> = T<C, W> o<W>)

• We use the (first) camera coordinate frame as the world reference frame unless we specify it.
• Besides above notations, we will follow MATLAB's notations.

Coordinate transformation

When you take a 3D point, ${\displaystyle X=[x_{1},x_{2},x_{3}]^{T}}$, each value x(j) of the point is represented with respect to a specific coordinate system (e.g., the sensor coordinate system). The point X can also be represented in a different coordinate system (e.g., a corner of the room), ${\displaystyle X'=[x'_{1},x'_{2},x'_{3}]^{T}}$ . To distinguish between two set of values, we can write the symbol along with the corresponding coordinate system: ${\displaystyle X^{s}=[x_{1},x_{2},x_{3}]^{T}}$ and ${\displaystyle X^{c}=[x'_{1},x'_{2},x'_{3}]^{T}}$, where ${\displaystyle X^{s}}$ is a point w.r.t. the sensor coordinate system and ${\displaystyle X^{c}}$ is the same point w.r.t. the corner coordinate system. There is a rigid motion between two coordinate systems and the motion is expressed by a (3×3) rotation matrix ${\displaystyle R}$ and a (3×1) translation vector ${\displaystyle t}$: ${\displaystyle X^{s}=RX^{c}+t}$. In homogeneous coordinate system, the rigid motion can be written as ${\displaystyle X^{s}={^{s}}T_{c}X^{c},{^{s}}T_{c}=[{^{s}}R_{c},{^{s}}t_{c};0,0,0,1]}$. Note that we put the symbol of the source coordinate system (corner) and the target system (sensor). A series of transformations is represented as ${\displaystyle X^{s}={^{s}}T_{s_{3}}{^{s_{3}}}T_{s_{2}}{^{s_{2}}}T_{s_{1}}c{^{s_{1}}}T_{c}X^{c}}$, where ${\displaystyle {^{s}}T_{c}={^{s}}T_{s_{3}}{^{s_{3}}}T_{s_{2}}{^{s_{2}}}T_{s_{1}}c{^{s_{1}}}T_{c}}$.

Calibration

Why calibration?

• Precision of the Kinect sensor

Overall process:

 a) RGB camera ---- RGB image --(1)-- Rectified RGB image --(3)-- Registered RGB image (w.r.t. IR image)

 b) IR camera ---- IR image --(2)--  Rectified IR image

 c) projector --(internal calibration + IR image)-- Depth image --(4)-- Rectified Depth image --(5)-- 3D points                         

• Prob. 1: (1), (2), (4) - radial distortion
• Method:

http://www.ros.org/wiki/kinect_calibration/technical

http://nicolas.burrus.name/index.php/Research/KinectCalibration

• Prob. 2: (3) - registration of RGB image to IR image
• Method: planar homography

• Prob. 3: (5) - 3D reconstruction
• Method: triangulation

• All in one

http://www.ee.oulu.fi/~dherrera/kinect/

Geometric primitives and uncertainty

Uncertainty and geometric constraint

Uncertainty of geometric entity

• Uncertainty of a point?
• Uncertainty of a line?
• Uncertainty of a direction?
• Uncertainty of a plane?

An example: triangle

3 angles of triangle

Imagine that you are measuring three angles of a triangle with a protractor. From the first trial, you got values [30.5, 60.0, 90.1] degrees. Which are not perfect because the summation of them is not equal to 180! So, you have decided to measure angles 10 times and average them. However, still you are not sure that the sum of the averages will be 180.

Figure XXX shows the relations between three angles in a triangle.

Method

1) Define a manifold M (e.g. Figure XXX)
2) Measure a set of values p = [a(i), b(i), c(i)] for N times
3) Compute the mean m and the covariance matrix S
4) Project p and S onto the manifold M: q, H

output:
  angles <-- q
  variance <-- f(H)

Averaging multiple observations

When you observed multiple rigid transformations for the extrinsic calibration of N sensors, you can use "Dual Quaternions" for to compute an "average" rigid transformation.

The Algorithm

* input: {R(i), t(i)} % a set of R, t
* output: R*, t*      % optimal R, t
Procedure:
* convert all R(i), t(i) to DQ(i) % Dual Quaternions
* find the mean DQ* from {DQ(i)}
* convert DQ* to R*, t*

Dual Quaternions

"Similar to the way that rotations in 3D space can be represented by quaternions of unit length, rigid motions in 3D space can be represented by dual quaternions of unit length. This fact is used in theoretical kinematics, and in applications to 3D computer graphics, robotics and computer vision." Wiki

L. Kavan, S. Collins, C. O'Sullivan, J. Zara (2006) Dual Quaternions for Rigid Transformation Blending, Technical report, Trinity College Dublin.

Matlab toolbox with a nice document

Working on the manifold

http://stat.fsu.edu/~anuj/CVPR_Tutorial/ShortCourse.htm Differential Geometric Methods for Shape Analysis and Activity Recognition

Registration of multiple views

* input: X, Y
* output: R, t

Solution = [];
For i=1:N
  randomly select 5 corresponding points from X, Y: s1 = {X(j), X(k), X(l)} and s2 = {Y(j), Y(k), Y(l)}
  estimate T' = (R', t') using s1, s2
  compute Y^ = T' X
  store the number of inliers, m(i), to Solution: S(i) = (R'(i), t'(i), m(i))
Endfor
Answer = arg max(i) (S(i).m(i))

Multiple Sensor Systems

Applications

Related topics

Appendices

openni::VideoFrameRef depthFrame;
const openni::DepthPixel* pDepth = (const openni::DepthPixel*) depthFrame.getData;
int width = depthFrame.getWidth();
int height = depthFrame.getHeight();
for (int i = 0; i < height; ++i)
    for (int j = 0; j < width; ++j, ++pDepth)
    if(i == height/2 && j == width/2 && *pDepth != 0 )//The if sentence depends on which point you want, this is the center point for example
    {} //assign the value to a 3D point structure

/* Display depth map */
// Matthias Hernandez: 07/27/2011 - University of Southern California
// m31hernandez@gmail.com
// display the depthMap from openNI input in the openCV image ‘out’
void displayDepthMap(IplImage *out, const XnDepthPixel* pDepthMap, XnDepthPixel max_depth, XnDepthPixel min_depth) {
    uchar *data = (uchar *)out->imageData;
    int step = out->widthStep;
    int channels = out->nChannels;
    float normalize = (float)(max_depth-min_depth)/255.0f;

    int index=0;
    for (int i=0; i<MAX_I; i++) {
   	 if (pDepthMap[i] < min_depth || pDepthMap[i] > max_depth) {
   		 for (int k=0; k<channels; k++)
   			 data[index++] = 0;
   	 } else {
   		 if (normalize != 0) {
   			 for (int k=0; k<channels; k++)
   				 data[index++] = (int)(255-(float)(pDepthMap[i]-min_depth)/normalize);
   		 }
   		 else
   			 for (int k=0; k<channels; k++)
   				 data[index++] = 255;
   	 }    
    }
    
}

/* Display image map */
// Matthias Hernandez: 07/27/2011 - University of Southern California
// m31hernandez@gmail.com
// display the image map from openNI input in the openCV image ‘out’
void displayImageMap(IplImage *out, const XnUInt8* pImageMap) {
    uchar *data = (uchar *)out->imageData;
    int step = out->widthStep;
    int channels = out->nChannels;    

    for (int i=0; i<MAX3_HR; i+=3) {
   	 data[i+2]    = pImageMap[i];
   	 data[i+1]    = pImageMap[i+1];
   	 data[i]   	 = pImageMap[i+2];
    }
}

#define XtoZ 1.114880018171494f
#define YtoZ 0.836160013628620f
#define MIN_DEPTH 450
#define MAX_DEPTH 800

void convertP2RW(float *pDepth, float *pReal, int x, int y, int w, int h) {
        int max_i = w * h;

	int i1 = (y * w + x),
		i2 = i1 + max_i,
		i3 = i2 + max_i;

	float Z = pDepth[i1];

    if (Z > MIN_DEPTH && Z < MAX_DEPTH){
        float X_rw = ( (float)x /(float)w -0.5f)*Z*XtoZ;
        float Y_rw = (0.5f-(float)y / (float)h)*Z*YtoZ;

		pReal[i1] = X_rw;
		pReal[i2] = Y_rw;
		pReal[i3] = Z;
	} else {
		pReal[i1] = 0.0f;
		pReal[i2] = 0.0f;
		pReal[i3] = 0.0f;
	}
}

// Of Course, the inverse function is 

void convertRW2P(float *pReal, float *pDepth,int x, int y, int w, int h) {
	 int max_i = w * h;

	int i1 = (y * w + x),
	i2 = i1 + max_i,
	i3 = i2 + max_i;

	float xR = pReal[i1];
	float yR = pReal[i2];
	float zR = pReal[i3];

	int ixi = (xR/zR/XtoZ + 0.5f)*(float)w + 0.5f;  // x
	int iyi = (-yR/zR/YtoZ + 0.5f)*(float)h + 0.5f; // y

	pDepth[i1] = zR;
}

  x = P (H inv(H)) X
  H: a projective transformation

 Q = H Q H’
 PAI is the null-vector of ADQ (Q a = 0)

References

• [Queirolo 2010] C. C. Queirolo, Luciano Silva, Olga R. P. Bellon, M. P. Segund, 3D Face Recognition using the Surface Interpenetration Measure: a Comparative Evaluation on the FRGC database, 3D Face Recognition Using Simulated Annealing and the Surface Interpenetration Measure. IEEE Trans. Pattern Anal. Mach. Intell. 32(2): 206-219 (2010).
• [CVonline] CVonline: The Evolving, Distributed, Non-Proprietary, On-Line Compendium of Computer Vision
• [Foerstner04] Uncertainty and Projective Geometry
• Computational Differential Geometry dmg.tuwien.ac.at
• Computational Geometry Algorithms Library

Depth sensors

PrimeSense www.primesense.com

Kinect - Xbox.com www.xbox.com/KINECT

Xtion PRO LIVE http://www.asus.com/Multimedia/Motion_Sensor/Xtion_PRO_LIVE/

softkinetic http://www.softkinetic.com/

velodynelidar.com http://velodynelidar.com/lidar/lidar.aspx

Introducing the Leap http://www.youtube.com/watch?v=_d6KuiuteIA

Mesa Imaging http://mesa-imaging.ch/

Tools

http://www.danielgm.net/cc/ CloudCompare - 3D point cloud and mesh processing software

http://www.cs.unc.edu/~isenburg/lastools/ LAStools: award-winning software for rapid converting, filtering, viewing, gridding, and compressing of LiDAR

Q&A

• {Put your questions here...}