diff --git a/include/vigra/delaunay.hxx b/include/vigra/delaunay.hxx new file mode 100644 index 000000000..19ff9160f --- /dev/null +++ b/include/vigra/delaunay.hxx @@ -0,0 +1,473 @@ +/************************************************************************/ +/* */ +/* Copyright 2013-2017 by Benjamin Seppke */ +/* */ +/* This file is part of the VIGRA computer vision library. */ +/* The VIGRA Website is */ +/* http://hci.iwr.uni-heidelberg.de/vigra/ */ +/* Please direct questions, bug reports, and contributions to */ +/* ullrich.koethe@iwr.uni-heidelberg.de or */ +/* vigra@informatik.uni-hamburg.de */ +/* */ +/* Permission is hereby granted, free of charge, to any person */ +/* obtaining a copy of this software and associated documentation */ +/* files (the "Software"), to deal in the Software without */ +/* restriction, including without limitation the rights to use, */ +/* copy, modify, merge, publish, distribute, sublicense, and/or */ +/* sell copies of the Software, and to permit persons to whom the */ +/* Software is furnished to do so, subject to the following */ +/* conditions: */ +/* */ +/* The above copyright notice and this permission notice shall be */ +/* included in all copies or substantial portions of the */ +/* Software. */ +/* */ +/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ +/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ +/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ +/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ +/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ +/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ +/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ +/* OTHER DEALINGS IN THE SOFTWARE. */ +/* */ +/************************************************************************/ + +/* + * Parts of this file, namely the first section including the delaunay + * function are taken from the LEMON distribution, lgf-gen.cc. Here, the + * following copyright applies: + * + * Copyright (C) 2003-2009 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef VIGRA_DELAUNAY_HXX +#define VIGRA_DELAUNAY_HXX + +#ifndef WITH_LEMON + #error "Should only be included with flag \"WITH_LEMON\"" +#endif + +#include +#include +#include +#include + +namespace vigra { + + namespace detail + { + using namespace lemon; + typedef dim2::Point Point; + + GRAPH_TYPEDEFS(ListGraph); + + struct Part + { + int prev, curr, next; + + Part(int p, int c, int n) : prev(p), curr(c), next(n) {} + }; + + inline std::ostream& operator<<(std::ostream& os, const Part& part) + { + os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')'; + return os; + } + + inline double circle_point(const Point& p, const Point& q, const Point& r) + { + double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y); + if (a == 0) return std::numeric_limits::quiet_NaN(); + + double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) + + (q.x * q.x + q.y * q.y) * (r.y - p.y) + + (r.x * r.x + r.y * r.y) * (p.y - q.y); + + double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) + + (q.x * q.x + q.y * q.y) * (r.x - p.x) + + (r.x * r.x + r.y * r.y) * (p.x - q.x); + + double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) + + (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) + + (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y); + + return d / (2 * a) + std::sqrt((d * d + e * e) / (4 * a * a) + f / a); + } + + inline bool circle_form(const Point& p, const Point& q, const Point& r) + { + return rot90(q - p) * (r - q) < 0.0; + } + + inline double intersection(const Point& p, const Point& q, double sx) + { + const double epsilon = 1e-8; + + if (p.x == q.x) return (p.y + q.y) / 2.0; + + if (sx < p.x + epsilon) return p.y; + if (sx < q.x + epsilon) return q.y; + + double a = q.x - p.x; + double b = (q.x - sx) * p.y - (p.x - sx) * q.y; + double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare(); + return (b - std::sqrt(d)) / a; + } + + struct YLess + { + YLess(const std::vector& points, double& sweep) + : _points(points), + _sweep(sweep) + { + } + + bool operator()(const Part& l, const Part& r) const + { + const double epsilon = 1e-8; + + // std::cerr << l << " vs " << r << std::endl; + double lbx = l.prev != -1 ? + intersection(_points[l.prev], _points[l.curr], _sweep) : + - std::numeric_limits::infinity(); + double rbx = r.prev != -1 ? + intersection(_points[r.prev], _points[r.curr], _sweep) : + - std::numeric_limits::infinity(); + double lex = l.next != -1 ? + intersection(_points[l.curr], _points[l.next], _sweep) : + std::numeric_limits::infinity(); + double rex = r.next != -1 ? + intersection(_points[r.curr], _points[r.next], _sweep) : + std::numeric_limits::infinity(); + + if (lbx > lex) std::swap(lbx, lex); + if (rbx > rex) std::swap(rbx, rex); + + if (lex < epsilon + rex && lbx + epsilon < rex) return true; + if (rex < epsilon + lex && rbx + epsilon < lex) return false; + return lex < rex; + } + + const std::vector& _points; + double& _sweep; + }; + + struct BeachIt; + + typedef std::multimap SpikeHeap; + + typedef std::multimap Beach; + + struct BeachIt + { + Beach::iterator it; + BeachIt(Beach::iterator iter) + : it(iter) + { + } + }; + } //end of namespace vigra::detail + +void delaunay(detail::ListGraph & g, detail::ListGraph::NodeMap& coords) +{ + using namespace detail; + + typedef std::vector > SiteHeap; + + std::vector points; + std::vector nodes; + + for (NodeIt it(g); it != INVALID; ++it) + { + nodes.push_back(it); + points.push_back(coords[it]); + } + + std::cerr << "Nodes: " << nodes.size() << "\n"; + + SiteHeap siteheap(points.size()); + + double sweep; + + for (int i = 0; i < int(siteheap.size()); ++i) + { + siteheap[i] = std::make_pair(points[i].x, i); + } + + std::sort(siteheap.begin(), siteheap.end()); + sweep = siteheap.front().first; + + YLess yless(points, sweep); + Beach beach(yless); + + SpikeHeap spikeheap; + + std::set > arcs; + + int siteindex = 0; + { + SiteHeap front; + + while (siteindex < int(siteheap.size()) && + siteheap[0].first == siteheap[siteindex].first) + { + front.push_back(std::make_pair(points[siteheap[siteindex].second].y, + siteheap[siteindex].second)); + ++siteindex; + } + + std::sort(front.begin(), front.end()); + + for (int i = 0; i < int(front.size()); ++i) + { + int prev = (i == 0 ? -1 : front[i - 1].second); + int curr = front[i].second; + int next = (i + 1 == int(front.size()) ? -1 : front[i + 1].second); + + beach.insert(std::make_pair(Part(prev, curr, next), + spikeheap.end())); + } + } + + while (siteindex < int(points.size()) || !spikeheap.empty()) + { + + SpikeHeap::iterator spit = spikeheap.begin(); + + if ( siteindex < int(points.size()) && + (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) + { + int site = siteheap[siteindex].second; + sweep = siteheap[siteindex].first; + + Beach::iterator bit = beach.upper_bound(Part(site, site, site)); + + if (bit->second != spikeheap.end()) + { + spikeheap.erase(bit->second); + } + + int prev = bit->first.prev; + int curr = bit->first.curr; + int next = bit->first.next; + + beach.erase(bit); + + SpikeHeap::iterator pit = spikeheap.end(); + if (prev != -1 && + circle_form(points[prev], points[curr], points[site])) + { + double x = circle_point(points[prev], points[curr], points[site]); + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); + pit->second.it = beach.insert(std::make_pair(Part(prev, curr, site), pit)); + } + else + { + beach.insert(std::make_pair(Part(prev, curr, site), pit)); + } + + beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end())); + + SpikeHeap::iterator nit = spikeheap.end(); + + if (next != -1 && + circle_form(points[site], points[curr],points[next])) + { + double x = circle_point(points[site], points[curr], points[next]); + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); + + nit->second.it = beach.insert(std::make_pair(Part(site, curr, next), nit)); + } + else + { + beach.insert(std::make_pair(Part(site, curr, next), nit)); + } + ++siteindex; + } + else + { + sweep = spit->first; + + Beach::iterator bit = spit->second.it; + + int prev = bit->first.prev; + int curr = bit->first.curr; + int next = bit->first.next; + + + // if(prev != -1 && next !=-1) + // std::cerr << "Triangle: (" << prev << ", " << curr << ", " << next << ")\n"; + { + std::pair arc; + + arc = prev < curr ? std::make_pair(prev, curr) : std::make_pair(curr, prev); + + if (arcs.find(arc) == arcs.end()) + { + arcs.insert(arc); + g.addEdge(nodes[prev], nodes[curr]); + } + + arc = curr < next ? std::make_pair(curr, next) : std::make_pair(next, curr); + + if (arcs.find(arc) == arcs.end()) + { + arcs.insert(arc); + g.addEdge(nodes[curr], nodes[next]); + } + } + + Beach::iterator pbit = bit; --pbit; + int ppv = pbit->first.prev; + Beach::iterator nbit = bit; ++nbit; + int nnt = nbit->first.next; + + if (bit->second != spikeheap.end()) + { + spikeheap.erase(bit->second); + } + if (pbit->second != spikeheap.end()) + { + spikeheap.erase(pbit->second); + } + if (nbit->second != spikeheap.end()) + { + spikeheap.erase(nbit->second); + } + + beach.erase(nbit); + beach.erase(bit); + beach.erase(pbit); + + SpikeHeap::iterator pit = spikeheap.end(); + + if (ppv != -1 && ppv != next && + circle_form(points[ppv], points[prev], points[next])) + { + double x = circle_point(points[ppv], points[prev], points[next]); + + if (x < sweep) + { + x = sweep; + } + pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); + pit->second.it = beach.insert(std::make_pair(Part(ppv, prev, next), pit)); + } + else + { + beach.insert(std::make_pair(Part(ppv, prev, next), pit)); + } + + SpikeHeap::iterator nit = spikeheap.end(); + if (nnt != -1 && prev != nnt && + circle_form(points[prev], points[next], points[nnt])) + { + double x = circle_point(points[prev], points[next], points[nnt]); + if (x < sweep) + { + x = sweep; + } + nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end()))); + nit->second.it = beach.insert(std::make_pair(Part(prev, next, nnt), nit)); + } + else + { + beach.insert(std::make_pair(Part(prev, next, nnt), nit)); + } + } + } + + //Add missing convex hull arcs to complete triangulation + for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) + { + int curr = it->first.curr; + int next = it->first.next; + + + if (next == -1) + { + continue; + } + + std::pair arc; + arc = curr < next ? std::make_pair(curr, next) : std::make_pair(next, curr); + + if (arcs.find(arc) == arcs.end()) + { + arcs.insert(arc); + g.addEdge(nodes[curr], nodes[next]); + } + } +} + +std::vector > trianglesFromDelaunay(detail::ListGraph & g, detail::ListGraph::NodeMap& nodeIDs) +{ + using namespace detail; + + std::vector > result; + + std::set > arcs; + + for(ArcIt aIt(g); aIt != INVALID; ++aIt) + { + int curr = nodeIDs[g.source(aIt)]; + int next = nodeIDs[g.target(aIt)]; + + if (curr < next) + arcs.insert(std::pair(curr, next)); + } + + for(ArcIt aIt(g); aIt != INVALID; ++aIt) + { + int curr = nodeIDs[g.source(aIt)]; + int next = nodeIDs[g.target(aIt)]; + //Order: curr < next + if (curr > next) + continue; + + int found=0; + + //Find third node of the triangle + for (NodeIt it(g); it != INVALID; ++it) + { + //Max. 2 triangles per edge + if(found == 2) + break; + + int cand = nodeIDs[it]; + + //Order: curr < next < cand + if (cand > next) + { + std::pair arc1; + std::pair arc2; + arc1 = curr < cand ? std::make_pair(curr, cand) : std::make_pair(cand, curr); + arc2 = next < cand ? std::make_pair(next, cand) : std::make_pair(cand, next); + + if (arcs.find(arc1) != arcs.end() && arcs.find(arc2) != arcs.end()) + { + //std::cerr << "Triangle: (" << curr << ", " << next << ", " << cand <<")\n"; + result.push_back(vigra::triple(curr, next, cand)); + found++; + } + } + } + } + return result; +} + +} //end of namespace vigra + +#endif //VIGRA_DELAUNAY_HXX diff --git a/include/vigra/piecewiseaffine_registration.hxx b/include/vigra/piecewiseaffine_registration.hxx new file mode 100644 index 000000000..c6edcfa13 --- /dev/null +++ b/include/vigra/piecewiseaffine_registration.hxx @@ -0,0 +1,275 @@ +/************************************************************************/ +/* */ +/* Copyright 2013-2017 by Benjamin Seppke */ +/* */ +/* This file is part of the VIGRA computer vision library. */ +/* The VIGRA Website is */ +/* http://hci.iwr.uni-heidelberg.de/vigra/ */ +/* Please direct questions, bug reports, and contributions to */ +/* ullrich.koethe@iwr.uni-heidelberg.de or */ +/* vigra@informatik.uni-hamburg.de */ +/* */ +/* Permission is hereby granted, free of charge, to any person */ +/* obtaining a copy of this software and associated documentation */ +/* files (the "Software"), to deal in the Software without */ +/* restriction, including without limitation the rights to use, */ +/* copy, modify, merge, publish, distribute, sublicense, and/or */ +/* sell copies of the Software, and to permit persons to whom the */ +/* Software is furnished to do so, subject to the following */ +/* conditions: */ +/* */ +/* The above copyright notice and this permission notice shall be */ +/* included in all copies or substantial portions of the */ +/* Software. */ +/* */ +/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ +/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ +/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ +/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ +/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ +/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ +/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ +/* OTHER DEALINGS IN THE SOFTWARE. */ +/* */ +/************************************************************************/ + +#ifndef VIGRA_PIECEWISEAFFINE_REGISTRATION_HXX +#define VIGRA_PIECEWISEAFFINE_REGISTRATION_HXX + +#ifndef WITH_LEMON + #error "Should only be included with flag \"WITH_LEMON\"" +#endif + +#include +#include + +#include "mathutil.hxx" +#include "matrix.hxx" +#include "linear_solve.hxx" +#include "tinyvector.hxx" +#include "splineimageview.hxx" +#include "affine_registration.hxx" +#include "delaunay.hxx" + +namespace vigra +{ + +/** \addtogroup Registration Image Registration + */ +//@{ + +//used vigra types for points and triangles +typedef vigra::TinyVector PointType; +typedef vigra::triple TriangleType; +typedef std::pair > TriangleTransformationType; + +/** + * Tiny helper to check whether a point is inside a triangl. + */ +inline +bool isInside(const PointType & point, const TriangleType & triangle, bool allowBorder=true) +{ + // Compute direction vectors + PointType v0 = triangle.third - triangle.first; + PointType v1 = triangle.second - triangle.first; + PointType v2 = point - triangle.first; + + // Compute dot products + double dot00 = dot(v0, v0), + dot01 = dot(v0, v1), + dot02 = dot(v0, v2), + dot11 = dot(v1, v1), + dot12 = dot(v1, v2); + + // Compute barycentric coordinates + double invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01), + u = (dot11 * dot02 - dot01 * dot12) * invDenom, + v = (dot00 * dot12 - dot01 * dot02) * invDenom; + + // Check if point is in triangle + if (allowBorder) + { + return (u >= 0) && (v >= 0) && (u + v <= 1); + } + else + { + return (u > 0) && (v > 0) && (u + v < 1); + } +} + +/********************************************************/ +/* */ +/* triangleTransformations2DFromCorrespondingPoints */ +/* */ +/********************************************************/ + +/** \brief Create a list of triangles of the dest image together with an affine transform for each. + + For use with \ref piecewiseAffineWarpImage(). + + This function performs a delaunay triangulation on the source image an computes the + affine transform from the dest image triangles to the first one. + + For the delaunay triangulation, LEMON is neccessary. +*/ +template +std::vector +triangleTransformations2DFromCorrespondingPoints(SrcPointIterator s, SrcPointIterator s_end, DestPointIterator d) +{ + using namespace lemon; + typedef dim2::Point Point; + + int point_count = s_end - s; + + //Create a new graph with all the points as nodes + ListGraph g; + ListGraph::NodeMap coords(g); + ListGraph::NodeMap idxs(g); + + //Add nodes and their geometric embedding + for (int i=0; i > triangles = trianglesFromDelaunay(g, idxs); + + //Fill the transformation matrixes + std::vector s_points(3), d_points(3); + std::vector result; + for (auto iter = triangles.begin(); iter!=triangles.end(); iter++) + { + int idx_a = iter->first, + idx_b = iter->second, + idx_c = iter->third; + + s_points[0] = s[idx_a]; s_points[1] = s[idx_b]; s_points[2] = s[idx_c]; + d_points[0] = d[idx_a]; d_points[1] = d[idx_b]; d_points[2] = d[idx_c]; + + result.push_back(TriangleTransformationType(TriangleType(d_points[0],d_points[1],d_points[2]), + affineMatrix2DFromCorrespondingPoints(d_points.begin(), d_points.end(), s_points.begin()))); + } + + return result; +} + +/********************************************************/ +/* */ +/* piecewiseAffineWarpImage */ +/* */ +/********************************************************/ + +/** \brief Warp an image according to a piecewise affine transformation. + + To get more information about the structure of the triangle+matrix + vector, see \ref polynomialMatrix2DFromCorrespondingPoints(). + + \#include \
+ Namespace: vigra + + pass 2D array views: + \code + namespace vigra { + template + void + piecewiseAffineWarpImage(SplineImageView const & src, + MultiArrayView<2, T2, S2> dest, + const std::vector & tri_trans); + } + \endcode + + \deprecatedAPI{piecewiseAffineWarpImage} + + pass arguments explicitly: + \code + namespace vigra { + template + void piecewiseAffineWarpImage(SplineImageView const & src, + DestIterator dul, DestIterator dlr, DestAccessor dest, + const std::vector & tri_trans); + } + \endcode + + use argument objects in conjunction with \ref ArgumentObjectFactories : + \code + namespace vigra { + template + void piecewiseAffineWarpImage(SplineImageView const & src, + triple dest, + const std::vector & tri_trans); + } + \endcode + \deprecatedEnd + */ +doxygen_overloaded_function(template <...> void piecewiseAffineWarpImage) + +template +void piecewiseAffineWarpImage(SplineImageView const & src, + DestIterator dul, DestIterator dlr, DestAccessor dest, + const std::vector & tri_trans) +{ + double w = dlr.x - dul.x; + double h = dlr.y - dul.y; + + for(double y = 0.0; y < h; ++y, ++dul.y) + { + typename DestIterator::row_iterator rd = dul.rowIterator(); + for(double x=0.0; x < w; ++x, ++rd) + { + for(std::vector::const_iterator iter=tri_trans.begin(); iter!=tri_trans.end(); ++iter) + { + if(isInside(PointType(x,y), iter->first)) + { + vigra::Matrix transformation = iter->second; + double sx = transformation(0,0)*x + transformation(0,1)*y + transformation(0,2); + double sy = transformation(1,0)*x + transformation(1,1)*y + transformation(1,2); + + if(src.isInside(sx, sy)) + dest.set(src(sx, sy), rd); + } + } + } + } +} + +template +inline +void piecewiseAffineWarpImage(SplineImageView const & src, + triple dest, + const std::vector & tri_trans) +{ + piecewiseAffineWarpImage(src, dest.first, dest.second, dest.third, tri_trans); +} + + +template +inline +void piecewiseAffineWarpImage(SplineImageView const & src, + MultiArrayView<2, T2, S2> dest, + const std::vector & tri_trans) +{ + piecewiseAffineWarpImage(src, destImageRange(dest), tri_trans); +} + + +//@} + +} // namespace vigra + + +#endif /* VIGRA_PIECEWISEAFFINE_REGISTRATION_HXX */ diff --git a/test/registration/CMakeLists.txt b/test/registration/CMakeLists.txt index 78bb20228..23b38cca4 100644 --- a/test/registration/CMakeLists.txt +++ b/test/registration/CMakeLists.txt @@ -4,11 +4,23 @@ if(FFTW3_FOUND) INCLUDE_DIRECTORIES(${SUPPRESS_WARNINGS} ${FFTW3_INCLUDE_DIR}) ADD_DEFINITIONS(-DHasFFTW3) - VIGRA_ADD_TEST(test_registration test.cxx LIBRARIES ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES} ${THREADING_LIBRARIES} vigraimpex) + if(LEMON_FOUND) + INCLUDE_DIRECTORIES(${SUPPRESS_WARNINGS} ${LEMON_INCLUDE_DIR}) + ADD_DEFINITIONS(-DWITH_LEMON) + VIGRA_ADD_TEST(test_registration test.cxx LIBRARIES ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES} ${THREADING_LIBRARIES} ${LEMON_LIBRARIES} vigraimpex) + else() + VIGRA_ADD_TEST(test_registration test.cxx LIBRARIES ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES} ${THREADING_LIBRARIES} vigraimpex) + endif() else() MESSAGE(STATUS "** WARNING: fftw not found. Fourier-domain registration tests will not be executed") - VIGRA_ADD_TEST(test_registration test.cxx LIBRARIES ${THREADING_LIBRARIES} vigraimpex) + if(LEMON_FOUND) + INCLUDE_DIRECTORIES(${SUPPRESS_WARNINGS} ${LEMON_INCLUDE_DIR}) + ADD_DEFINITIONS(-DWITH_LEMON) + VIGRA_ADD_TEST(test_registration test.cxx LIBRARIES ${THREADING_LIBRARIES} ${LEMON_LIBRARIES} vigraimpex) + else() + VIGRA_ADD_TEST(test_registration test.cxx LIBRARIES ${THREADING_LIBRARIES} vigraimpex) + endif() endif() VIGRA_COPY_TEST_DATA(nuernberg-1991.png nuernberg-1995.png) diff --git a/test/registration/test.cxx b/test/registration/test.cxx index 27fa12029..7610a67bd 100755 --- a/test/registration/test.cxx +++ b/test/registration/test.cxx @@ -47,6 +47,10 @@ #include #include +#ifdef WITH_LEMON +#include +#endif + using namespace vigra; static int pointdata[] = { @@ -809,6 +813,127 @@ struct RadialBasisRegistrationTestSuite } }; +#ifdef WITH_LEMON +struct PiecewiseAffineIdentityTest +{ + std::vector > s_points; + + PiecewiseAffineIdentityTest() + : s_points(srcPoints()) + { + } + + void testInit() + { + /** + * First test: If point sets are equal -> a number of triangles each assigned + * with an identity matrix representation should be the result! + */ + typedef std::vector TTVectorType; + + TTVectorType identities = triangleTransformations2DFromCorrespondingPoints(s_points.begin(), s_points.end(), s_points.begin()); + + /** + * Estimated result: + * + * Triangles... and each with an identity matrix + */ + Matrix reference = linalg::identityMatrix(3); + + for (TTVectorType::const_iterator iter=identities.begin(); iter!=identities.end(); ++iter) + { + shouldEqualToleranceMatrices(iter->second, reference, test_epsilon); + } + } +}; + + + + +struct PiecewiseAffineRegistrationTest +{ + vigra::BImage s_img; + vigra::BImage d_img; + + std::vector > s_points; + std::vector > d_points; + + PiecewiseAffineRegistrationTest() + : s_points(srcPoints()), + d_points(destPoints()) + { + vigra::ImageImportInfo info1("nuernberg-1991.png"); + s_img.resize(info1.width(), info1.height()); + vigra::importImage(info1, destImage(s_img)); + + vigra::ImageImportInfo info2("nuernberg-1995.png"); + d_img.resize(info2.width(), info2.height()); + vigra::importImage(info2, destImage(d_img)); + } + + void testInit() + { + /** + * First test: If point sets are equal -> a number of triangles each assigned + * with an identity matrix representation should be the result! + */ + + typedef std::vector TTVectorType; + + TTVectorType tri_trans = triangleTransformations2DFromCorrespondingPoints(s_points.begin(), s_points.end(), d_points.begin()); + + /** + * Estimated result: + * + * Triangles... and each with a matrix that transforms dest Corners to src corners! + */ + std::vector d_test_points(3); + + for (TTVectorType::const_iterator iter=tri_trans.begin(); iter!=tri_trans.end(); ++iter) + { + d_test_points[0] = iter->first.first; + d_test_points[1] = iter->first.second; + d_test_points[2] = iter->first.third; + + vigra::Matrix transformation = iter->second; + + for(int p=0; p<3; p++) + { + for(size_t i=0; i< d_points.size(); i++) + { + //if current d_point is located at position i of all d_points: + //test if it transforms to s_point of same index + if (d_test_points[p] == d_points[i]) + { + double sx = transformation(0,0)*d_test_points[p][0] + transformation(0,1)*d_test_points[p][1] + transformation(0,2); + double sy = transformation(1,0)*d_test_points[p][0] + transformation(1,1)*d_test_points[p][1] + transformation(1,2); + + shouldEqualTolerance(sx, s_points[i][0], test_epsilon); + shouldEqualTolerance(sy, s_points[i][1], test_epsilon); + } + } + } + } + /** + * visual interpretation by means of the warped image: + */ + BImage temp = d_img; + vigra::piecewiseAffineWarpImage(SplineImageView<2,unsigned char>(srcImageRange(s_img)), destImageRange(temp), tri_trans); + vigra::exportImage(srcImageRange(temp), vigra::ImageExportInfo("res-piecewise_affine.png")); + } +}; + +struct PiecewiseAffineRegistrationTestSuite +: public vigra::test_suite +{ + PiecewiseAffineRegistrationTestSuite() + : vigra::test_suite("PiecewiseAffineRegistrationTestSuite") + { + add( testCase( &PiecewiseAffineIdentityTest::testInit)); + add( testCase( &PiecewiseAffineRegistrationTest::testInit)); + } +}; +#endif /* WITH_LEMON */ struct RegistrationTestCollection @@ -823,6 +948,9 @@ struct RegistrationTestCollection add( new ProjectiveRegistrationTestSuite); add( new PolynomialRegistrationTestSuite); add( new RadialBasisRegistrationTestSuite); +#ifdef WITH_LEMON + add( new PiecewiseAffineRegistrationTestSuite); +#endif } };