Skip to content

Commit

Permalink
Remove 'thresh' option in fitwarp2d.
Browse files Browse the repository at this point in the history 
Singular-value decomposition of the basis vectors can return
singular values that differ by 8--10 orders of magnitude. The
default value of thresh (1e-5) would discard several of these smaller
values.  That's a problem because the _inverses_ of the singular
values are actually used to produce the solution in _svd. Discarding
the largest singular values didn't seem to help either. Since it
wasn't entirely clear why some of the singular values were being
discarded, I commented out that line of code, and the associated
documentation, so now we keep all of the singular values.

I also added a detailed explanation of what fitwarp2d, _svd, and
_mkbasis actually do.  And of course, some tests.

If somebody with a lot of linear algebra experience is convinced
of a a good reason to discard some of these singular values, we
can revisit this (hence why I haven't removed the code, just
commented it out for now).
  • Loading branch information
@d-lamb
d-lamb committed Sep 29, 2017
1 parent dc51393 commit 0fc6c84
Show file tree
Hide file tree
Showing 2 changed files with 201 additions and 17 deletions.
161 changes: 145 additions & 16 deletions Lib/Image2D/image2d.pd
View file
Expand Up @@ -1721,7 +1721,7 @@ a coordinate transformation.
=for usage
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf. { options } )
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf, { options } )
Given a set of points in the output plane (C<$u,$v>), find
the best-fit (using singular-value decomposition) 2D polynomial
Expand All @@ -1748,9 +1748,15 @@ Options:
=for options
FIT - which terms to fit? default ones(byte,$nf,$nf)
=begin comment
old option, caused trouble
THRESH - in svd, remove terms smaller than THRESH * max value
default is 1.0e-5
=end comment
=over 4
=item FIT
Expand All @@ -1762,11 +1768,15 @@ used for the x and y polynomials; otherwise
C<$fit-E<gt>slice(":,:,(0)")> will be used for C<$px> and
C<$fit-E<gt>slice(":,:,(1)")> will be used for C<$py>.
=begin comment
=item THRESH
Remove all singular values whose valus is less than C<THRESH>
times the largest singular value.
=end comment
=back
The number of points must be at least equal to the number of
Expand Down Expand Up @@ -1810,6 +1820,40 @@ terms to fit (C<$nf*$nf> points for the default value of C<FIT>).
[ 1.25 -5.8546917e-18]
]
# A higher-degree polynomial should not affect the answer much, but
# will require more control points
$x = $x->glue(0,pdl(50,12.5, 37.5, 12.5, 37.5));
$y = $y->glue(0,pdl(50,12.5, 37.5, 37.5, 12.5));
$u = $u->glue(0,pdl(73,20,40,20,40));
$v = $v->glue(0,pdl(29,20,40,40,20));
( $px3, $py3 ) = fitwarp2d( $x, $y, $u, $v, 3 );
print "px3 =${px3}py3 =$py3";
px3 =
[
[-6.4981162e+08 71034917 -726498.95]
[ 49902244 -5415096.7 55945.388]
[ -807778.46 88457.191 -902.51612]
]
py3 =
[
[-6.2732159e+08 68576392 -701354.77]
[ 48175125 -5227679.8 54009.114]
[ -779821.18 85395.681 -871.27997]
]
#This illustrates an important point about singular value
#decompositions that are used in fitwarp2d: like all SVDs, the
#rotation matrices are not unique, and so the $px and $py returned
#by fitwarp2d are not guaranteed to be the "simplest" solution.
#They do still work, though:
($x3,$y3) = applywarp2d($px3,$py3,$u,$v);
print approx $x3,$x,1e-4;
[1 1 1 1 1 1 1 1 1]
print approx $y3,$y;
[1 1 1 1 1 1 1 1 1]
=head2 applywarp2d
=for ref
Expand All @@ -1828,23 +1872,95 @@ for more information on the format of C<$px> and C<$py>.
=cut
# use SVD to fit data. Assuming no errors.
sub _svd ($$$) {
=pod
=begin comment
Some explanation of the following three subroutines (_svd, _mkbasis,
and fitwarp2d): See Wolberg 1990 (full ref elsewhere in this
documentation), Chapter 3.6 "Polynomial Transformations". The idea is
that, given a set of control points in the input and output images
denoted by coordinates (x,y) and (u,v), we want to create a polynomial
transformation that relates u to linear combinations of powers of x
and y, and another that relates v to powers of x and y.
The transformations used here and by Wolberg differ slightly, but the
basic idea is something like this. For each u and each v, define a
transform:
u = (sum over j) (sum over i) a_{ij} x**i * y**j , (eqn 1)
v = (sum over j) (sum over i) b_{ij} x**i * y**j . (eqn 2)
We can write this in matrix notation. Given that there are M control
points, U is a column vector with M rows. A and B are vectors containing
the N coefficients (related to the degree of the polynomial fit). W
is an MxN matrix of the basis row-vectors (the x**i y**j).
The matrix equations we are trying to solve are
U = W A , (eqn 3)
V = W B . (eqn 4)
We need to find the A and B column vectors, those are the coefficients
of the polynomial terms in W. W is not square, so it has no inverse.
But is has a pseudo-inverse W^+ that is NxM. We are going to use that
pseudo-inverse to isolate A and B, like so:
W^+ U = W^+ W A = A (eqn 5)
W^+ V = W^+ W B = B (eqn 6)
We are going to get W^+ by performing a singular value decomposition
of W (to use some of the variable names below):
W = $svd_u x SIGMA x $svd_v->transpose (eqn 7)
W^+ = $svd_v x SIGMA^+ x $svd_u->transpose . (eqn 8)
Here SIGMA is a square diagonal matrix that contains the singular
values of W that are in the variable $svd_w. SIGMA^+ is the
pseudo-inverse of SIGMA, which is calculated by replacing the non-zero
singular values with their reciprocals, and then taking the transpose
of the matrix (which is a no-op since the matrix is square and
diagonal).
So the code below does this:
_mkbasis computes the matrix W, given control coordinates u and v and
the maximum degree of the polynomial (and the terms to use).
_svd takes the svd of W, computes the pseudo-inverse of W, and then
multiplies that with the U vector (there called $y). The output of
_svd is the A or B vector in eqns 5 & 6 above. Rarely is the matrix
multiplication explicit, unfortunately, so I have added EXPLANATIONs
below.
=end comment
=cut
sub _svd ($$) {
my $basis = shift;
my $y = shift;
my $thresh = shift;
# my $thresh = shift;
# if we had errors for these points, would normalise the
# basis functions, and the output array, by these errors here
# perform the SVD
my ( $svd_u, $svd_w, $svd_v ) = svd( $basis );
# remove any singular values
$svd_w *= ( $svd_w >= ($svd_w->max * $thresh ) );
# DAL, 09/2017: the reason for removing ANY singular values, much less
#the smallest ones, is not clear. I commented the line below since this
#actually removes the LARGEST values in SIGMA^+.
# $svd_w *= ( $svd_w >= ($svd_w->max * $thresh ) );
# The line below would instead remove the SMALLEST values in SIGMA^+, but I can see no reason to include it either.
# $svd_w *= ( $svd_w <= ($svd_w->min / $thresh ) );
# perform the back substitution
#
# EXPLANATION: same thing as $svd_u->transpose x $y->transpose.
my $tmp = $y x $svd_u;
#EXPLANATION: the division by (the non-zero elements of) $svd_w (the singular values) is
#equivalent to $sigma_plus x $tmp, so $tmp is now SIGMA^+ x $svd_u x $y
if ( $PDL::Bad::Status ) {
$tmp /= $svd_w->setvaltobad(0.0);
$tmp->inplace->setbadtoval(0.0);
Expand All @@ -1855,12 +1971,23 @@ sub _svd ($$$) {
$tmp *= ( 1 - $mask );
}
my $ans = sumover( $svd_v * $tmp );
return $ans;
#EXPLANATION: $svd_v x SIGMA^+ x $svd_u x $y
return sumover( $svd_v * $tmp );
} # sub: _svd()
#_mkbasis returns a piddle in which the k(=j*n+i)_th column is v**j * u**i
#k=0 j=0 i=0
#k=1 j=0 i=1
#k=2 j=0 i=2
#k=3 j=1 i=0
# ...
#each row corresponds to a control point
#and the rows for each of the control points look like this, e.g.:
# (1) (u) (u**2) (v) (vu) (v(u**2)) (v**2) ((v**2)u) ((v**2)(u**2))
# and so on for the next control point.
sub _mkbasis ($$$$) {
my $fit = shift;
my $npts = shift;
Expand Down Expand Up @@ -1896,7 +2023,7 @@ sub PDL::fitwarp2d {
my $v = shift;
my $nf = shift;
my $opts = PDL::Options->new( { FIT => ones(byte,$nf,$nf), THRESH => 1.0e-5 } );
my $opts = PDL::Options->new( { FIT => ones(byte,$nf,$nf) } ); #, THRESH => 1.0e-5 } );
$opts->options( $_[0] ) if $#_ > -1;
my $oref = $opts->current();
Expand All @@ -1906,9 +2033,9 @@ sub PDL::fitwarp2d {
unless $npts == $y->nelem and $npts == $u->nelem and $npts == $v->nelem
and $x->getndims == 1 and $y->getndims == 1 and $u->getndims == 1 and $v->getndims == 1;
my $svd_thresh = $$oref{THRESH};
croak "fitwarp2d: THRESH option must be >= 0."
if $svd_thresh < 0;
# my $svd_thresh = $$oref{THRESH};
# croak "fitwarp2d: THRESH option must be >= 0."
# if $svd_thresh < 0;
my $fit = $$oref{FIT};
my $fit_ndim = $fit->getndims();
Expand All @@ -1933,20 +2060,22 @@ sub PDL::fitwarp2d {
$ncoeff = $ncoeffx > $ncoeffy ? $ncoeffx : $ncoeffy;
}
croak "fitwarp2d: number of points must be >= \$ncoeff"
croak "fitwarp2d: number of points ($npts) must be >= \$ncoeff ($ncoeff)"
unless $npts >= $ncoeff;
# create the basis functions for the SVD fitting
my ( $basisx, $basisy );
$basisx = _mkbasis( $fitx, $npts, $u, $v );
if ( $fit_ndim == 2 ) {
$basisy = $basisx;
} else {
$basisy = _mkbasis( $fity, $npts, $u, $v );
}
my $px = _svd( $basisx, $x, $svd_thresh );
my $py = _svd( $basisy, $y, $svd_thresh );
my $px = _svd( $basisx, $x ); # $svd_thresh);
my $py = _svd( $basisy, $y ); # $svd_thresh);
# convert into $nf x $nf element piddles, if necessary
my $nf2 = $nf * $nf;
Expand Down
57 changes: 56 additions & 1 deletion t/image2d.t
View file
@@ -1,7 +1,7 @@
# -*-perl-*-
#

use Test::More tests => 29;
use Test::More tests => 47;
use Test::Exception;

use PDL;
Expand Down Expand Up @@ -210,3 +210,58 @@ polyfillv($im2,$ps) .= 25;
polyfill($im,$ps,25);
ok(all($im == $im2), "polyfill using default algorithm");
}

#######################
#warp2d and friends
{#this just runs the example in the documentation

my $x = pdl( 0, 0, 100, 100 );
my $y = pdl( 0, 100, 100, 0 );
# get warped to these positions
my $u = pdl( 10, 10, 90, 90 );
my $v = pdl( 10, 90, 90, 10 );
#
# shift of origin + scale x/y axis only
my $fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] );
my ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } );
ok(all(approx($px,pdl([-12.5,1.25],[0,0]),1e-13)),'px fitwarp2d linear restricted');
ok(all(approx($py,pdl([-12.5,0],[1.25,0]))),'py fitwarp2d linear restricted');
# Compared to allowing all 4 terms
( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 );
ok(all(approx($px, pdl([-12.5,1.25],[0,0]))),'px fitwarp2d linear unrestricted');
ok(all(approx($py, pdl([-12.5,0],[1.25,0]))),'py fitwarp2d linear unrestricted');
# A higher-degree polynomial should not affect the answer much, but
# will require more control points
$x = $x->glue(0,pdl(50,12.5, 37.5, 12.5, 37.5));
$y = $y->glue(0,pdl(50,12.5, 37.5, 37.5, 12.5));
$u = $u->glue(0,pdl(73,20,40,20,40));
$v = $v->glue(0,pdl(29,20,40,40,20));
my ( $px3, $py3 ) = fitwarp2d( $x, $y, $u, $v, 3 );
my ($x3,$y3) = applywarp2d($px3,$py3,$u,$v);
ok(all(approx($x3,$x,1e-4)),'px fitwarp2d quadratic unrestricted');
ok(all(approx($y3,$y,1e-4)),'py fitwarp2d quadratic unrestricted');

#define a simple grid image
my $img = (xvals(50,50) % 8 == 5 ) * (yvals(50,50) % 9 == 6); #stretch the y control points out a bit, and offset them too.
#get the control points
($u,$v) = whichND($img)->mv(0,-1)->double->dog;
#and shift it by 1 in horizontal and vertical directions
my $shift = $img->range([1,1],[$img->dims],'p');
#get the control points of the shifted image
($x,$y) = whichND($shift)->mv(0,-1)->double->dog;

use PDL::NiceSlice;
# we try 1st-, 2nd-, and 3rd-order fits, with and without restrictions to shift-and-scale-only
foreach my $deg(2,3,4){
my $fit = zeroes(byte,$deg,$deg,2);
$fit(:,(0),(0)).=1;
$fit((0),:,(1)).=1;
foreach my $restrict_fit(1,0){
my ($pxn,$pyn) = fitwarp2d($x,$y,$u,$v,$deg,$restrict_fit?{FIT=>$fit}:{});
my $out = warp2d($shift,$pxn,$pyn);
ok(all(approx($out,$img,1e-3)),'warp2d ' . ($restrict_fit?'':'un') . "restricted deg $deg values approx");
ok(all($out->rint==$img),'warp2d ' . ($restrict_fit?'':'un') . "restricted deg $deg rint exact");
}
}

}

0 comments on commit 0fc6c84

Please sign in to comment.