1 | %vgg_poly3d_orthorectify Frontoparallel rectification of image of 3D polygon.
|
---|
2 | % [H,imsize] = vgg_poly3d_orthorectify(Q,u) finds homography H that
|
---|
3 | % removes perspective distortion of image Q*hom(u) of a 3D polygon. Parameters:
|
---|
4 | % u ... double(2,N), inhomog. coordinates of polygon vertices measured in an orthonormal
|
---|
5 | % coordinate frame in the 3D polygon plane.
|
---|
6 | % Q ... double(3,3), homography that maps u to image plane.
|
---|
7 | % Image of the polygon is thus x = Q*hom(u) (homog. coords.).
|
---|
8 | % H ... double(3,3), rectifying homography
|
---|
9 | % imsize ... double(2,1), size of the target image
|
---|
10 | % The rectification is done as follows:
|
---|
11 | % Irect = vgg_warp_H(I{k},inv(H),'linear',[1 imsize(2) 1 imsize(1)]).
|
---|
12 | %
|
---|
13 | % [H,imsize] = vgg_poly3d_orthorectify(Q,u,smax) allows to specify maximum ratio of
|
---|
14 | % area element in output and input image (if omitted, smax=5). This is to prevent
|
---|
15 | % the output image from becoming huge if the perspective distortion is large.
|
---|
16 | %
|
---|
17 | % If Q is K-cell of double(3,3), rectification is done simultaneously for K images
|
---|
18 | % of the 3D polygon, so that pixels with the same coordinates correspond in the
|
---|
19 | % rectified images. Then H is K-cell of double(3,3) and imsize is double(2,K).
|
---|
20 | %
|
---|
21 | % [H,imsize] = vgg_poly3d_orthorectify(Q,u,smax,in_imsize) specifies sizes of input
|
---|
22 | % images. This is needed if some vertices project outside some input images.
|
---|
23 |
|
---|
24 | function [H,imsize] = vgg_poly3d_frontorectify(H,u,smax,imsize)
|
---|
25 |
|
---|
26 | if ~iscell(H)
|
---|
27 | H = {H};
|
---|
28 | end
|
---|
29 |
|
---|
30 | if nargin<3
|
---|
31 | smax = 5;
|
---|
32 | end
|
---|
33 |
|
---|
34 | % find scales of H in all u
|
---|
35 | for k = 1:length(H)
|
---|
36 | if nargin>3 % consider only vertices inside input image
|
---|
37 | i = boxclipu( [[1;1] imsize(:,k)], vgg_get_nonhomg(H{k}*vgg_get_homg(u)) );
|
---|
38 | else
|
---|
39 | i = logical(ones(1,size(u,2)));
|
---|
40 | end
|
---|
41 | s(k,~i) = nan;
|
---|
42 | if nnz(i)==0, continue, end
|
---|
43 | s(k,i) = sqrt(abs( dethomog(H{k},u(:,i)) ));
|
---|
44 | end
|
---|
45 |
|
---|
46 | % determine the scale;
|
---|
47 | % forbid scaling polygons more than smax-times
|
---|
48 | s = s(:);
|
---|
49 | s(isnan(s)) = [];
|
---|
50 | if isempty(s) % polygon is not in the image
|
---|
51 | imsize = [0;0];
|
---|
52 | return
|
---|
53 | end
|
---|
54 | if nargin < 5
|
---|
55 | smax = inf;
|
---|
56 | end
|
---|
57 | s = min(s) * min(smax,max(s)/min(s));
|
---|
58 |
|
---|
59 | % scale H and u consistently by s
|
---|
60 | u = u*s;
|
---|
61 | Hs = diag([1/s 1/s 1]);
|
---|
62 | for k = 1:length(H)
|
---|
63 | H{k} = H{k}*Hs;
|
---|
64 | end
|
---|
65 |
|
---|
66 | % translate the polygon to the origin
|
---|
67 | b = bboxu(u);
|
---|
68 | imsize = ceil( b(:,2) - b(:,1) );
|
---|
69 | Ht = eye(3);
|
---|
70 | Ht(1:2,3) = b(:,1) + 1;
|
---|
71 |
|
---|
72 | for k = 1:length(H)
|
---|
73 | H{k} = H{k}*Ht;
|
---|
74 | end
|
---|
75 |
|
---|
76 | if length(H)==1
|
---|
77 | H = H{1};
|
---|
78 | end
|
---|
79 |
|
---|
80 | return
|
---|
81 |
|
---|
82 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
|
---|
83 |
|
---|
84 | % D = dethomog(H,q) Determinant of Jacobian of homography transformation r = nhom(H*hom(q)), D = det(dr/dq). Works vector-wise for size(q)=[2 ?].
|
---|
85 | %
|
---|
86 | function D = dethomog(H,q)
|
---|
87 | % Derivative of homography r = h(q) = nhom(H*[q;1]) is:
|
---|
88 | % dh = [eye(2) -r]/y(3)*H(:,1:2), where y = H*[q;1].
|
---|
89 | q(3,:) = 1;
|
---|
90 | r = vgg_get_nonhomg(H*q);
|
---|
91 | D = ((H(1,1)-r(1,:)*H(3,1)).*(H(2,2)-r(2,:)*H(3,2)) - (H(1,2)-r(1,:)*H(3,2)).*(H(2,1)-r(2,:)*H(3,1))) ./ (H(3,:)*q).^2;
|
---|
92 | return
|
---|
93 |
|
---|
94 |
|
---|
95 | function i = boxclipu(b,u)
|
---|
96 | % i = boxclipu(b,u) Returns logical indices of points inside box b.
|
---|
97 | i = all( u >= b(:,row1(u)) & u <= b(:,2*row1(u)), 1 );
|
---|
98 | return
|
---|
99 |
|
---|
100 |
|
---|
101 | function ub = bboxu(u)
|
---|
102 | % ub = bboxu(u) Bounding box of point cloud u.
|
---|
103 | ub = [min(u')' max(u')'];
|
---|
104 | return |
---|