code

H_from_3xP.m

@@ -0,0 +1,71 @@
+function [H, X] = H_from_3xP(c, P, CP)
+
+% get the homography matrix from two matched triangles, camera matrixes and 
+% the motion between cameras
+% c  -   the matched triangles
+% P  -   the camera matrixes
+% CP -   the first camera's position
+% H  -   the homography matrix (output)
+% see-   "Image registration for image-based rendering."IEEE TIP 2005.
+
+I = eye(3);
+
+%{
+if all(all(c(:,1,:) - c(:,2,:)))
+    H = I;
+    return;
+end
+%}
+
+% get the 3D point of the plane in world coordinate
+X1 = X_from_xP(c(:,:,1), P);
+X2 = X_from_xP(c(:,:,2), P);
+X3 = X_from_xP(c(:,:,3), P);
+
+% the parameters of the plane
+PI = [cross(X1(1:3) - X3(1:3), X2(1:3) - X3(1:3)); ...
+    -X3(1:3)'*cross(X1(1:3), X2(1:3))];
+PI = PI./PI(4);
+
+% projective matrix
+N = inv(P(1:3, 1:3, 1));
+
+% H1 can transform the points in image 1 to the 3D points.
+H1 = [( I - ( CP * PI(1:3)' ) ./ ( PI(1:3)' * CP + 1) ) * N;...
+    -(PI(1:3)' * N) ./ ( PI(1:3)' * CP + 1)]; 
+
+% H2 can transform the 3D points to the points in image 2.
+H2 = P(:,:,2);
+
+% the final homography matrix
+H = H2 * H1;
+
+% the 3D triangle
+X = [X1(1:3), X2(1:3), X3(1:3)];
+
+% affine matrix estimation
+x1 = H*[c(:,1,1);1]; x1 = x1(1:2)/x1(3);
+x2 = H*[c(:,1,2);1]; x2 = x2(1:2)/x2(3);
+x3 = H*[c(:,1,3);1]; x3 = x3(1:2)/x3(3);
+
+xp1 = c(:,2,1);
+xp2 = c(:,2,2);
+xp3 = c(:,2,3);
+
+A = [0 0 0 -x1(1) -x1(2) -1
+    x1(1) x1(2) 1 0 0 0
+    0 0 0 -x2(1) -x2(2) -1
+    x2(1) x2(2) 1 0 0 0
+    0 0 0 -x3(1) -x3(2) -1
+    x3(1) x3(2) 1 0 0 0];
+
+b = [-xp1(2); xp1(1); -xp2(2); xp2(1); -xp3(2); xp3(1)];
+affine = inv(A)*b;
+affine = [affine(1) affine(2) affine(3);
+    affine(4) affine(5) affine(6);
+    0 0 1];
+H = affine*H;
+
+
+
+

PointInSphericalTriangle2.m

@@ -0,0 +1,27 @@
+function result = PointInSphericalTriangle2(tri, pt, lHand)
+
+% test whether a point on unit sphere falls on a spherical triangle
+% tri:      3*3 matrix, each column represents a vertex of the spherical triangle
+% pt:       the point on the unit sphere to test
+% lHand:    weather the coordinate system is left hand
+% see       http://www.cs.brown.edu/~scd/facts.html
+
+% test whether the point is on the inside of all three triangle edges
+% be care of the orientation of the triangle (clockwise)
+
+if lHand
+    t1 = det([pt tri(:,1) tri(:,2)]);
+    t2 = det([pt tri(:,2) tri(:,3)]);
+    t3 = det([pt tri(:,3) tri(:,1)]);
+else
+    t1 = det([pt tri(:,2) tri(:,1)]);
+    t2 = det([pt tri(:,3) tri(:,2)]);
+    t3 = det([pt tri(:,1) tri(:,3)]);
+end
+
+% use 'eps' instead of 0 in case of floating point relative acuracy
+if(t1>eps && t2>eps && t3>eps)
+    result = 1;
+else
+    result = 0;
+end

SphTriOverlap.m

@@ -0,0 +1,82 @@
+function result = SphTriOverlap(tri1, tri2)
+
+% test whether there are overlap regions between two spherical triangles
+% tri1: 3*3 matrix, each column represents a vertex of the triangle
+% tri2: the second spherical triangle
+
+
+%% intersection test
+result = greatCircleArcIntersect(tri1(:,1), tri1(:,2), tri2(:,1), tri2(:,2));
+if result == 1
+    return;
+end
+
+result = greatCircleArcIntersect(tri1(:,1), tri1(:,2), tri2(:,2), tri2(:,3));
+if result == 1
+    return;
+end
+
+result = greatCircleArcIntersect(tri1(:,1), tri1(:,2), tri2(:,3), tri2(:,1));
+if result == 1
+    return;
+end
+
+result = greatCircleArcIntersect(tri1(:,2), tri1(:,3), tri2(:,1), tri2(:,2));
+if result == 1
+    return;
+end
+
+result = greatCircleArcIntersect(tri1(:,2), tri1(:,3), tri2(:,2), tri2(:,3));
+if result == 1
+    return;
+end
+
+result = greatCircleArcIntersect(tri1(:,2), tri1(:,3), tri2(:,3), tri2(:,1));
+if result == 1
+    return;
+end
+result = greatCircleArcIntersect(tri1(:,3), tri1(:,1), tri2(:,1), tri2(:,2));
+if result == 1
+    return;
+end
+
+result = greatCircleArcIntersect(tri1(:,3), tri1(:,1), tri2(:,2), tri2(:,3));
+if result == 1
+    return;
+end
+
+result = greatCircleArcIntersect(tri1(:,3), tri1(:,1), tri2(:,3), tri2(:,1));
+if result == 1
+    return;
+end
+
+%% inside test
+result = PointInSphericalTriangle2(tri1, tri2(:,1), 1);
+if result == 1
+    return;
+end
+
+result = PointInSphericalTriangle2(tri1, tri2(:,2), 1);
+if result == 1
+    return;
+end
+
+result = PointInSphericalTriangle2(tri1, tri2(:,3), 1);
+if result == 1
+    return;
+end
+
+result = PointInSphericalTriangle2(tri2, tri1(:,1), 1);
+if result == 1
+    return;
+end
+
+result = PointInSphericalTriangle2(tri2, tri1(:,2), 1);
+if result == 1
+    return;
+end
+
+result = PointInSphericalTriangle2(tri2, tri1(:,3), 1);
+if result == 1
+    return;
+end

SphTriTopoConsist.m

@@ -0,0 +1,15 @@
+function bConsist = SphTriTopoConsist(TRIS, spherePts, sample, srcTri)
+bConsist = 1;
+
+if TriOriReverse(srcTri, sample)
+    bConsist = 0;
+    return;
+end
+
+for i=1:size(TRIS,1)
+    tri = spherePts(:,TRIS(i,:));
+    if SphTriOverlap(tri, sample)
+        bConsist = 0;
+        return;
+    end
+end

TriOriReverse.m

@@ -0,0 +1,30 @@
+function bReverse = TriOriReverse(tri1,tri2)
+
+% test whether two triangles are orientation reversed
+% tri1: n*3 matrix, each column represents a vertex of the triangle
+% tri2: the second spherical triangle
+
+if ~all(size(tri1)==size(tri2))
+    error('Triangle must have the same dimension!!');
+end
+
+[rows, npts] = size(tri1);
+
+if npts ~=3
+    error('The number of the vertexs must be 3!!');
+end
+
+ if rows~=2 && rows~=3
+     error('The vertex must be 2D or 3D!!');
+ end
+
+if rows == 2
+    tri1 = [tri1; ones(1, 3)];
+    tri2 = [tri2; ones(1 ,3)];
+end
+
+
+det1 = det(tri1);
+det2 = det(tri2);
+
+bReverse = (det1 * det2) <= eps;

X_from_xP.m

@@ -0,0 +1,31 @@
+function X = X_from_xP(c,P)
+
+% Estimation of 3D point from image matches and camera matrices, linear.
+% c -   the correspondence in 2D or 3D (homogeneous)
+% P -   the camera matrices
+% X -   (output) estimated 3D point
+
+K = size(P,3);
+%A = zeros(K*3, 4);
+A = zeros(K*2, 4);
+
+if size(c,1) ~= 2 && size(c,1) ~=3
+    error('The points must be 2 or 3 dimensions!\n');
+end
+
+if size(c,1) == 2
+    c(end+1,:) = 1;
+end
+
+for k = 1:K
+%    crossM = [0         -c(3,k)     c(2,k)
+%              c(3,k)    0           -c(1,k)
+%              -c(2,k)   c(1,k)      0];
+%    A((k-1)*3+1:(k-1)*3+3,:) = crossM * P(:,:,k);
+    A((k-1)*2+1:(k-1)*2+2,:) = [c(1,k)*P(3,:,k) - c(3,k)*P(1,:,k);
+                                c(2,k)*P(3,:,k) - c(3,k)*P(2,:,k);];
+end
+
+[dummy,dummy,X] = svd(A,0);
+X = X(:,end);
+X = X./X(4);

bc_from_tri.m

@@ -0,0 +1,36 @@
+function bc = bc_from_tri(tri,x)
+
+% get the barycentric coordinates of a given point
+% tri - the three points, which are represented by column vectors
+% x   - the given point
+% bc  - (output) the barycentric coordinates
+
+if size(tri, 1) ~= size(x, 1)
+    error('different dimensions!');
+end
+
+if size(tri, 1) ~= 2 && size(tri, 1) ~= 3
+    error('The dimension must be 2 or 3!');
+end
+
+if size(tri, 1) == 2
+    tri(end+1,:) = 1;
+end
+
+if size(x, 1) == 2
+    x(end+1,:) = 1;
+end
+
+BA = tri(:, 2) - tri(:, 1);
+CA = tri(:, 3) - tri(:, 1);
+
+N = cross(BA, CA);
+SNN = sum(N.^2);
+
+NA = cross(tri(:, 3) - tri(:, 2), x - tri(:, 2));
+NB = cross(tri(:, 1) - tri(:, 3), x - tri(:, 3));
+% NC = cross(tri(:, 2) - tri(:, 1), x - tri(:, 1));
+alpha = N' * NA / SNN;
+beta = N' * NB / SNN;
+gamma = 1 - alpha - beta;
+bc = [alpha; beta; gamma];

cube2sphere.m

@@ -0,0 +1,57 @@
+function sphPos = cube2sphere(faceIdx, d, x, y)
+
+% transform a point on the cube map to a unit 3D vector
+%
+% faceIdx:  the face index
+%       d:  focal length
+%   (x,y):  the coordinate
+%
+%  sphPos:  the output unit vector
+
+switch faceIdx
+    case 1
+        if x <= d
+            if y <= d
+                phi = atan2(d-x, d-y);	
+            else
+                phi = atan2(y-d, d-x) + 0.5 * pi;
+            end
+        else
+            if y <= d
+                phi = atan2(d-y, x-d) + 1.5 * pi;
+            else
+                phi = atan2(x-d, y-d) + pi;
+            end
+        end
+        theta = atan2(hypot(x-d, y-d), d);
+    case 2
+        phi = atan2(x-d, d) + 0.5 * pi;
+        theta = atan2(y-d, d / abs(sin(phi))) + 0.5 * pi;
+    case 3
+        phi = atan2(x-d, d) + pi;
+        theta = atan2(y-d, d / abs(cos(phi))) + 0.5 * pi;
+    case 4
+        phi = atan2(x-d, d) + 1.5 * pi;
+        theta = atan2(y-d, d / abs(sin(phi))) + 0.5 * pi;
+    case 5
+        phi = atan2(x-d,d);
+        theta = atan2(y-d, d / abs(cos(phi))) + 0.5 * pi;
+    case 6
+        if x <= d
+            if y <= d
+                phi = atan2(d-y, d-x) + 0.5 * pi;
+            else
+                phi = atan2(d-x, y-d);
+            end
+        else
+            if y <= d
+                phi = atan2(x-d, d-y) + pi;
+            else
+                phi = atan2(y-d, x-d) + 1.5 * pi;
+            end
+        end
+        theta = pi - atan2(hypot(x-d, y-d), d);
+end
+phi = phi + pi;
+sphPos = [sin(theta)*cos(phi); sin(theta)*sin(phi); cos(theta)];
+