2d fit problem ?
jjm...@my-dejanews.com
Hi, I need some help with the following problem:
I have Accurate map coordinates(x,y) that I would like to
overlay on a distorted image(u,v). From the research I've
done a common way of handling this problem (rectification)
is with known control points and transformation functions:
u = f(x,y)
v = g(x,y)
I would like to use second order polynomials for the transform
function:
u = a0 + a1x + a2y + a3xy + a4x2 + a5x2
v = b0 + b1x + b2y + b3xy + b4x2 + b5x2
I can come up with a number of control points (where both x,y and u,v are
known). My problem is I don't know how to find the values for the
constants a0..a6, b0..b5. It looks like a least squares fit problem
but in 2 dimensions. The solution I'm looking for would be a function
I can hand a set of control points to and get the constants back, The
more controls the better the fit.
My programming skills are strong but I'm lacking the math. Source code
or a good description of how to (pointers to web sites, books etc.)
would be greatly appreciated.
Thanks,
jj
-----== Posted via Deja News, The Leader in Internet Discussion ==-----
http://www.dejanews.com/ Now offering spam-free web-based newsreading
David Ziskind
jjm...@my-dejanews.com wrote in article
<6kv545$1cv$1...@nnrp1.dejanews.com> regarding a problem class that
includes second degree, two dimensional least squares regression
(LSR). Basically, LSR for any degree and any dimension is computed
following the same approach used for the classic case of one degree
and one dimension. Second and higher degree introduces a certain
complexity as to linearly ordering the parameters.
In the application discussed, introduce the following definitions:
P(n) is the label for the "n th point". x(n,1), and x(n,2) are the
first and second coordinates of P(n) using the Accurate Map. x(n) is
a two position vector. In this post, arguments that *might* be
thought of as subscripts will be written as arguments, rather than as
subscripts. Commas and parentheses are omitted as possible.
Similarly, u(n,1) and u(n,2) are the first and second coordinates of
P(n) using the Distorted Map.
F, a 2 dim into 2 dim function, maps an x-space vector into a u-space
vector. F is the approximating function and the objective is to
design it so that F(x(n)) is as close to u(n) as possible. To allow
this design/construction, F actually admits as arguments 2 x 6 = 12
additional parameters, so that technically the approximating function
is a mapping from 14 space into 2 space. However, we will think of F
as the 2 dim into 2 dim function obtained when the first 12 arguments
are fixed. The polynomial form for F will be deferred for the
moment.
Following the traditional least squares plan, the error function E is
defined by: E(a) =def SS{k=1;2: n=1;T: (u(nk)-(F(x(n))k))^2} (1)
- S{n=a;b: f(n)} means the sum of f(n) taken over n = a through b.
SS{n=a;b: m=c;d: f(n,m)} is the analogous double sum.
- E(a) is the evaluation of the error function when the parameter
vector has value a.
- E(a) has a-type dependency because F(x(n)) on the right hand side
has the above mentioned concealed dependency on the 12 parameters.
- The first 6 of the 12 parameters are the 6 numbers obtained when k
in the following list is set = to 1:
a(k00) a(k10) a(k11) a(k20) a(k21) a(k22) . (See more below)
- (F(x(n))k means the k th element of F evaluated at x(n).
- T is the number of control points.
Property LSM (Least Squares Minimization) is defined by:
a* has the LSM property if and only if (2)
E(a*) <= E(a) for all nearby a
The computational objective is find a 12 component vector which has
the LSM property. The LSM property is the same as the local minimum
property.
As is traditional in these problems, the solution focuses on
computing a "critical point". Once such a critical point has been
identified, there are various tests that can be applied to determine
whether the critical point is a local minimum, is not (perhaps is a
saddle point), etc. A good discussion can be found in Bronshtein and
Semendyayev (Handbook of Mathematics, Van Nostrand Reinhold, New
York, 1985.) In this post, the end issue as to whether such a
critical point is indeed a local minimum, is left open.
Finding a* is done by a piecemeal approach -- the first six elements
of the 12 component vector "a" are computed in one step (k=1 related
elements), followed by a second step for the last six (k=2 related
elements).
We now are ready to define the form of F. Let w be a generic two
position vector which represents one of: x(1), x(2), x(3), .... etc.
(F(w))k =def SS{j=0;2: m=0;j: a(kjm)*w(j)*w(m)} (3)
with k = 1, 2;
and the understanding that on the right-hand side, expressions
obtained by making a replacement for w *which end in (0)* are
non-literal and should be intepreted as unity (1). This device
places the constant and linear terms into a quadratic form
simplifying the notation.
Example: (F(x(3))1 produces on the right-hand side various terms
one of which is: a(1,1,0)*(x(3))(1)*(x(3))(0). The right term
"(x((3))(0)" is interpreted as 1. Note also that x(3,0) is equal
to (x(3))(0) by definition and so likewise is interpreted as 1,
etc.
By the calculus, IF a* has the LSM property, THEN for the appropriate
12 (kjm):
(derivative of E(a) with respect to a(kjm))(evaluated at a=a*)
= 0 (4)
In this prescription, differentiating w/r to "a(1,0,0)" means 1st
partial, w/r to "a(1,1,0)" means 2nd partial, ..., w/r to "a(2,0,0)"
means 7th partial, etc.
******************************************************************
In what follows, the use of a* to suggest a distinguished value of
"a" with the LSM property, will be dropped. From this point on, for
brevity, "a" (without the star) will be used represent a value with
(hopefully) the LSM property (ie, a local minimum). There is no
further need to refer to generic "a" so there should be no confusion.
******************************************************************
Differentiating E(a) is very simple due to the polynomial form of F.
Before rearrangment and relabeling, the result FOR ONE FIXED VALUE OF
K, SAY FOR EXAMPLE K=1 but left general below, is:
S{n=1;T: u(nk)*x(nj)*x(nm)} = (5)
SSS{n=1;T: p=0;2: q=0;p: a(kpq)*x(np)*x(nq)*x(nj)*x(nm)}
with j = 0, 1, 2
and m = 0, ..., j
Basically, the above relation gives the wanted recipe for calculating
the critical point. The main problem is that the six numbers
a(k,0,0), a(k,1,0), ... are enumerated with double subscripts and to
solve the system a linear ordering is wanted. It is time to
introduce the L(eft) and R(ight) functions.
We start out by writing down a linear ordering of the "a" elements as
they might appear in matrix form:
a(k,0,0) 1
a(k,1,0) 2 a(k,1,1) 3
a(k,2,0) 4 a(k,2,1) 5 a(k,2,2) 6
The linear ordering labels (the numbers 1, 2, 3, ... to the right)
are written in as the a's appear traversing the lower triangular
matrix.
The L and R values for argument = 1, 2, 3, ... are determined by:
going to the argument, taking the L value as the left hand subscript
(forget about k), and taking the R value as the right hand subscript.
Thus:
If z = 1, 2, 3, 4, 5, 6 (6)
then L(z) = 0, 1, 1, 2, 2, 2
and R(z) = 0, 0, 1, 0, 1, 2
In the system of equations (5), we substitute (Lz, Rz) for (j, m) and
"substitute" (Ly, Ry) for (p, q). (The latter is basically a
formalism.) The validity of both operations is provable. We obtain
the logically equivalent system:
S{n=1;T: u(n,k)*x(n,Lz)*x(n,Rz)} = (7)
S{y=1;6: S{n=1;T: x(n,Ly)*x(n,Ry)*x(n,Lz)*x(n,Rz)}
* a(k,Ly,Ry)
}
with z = 1, 2, ..., 6
As should be emphasized, k IS FIXED. We have six equations in six
unknowns, namely: a(k,L1,R1), a(k,L2,R2), ..., a(k,L6,R6).
The first product term on the right is the (z,y) th term of a certain
matrix; call it M. The term on the left is the z th term of a
certain vector; call it Q.
1. Compute and store the components M(z,y) by the formula
M(z,y) = S{n=1;T: x(n,Ly)*x(n,Ry)*x(n,Lz)*x(n,Rz)}
2. Compute and store the components Q(z) by the formula
Q(z) = S{n=1;T: u(n,k)*x(n,Lz)*x(n,Rz)}
3. Invert the 6 by 6 matrix M.
4. Multiply Inverse-of-M onto Q.
5. The result is the six component "a" vector, the six unknowns
mentioned above. Note that element one of the "a" vector is in
the triple indexing scheme, a(k,0,0), etc.
6. Do 1 through 5 for k=1 and then for k=2. (Skip 1 and 3 after the
the first pass -- M is free of k.) The result will be a set
of a(k,j,m) entries for k = 1 and 2, and (j,m) = lower triangular
matrix indices. These a(k,j,m) entries are just the values wanted
for the (F(w))k approximating function defined above in (3).
These "a" entries determine F for a given set of control points and
thus the problem is solved (provided however that the critical point
is in fact a local minimum). In the F definition, for exposition, we
said w represented a control point; in the application, w can be
anything.
For higher degree (>2) approximating polynomials (say degree d), I
suggest looking at a d level tree approach as a way of linearly
ordering the (d+1) indexed "a" elements.
David Ziskind
zis...@ntplx.net