Beefy Boxes and Bandwidth Generously Provided by pair Networks
Welcome to the Monastery
 
PerlMonks  

Abstract image registration or feature detection [UPDATED w examples]

by kikuchiyo (Friar)
on Jul 01, 2022 at 14:55 UTC ( #11145220=perlquestion: print w/replies, xml ) Need Help??

kikuchiyo has asked for the wisdom of the Perl Monks concerning the following question:

I have a task that can be best summarized by the keywords in the title, and I wonder if there is a somewhat ready-made solution, preferably in Perl, that I've overlooked.

I have a set of points in a plane (originally coordinates of known features in an image), in two versions: one from a reference version of the image, the other from a distorted and warped version of the same image. The points themselves belong to two subsets: for the first subset, let's call them "known" points, I know the coordinates from both images, and for the second subset, "target" points, I know their coordinates only from the reference image. What I want is to determine the coordinates of these "target" points, based on the transformation determined by the corresponding "known" points.

Maybe I'm not looking right, but I haven't found anything besides Imagemagick, dodgy Matlab recipes and a bunch of research articles.

UPDATE:

Here are some example datasets. Columns are x, y, name. Points named p01..p12 are the "known" points, or the red points from the image downthread, q01..q05 are the "target" or blue points.

And here is an explanatory drawing, repeated for visibility: https://imgur.com/KxH62Sa

Each labeled red point from the reference image corresponds to the same labeled red point on the distorted image. And I know the coordinates for both. Similarly, one blue point on the reference image corresponds to one blue point on the distorted image, but I only know their coordinates on the reference, and I want to find them on the other.

Replies are listed 'Best First'.
Re: Abstract image registration or feature detection
by Corion (Patriarch) on Jul 01, 2022 at 15:00 UTC

    The first step would be to identify the algorithm suitable to your set of features.

    I'm not sure what kind of input data and training data you have. Maybe Scale-invariant feature transform ("SIFT") is already enough? I have a smallish implementation in Image::CCV, but most likely once you have the algorithm, you can find Perl bindings for whatever library implements it.

Re: Abstract image registration or feature detection
by pryrt (Abbot) on Jul 02, 2022 at 18:14 UTC
    I am hoping that etj will come by and make this more idiomatic PDL. But given that I'm a PDL hack, this is what I can come up with, by taking the approximate coordinates from the image supplied and doing some manual math using PDL:

    #!perl -l use 5.012; # //, strict, say use warnings; use PDL; use PDL::Matrix; print my $red_in = (mpdl [[58,48], [108,155], [186,80], [255,191], [3 +31,48]])->transpose->append(ones(1)); print my $red_out = (mpdl [[471,15], [531,141], [603,90], [682,227], [ +747,107]])->transpose->append(ones(1)); print my $blue_in = (mpdl [[125,73], [197,158], [282,94]])->transpose- +>append(ones(1)); print my $blue_out_for_checking = (mpdl [[542,64], [622,175], [701,138 +]])->transpose->append(ones(1)); print for my $sx = sum( $red_in->slice(0) ), my $sy = sum( $red_in->slice(1) ), my $s1 = sum( $red_in->slice(2) ), my $sxx = sum( $red_in->slice(0)**2 ), my $syy = sum( $red_in->slice(1)**2 ), my $sxy = sum( $red_in->slice(0)*$red_in->slice(1) ), my $sv = sum( $red_out->slice(0) ), my $sw = sum( $red_out->slice(1) ), my $sxv = sum( $red_in->slice(0)*$red_out->slice(0) ), my $sxw = sum( $red_in->slice(0)*$red_out->slice(1) ), my $syv = sum( $red_in->slice(1)*$red_out->slice(0) ), my $syw = sum( $red_in->slice(1)*$red_out->slice(1) ), ; print my $G = mpdl [ [ $sxx, $sxy, $sx ], [ $sxy, $syy, $sy ], [ $sx, $sy, $s1 ], ]; print my $Q = mpdl [ [ $sxv, $sxw ], [ $syv, $syw ], [ $sv, $sw ], ]; print my $P = inv($G) x $Q; print my $abcdef001 = ($P->transpose)->append(zeros(1))->set(2,2,1); print my $blue_out_calc = $abcdef001 x $blue_in; print $blue_out_calc - $blue_out_for_checking; # no more than one pixe +l off
    __END__ I am not a PDL expert, so I am sure there's a way with its linear alge +bra to skip doing all the math. But I enjoy the math derivation. You have a from IN to OUT, IN and OUT are of the form where each point #i is a column vector of t +he form [ x_i ] [ y_i ] So for each point, you can think of it as OUT = M * IN + T OUT = [ a b ] . IN + [ c ] [ d e ] [ f ] But if you want to work with all of IN and OUT at once, rather than in +dividually (embrace the power of the matrix), you can pad it out into a single multiplication as OUT = [ a b c ] [ x1 x2 x3 ... xN ] = [ v1 v2 v3 ... vN ] [ d e f ] * [ y1 y2 y3 ... yN ] [ w1 w2 w3 ... wN ] [ 0 0 1 ] [ 1 1 1 ... 1 ] [ 1 1 1 ... 1 ] So if you want to figure out the expanded M, you want to best fit on O +UT = M * IN (with the extra rows) To figure out the best fit, you basically have the equations v_i = a*Xi + b*Yi + c*1 w_i = d*Xi + e*Yi + f*1 Best fit is least-sum-squared-error, so err_v_i = a*Xi + b*Yi + c*1 - Vi err_w_i = d*Xi + e*Yi + f*1 - Wi And the sq error esq_Vi = a*Xi + b*Yi + c*1 + Vi + 2ab*Xi*Yi + 2ac*Xi - 2a*Xi +*Vi + 2bc*Yi - 2b*Yi*Vi - 2c*Vi esq_Wi = exact same form (I will just show derivation on the Vi te +rms for now) Sum up the total error: SSEv = sum(a*Xi + b*Yi + c*1 + Vi + 2ab*Xi*Yi + 2ac*Xi - 2a* +Xi*Vi + 2bc*Yi - 2b*Yi*Vi - 2c*Vi) = sum(a*Xi) + sum(b*Yi) + sum(c*1) + sum(Vi) + sum(2ab* +Xi*Yi) + sum(2ac*Xi) - sum(2a*Xi*Vi) + sum(2bc*Yi) - sum(2b*Yi*Vi) - +sum(2c*Vi) = a*sum(Xi) + b*sum(Yi) + c*sum(1) + sum(Vi) + 2ab*sum( +Xi*Yi) + 2ac*sum(Xi) - 2a*sum(Xi*Vi) + 2bc*sum(Yi) - 2b*sum(Yi*Vi) - +2c*sum(Vi) Need to find the minimum with respect to each of a, b, c (and equivale +ntly d,e,f) by setting each partial derivative to 0 dSSEv/da = 2a*sum(Xi) + 0 + 0 + 0 + 2b*sum(Xi*Yi) + 2c*sum(Xi) - +2*sum(Xi*Vi) + 0 - 0 - 0 = 0 dSSEv/db = 0 + 2b*sum(Yi) + 0 + 0 + 2a*sum(Xi*Yi) + 0 - 0 + 2c*su +m(Yi) - 2*sum(Yi*Vi) - 0 = 0 dSSEv/dc = 0 + 0 + 2c*sum(1) + 0 + 0 + 2a*sum(Xi) - 0 + 2b*sum(Yi) + - 0 - 2*sum(Vi) = 0 Divide by two and rearrange those three, grouping into a,b,c on one si +de and coefficientless on the other dSSEv/da: a*sum(Xi) + b*sum(Xi*Yi) + c*sum(Xi) = sum(Xi*Vi) dSSEv/db: a*sum(Xi*Yi) + b*sum(Yi) + c*sum(Yi) = sum(Yi*Vi) dSSEv/dc: a*sum(Xi) + b*sum(Yi) + c*sum(1) = sum(Vi) But that is a matrix multiplication: [ sum(Xi) sum(Xi*Yi) sum(Xi) ] [ a ] = [ sum(Xi*Vi) + ] [ sum(Xi*Yi) sum(Yi) sum(Yi) ] [ b ] = [ sum(Yi*Vi) + ] [ sum(Xi) sum(Yi) sum(1) ] [ c ] = [ sum(Vi) + ] And similarly from the esq_Wi, you get: [ sum(Xi) sum(Xi*Yi) sum(Xi) ] [ d ] = [ sum(Xi*Wi) + ] [ sum(Xi*Yi) sum(Yi) sum(Yi) ] [ e ] = [ sum(Yi*Wi) + ] [ sum(Xi) sum(Yi) sum(1) ] [ f ] = [ sum(Wi) + ] Note that the matrix on the left is the same for both. So we can merg +e those two into a single matrix multiplication [ sum(Xi) sum(Xi*Yi) sum(Xi) ] [ a d ] = [ sum(Xi*V +i) sum(Xi*Wi) ] [ sum(Xi*Yi) sum(Yi) sum(Yi) ] [ b e ] = [ sum(Yi*V +i) sum(Yi*Wi) ] [ sum(Xi) sum(Yi) sum(1) ] [ c f ] = [ sum(Vi) + sum(Wi) ] If we call those G x P = Q, then we can solve for the parameter matrix + P by P = inv(G) x Q But you can get the original abcdef001 matrix by transposing P and add +ing a row of (001) Then to find the blue outputs, just do BLUE_OUT = abcdef001 * BLUE_IN


    Edit: when I said " # no more than one pixel off" , that was not meant as a guarantee; just that this example was only approximately 1 pixel from the supplied locations
      Ironically given your kind words, and as proved by some emails just today on the PDL mailing lists, I am not at all a linear-algebra expert. However, if you're trying to solve linear least-squares systems, then check out mlls (using QR or LQ factorisation) and mllss (using SVD) in PDL::LinearAlgebra. If you want to try Levenberg-Marquardt, check out PDL::Fit::Levmar.

      One other observation for more idiomatic PDL usage is that if you're taking several slices (e.g. your sx, sy, s1) and doing the same to each, it'll be quicker to do something like:

      my ($sx, $sy, $s1) = $red_in->slice('0:2')->mv(0,-1)->sumover->dog;

      Thanks!

      This looks interesting, but unfortunately it dies for me with Dim mismatch in matmult of 3x3 x 3x2: 3 != 2 at the print my $P = inv($G) x $Q; line. PDL version 2.057 on Perl 5.34.1, if that matters.

        Weird, perl 5.30.0 with PDL 2.019 (which is the combo that shipped with Strawberry Perl's PDL edition) doesn't give that error.

        These are the outputs for my run (adding a print "perl version $]";print "PDL version $PDL::VERSION"; at the beginning to print the versions):

        Mathematically speaking, a [3x3] x [3x2] is the right balance of dimensions for a successful matrix multiplication (resulting in a [3x2] matrix). The use of PDL::Matrix used to be the way to make PDL match mathematical order for the dimensions, but maybe they've changed that between 2.019 and 2.057. (I don't think I can test, because the Strawberry Perl PDL setup is rather dependent on the libraries that Strawberry bundles, and I don't know if upgrading PDL will work... maybe I'll try in a sandbox to avoid ruining my "main" installation).

        I tried getting rid of the use PDL::Matrix, but then that gets rid of the mpdl command, causing errors during perl's compile stage, and I'm not sure what hoops I would need to jump through to make my script compatible with both 2.019 and 2.057.

Re: Abstract image registration or feature detection
by vr (Curate) on Jul 03, 2022 at 23:39 UTC

    Task as stated looks pretty much as case for interpolation to me. One trick pony, etc., but I even found my answer (Re: which function in PDL can do the same thing as matlab pcolor?) through SS for "Delaunay", see links in last paragraph there, especially ugly formulae for barycentric weights (didn't find ready-made solution on CPAN). Code is adjusted version of what I did a few years ago for some image hacking, unrelated to feature recognition.

    No idea what data pryrt used for a sample (same "red-blue" terms used here), and whether they are supposed to represent a trivial demo, or, indeed, an "interesting" non-affine distortion was applied. Anyway, results below seem to match well with his, even if it's just a hack -- X and Y of target are treated as separate functions to interpolate, with, furthermore, simple linear interpolation.

    use strict; use warnings; use feature qw/ say /; use POSIX qw/ round /; # use 5.022; use List::Util qw/ sum /; use Math::Geometry::Delaunay qw/ TRI_CONSTRAINED /; use constant EPSILON_NEGATIVE => -1e-6; sub key { pack 'II', @{ $_[ 0 ]}} my @red_in = ( [58,48], [108,155], [186,80], [255,191], [331,48] ); my @red_out = ( [471,15], [531,141], [603,90], [682,227], [747,107] ); my %mapped = map { key( $red_in[ $_ ]), $red_out[ $_ ]} 0 .. $#red_in +; my $tri = Math::Geometry::Delaunay-> new( TRI_CONSTRAINED ); $tri-> addPoints( \@red_in ); $tri-> triangulate; my @blue_in = ( [125,73], [197,158], [282,94] ); BLUE: for my $blue ( @blue_in ) { TRI: for my $elem ( @{ $tri-> elements }) { my $y23 = $elem-> [ 1 ][ 1 ] - $elem-> [ 2 ][ 1 ]; my $x32 = $elem-> [ 2 ][ 0 ] - $elem-> [ 1 ][ 0 ]; my $x13 = $elem-> [ 0 ][ 0 ] - $elem-> [ 2 ][ 0 ]; my $y13 = $elem-> [ 0 ][ 1 ] - $elem-> [ 2 ][ 1 ]; my $denominator = $y23 * $x13 + $x32 * $y13; my $xx3 = $blue-> [ 0 ] - $elem-> [ 2 ][ 0 ]; my $yy3 = $blue-> [ 1 ] - $elem-> [ 2 ][ 1 ]; my @weights; next TRI if EPSILON_NEGATIVE > ( $weights[ 0 ] = ( $y23 * $xx3 + $x32 * $yy3 ) / $denominator ); next TRI if EPSILON_NEGATIVE > ( $weights[ 1 ] = ( -$y13 * $xx3 + $x13 * $yy3 ) / $denominator ); next TRI if EPSILON_NEGATIVE > ( $weights[ 2 ] = 1 - $weights[ 0 ] - $weights[ 1 ]); printf "[%3d,%3d] => [%3d,%3d]\n", @$blue, map { my $coord = $_; round sum map { $weights[ $_ ] * $mapped{ key $elem-> [ $_ ]}[ $co +ord ] } 0 .. 2 # 3 vertices } 0, 1; # x, y next BLUE } die "point @$blue outside convex hull" } __END__ [125, 73] => [541, 63] [197,158] => [621,174] [282, 94] => [701,137]

      Thanks for this!

      I've updated the main node with real-ish example data, and while this algorithm is not perfect (it gives errors up to 7 pixels for some points), it might be good enough for my use case.

Re: Abstract image registration or feature detection
by LanX (Sage) on Jul 01, 2022 at 15:07 UTC
    > What I want is to determine the coordinates of these "target" points, based on the transformation determined by the corresponding "known" points.

    could you please elaborate? Do you know the corresponding point pairs (known_i,target_i)?

    There are many "transformations", I remember seeing a projection model based on a 4x4 matrix in uni.

    Determining the entries would mean to solve the resulting linear equations.

    Again it depends on the class of allowed transformations.

    Cheers Rolf
    (addicted to the Perl Programming Language :)
    Wikisyntax for the Monastery

      > could you please elaborate? Do you know the corresponding point pairs (known_i,target_i)?

      Each known_i and target_i point represents an unique named feature of the original image, and their locations (thus, the vectors between any two known_i and target_j point) are fixed on the reference image.

      > There are many "transformations", I remember seeing a projection model based on a 4x4 matrix in uni.

      The list of possible transformations include translation, rotation, skewing (so all affine transformations) plus barrel distortion, lens distortion etc. Think of a reference version of a poster in an image authoring software vs. the same poster photographed with a potato camera, not necessarily from a head-on orientation.

        > the vectors between any two known_i and target_j point

        I'm confused. Did you answer if you know which target point belongs to which source point or if you need to determine those "vectors"?

        > Think of a reference version of a poster in an image authoring software vs. the same poster photographed with a potato camera, not necessarily from a head-on orientation.

        Please take this for a start Camera Matrix

        It's a 3x4 matrix I remembered 4x4, so there might be more.

        Cheers Rolf
        (addicted to the Perl Programming Language :)
        Wikisyntax for the Monastery

        translation, rotation, skewing (so all affine transformations)
        You'll probably have more points than parameters in the affine transformation matrix, making this an overdetermined linear system. This is good: you can solve it in the least squares sense and handle a bit of the noise in the coordinates this way.
        plus barrel distortion, lens distortion
        This could be harder, because it's definitely not linear. Panorama stitching software like Hugin can let you optimise the radial distortion parameters at the same time as optimising the affine transform matrices; maybe you could look at its source code for some inspiration. Wikipedia lists a formula for a generalised distortion model, which you could fit the parameters for using Levenberg-Marquardt.
Re: Abstract image registration or feature detection
by Anonymous Monk on Jul 01, 2022 at 15:32 UTC
    The points themselves belong to two subsets: for the first subset, let's call them "known" points, I know the coordinates from both images
    Use regression to fit the parameters of the function transforming your points from the first subset to the second subset. Depending on the form of the distorsion, it could be multivariate linear regression, nonlinear regression you would need to provide a function (and starting values) for, or some sort of data-driven approach (P-splines, PLS, support vector machines with nonlinear basis functions, XGBoost...) that doesn't need an explicit function form and uses its many degrees of freedom to fit arbitrary function shapes. The PDL module family probably contains enough Levenberg-Marquardt and linear algebra to get you started.
    and for the second subset, "target" points, I know their coordinates only from the reference image.
    Use the parameters obtained above to perform the transformation.

    If you didn't know the corresponding points in both images, the inner part of the problem would stay the same, and the outer part of the problem would use Iterative closest point to try to estimate and update those correspondences between points.

Re: Abstract image registration or feature detection
by tybalt89 (Prior) on Jul 02, 2022 at 00:07 UTC

    How about providing us with a test data set or two...

    ( I may or may not be working or playing with something which may or may not work. (Hope i qualified that enough :) )

      > ( I may or may not be working or playing with something which may or may not work. (Hope i qualified that enough :) )

      maybe or maybe not ... :)

      Cheers Rolf
      (addicted to the Perl Programming Language :)
      Wikisyntax for the Monastery

Re: Abstract image registration or feature detection
by perlfan (Vicar) on Jul 03, 2022 at 07:22 UTC
    Image::PHash looks really great for duplicate and mirrored image detection. The best I've seen so far is unfortunately not accessible via Perl. YOLO/darknet is pretty neat, but the rub comes down to creating training data. I found that "roboflow dot com" makes it about as straightforward as one can get. Problem is, it's very expensive. Maybe with enough experience with it there can be a "ready made" Perl solution. Alien::OpenCV exists but I don't think that's gonna get you what you want. YMMV.
Re: Abstract image registration or feature detection [UPDATED w examples]
by kikuchiyo (Friar) on Jul 05, 2022 at 16:07 UTC

    It turns out that PDL has a set of routines for precisely this task: fitwarp2d, applywarp2d.

    The results are not that great, but it works out of the box.

      How are you measuring "results are not that great" and what are the results against the two test cases you have posted ?

      I'm curious how the algorithm I was playing with compares.

        How are you measuring "results are not that great"

        Not very scientifically, I've just looked at the largest difference in the coordinates.

        I've also plotted the results: https://imgur.com/a/W7j3yC9

        what are the results against the two test cases you have posted

        target_41:

        339 508 345.92214 512.55747 314 518 318.48381 516.9291 265 533 274.68006 530.37545 245 537 248.82152 537.19123 228 544 227.20604 544.26418

        target_996:

        207 550 213.15752 551.70605 206 569 209.79236 571.16027 205 602 208.69494 603.29318 203 621 207.30388 622.07696 201 638 207.08208 637.97918

        Perhaps "not that great" was not that great a description; this algorithm in fact produces comparable results to those in other replies in the thread.

        Note that I've messed up the direction of the transformation in my reply, the lines at the end should be

        my ($px, $py) = fitwarp2d($red_out_x, $red_out_y, $red_in_x, $red_in_y +, 2); my ($blue_new_x, $blue_new_y) = applywarp2d($px, $py, $blue_in_x, $blu +e_in_y); #say for $blue_out_x, $blue_out_y, $blue_new_x, $blue_new_y, $blue_out +_x - $blue_new_x, $blue_out_y - $blue_new_y;

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://11145220]
Approved by marto
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others about the Monastery: (8)
As of 2022-08-19 08:27 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found

    Notices?