Algorithm for obtaining the surface of a mirror
My colleague and I have been trying to implement an algorithm described in the paper "Recovering local shape of a mirror surface from reflection of a regular grid", primary author of which being Silvio Savarese. I believe we have come a way in our understanding of it, but I wondered if any one else has attempted it's implementation, or succeeded, and would be willing to provide some insight?
More specifically, the algorithm is suppose to be able to obtain the position, normal, and curvature of a reflective surface at a point using only 2 intersecting lines at that point, and their vertices. I have already successfully implemented an algorithm similar to that described in "3D Acquisition of Mirroring Objects" in c++ and Objective-C using OpenCV, but that unfortunately requires an initial depth value (at least it does for us because the solution they provide for obtaining the initial depth in sec 5.4 doesn't seem to work, probably because we have much lower resolution as we can't use calibrated color as a cue as they do). Which, since having to know an initial point is not very ideal in our application, we're exploring Silvio's method as a possible alternative.
Heres what we have gleaned of his method so far (keeping in mind that all I care about, at the moment, is the value 's' and it's corresponding normal 'nr' which can be easily obtained once s is known):
Calibrate the camera and scene (already done and tested and used in our prior algorithm)
Obtain qDot for the 2 lines by Eq 32 and project them onto the camera sensor plane 'l'
units along the view direction 'v'
Obtain the corresponding deltaP's (basic math since the scene is already calibrated)
Calculate the tanPhi's by normalizing qDot and taking the tangent of the angle between
that and the image plane's x axis.
for(double s=startValue, s<=endValue; s+=stepSize)
{
Solve for the vector position 'r' on the mirror surface (Eq 1)
Solve for the matrix 'R' via solving for 'np' and 'nr' to obtain conversion between XYZ
and UVW space (Page 4of14, second paragraph)
Solve cosTheta (shown with Eq3) and sinTheta (obtained by equation for cosTheta
via trig rules)
Solve Bu and Bv (shown with Eq3)
Solve for matrix 'T' (Eq 12)
Solve for h1, h2, and h3 (described right below Eq 4 on page 6of14)
Solve for V (Described half way down on page 8of14)
Check to see if Eq 17 is true (or within some epsilon), if is, return s,
else continue.
}
I realize that is not the most effective numerical solver for s (we're actually using a bit better solution), but understand that this is just a example and I've already graphed all chosen values for s (1in to 1000in with a 0.1in step, 40in or so being the correct value) at several regions on our part and it only approaches solving for Eq 17 one time, and diverges all others. So something is wrong, because it never actually reaches a solution. Furthermore, if you run this algorithm on all the points we have on our part, and project their corresponding s values (the ones that came closest to solving E17) into 3 space, it looks like noise and has no resemblance to our part.
I've tried several different methods for calculating tanPhi, and am not totally sure what is the correct way. Calculating it in image space, or calculating it by the angle between rDot and the U axis (obtained by Eq11 and the cross of np and nr)? The paper doesn't seem totally clear on that, or how to do it. We don't actually have intersecting lines, just dots, so I know tanPhi will be a little off just because it's changing as deltaT varies (which our deltaT is relatively large) but it's not varying that much, the surface we're looking at is only about 1200mm-1400mm radius, and I would think if we had it just a little off, the part would look noisy but still generally correct... instead it is totally off.
'l' I am still a little confused on as well, in another stackoverflow post 'Camera Properties for 3d using specular reflection', 'l' was stated to be the focal length, but there is an x and y component to the focal length in a pinhole camera model intrinsic matrix... Can we just average the 2? They are fairly close in our case, but what if they weren't would that still be sufficient?
I'm also not totally sure what needs to be converted to UVW space ("the principle reference system"), and what doesn't. I know it states on page 4, just before sec 2.1, that all calculations from here are in UVW (unless otherwise stated), but some things seem to make more sense to calculate in XYZ (like cosTheta and 'T') and he never mentions they should be. In Silvios phd thesis, for example, he doesn't say that until after obtaining cosTheta. So I'm assuming cosTeta can be obtained in XYZ, as it's just the cos of an angle that can be obtained in any space? But can't 'T' be obtained in any space as well, by the same logic, since it's just a rotation matrix?
And, with the converting (XYZ to UVW) does it really make sense to subtract r in every instance? I know for a point like 'p0' it makes sense... your just centering it on r, but for directions does it really make any sense (for calculating something like 'T')? For that you'd have to convert 'd' to UVW as well (if 'T' is in fact suppose to be calculated in UVW), but 'd' is a direction... centering it on r makes no sense. He's not very explicit on that.
Last thing I can think of, and this probably can wait until the math is all right, but can s be solved for in a closed form of the sort "s=f(measured and calibrated values)"? Silvio eludes to the fact that rDot can be solved for in closed form, after knowing s, but never really seems to say where s can be or not. I just wonder if this is going to forever require a numerical solver or not? The one I quick wrote up in c++ seems to work well, and converges quickly, but it'd obviously be nice if this was a closed function, plug and chug sort of speak!
My colleague has attempted to contact Silvio, as well as a few of the co-authors, but so far no luck. I would be greatly appreciative if someone could provide some insight here! It would be nice to see an algorithm like this out and in peoples hands, perhaps even implemented in OpenCV one day! I know obtaining the surface of a specular reflector is still very much an open problem, but it would be nice if this algorithm, and others like it, were available for testing!
The reconstruction of the mirror surface will result in a couple of parameters which describe that surface at position vector $r$. One of these, ranked first order, is the distance $s$ from that point at $r$ to the center of the camera at $c$.
The given paper is not too clear about how to recover $s$ from the data. Anyway, as the authors point out, the calculation of that distance $s$ seems to be the main subject of the earlier paper [2]:
S. Savarese and P. Perona: Local Analysis for 3D Reconstruction of Specular Surfaces - Part ii. ECCV, II 759–774 (2002).
In the end $s$ is characterized there as parameter at which some matrix $H=H(s)$ turns ambigious, see figure 4 on p. 769:
So the authors use numerical methods to solve $$ \mbox{det } H(s) = 0 $$ for $s$, see [2], eq. $(12)$, which seems different from your method, if I understood it correctly.