From a72bdd9795419e5252d1e7284ef68c71e07ec020 Mon Sep 17 00:00:00 2001 From: Patrick Faion Date: Wed, 19 Feb 2020 11:22:39 +0100 Subject: [PATCH] Remove old cython/cpp references throughout the code base --- .gitignore | 4 - deployment/bundle.bat | 3 - deployment/bundle.sh | 6 - pupil_src/launchables/service.py | 5 - pupil_src/launchables/world.py | 5 - .../shared_cpp/include/DistancePointEllipse.h | 185 --------- .../shared_cpp/include/ceres/CeresUtils.h | 110 ------ .../ceres/EigenQuaternionParameterization.h | 63 ---- .../ceres/Fixed3DNormParametrization.h | 160 -------- pupil_src/shared_cpp/include/common/colors.h | 18 - .../shared_cpp/include/common/constants.h | 19 - pupil_src/shared_cpp/include/common/traits.h | 36 -- pupil_src/shared_cpp/include/common/types.h | 133 ------- .../shared_cpp/include/geometry/Circle.h | 57 --- pupil_src/shared_cpp/include/geometry/Conic.h | 93 ----- .../shared_cpp/include/geometry/Conicoid.h | 119 ------ .../shared_cpp/include/geometry/Ellipse.h | 156 -------- .../shared_cpp/include/geometry/Sphere.h | 56 --- pupil_src/shared_cpp/include/math/distance.h | 185 --------- pupil_src/shared_cpp/include/math/intersect.h | 127 ------- pupil_src/shared_cpp/include/mathHelper.h | 351 ------------------ pupil_src/shared_cpp/include/projection.h | 293 --------------- pupil_src/shared_cpp/include/solve.h | 114 ------ .../shared_modules/cython_methods/.gitignore | 1 - .../shared_modules/cython_methods/__init__.py | 18 - .../shared_modules/cython_methods/build.py | 39 -- .../shared_modules/cython_methods/methods.pyx | 46 --- .../shared_modules/cython_methods/setup.py | 42 --- update_license_header.py | 1 - 29 files changed, 2445 deletions(-) delete mode 100644 pupil_src/shared_cpp/include/DistancePointEllipse.h delete mode 100644 pupil_src/shared_cpp/include/ceres/CeresUtils.h delete mode 100644 pupil_src/shared_cpp/include/ceres/EigenQuaternionParameterization.h delete mode 100644 pupil_src/shared_cpp/include/ceres/Fixed3DNormParametrization.h delete mode 100644 pupil_src/shared_cpp/include/common/colors.h delete mode 100644 pupil_src/shared_cpp/include/common/constants.h delete mode 100644 pupil_src/shared_cpp/include/common/traits.h delete mode 100644 pupil_src/shared_cpp/include/common/types.h delete mode 100644 pupil_src/shared_cpp/include/geometry/Circle.h delete mode 100644 pupil_src/shared_cpp/include/geometry/Conic.h delete mode 100644 pupil_src/shared_cpp/include/geometry/Conicoid.h delete mode 100644 pupil_src/shared_cpp/include/geometry/Ellipse.h delete mode 100644 pupil_src/shared_cpp/include/geometry/Sphere.h delete mode 100644 pupil_src/shared_cpp/include/math/distance.h delete mode 100644 pupil_src/shared_cpp/include/math/intersect.h delete mode 100644 pupil_src/shared_cpp/include/mathHelper.h delete mode 100644 pupil_src/shared_cpp/include/projection.h delete mode 100644 pupil_src/shared_cpp/include/solve.h delete mode 100644 pupil_src/shared_modules/cython_methods/.gitignore delete mode 100644 pupil_src/shared_modules/cython_methods/__init__.py delete mode 100644 pupil_src/shared_modules/cython_methods/build.py delete mode 100644 pupil_src/shared_modules/cython_methods/methods.pyx delete mode 100644 pupil_src/shared_modules/cython_methods/setup.py diff --git a/.gitignore b/.gitignore index b29aefeb60..4a6465142b 100644 --- a/.gitignore +++ b/.gitignore @@ -39,10 +39,6 @@ settings/ .idea #directory for test file will be ignored Pupil_Test_Files -detector_2d.cpp -detector_3d.cpp -CircleGoodnessTest.cpp -SphereCircleTest.cpp *.mkv *.egg-info hmd_cal_data diff --git a/deployment/bundle.bat b/deployment/bundle.bat index 1f7d3b5db5..e2690e65c8 100644 --- a/deployment/bundle.bat +++ b/deployment/bundle.bat @@ -35,9 +35,6 @@ if not exist %release_dir% ( set PATH=%PATH%;C:\Python36\Lib\site-packages\scipy\.libs set PATH=%PATH%;C:\Python36\Lib\site-packages\zmq -python ..\pupil_src\shared_modules\cython_methods\build.py -python ..\pupil_src\shared_modules\calibration_routines\optimization_calibration\build.py - call :Bundle capture %current_tag% call :Bundle service %current_tag% call :Bundle player %current_tag% diff --git a/deployment/bundle.sh b/deployment/bundle.sh index c94fb97985..20fa5b87d3 100755 --- a/deployment/bundle.sh +++ b/deployment/bundle.sh @@ -13,12 +13,6 @@ fi echo "release_dir: ${release_dir}" mkdir ${release_dir} -printf "\n##########\nBuilding cython modules\n##########\n\n" -python3 ../pupil_src/shared_modules/cython_methods/build.py - -printf "\n##########\nBuilding calibration methods\n##########\n\n" -python3 ../pupil_src/shared_modules/calibration_routines/optimization_calibration/build.py - # bundle Pupil Capture printf "\n##########\nBundling Pupil Capture\n##########\n\n" cd deploy_capture diff --git a/pupil_src/launchables/service.py b/pupil_src/launchables/service.py index 8df3181b20..6d9d7e386e 100644 --- a/pupil_src/launchables/service.py +++ b/pupil_src/launchables/service.py @@ -93,11 +93,6 @@ def stop_eye_process(eye_id): import audio from uvc import get_time_monotonic - # trigger pupil detector cpp build: - import pupil_detectors - - del pupil_detectors - # Plug-ins from plugin import Plugin, Plugin_List, import_runtime_plugins from calibration_routines import calibration_plugins, gaze_mapping_plugins diff --git a/pupil_src/launchables/world.py b/pupil_src/launchables/world.py index b56407cee1..e7fb49779b 100644 --- a/pupil_src/launchables/world.py +++ b/pupil_src/launchables/world.py @@ -134,11 +134,6 @@ def set_detection_mapping_mode(new_mode): import audio - # trigger pupil detector cpp build: - import pupil_detectors - - del pupil_detectors - # Plug-ins from plugin import ( Plugin, diff --git a/pupil_src/shared_cpp/include/DistancePointEllipse.h b/pupil_src/shared_cpp/include/DistancePointEllipse.h deleted file mode 100644 index fe75babd94..0000000000 --- a/pupil_src/shared_cpp/include/DistancePointEllipse.h +++ /dev/null @@ -1,185 +0,0 @@ -// Geometric Tools, LLC -// Copyright (c) 1998-2014 -// Distributed under the Boost Software License, Version 1.0. -// http://www.boost.org/LICENSE_1_0.txt -// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt -// -// File Version: 5.0.3 (2013/01/03) -// -// Modified by Lech Swirski 2013 - -#ifndef DistancePointEllipse_h__ -#define DistancePointEllipse_h__ - -#include -#include "geometry/Ellipse.h" - -namespace singleeyefitter { - - //---------------------------------------------------------------------------- - // The ellipse is (x0/a)^2 + (x1/b)^2 = 1 with a >= b. The query point is - // (p0,p1) with p0 >= 0 and p1 >= 0. The function returns the distance from - // the query point to the ellipse. It also computes the ellipse point (x0,x1) - // in the first quadrant that is closest to (p0,p1). - //---------------------------------------------------------------------------- - template - Real DistancePointEllipseSpecial(Real a, Real b, const Array& p, Eigen::Matrix& x) - { - Real distance; - - if (p.y() > Real(0)) { - if (p.x() > Real(0)) { - // Bisect to compute the root of F(t) for t >= -e1*e1. - Eigen::Array esqr(a * a, b * b); - Eigen::Array ep(a * p.x(), b * p.y()); - Real t0 = -esqr.y() + ep.y(); - Real t1 = -esqr.y() + ep.matrix().norm(); - Real t = t0; - const int imax = 2 * std::numeric_limits::max_exponent; - - for (int i = 0; i < imax; ++i) { - t = Real(0.5) * (t0 + t1); - - if (t == t0 || t == t1) { - break; - } - - Real r[2] = { ep.x() / (t + esqr[0]), ep.y() / (t + esqr[1]) }; - Real f = r[0] * r[0] + r[1] * r[1] - Real(1); - - if (f > Real(0)) { - t0 = t; - - } else if (f < Real(0)) { - t1 = t; - - } else { - break; - } - } - - x = esqr * p / (t + esqr); - distance = (x - p.matrix()).norm(); - - } else { // y0 == 0 - x[0] = (Real) 0; - x[1] = b; - distance = fabs(p.y() - b); - } - - } else { // y1 == 0 - Real denom0 = a * a - b * b; - Real e0y0 = a * p.x(); - - if (e0y0 < denom0) { - // y0 is inside the subinterval. - Real x0de0 = e0y0 / denom0; - Real x0de0sqr = x0de0 * x0de0; - x[0] = a * x0de0; - x[1] = b * sqrt(fabs(Real(1) - x0de0sqr)); - Real d0 = x[0] - p.x(); - distance = sqrt(d0 * d0 + x[1] * x[1]); - - } else { - // y0 is outside the subinterval. The closest ellipse point has - // x1 == 0 and is on the domain-boundary interval (x0/e0)^2 = 1. - x[0] = a; - x[1] = Real(0); - distance = fabs(p.x() - a); - } - } - - return distance; - } - //---------------------------------------------------------------------------- - // The ellipse is (x0/e0)^2 + (x1/e1)^2 = 1. The query point is (y0,y1). - // The function returns the distance from the query point to the ellipse. - // It also computes the ellipse point (x0,x1) that is closest to (y0,y1). - //---------------------------------------------------------------------------- - template - Real DistancePointEllipse(const Real e[2], const Real y[2], Real x[2]) - { - // Determine reflections for y to the first quadrant. - bool reflect[2]; - int i, j; - - for (i = 0; i < 2; ++i) { - reflect[i] = (y[i] < (Real) 0); - } - - // Determine the axis order for decreasing extents. - int permute[2]; - - if (e[0] < e[1]) { - permute[0] = 1; permute[1] = 0; - - } else { - permute[0] = 0; permute[1] = 1; - } - - int invpermute[2]; - - for (i = 0; i < 2; ++i) { - invpermute[permute[i]] = i; - } - - Real locE[2], locY[2]; - - for (i = 0; i < 2; ++i) { - j = permute[i]; - locE[i] = e[j]; - locY[i] = y[j]; - - if (reflect[j]) { - locY[i] = -locY[i]; - } - } - - Real locX[2]; - Real distance = DistancePointEllipseSpecial(locE, locY, locX); - - // Restore the axis order and reflections. - for (i = 0; i < 2; ++i) { - j = invpermute[i]; - - if (reflect[j]) { - locX[j] = -locX[j]; - } - - x[i] = locX[j]; - } - - return distance; - } - //---------------------------------------------------------------------------- - - template - Scalar DistancePointEllipse(const singleeyefitter::Ellipse2D& ellipse, Scalar x, Scalar y) - { - Eigen::Matrix A; - A << cos(ellipse.angle), sin(ellipse.angle), - -sin(ellipse.angle), cos(ellipse.angle); - Eigen::Matrix p(x - ellipse.center.x(), y - ellipse.center.y()); - Eigen::Matrix Ap = A * p; - // Flip signs to make sure Ap is in the positive quadrant - Eigen::Matrix Ap_pos = Ap; - - for (int i = 0; i < 2; ++i) { - if (Ap[i] < 0) { Ap_pos[i] = -Ap_pos[i]; } - } - - assert(ellipse.major_radius > ellipse.minor_radius); - Eigen::Matrix el_x; - auto distance = DistancePointEllipseSpecial(ellipse.major_radius, ellipse.minor_radius, Ap_pos.array(), el_x); - - // Flip signs back - for (int i = 0; i < 2; ++i) { - if (Ap[i] < 0) { el_x[i] = -el_x[i]; } - } - - return distance; - } - -} - -#endif // DistancePointEllipse_h__ diff --git a/pupil_src/shared_cpp/include/ceres/CeresUtils.h b/pupil_src/shared_cpp/include/ceres/CeresUtils.h deleted file mode 100644 index 5928a50660..0000000000 --- a/pupil_src/shared_cpp/include/ceres/CeresUtils.h +++ /dev/null @@ -1,110 +0,0 @@ - -#ifndef CERESUTILS_H__ -#define CERESUTILS_H__ - -#include -#include - -namespace pupillabs { - - -template inline -void EigenQuaternionToScaledRotation(const T q[4], T R[3 * 3]) { - EigenQuaternionToScaledRotation(q, RowMajorAdapter3x3(R)); -} - -template inline -void EigenQuaternionToScaledRotation(const T q[4], - const ceres::MatrixAdapter& R) { - // Make convenient names for elements of q. - T a = q[3]; - T b = q[0]; - T c = q[1]; - T d = q[2]; - // This is not to eliminate common sub-expression, but to - // make the lines shorter so that they fit in 80 columns! - T aa = a * a; - T ab = a * b; - T ac = a * c; - T ad = a * d; - T bb = b * b; - T bc = b * c; - T bd = b * d; - T cc = c * c; - T cd = c * d; - T dd = d * d; - - R(0, 0) = aa + bb - cc - dd; R(0, 1) = T(2) * (bc - ad); R(0, 2) = T(2) * (ac + bd); // NOLINT - R(1, 0) = T(2) * (ad + bc); R(1, 1) = aa - bb + cc - dd; R(1, 2) = T(2) * (cd - ab); // NOLINT - R(2, 0) = T(2) * (bd - ac); R(2, 1) = T(2) * (ab + cd); R(2, 2) = aa - bb - cc + dd; // NOLINT -} - -template inline -void EigenQuaternionToRotation(const T q[4], T R[3 * 3]) { - EigenQuaternionToRotation(q, RowMajorAdapter3x3(R)); -} - -template inline -void EigenQuaternionToRotation(const T q[4], - const ceres::MatrixAdapter& R) { - EigenQuaternionToScaledRotation(q, R); - - T normalizer = q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]; - CHECK_NE(normalizer, T(0)); - normalizer = T(1) / normalizer; - - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - R(i, j) *= normalizer; - } - } -} - -template inline -void EigenUnitQuaternionRotatePoint(const T q[4], const T pt[3], T result[3]) { - const T t2 = q[3] * q[0]; - const T t3 = q[3] * q[1]; - const T t4 = q[3] * q[2]; - const T t5 = -q[0] * q[0]; - const T t6 = q[0] * q[1]; - const T t7 = q[0] * q[2]; - const T t8 = -q[1] * q[1]; - const T t9 = q[1] * q[2]; - const T t1 = -q[2] * q[2]; - result[0] = T(2) * ((t8 + t1) * pt[0] + (t6 - t4) * pt[1] + (t3 + t7) * pt[2]) + pt[0]; // NOLINT - result[1] = T(2) * ((t4 + t6) * pt[0] + (t5 + t1) * pt[1] + (t9 - t2) * pt[2]) + pt[1]; // NOLINT - result[2] = T(2) * ((t7 - t3) * pt[0] + (t2 + t9) * pt[1] + (t5 + t8) * pt[2]) + pt[2]; // NOLINT -} - -template inline -void EigenQuaternionRotatePoint(const T q[4], const T pt[3], T result[3]) { - // 'scale' is 1 / norm(q). - const T scale = T(1) / sqrt(q[0] * q[0] + - q[1] * q[1] + - q[2] * q[2] + - q[3] * q[3]); - - // Make unit-norm version of q. - const T unit[4] = { - scale * q[0], - scale * q[1], - scale * q[2], - scale * q[3], - }; - - EigenUnitQuaternionRotatePoint(unit, pt, result); -} - -template inline -void EigenQuaternionProduct(const T z[4], const T w[4], T zw[4]) { - zw[0] = z[0] * w[3] + z[1] * w[2] - z[2] * w[1] + z[3] * w[0]; - zw[1] = - z[0] * w[2] + z[1] * w[3] + z[2] * w[0] + z[3] * w[1]; - zw[2] = z[0] * w[1] - z[1] * w[0] + z[2] * w[3] + z[3] * w[2]; - zw[3] = - z[0] * w[0] - z[1] * w[1] - z[2] * w[2] + z[3] * w[3]; -} - - -} // pupillabs - - -#endif /* end of include guard: CERESUTILS_H__ */ diff --git a/pupil_src/shared_cpp/include/ceres/EigenQuaternionParameterization.h b/pupil_src/shared_cpp/include/ceres/EigenQuaternionParameterization.h deleted file mode 100644 index f390a30f42..0000000000 --- a/pupil_src/shared_cpp/include/ceres/EigenQuaternionParameterization.h +++ /dev/null @@ -1,63 +0,0 @@ - -/* - Created by Lloyd Hughes on 2014/04/11. - Copyright (c) 2014 Lloyd Hughes. All rights reserved. - hughes.lloyd@gmail.com -*/ - - -#ifndef EIGENQUATERNIONPARAMETERIZATION_H__ -#define EIGENQUATERNIONPARAMETERIZATION_H__ - - -#include -#include - - -namespace pupillabs { - -// Plus(x, delta) = [cos(|delta|), sin(|delta|) delta / |delta|] * x -// with * being the quaternion multiplication operator. Here we assume -// that the first element of the quaternion vector is the real (cos -// theta) part. -class EigenQuaternionParameterization : public ceres::LocalParameterization { -public: - virtual ~EigenQuaternionParameterization() {} - - virtual bool Plus(const double* x_raw, const double* delta_raw, double* x_plus_delta_raw) const { - const Eigen::Map x(x_raw); - const Eigen::Map delta(delta_raw); - - Eigen::Map x_plus_delta(x_plus_delta_raw); - - const double delta_norm = delta.norm(); - if ( delta_norm > 0.0 ){ - const double sin_delta_by_delta = sin(delta_norm) / delta_norm; - Eigen::Quaterniond tmp( cos(delta_norm), sin_delta_by_delta*delta[0], sin_delta_by_delta*delta[1], sin_delta_by_delta*delta[2] ); - - x_plus_delta = tmp*x; - } - else { - x_plus_delta = x; - } - return true; - } - - virtual bool ComputeJacobian(const double* x, double* jacobian) const { - jacobian[0] = x[3]; jacobian[1] = x[2]; jacobian[2] = -x[1]; // NOLINT x - jacobian[3] = -x[2]; jacobian[4] = x[3]; jacobian[5] = x[0]; // NOLINT y - jacobian[6] = x[1]; jacobian[7] = -x[0]; jacobian[8] = x[3]; // NOLINT z - jacobian[9] = -x[0]; jacobian[10] = -x[1]; jacobian[11] = -x[2]; // NOLINT w - return true; - } - - virtual int GlobalSize() const { return 4; } - virtual int LocalSize() const { return 3; } - -}; - - -} // pupillabs - - -#endif /* end of include guard: EIGENQUATERNIONPARAMETERIZATION_H__ */ diff --git a/pupil_src/shared_cpp/include/ceres/Fixed3DNormParametrization.h b/pupil_src/shared_cpp/include/ceres/Fixed3DNormParametrization.h deleted file mode 100644 index 3b3ac502ba..0000000000 --- a/pupil_src/shared_cpp/include/ceres/Fixed3DNormParametrization.h +++ /dev/null @@ -1,160 +0,0 @@ -/* - -Copyright (C) 2014, University of Oulu, all rights reserved. -Copyright (C) 2014, NVIDIA Corporation, all rights reserved. -Third party copyrights are property of their respective owners. - -Redistribution and use in source and binary forms, with or without modification, -are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - * Neither the name of UNIVERSITY OF OULU, NVIDIA CORPORATION nor the names of its - contributors may be used to endorse or promote products derived - from this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY -EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR -CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, -EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, -PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR -PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY -OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -*/ - -#ifndef FIXED3DNORMPARAMETRIZATION_H__ -#define FIXED3DNORMPARAMETRIZATION_H__ - - -#include -#include - -namespace pupillabs -{ - -/** - * @brief A parameterization class that is used for CERES solver. It parametrizes a 3D vector (like a translation) with two components, keeping the L2 norm fixed - */ -class Fixed3DNormParametrization: public ceres::LocalParameterization -{ -public: - Fixed3DNormParametrization(double norm) - : mFixedNorm(norm) - { - } - virtual ~Fixed3DNormParametrization() - { - } - - virtual int GlobalSize() const - { - return 3; - } - virtual int LocalSize() const - { - return 2; - } - - // Calculates two vectors that are orthogonal to X. - // It first picks a non-colinear point C then basis1=(X-C) x C and basis2=X x basis1 - static void GetBasis(const double *x, double *basis1, double *basis2) - { - const double kThreshold = 0.1; - - //Check that the point we use is not colinear with x - if (x[1] > kThreshold || x[1] < -kThreshold || x[2] > kThreshold || x[2] < -kThreshold) - { - //Use C=[1,0,0] - basis1[0] = 0; - basis1[1] = x[2]; - basis1[2] = -x[1]; - - basis2[0] = -(x[1] * x[1] + x[2] * x[2]); - basis2[1] = x[0] * x[1]; - basis2[2] = x[0] * x[2]; - } - else - { - //Use C=[0,1,0] - basis1[0] = -x[2]; - basis1[1] = 0; - basis1[2] = x[0]; - - basis2[0] = x[0] * x[1]; - basis2[1] = -(x[0] * x[0] + x[2] * x[2]); - basis2[2] = x[1] * x[2]; - } - double norm; - norm = sqrt(basis1[0] * basis1[0] + basis1[1] * basis1[1] + basis1[2] * basis1[2]); - basis1[0] /= norm; - basis1[1] /= norm; - basis1[2] /= norm; - - norm = sqrt(basis2[0] * basis2[0] + basis2[1] * basis2[1] + basis2[2] * basis2[2]); - basis2[0] /= norm; - basis2[1] /= norm; - basis2[2] /= norm; - - } - - virtual bool Plus(const double *x, const double *delta, double *x_plus_delta) const - { - double basis1[3]; - double basis2[3]; - - //Translation is constrained - GetBasis(x, basis1, basis2); - - x_plus_delta[0] = x[0] + delta[0] * basis1[0] + delta[1] * basis2[0]; - x_plus_delta[1] = x[1] + delta[0] * basis1[1] + delta[1] * basis2[1]; - x_plus_delta[2] = x[2] + delta[0] * basis1[2] + delta[1] * basis2[2]; - - double norm = sqrt( - x_plus_delta[0] * x_plus_delta[0] + x_plus_delta[1] * x_plus_delta[1] + x_plus_delta[2] * x_plus_delta[2]); - double factor = mFixedNorm / norm; - x_plus_delta[0] *= factor; - x_plus_delta[1] *= factor; - x_plus_delta[2] *= factor; - - return true; - } - - virtual bool ComputeJacobian(const double *x, double *jacobian) const - { - typedef Eigen::Matrix Matrix32d; - Matrix32d &jacobian_ = *(Matrix32d *)jacobian; - double basis1[3]; - double basis2[3]; - - //Translation is special - GetBasis(x, basis1, basis2); - - jacobian_(0, 0) = basis1[0]; - jacobian_(1, 0) = basis1[1]; - jacobian_(2, 0) = basis1[2]; - - jacobian_(0, 1) = basis2[0]; - jacobian_(1, 1) = basis2[1]; - jacobian_(2, 1) = basis2[2]; - return true; - } - - -protected: - const double mFixedNorm; -}; - -} //namespace pupillabs - - -#endif /* end of include guard: FIXED3DNORMPARAMETRIZATION_H__ */ diff --git a/pupil_src/shared_cpp/include/common/colors.h b/pupil_src/shared_cpp/include/common/colors.h deleted file mode 100644 index e00d760234..0000000000 --- a/pupil_src/shared_cpp/include/common/colors.h +++ /dev/null @@ -1,18 +0,0 @@ -#ifndef singleeyefitter_colors_h__ -#define singleeyefitter_colors_h__ - -#include - - -namespace singleeyefitter { - - const cv::Scalar_ mRed_color = {0, 0, 255}; - const cv::Scalar_ mGreen_color = {0, 255, 0}; - const cv::Scalar_ mBlue_color = {255, 0, 0}; - const cv::Scalar_ mRoyalBlue_color = {255, 100, 100}; - const cv::Scalar_ mYellow_color = {255, 255, 0}; - const cv::Scalar_ mWhite_color = {255, 255, 255}; - -} // singleeyefitter namespace - -#endif //singleeyefitter_types_h__ diff --git a/pupil_src/shared_cpp/include/common/constants.h b/pupil_src/shared_cpp/include/common/constants.h deleted file mode 100644 index 54bd405097..0000000000 --- a/pupil_src/shared_cpp/include/common/constants.h +++ /dev/null @@ -1,19 +0,0 @@ - -#ifndef CONSTANTS_H__ -#define CONSTANTS_H__ - -#include -namespace singleeyefitter { - - namespace constants { - const double PI = std::acos(-1); - const double TWO_PI = 2.0 * std::acos(-1); - - } // constants - - -} // singleeyefitter - - -#endif /* end of include guard: CONSTANTS_H__ */ - diff --git a/pupil_src/shared_cpp/include/common/traits.h b/pupil_src/shared_cpp/include/common/traits.h deleted file mode 100644 index 672ed9fbe2..0000000000 --- a/pupil_src/shared_cpp/include/common/traits.h +++ /dev/null @@ -1,36 +0,0 @@ -#ifndef singleeyefitter_traits_h__ -#define singleeyefitter_traits_h__ - - -#include - -namespace singleeyefitter { - - struct scalar_tag {}; - struct ceres_jet_tag {}; - - template - struct ad_traits; - - template - struct ad_traits::value >::type > { - typedef scalar_tag ad_tag; - typedef T scalar; - static inline scalar value(const T& x) { return x; } - }; - - template - struct ad_traits<::ceres::Jet> { - typedef ceres_jet_tag ad_tag; - typedef T scalar; - static inline scalar get(const ::ceres::Jet& x) { return x.a; } - }; - - template - struct ad_traits < T, typename std::enable_if < !std::is_same::type>::value >::type > - : public ad_traits::type> { - }; -} - - -#endif //singleeyefitter_traits_h__ diff --git a/pupil_src/shared_cpp/include/common/types.h b/pupil_src/shared_cpp/include/common/types.h deleted file mode 100644 index 2e74c1f074..0000000000 --- a/pupil_src/shared_cpp/include/common/types.h +++ /dev/null @@ -1,133 +0,0 @@ - -#ifndef singleeyefitter_types_h__ -#define singleeyefitter_types_h__ - -#include "geometry/Ellipse.h" -#include "geometry/Circle.h" -#include "geometry/Sphere.h" -#include "projection.h" - -#include -#include -#include - -#include - - -namespace singleeyefitter { - - - //######## 2D Detector ############ - typedef std::vector > Contours_2D; - typedef std::vector Contour_2D; - typedef std::vector Edges2D; - typedef std::vector ContourIndices; - typedef Ellipse2D Ellipse; - - //######## 3D Detector ############ - - typedef Eigen::Matrix Vector2; - typedef Eigen::Matrix Vector3; - typedef Eigen::ParametrizedLine Line; - typedef Eigen::ParametrizedLine Line3; - typedef Circle3D Circle; - typedef size_t Index; - - typedef std::vector Contour3D; - typedef std::vector Edges3D; - typedef std::vector> Contours3D; - - struct ConfidenceValue{ - ConfidenceValue(double v,double c) - { - value = v; - confidence = c; - }; - ConfidenceValue() - { - value = 0; - confidence = 0; - }; - double value; - double confidence; - }; - - // general time - typedef std::chrono::steady_clock Clock; - - - // every coordinates are relative to the roi - struct Detector2DResult { - double confidence = 0.0 ; - Ellipse ellipse = Ellipse::Null; - Edges2D final_edges; // edges used to fit the final ellipse in 2D - Edges2D raw_edges; - cv::Rect current_roi; // contains the roi for this results - double timestamp = 0.0; - int image_width = 0; - int image_height = 0; - - }; - - struct ModelDebugProperties{ - Sphere sphere; - Sphere initialSphere; - std::vector binPositions; - double maturity; - double solverFit; - double confidence; - double performance; - double performanceGradient; - int modelID; - double birthTimestamp; - }; - - struct Detector3DResult { - double confidence = 0.0 ; - Circle circle = Circle::Null; - Ellipse ellipse = Ellipse::Null; // the circle projected back to 2D - Sphere sphere = Sphere::Null; - Ellipse projectedSphere = Ellipse::Null; // the sphere projected back to 2D - double timestamp; - int modelID = 0; - double modelBirthTimestamp = 0.0; - double modelConfidence = 0.0; - //-------- For visualization ---------------- - // just valid if we want it for visualization - Edges3D edges; - Circle predictedCircle = Circle::Null; - std::vector models; - }; - - // use a struct for all properties and pass it to detect method every time we call it. - // Thus we don't need to keep track if GUI is updated and cython handles conversion from Dict to struct - struct Detector2DProperties { - int intensity_range; - int blur_size; - float canny_treshold; - float canny_ration; - int canny_aperture; - int pupil_size_max; - int pupil_size_min; - float strong_perimeter_ratio_range_min; - float strong_perimeter_ratio_range_max; - float strong_area_ratio_range_min; - float strong_area_ratio_range_max; - int contour_size_min; - float ellipse_roundness_ratio; - float initial_ellipse_fit_treshhold; - float final_perimeter_ratio_range_min; - float final_perimeter_ratio_range_max; - float ellipse_true_support_min_dist; - float support_pixel_ratio_exponent; - - }; - - struct Detector3DProperties { - float model_sensitivity; - bool model_is_frozen; - }; - -} // singleeyefitter namespace - -#endif //singleeyefitter_types_h__ diff --git a/pupil_src/shared_cpp/include/geometry/Circle.h b/pupil_src/shared_cpp/include/geometry/Circle.h deleted file mode 100644 index 41da5e9e5a..0000000000 --- a/pupil_src/shared_cpp/include/geometry/Circle.h +++ /dev/null @@ -1,57 +0,0 @@ -#ifndef _CIRCLE_H_ -#define _CIRCLE_H_ - -#include - -namespace singleeyefitter { - - template - class Circle3D { - public: - typedef T Scalar; - typedef Eigen::Matrix Vector; - - Vector center, normal; - Scalar radius; - - Circle3D() : center(0, 0, 0), normal(0, 0, 0), radius(0) - { - } - Circle3D(Vector center, Vector normal, Scalar radius) - : center(std::move(center)), normal(std::move(normal)), radius(std::move(radius)) - { - } - - static const Circle3D Null; - - }; - - template - const Circle3D Circle3D::Null = Circle3D(); - - template - bool operator== (const Circle3D& s1, const Circle3D& s2) - { - return s1.center == s2.center - && s1.normal == s2.normal - && s1.radius == s2.radius; - } - template - bool operator!= (const Circle3D& s1, const Circle3D& s2) - { - return s1.center != s2.center - || s1.normal != s2.normal - || s1.radius != s2.radius; - } - - template - std::ostream& operator<< (std::ostream& os, const Circle3D& circle) - { - return os << "Circle { center: (" << circle.center[0] << "," << circle.center[1] << "," << circle.center[2] << "), " - "normal: (" << circle.normal[0] << "," << circle.normal[1] << "," << circle.normal[2] << "), " - "radius: " << circle.radius << " }"; - } - -} - -#endif//_CIRCLE_H_ diff --git a/pupil_src/shared_cpp/include/geometry/Conic.h b/pupil_src/shared_cpp/include/geometry/Conic.h deleted file mode 100644 index ee9d8f2b52..0000000000 --- a/pupil_src/shared_cpp/include/geometry/Conic.h +++ /dev/null @@ -1,93 +0,0 @@ -#ifndef _CONIC_H_ -#define _CONIC_H_ - -#include -#include "mathHelper.h" - -namespace singleeyefitter { - - template - class Ellipse2D; - - template - class Conic { - public: - typedef T Scalar; - - Scalar A, B, C, D, E, F; - - Conic(Scalar A, Scalar B, Scalar C, Scalar D, Scalar E, Scalar F) - : A(A), B(B), C(C), D(D), E(E), F(F) - { - } - - template - explicit Conic(const Ellipse2D& ellipse) - { - using std::sin; - using std::cos; - using singleeyefitter::math::sq; - auto ax = cos(ellipse.angle); - auto ay = sin(ellipse.angle); - auto a2 = sq(ellipse.major_radius); - auto b2 = sq(ellipse.minor_radius); - A = ax * ax / a2 + ay * ay / b2; - B = 2 * ax * ay / a2 - 2 * ax * ay / b2; - C = ay * ay / a2 + ax * ax / b2; - D = (-2 * ax * ay * ellipse.center[1] - 2 * ax * ax * ellipse.center[0]) / a2 - + (2 * ax * ay * ellipse.center[1] - 2 * ay * ay * ellipse.center[0]) / b2; - E = (-2 * ax * ay * ellipse.center[0] - 2 * ay * ay * ellipse.center[1]) / a2 - + (2 * ax * ay * ellipse.center[0] - 2 * ax * ax * ellipse.center[1]) / b2; - F = (2 * ax * ay * ellipse.center[0] * ellipse.center[1] + ax * ax * ellipse.center[0] * ellipse.center[0] + ay * ay * ellipse.center[1] * ellipse.center[1]) / a2 - + (-2 * ax * ay * ellipse.center[0] * ellipse.center[1] + ay * ay * ellipse.center[0] * ellipse.center[0] + ax * ax * ellipse.center[1] * ellipse.center[1]) / b2 - - 1; - } - - Scalar operator()(Scalar x, Scalar y) const - { - return A * x * x + B * x * y + C * y * y + D * x + E * y + F; - } - - template - Conic transformed(const Eigen::MatrixBase& a, const Eigen::MatrixBase& t) const - { - static_assert(ADerived::RowsAtCompileTime == 2 && ADerived::ColsAtCompileTime == 2, "Affine transform must be 2x2 matrix"); - static_assert(TDerived::IsVectorAtCompileTime && TDerived::SizeAtCompileTime == 2, "Translation must be 2 element vector"); - // We map x,y to a new space using - // [x y] -> affine*[x y] + translation - // - // Using a for affine and t for translation: - // x -> a_00*x + a01*y + t0 - // y -> a_10*x + a11*y + t1 - // - // So - // Ax^2 + Bxy + Cy^2 + Dx + Ey + F - // becomes - // A(a_00*x + a01*y + t0)(a_00*x + a01*y + t0) - // + B(a_00*x + a01*y + t0)(a_10*x + a11*y + t1) - // + C(a_10*x + a11*y + t1)(a_10*x + a11*y + t1) - // + D(a_00*x + a01*y + t0) - // + E(a_10*x + a11*y + t1) - // + F - // - // Collecting terms gives: - return Conic( - A * sq(a(0, 0)) + B * a(0, 0) * a(1, 0) + C * sq(a(1, 0)), - 2 * A * a(0, 0) * a(0, 1) + B * a(0, 0) * a(1, 1) + B * a(0, 1) * a(1, 0) + 2 * C * a(1, 0) * a(1, 1), - A * sq(a(0, 1)) + B * a(0, 1) * a(1, 1) + C * sq(a(1, 1)), - 2 * A * a(0, 0) * t(0) + B * a(0, 0) * t(1) + B * a(1, 0) * t(0) + 2 * C * a(1, 0) * t(1) + D * a(0, 0) + E * a(1, 0), - 2 * A * a(0, 1) * t(0) + B * a(0, 1) * t(1) + B * a(1, 1) * t(0) + 2 * C * a(1, 1) * t(1) + D * a(0, 1) + E * a(1, 1), - A * sq(t(0)) + B * t(0) * t(1) + C * sq(t(1)) + D * t(0) + E * t(1) + F - ); - } - }; - - template - std::ostream& operator<< (std::ostream& os, const Conic& conic) - { - return os << "Conic { " << conic.A << "x^2 + " << conic.B << "xy + " << conic.C << "y^2 + " << conic.D << "x + " << conic.E << "y + " << conic.F << " = 0 } "; - } - -} - -#endif diff --git a/pupil_src/shared_cpp/include/geometry/Conicoid.h b/pupil_src/shared_cpp/include/geometry/Conicoid.h deleted file mode 100644 index 6a906f52af..0000000000 --- a/pupil_src/shared_cpp/include/geometry/Conicoid.h +++ /dev/null @@ -1,119 +0,0 @@ -#ifndef _CONICOID_H_ -#define _CONICOID_H_ - -#include - -namespace singleeyefitter { - - template - class Conic; - - // Conicoid (quartic surface) of the form: - // Ax^2 + By^2 + Cz^2 + 2Fyz + 2Gzx + 2Hxy + 2Ux + 2Vy + 2Wz + D = 0 - template - class Conicoid { - public: - typedef T Scalar; - - Scalar A, B, C, F, G, H, U, V, W, D; - - Conicoid(Scalar A, Scalar B, Scalar C, Scalar F, Scalar G, Scalar H, Scalar D) - : A(A), B(B), C(C), F(F), G(G), H(H), U(U), V(V), W(W), D(D) - { - } - - template - explicit Conicoid(const Conic& conic, const Eigen::MatrixBase& vertex) - { - static_assert(Derived::IsVectorAtCompileTime && Derived::SizeAtCompileTime == 3, "Cone vertex requires 3 element vector as vector type"); - using math::sq; - // Finds conicoid with given conic base and vertex - // Assumes conic is on the plane z = 0 - auto alpha = vertex[0]; - auto beta = vertex[1]; - auto gamma = vertex[2]; - A = sq(gamma) * conic.A; - B = sq(gamma) * conic.C; - C = conic.A * sq(alpha) + conic.B * alpha * beta + conic.C * sq(beta) + conic.D * alpha + conic.E * beta + conic.F; - F = -gamma * (conic.C * beta + conic.B / 2 * alpha + conic.E / 2); - G = -gamma * (conic.B / 2 * beta + conic.A * alpha + conic.D / 2); - H = sq(gamma) * conic.B / 2; - U = sq(gamma) * conic.D / 2; - V = sq(gamma) * conic.E / 2; - W = -gamma * (conic.E / 2 * beta + conic.D / 2 * alpha + conic.F); - D = sq(gamma) * conic.F; - } - - Scalar operator()(Scalar x, Scalar y, Scalar z) const - { - return A * sq(x) + B * sq(y) + C * sq(z) + 2 * F * y * z + 2 * G * x * z + 2 * H * x * y + 2 * U * x + 2 * V * y + 2 * W * z + D; - } - - Conic intersectZ(Scalar z = Scalar(0)) const - { - // Finds conic at given z intersection - // Ax^2 + By^2 + Cz^2 + 2Fyz + 2Gzx + 2Hxy + 2Ux + 2Vy + 2Wz + D = 0 - // becomes - // Ax^2 + Bxy + Cy^2 + Fx + Ey + D = 0 - return Conic(A, - 2 * H, - B, - 2 * G * z + 2 * U, - 2 * F * z + 2 * V, - C * sq(z) + 2 * W * z + D); - } - - template - Conicoid transformed(const Eigen::MatrixBase& a, const Eigen::MatrixBase& t) const - { - static_assert(ADerived::RowsAtCompileTime == 3 && ADerived::ColsAtCompileTime == 3, "Affine transform must be 3x3 matrix"); - static_assert(TDerived::IsVectorAtCompileTime && TDerived::SizeAtCompileTime == 3, "Translation must be 3 element vector"); - // We map x,y,z to a new space using - // [x y z] -> affine*[x y z] + translation - // - // Using a for affine and t for translation: - // x -> a_00*x + a01*y + a02*z + t0 - // y -> a_10*x + a11*y + a12*z + t1 - // z -> a_20*x + a21*y + a22*z + t2 - // - // So - // Ax^2 + By^2 + Cz^2 + 2Fyz + 2Gzx + 2Hxy + 2Ux + 2Vy + 2Wz + D - // becomes - // A(a_00*x + a01*y + a02*z + t0)(a_00*x + a01*y + a02*z + t0) - // + B(a_10*x + a11*y + a12*z + t1)(a_10*x + a11*y + a12*z + t1) - // + C(a_20*x + a21*y + a22*z + t2)(a_20*x + a21*y + a22*z + t2) - // + 2F(a_10*x + a11*y + a12*z + t1)(a_20*x + a21*y + a22*z + t2) - // + 2G(a_20*x + a21*y + a22*z + t2)(a_00*x + a01*y + a02*z + t0) - // + 2H(a_00*x + a01*y + a02*z + t0)(a_10*x + a11*y + a12*z + t1) - // + 2U(a_00*x + a01*y + a02*z + t0) - // + 2V(a_10*x + a11*y + a12*z + t1) - // + 2W(a_20*x + a21*y + a22*z + t2) - // + D - // - // Collecting terms gives: - return Conicoid( - A * sq(a(0, 0)) + B * sq(a(1, 0)) + C * sq(a(2, 0)) + Scalar(2) * F * a(1, 0) * a(2, 0) + Scalar(2) * G * a(0, 0) * a(2, 0) + Scalar(2) * H * a(0, 0) * a(1, 0), - A * sq(a(0, 1)) + B * sq(a(1, 1)) + C * sq(a(2, 1)) + Scalar(2) * F * a(1, 1) * a(2, 1) + Scalar(2) * G * a(0, 1) * a(2, 1) + Scalar(2) * H * a(0, 1) * a(1, 1), - A * sq(a(0, 2)) + B * sq(a(1, 2)) + C * sq(a(2, 2)) + Scalar(2) * F * a(1, 2) * a(2, 2) + Scalar(2) * G * a(0, 2) * a(2, 2) + Scalar(2) * H * a(0, 2) * a(1, 2), - A * a(0, 1) * a(0, 2) + B * a(1, 1) * a(1, 2) + C * a(2, 1) * a(2, 2) + F * a(1, 1) * a(2, 2) + F * a(1, 2) * a(2, 1) + G * a(0, 1) * a(2, 2) + G * a(0, 2) * a(2, 1) + H * a(0, 1) * a(1, 2) + H * a(0, 2) * a(1, 1), - A * a(0, 2) * a(0, 0) + B * a(1, 2) * a(1, 0) + C * a(2, 2) * a(2, 0) + F * a(1, 2) * a(2, 0) + F * a(1, 0) * a(2, 2) + G * a(0, 2) * a(2, 0) + G * a(0, 0) * a(2, 2) + H * a(0, 2) * a(1, 0) + H * a(0, 0) * a(1, 2), - A * a(0, 0) * a(0, 1) + B * a(1, 0) * a(1, 1) + C * a(2, 0) * a(2, 1) + F * a(1, 0) * a(2, 1) + F * a(1, 1) * a(2, 0) + G * a(0, 0) * a(2, 1) + G * a(0, 1) * a(2, 0) + H * a(0, 0) * a(1, 1) + H * a(0, 1) * a(1, 0), - A * a(0, 0) * t(0) + B * a(1, 0) * t(1) + C * a(2, 0) * t(2) + F * a(1, 0) * t(2) + F * a(2, 0) * t(1) + G * a(0, 0) * t(2) + G * a(2, 0) * t(0) + H * a(0, 0) * t(1) + H * a(1, 0) * t(0) + U * a(0, 0) + V * a(1, 0) + W * a(2, 0), - A * a(0, 1) * t(0) + B * a(1, 1) * t(1) + C * a(2, 1) * t(2) + F * a(1, 1) * t(2) + F * a(2, 1) * t(1) + G * a(0, 1) * t(2) + G * a(2, 1) * t(0) + H * a(0, 1) * t(1) + H * a(1, 1) * t(0) + U * a(0, 1) + V * a(1, 1) + W * a(2, 1), - A * a(0, 2) * t(0) + B * a(1, 2) * t(1) + C * a(2, 2) * t(2) + F * a(1, 2) * t(2) + F * a(2, 2) * t(1) + G * a(0, 2) * t(2) + G * a(2, 2) * t(0) + H * a(0, 2) * t(1) + H * a(1, 2) * t(0) + U * a(0, 2) + V * a(1, 2) + W * a(2, 2), - A * sq(t(0)) + B * sq(t(1)) + C * sq(t(2)) + Scalar(2) * F * t(1) * t(2) + Scalar(2) * G * t(0) * t(2) + Scalar(2) * H * t(0) * t(1) + Scalar(2) * U * t(0) + Scalar(2) * V * t(1) + Scalar(2) * W * t(2) + D - ); - } - }; - - template - std::ostream& operator<< (std::ostream& os, const Conicoid& conicoid) - { - return os << "Conicoid { " << conicoid.A << "x^2 + " << conicoid.B << "y^2 + " << conicoid.C << "z^2 + " - "2*" << 2 * conicoid.F << "yz + 2*" << 2 * conicoid.G << "zx + 2*" << 2 * conicoid.H << "xy + " - "2*" << 2 * conicoid.U << "x + 2*" << 2 * conicoid.V << "y + 2*" << 2 * conicoid.W << "z + " << conicoid.D << " = 0 }"; - } - -} - -#endif diff --git a/pupil_src/shared_cpp/include/geometry/Ellipse.h b/pupil_src/shared_cpp/include/geometry/Ellipse.h deleted file mode 100644 index 244d0306cf..0000000000 --- a/pupil_src/shared_cpp/include/geometry/Ellipse.h +++ /dev/null @@ -1,156 +0,0 @@ -#ifndef _ELLIPSE_H_ -#define _ELLIPSE_H_ - -#include "common/constants.h" - -#include - -namespace singleeyefitter { - - template - class Conic; - - template - class Ellipse2D { - public: - typedef T Scalar; - typedef Eigen::Matrix Vector; - Vector center; - Scalar major_radius; - Scalar minor_radius; - Scalar angle; - - Ellipse2D() - : center(0, 0), major_radius(0), minor_radius(0), angle(0) - { - } - template - Ellipse2D(const Eigen::EigenBase& center, Scalar major_radius, Scalar minor_radius, Scalar angle) - : center(center), major_radius(major_radius), minor_radius(minor_radius), angle(angle) - { - } - Ellipse2D(Scalar x, Scalar y, Scalar major_radius, Scalar minor_radius, Scalar angle) - : center(x, y), major_radius(major_radius), minor_radius(minor_radius), angle(angle) - { - } - template - explicit Ellipse2D(const Conic& conic) - { - using std::atan2; - using std::sin; - using std::cos; - using std::sqrt; - using std::abs; - angle = 0.5 * atan2(conic.B, conic.A - conic.C); - auto cost = cos(angle); - auto sint = sin(angle); - auto sin_squared = sint * sint; - auto cos_squared = cost * cost; - auto Ao = conic.F; - auto Au = conic.D * cost + conic.E * sint; - auto Av = -conic.D * sint + conic.E * cost; - auto Auu = conic.A * cos_squared + conic.C * sin_squared + conic.B * sint * cost; - auto Avv = conic.A * sin_squared + conic.C * cos_squared - conic.B * sint * cost; - // ROTATED = [Ao Au Av Auu Avv] - auto tucenter = -Au / (2.0 * Auu); - auto tvcenter = -Av / (2.0 * Avv); - auto wcenter = Ao - Auu * tucenter * tucenter - Avv * tvcenter * tvcenter; - center[0] = tucenter * cost - tvcenter * sint; - center[1] = tucenter * sint + tvcenter * cost; - major_radius = sqrt(abs(-wcenter / Auu)); - minor_radius = sqrt(abs(-wcenter / Avv)); - - if (major_radius < minor_radius) { - std::swap(major_radius, minor_radius); - angle = angle + constants::PI / 2; - } - - if (angle > constants::PI ) - angle = angle - constants::PI ; - } - - EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF((sizeof(Vector) % 16) == 0) - - Vector major_axis() const - { - using std::sin; - using std::cos; - return Vector(major_radius * sin(angle), major_radius * cos(angle)); - } - Vector minor_axis() const - { - using std::sin; - using std::cos; - return Vector(-minor_radius * cos(angle), minor_radius * sin(angle)); - } - - Scalar circumference() const - { - using std::abs; - using std::sqrt; - using std::pow; - return constants::PI * abs(3.0 * (major_radius + minor_radius) - - sqrt(10.0 * major_radius * minor_radius + 3.0 * - (pow(major_radius, 2) + pow(minor_radius, 2)))); - } - Scalar area() const - { - return constants::PI * major_radius * minor_radius; - } - - - static const Ellipse2D Null; - - private: - - }; - - template - const Ellipse2D Ellipse2D::Null = Ellipse2D(); - - template - bool operator==(const Ellipse2D& el1, const Ellipse2D& el2) - { - return el1.center[0] == el2.center[0] && - el1.center[1] == el2.center[1] && - el1.major_radius == el2.major_radius && - el1.minor_radius == el2.minor_radius && - el1.angle == el2.angle; - } - template - bool operator!=(const Ellipse2D& el1, const Ellipse2D& el2) - { - return !(el1 == el2); - } - - template - std::ostream& operator<< (std::ostream& os, const Ellipse2D& ellipse) - { - return os << "Ellipse { center: (" << ellipse.center[0] << "," << ellipse.center[1] << "), a: " << - ellipse.major_radius << ", b: " << ellipse.minor_radius << ", theta: " << (ellipse.angle / constants::PI) << "pi }"; - } - - template - Ellipse2D scaled(const Ellipse2D& ellipse, U scale) - { - return Ellipse2D( - ellipse.center[0].a, - ellipse.center[1].a, - ellipse.major_radius.a, - ellipse.minor_radius.a, - ellipse.angle.a); - } - - template - inline Eigen::Matrix::type, 2, 1> pointAlongEllipse(const Ellipse2D& el, Scalar2 t) - { - using std::sin; - using std::cos; - auto xt = el.center.x() + el.major_radius * cos(el.angle) * cos(t) - el.minor_radius * sin(el.angle) * sin(t); - auto yt = el.center.y() + el.major_radius * sin(el.angle) * cos(t) + el.minor_radius * cos(el.angle) * sin(t); - return Eigen::Matrix::type, 2, 1>(xt, yt); - } - -} - -#endif diff --git a/pupil_src/shared_cpp/include/geometry/Sphere.h b/pupil_src/shared_cpp/include/geometry/Sphere.h deleted file mode 100644 index 70432341f7..0000000000 --- a/pupil_src/shared_cpp/include/geometry/Sphere.h +++ /dev/null @@ -1,56 +0,0 @@ -#ifndef _SPHERE_H_ -#define _SPHERE_H_ - -#include - -namespace singleeyefitter { - - template - class Sphere { - public: - typedef T Scalar; - typedef Eigen::Matrix Vector; - - Vector center; - Scalar radius; - - Sphere() : center(0, 0, 0), radius(0) - { - } - Sphere(Vector center, Scalar radius) - : center(std::move(center)), radius(std::move(radius)) - { - } - - static const Sphere Null; - - private: - - }; - - template - const Sphere Sphere::Null = Sphere(); - - template - bool operator== (const Sphere& s1, const Sphere& s2) - { - return s1.center == s2.center - && s1.radius == s2.radius; - } - template - bool operator!= (const Sphere& s1, const Sphere& s2) - { - return s1.center != s2.center - || s1.radius != s2.radius; - } - - template - std::ostream& operator<< (std::ostream& os, const Sphere& circle) - { - return os << "Sphere { center: (" << circle.center[0] << "," << circle.center[1] << "," << circle.center[2] << "), " - "radius: " << circle.radius << " }"; - } - -} - -#endif//_SPHERE_H_ diff --git a/pupil_src/shared_cpp/include/math/distance.h b/pupil_src/shared_cpp/include/math/distance.h deleted file mode 100644 index 4a029d1790..0000000000 --- a/pupil_src/shared_cpp/include/math/distance.h +++ /dev/null @@ -1,185 +0,0 @@ -#ifndef distance_h__ -#define distance_h__ - -#include -#include -#include "geometry/Ellipse.h" -#include "DistancePointEllipse.h" - -using namespace singleeyefitter; - -template -auto euclidean_distance(const Eigen::MatrixBase& p1, const Eigen::MatrixBase& p2) -> decltype((p1 - p2).norm()) -{ - return (p1 - p2).norm(); -} -template -auto euclidean_distance_squared(const Eigen::MatrixBase& p1, const Eigen::MatrixBase& p2) -> decltype((p1 - p2).squaredNorm()) -{ - return (p1 - p2).squaredNorm(); -} -template -auto euclidean_distance(const Eigen::MatrixBase& point, - const Eigen::ParametrizedLine& line) -> decltype(point.norm()) -{ - return ((line.origin() - point) - ((line.origin() - point).dot(line.direction())) * line.direction()).norm(); -} - -template -auto euclidean_distance_squared(const Eigen::MatrixBase& point, - const Eigen::ParametrizedLine& line) -> decltype(point.norm()) -{ - return ((line.origin() - point) - ((line.origin() - point).dot(line.direction())) * line.direction()).squaredNorm(); -} - - -template -Scalar euclidean_distance(const Eigen::ParametrizedLine& line1, - const Eigen::ParametrizedLine& line2) -{ - return std::sqrt( euclidean_distance_squared(line1, line2) ) ; -} - -template -Scalar euclidean_distance_squared(const Eigen::ParametrizedLine& line1, - const Eigen::ParametrizedLine& line2) -{ - - auto closestPoints = closest_points_on_line(line1, line2); - auto diff = closestPoints.first - closestPoints.second; - return diff.dot(diff); -} - -template -std::pair< typename Eigen::ParametrizedLine::VectorType , typename Eigen::ParametrizedLine::VectorType > -closest_points_on_line(const Eigen::ParametrizedLine& line1, - const Eigen::ParametrizedLine& line2) -{ - typedef typename Eigen::ParametrizedLine::VectorType Vector; - Vector diff = line1.origin() - line2.origin(); - Scalar a01 = -line1.direction().dot(line2.direction()); - Scalar b0 = diff.dot(line1.direction()); - Scalar s0, s1; - - if (std::abs(a01) < Scalar(1) ) - { - // Lines are not parallel. - Scalar det = Scalar(1) - a01 * a01; - Scalar b1 = -diff.dot(line2.direction()); - s0 = (a01 * b1 - b0) / det; - s1 = (a01 * b0 - b1) / det; - } - else - { - // Lines are parallel, select any pair of closest points. - s0 = -b0; - s1 = Scalar(0); - } - - Vector closestPoint1 = line1.origin() + s0 * line1.direction(); - Vector closestPoint2 = line2.origin() + s1 * line2.direction(); - return std::pair(closestPoint1, closestPoint2); -} - -template -Scalar euclidean_distance(const Eigen::Matrix& p, const Eigen::Matrix& v, const Eigen::Matrix& w) -{ - // Return minimum distance between line segment vw and point p - auto l2 = (v - w).squaredNorm(); - - if (l2 == 0.0) - return euclidean_distance(p, v); // v == w case - - // Consider the line extending the segment, parameterized as v + t (w - v). - // We find projection of point p onto the line. - // It falls where t = [(p-v) . (w-v)] / |w-v|^2 - auto t = (p - v).dot(w - v) / l2; - - if (t < 0.0) - return euclidean_distance(p, v); // Beyond the 'v' end of the segment - else if (t > 1.0) - return euclidean_distance(p, w); // Beyond the 'w' end of the segment - - auto projection = v + t * (w - v); // Projection falls on the segment - return euclidean_distance(p, projection); -} - -template -Scalar euclidean_distance(const Eigen::Matrix& point, const std::vector>& polygon) -{ - auto from = polygon.back(); - Scalar min_distance = std::numeric_limits::infinity(); - - for (const auto& to : polygon) { - min_distance = std::min(min_distance, euclidean_distance(point, from, to)); - from = to; - } - - return min_distance; -} -template -Scalar euclidean_distance( const std::vector>& polygon) -{ - auto from = polygon.back(); - Scalar distance = 0.0; - - for (const auto& to : polygon) { - distance += euclidean_distance(from, to); - from = to; - } - - return distance; -} - - -template -Scalar euclidean_distance(const Eigen::Matrix& point, const Ellipse2D& ellipse) -{ - return DistancePointEllipse(ellipse, point[0], point[1]); -} -template -Scalar euclidean_distance(const Scalar x, const Scalar y, const Ellipse2D& ellipse) -{ - return DistancePointEllipse(ellipse, x, y); -} - -template -Scalar oneway_hausdorff_distance(const Ellipse2D& ellipse, const TOther& other) -{ - Scalar max_dist = -1; - - for (Scalar i = 0; i < 100; ++i) { - Scalar t = i * 2 * M_PI / 100; - auto pt = pointAlongEllipse(ellipse, t); - Scalar i_dist = euclidean_distance(pt, other); - max_dist = std::max(max_dist, i_dist); - } - - return max_dist; -} -template -Scalar oneway_hausdorff_distance(const std::vector>& polygon, const TOther& other) -{ - Scalar max_dist = -1; - - for (const auto& pt : polygon) { - Scalar pt_dist = euclidean_distance(pt, other); - max_dist = std::max(max_dist, pt_dist); - } - - return max_dist; -} - -template -Scalar hausdorff_distance(const Ellipse2D& ellipse, const TOther& other) -{ - return std::max(oneway_hausdorff_distance(ellipse, other), oneway_hausdorff_distance(ellipse, other)); -} -template -typename std::enable_if < !std::is_same>::value, Scalar >::type - hausdorff_distance(const TOther& other, const Ellipse2D& ellipse) -{ - return hausdorff_distance(ellipse, other); -} - -#endif // distance_h__ diff --git a/pupil_src/shared_cpp/include/math/intersect.h b/pupil_src/shared_cpp/include/math/intersect.h deleted file mode 100644 index c73c3ec076..0000000000 --- a/pupil_src/shared_cpp/include/math/intersect.h +++ /dev/null @@ -1,127 +0,0 @@ -#ifndef __INTERSECT_H__ -#define __INTERSECT_H__ - -#include -#include -#include "geometry/Sphere.h" -#include "mathHelper.h" - -namespace singleeyefitter { - -// Finds the intersection of 2 lines in 2D - template - Eigen::Matrix intersect(const Eigen::ParametrizedLine& line1, const Eigen::ParametrizedLine& line2) - { - Scalar x1 = line1.origin().x(); - Scalar y1 = line1.origin().y(); - Scalar x2 = (line1.origin() + line1.direction()).x(); - Scalar y2 = (line1.origin() + line1.direction()).y(); - Scalar x3 = line2.origin().x(); - Scalar y3 = line2.origin().y(); - Scalar x4 = (line2.origin() + line2.direction()).x(); - Scalar y4 = (line2.origin() + line2.direction()).y(); - Scalar denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); - Scalar px = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * x4 - y3 * x4)) / denom; - Scalar py = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * x4 - y3 * x4)) / denom; - return Eigen::Matrix(px, py); - } - - template - typename Range::value_type::VectorType nearest_intersect(const Range& lines); - - namespace detail { - template - struct intersect_helper { - static Eigen::Matrix nearest(const Eigen::ParametrizedLine& line1, const Eigen::ParametrizedLine& line2) - { - std::vector> lines; - lines.push_back(line1); - lines.push_back(line2); - return nearest_intersect(lines); - } - }; - - // Specialise for 2D - template - struct intersect_helper { - static Eigen::Matrix nearest(const Eigen::ParametrizedLine& line1, const Eigen::ParametrizedLine& line2) - { - return intersect(line1, line2); - } - }; - } - -// Finds the intersection (in a least-squares sense, or exact in 2D) of 2 lines - template - Eigen::Matrix nearest_intersect(const Eigen::ParametrizedLine& line1, const Eigen::ParametrizedLine& line2) - { - return detail::intersect_helper::nearest(line1, line2); - } - -// Finds the intersection (in a least-squares sense) of multiple lines - template - typename Range::value_type::VectorType nearest_intersect(const Range& lines) - { - // NOTE: The type 'Range' will be of type std::vector> - typedef typename Range::value_type::VectorType Vector; - typedef typename Eigen::Matrix Matrix; - static_assert(Vector::ColsAtCompileTime == 1, "Requires column vector"); - //size_t N = lines.size(); - std::vector Ivv; - Matrix A = Matrix::Zero(); - Vector b = Vector::Zero(); - size_t i = 0; - - for (auto& line : lines) { - Vector vi = line.direction(); - Vector pi = line.origin(); - Matrix Ivivi = Matrix::Identity() - vi * vi.transpose(); - Ivv.push_back(Ivivi); - A += Ivivi; - b += (Ivivi * pi); - i++; - } - - Vector x = A.partialPivLu().solve(b); - /*if nargout == 2 - % Calculate error term - for i = 1:n - pi = p(i,:)'; - e = e + ( (x-pi)' * Ivv(:,:,i) * (x-pi) ); - end - end*/ - return x; - } - -// Finds the intersection of a line and a sphere - template - bool intersect(const Eigen::ParametrizedLine& line, const Sphere& sphere, std::pair, Eigen::Matrix>& intersected_points ) - { - using std::sqrt; - using math::sq; - typedef typename Eigen::ParametrizedLine::VectorType Vector; - assert(std::abs(line.direction().norm() - 1) < 0.0001); - Vector v = line.direction(); - // Put p at origin - Vector p = line.origin(); - Vector c = sphere.center - p; - Scalar r = sphere.radius; - // From Wikipedia :) - Scalar vcvc_cc_rr = sq(v.dot(c)) - c.dot(c) + sq(r); - - if (vcvc_cc_rr < 0) { - return false; - } - - Scalar s1 = v.dot(c) - sqrt(vcvc_cc_rr); - Scalar s2 = v.dot(c) + sqrt(vcvc_cc_rr); - Vector p1 = p + s1 * v; - Vector p2 = p + s2 * v; - intersected_points.first = p1; - intersected_points.second = p2; - return true; - } - -} - -#endif//__INTERSECT_H__ diff --git a/pupil_src/shared_cpp/include/mathHelper.h b/pupil_src/shared_cpp/include/mathHelper.h deleted file mode 100644 index 8966dcc89d..0000000000 --- a/pupil_src/shared_cpp/include/mathHelper.h +++ /dev/null @@ -1,351 +0,0 @@ -#ifndef singleeyefitter_math_h__ -#define singleeyefitter_math_h__ - -#include -#include - -#include "common/traits.h" - -namespace singleeyefitter { - - namespace math { - -#define MAKE_SQ(TYPE) \ - inline auto sq(TYPE val) -> decltype(val*val) { return val * val; } - - MAKE_SQ(float) - MAKE_SQ(double) - MAKE_SQ(long double) - MAKE_SQ(char) - MAKE_SQ(short) - MAKE_SQ(int) - MAKE_SQ(long) - MAKE_SQ(long long) - MAKE_SQ(unsigned char) - MAKE_SQ(unsigned short) - MAKE_SQ(unsigned int) - MAKE_SQ(unsigned long) - MAKE_SQ(unsigned long long) - -#undef MAKE_STD_SQ - - template - inline T round(T value, T precision) { - T factor = T(1) / precision; - return floor( value * factor + 0.5 ) / factor; - } - - template - inline T clamp(T val, TMin min_val = std::numeric_limits::min(), TMax max_val = std::numeric_limits::max()) - { - if (min_val > max_val) - return clamp(val, max_val, min_val); - - if (val <= min_val) - return min_val; - - if (val >= max_val) - return max_val; - - return val; - } - - template - inline T lerp(const T& val1, const T& val2, const T& alpha) - { - return val1 * (1.0 - alpha) + val2 * alpha; - } - - template - Scalar getAngleABC(const T& a, const T& b, const T& c) - { - T ab = { b.x - a.x, b.y - a.y }; - T cb = { b.x - c.x, b.y - c.y }; - Scalar dot = ab.dot(cb); // dot product - Scalar cross = ab.cross(cb); // cross product - Scalar alpha = atan2(cross, dot); - return alpha * Scalar(180.0) / M_PI; - } - - template - inline T smootherstep(T edge0, T edge1, T x, scalar_tag) - { - if (x >= edge1) - return T(1); - else if (x <= edge0) - return T(0); - else { - x = (x - edge0) / (edge1 - edge0); - return x * x * x * (x * (x * T(6) - T(15)) + T(10)); - } - } - template - inline ::ceres::Jet smootherstep(T edge0, T edge1, const ::ceres::Jet& f, ceres_jet_tag) - { - if (f.a >= edge1) - return ::ceres::Jet(1); - else if (f.a <= edge0) - return ::ceres::Jet(0); - else { - T x = (f.a - edge0) / (edge1 - edge0); - // f is referenced by this function, so create new value for return. - ::ceres::Jet g; - g.a = x * x * x * (x * (x * T(6) - T(15)) + T(10)); - g.v = f.v * (x * x * (x * (x * T(30) - T(60)) + T(30)) / (edge1 - edge0)); - return g; - } - } - template - inline ::ceres::Jet smootherstep(T edge0, T edge1, ::ceres::Jet&& f, ceres_jet_tag) - { - if (f.a >= edge1) - return ::ceres::Jet(1); - else if (f.a <= edge0) - return ::ceres::Jet(0); - else { - T x = (f.a - edge0) / (edge1 - edge0); - // f is moved into this function, so reuse it. - f.a = x * x * x * (x * (x * T(6) - T(15)) + T(10)); - f.v *= (x * x * (x * (x * T(30) - T(60)) + T(30)) / (edge1 - edge0)); - return f; - } - } - template - inline auto smootherstep(typename ad_traits::scalar edge0, typename ad_traits::scalar edge1, T&& val) - -> decltype(smootherstep(edge0, edge1, std::forward(val), typename ad_traits::ad_tag())) - { - return smootherstep(edge0, edge1, std::forward(val), typename ad_traits::ad_tag()); - } - - template - inline T norm(T x, T y, scalar_tag) - { - using std::sqrt; - using math::sq; - return sqrt(sq(x) + sq(y)); - } - template - inline ::ceres::Jet norm(const ::ceres::Jet& x, const ::ceres::Jet& y, ceres_jet_tag) - { - T anorm = norm(x.a, y.a, scalar_tag()); - ::ceres::Jet g; - g.a = anorm; - g.v = (x.a / anorm) * x.v + (y.a / anorm) * y.v; - return g; - } - template - inline typename std::decay::type norm(T&& x, T&& y) - { - return norm(std::forward(x), std::forward(y), typename ad_traits::ad_tag()); - } - - template - inline auto Heaviside(T&& val, typename ad_traits::scalar epsilon) -> decltype(smootherstep(-epsilon, epsilon, std::forward(val))) - { - return smootherstep(-epsilon, epsilon, std::forward(val)); - } - - template - Eigen::Matrix sph2cart(T r, T theta, T psi) - { - using std::sin; - using std::cos; - return r * Eigen::Matrix(sin(theta) * cos(psi), cos(theta), sin(theta) * sin(psi)); - } - - template - Eigen::Matrix cart2sph(T x, T y, T z) - { - using std::sin; - using std::cos; - using std::sqrt; - double r = sqrt( x*x + y*y + z*z); - double theta = acos( y / r); - double psi = atan2(z, x ); - return Eigen::Matrix(theta,psi); - } - template - Eigen::Matrix cart2sph(const Eigen::Matrix& m ) - { - return cart2sph( m.x(), m.y(), m.z()); - } - - template - inline ::ceres::Jet sq(::ceres::Jet val) - { - val.v *= 2 * val.a; - val.a *= val.a; - return val; - } - template - inline int sign(const T& z) - { - return (z == 0) ? 0 : (z < 0) ? -1 : 1; - } - - template - T haversine(T theta1, T psi1, T theta2, T psi2 ) - { - using std::sin; - using std::cos; - using std::acos; - using std::asin; - using std::atan2; - using std::sqrt; - using singleeyefitter::math::sq; - - if (theta1 == theta2 && psi1 == psi2) { - return T(0); - } - // Haversine distance - auto dist = T(2) * asin(sqrt( (sin((theta2 - theta1) / T(2))*sin((theta2 - theta1) / T(2))) + cos(theta1) * cos(theta2) * (sin((psi2 - psi1) / T(2))*sin((psi2 - psi1) / T(2))) )); - return dist; - - } - - - - // Hash function for Eigen matrix and vector. - // The code is from `hash_combine` function of the Boost library. See - // http://www.boost.org/doc/libs/1_55_0/doc/html/hash/reference.html#boost.hash_combine . - template - struct matrix_hash : std::unary_function { - std::size_t operator()(T const& matrix) const - { - // Note that it is oblivious to the storage order of Eigen matrix (column- or - // row-major). It will give you the same hash value for two different matrices if they - // are the transpose of each other in different storage order. - size_t seed = 0; - for (size_t i = 0; i < matrix.size(); ++i) { - auto elem = *(matrix.data() + i); - seed ^= std::hash()(elem) + 0x9e3779b9 + (seed << 6) + (seed >> 2); - } - return seed; - } - }; - - template - class SMA{ //simple moving average - - public: - - SMA( int windowSize ) : mWindowSize(windowSize), mAverage(0.0), mAverageDirty(true) - {}; - - void addValue( T value ){ - mValues.push_back( value ); - // calculate moving average of value - if( mValues.size() <= mWindowSize || mAverageDirty ){ - mAverageDirty = false; - mAverage = 0.0; - for(auto& element : mValues){ - mAverage += element; - } - mAverage /= mValues.size(); - }else{ - // we can optimize if the wanted window size is reached - T first = mValues.front(); - mValues.pop_front(); - mAverage += value/mWindowSize - first/mWindowSize; - } - } - - double getAverage() const { return mAverage; }; - int getWindowSize() const { return mWindowSize; }; - - void changeWindowSize( int windowSize){ - - if( windowSize < mWindowSize){ - - if( mValues.size() > windowSize ) - mAverageDirty = true; - while( mValues.size() > windowSize){ - mValues.pop_front(); - } - - } - mWindowSize = windowSize; - - } - - private: - - SMA(){}; - - std::list mValues; - int mWindowSize; - T mAverage; - bool mAverageDirty; // when we change the window size we need to recalculate from ground up - }; - - template - class WMA{ //weighted moving average - - public: - - WMA( int windowSize ) : mWindowSize(windowSize) , mDenominator(1.0), mNumerator(0.0), mAverage(0.0), mAverageDirty(true) - {}; - - void addValue( T value , T weight ){ - mValues.emplace_back( value, weight ); - // calculate weighted moving average of value - - if( mValues.size() <= mWindowSize || mAverageDirty){ - mAverageDirty = false; - mDenominator = 0.0; - mNumerator = 0.0; - for(auto& element : mValues){ - mNumerator += element.first * element.second; - mDenominator += element.second; - } - mAverage = mNumerator / mDenominator; - }else{ - // we can optimize if the wanted window size is reached - auto observation = mValues.front(); - mValues.pop_front(); - mDenominator -= observation.second; - mDenominator += weight; - - mNumerator -= observation.first * observation.second; - mNumerator += value * weight; - mAverage = mNumerator / mDenominator; - } - - } - - double getAverage() const { return mAverage; }; - int getWindowSize() const { return mWindowSize; }; - - void changeWindowSize( int windowSize){ - - if( windowSize < mWindowSize){ - - if( mValues.size() > windowSize ) - mAverageDirty = true; - while( mValues.size() > windowSize){ - mValues.pop_front(); - } - - } - mWindowSize = windowSize; - - } - - private: - - WMA(){}; - - std::list> mValues; - int mWindowSize; - T mDenominator; - T mNumerator; - T mAverage; - bool mAverageDirty; - }; - - } // math namespace - -} // singleeyefitter namespace - - -#endif // singleeyefitter_math_h__ diff --git a/pupil_src/shared_cpp/include/projection.h b/pupil_src/shared_cpp/include/projection.h deleted file mode 100644 index 4b0f0ba58b..0000000000 --- a/pupil_src/shared_cpp/include/projection.h +++ /dev/null @@ -1,293 +0,0 @@ -#ifndef singleeyefitter_project_h__ -#define singleeyefitter_project_h__ - -#include -#include -#include "geometry/Ellipse.h" -#include "geometry/Circle.h" -#include "geometry/Conic.h" -#include "geometry/Sphere.h" -#include "geometry/Conicoid.h" -#include "mathHelper.h" -#include "solve.h" - -namespace singleeyefitter { - - template - Conic project(const Circle3D& circle, Scalar focal_length) - { - typedef typename Circle3D::Vector Vector; - using math::sq; - Vector c = circle.center; - Vector n = circle.normal; - Scalar r = circle.radius; - Scalar f = focal_length; - // Construct cone with circle as base and vertex v = (0,0,0). - // - // For the circle, - // |p - c|^2 = r^2 where (p-c).n = 0 (i.e. on the circle plane) - // - // A cone is basically concentric circles, with center on the line c->v. - // For any point p, the corresponding circle center c' is the intersection - // of the line c->v and the plane through p normal to n. So, - // - // d = ((p - v).n)/(c.n) - // c' = d c + v - // - // The radius of these circles decreases linearly as you approach 0, so - // - // |p - c'|^2 = (r*|c' - v|/|c - v|)^2 - // - // Since v = (0,0,0), this simplifies to - // - // |p - (p.n/c.n)c|^2 = (r*|(p.n/c.n)c|/|c|)^2 - // - // |(c.n)p - (p.n)c|^2 / p.n \^2 - // ------------------- = r^2 * | --- | - // (c.n)^2 \ c.n / - // - // |(c.n)p - (p.n)c|^2 - r^2 * (p.n)^2 = 0 - // - // Expanding out p, c and n gives - // - // |(c.n)x - (x*n_x + y*n_y + z*n_z)c_x|^2 - // |(c.n)y - (x*n_x + y*n_y + z*n_z)c_y| - r^2 * (x*n_x + y*n_y + z*n_z)^2 = 0 - // |(c.n)z - (x*n_x + y*n_y + z*n_z)c_z| - // - // ((c.n)x - (x*n_x + y*n_y + z*n_z)c_x)^2 - // + ((c.n)y - (x*n_x + y*n_y + z*n_z)c_y)^2 - // + ((c.n)z - (x*n_x + y*n_y + z*n_z)c_z)^2 - // - r^2 * (x*n_x + y*n_y + z*n_z)^2 = 0 - // - // (c.n)^2 x^2 - 2*(c.n)*(x*n_x + y*n_y + z*n_z)*x*c_x + (x*n_x + y*n_y + z*n_z)^2 c_x^2 - // + (c.n)^2 y^2 - 2*(c.n)*(x*n_x + y*n_y + z*n_z)*y*c_y + (x*n_x + y*n_y + z*n_z)^2 c_y^2 - // + (c.n)^2 z^2 - 2*(c.n)*(x*n_x + y*n_y + z*n_z)*z*c_z + (x*n_x + y*n_y + z*n_z)^2 c_z^2 - // - r^2 * (x*n_x + y*n_y + z*n_z)^2 = 0 - // - // (c.n)^2 x^2 - 2*(c.n)*c_x*(x*n_x + y*n_y + z*n_z)*x - // + (c.n)^2 y^2 - 2*(c.n)*c_y*(x*n_x + y*n_y + z*n_z)*y - // + (c.n)^2 z^2 - 2*(c.n)*c_z*(x*n_x + y*n_y + z*n_z)*z - // + (x*n_x + y*n_y + z*n_z)^2 * (c_x^2 + c_y^2 + c_z^2 - r^2) - // - // (c.n)^2 x^2 - 2*(c.n)*c_x*(x*n_x + y*n_y + z*n_z)*x - // + (c.n)^2 y^2 - 2*(c.n)*c_y*(x*n_x + y*n_y + z*n_z)*y - // + (c.n)^2 z^2 - 2*(c.n)*c_z*(x*n_x + y*n_y + z*n_z)*z - // + (|c|^2 - r^2) * (n_x^2*x^2 + n_y^2*y^2 + n_z^2*z^2 + 2*n_x*n_y*x*y + 2*n_x*n_z*x*z + 2*n_y*n_z*y*z) - // - // Collecting conicoid terms gives - // - // [xyz]^2 : (c.n)^2 - 2*(c.n)*c_[xyz]*n_[xyz] + (|c|^2 - r^2)*n_[xyz]^2 - // [yzx][zxy] : - 2*(c.n)*c_[yzx]*n_[zxy] - 2*(c.n)*c_[zxy]*n_[yzx] + (|c|^2 - r^2)*2*n_[yzx]*n_[zxy] - // : 2*((|c|^2 - r^2)*n_[yzx]*n_[zxy] - (c,n)*(c_[yzx]*n_[zxy] + c_[zxy]*n_[yzx])) - // [xyz] : 0 - // 1 : 0 - Scalar cn = c.dot(n); - Scalar c2r2 = (c.dot(c) - sq(r)); - Vector ABC = (sq(cn) - 2.0 * cn * c.array() * n.array() + c2r2 * n.array().square()); - Scalar F = 2.0 * (c2r2 * n(1) * n(2) - cn * (n(1) * c(2) + n(2) * c(1))); - Scalar G = 2.0 * (c2r2 * n(2) * n(0) - cn * (n(2) * c(0) + n(0) * c(2))); - Scalar H = 2.0 * (c2r2 * n(0) * n(1) - cn * (n(0) * c(1) + n(1) * c(0))); - // Then set z=f to get conic which is the result of intersecting the cone with the focal plane - return Conic( - ABC(0), // x^2 (Ax^2) - H, // xy (Hxy) - ABC(1), // y^2 (By^2) - G * f /*+ Const(0)*/, // x (Gxz + Ux, z = f) - F * f /*+ Const(0)*/, // y (Fyz + Vy, z = f) - ABC(2) * sq(f) /*+ Const(0)*f + Const(0)*/ // 1 (Cz^2 + Wz + D, z = f) - ); - } - - /*template - typename Conic project(Conic conic, const Eigen::DenseBase& point, const Eigen::DenseBase& normal, Scalar focal_length) - { - // Consider two coordinate systems: - // camera (camera at 0, x,y aligned with image plane, z going away from camera) - // conic (conic on xy-plane, with plane normal = (0,0,1) and plane point = (0,0,0) ) - // - // To project conic lying on plane defined by point and normal (point corresponding to (0,0) in conic's 2D space), do: - // - // Input as in camera space, - // Transform to conic space, - // Form conicoid with conic as base and camera center as vertex - // Transform back to camera space - // Intersect conicoid with image plane (z=f) - - Eigen::Matrix camera_center(0,0,0); - }*/ - - - template - Ellipse2D project(const Sphere& sphere, Scalar focal_length) - { - return Ellipse2D( - focal_length * sphere.center.template head<2>() / sphere.center[2], - focal_length * sphere.radius / sphere.center[2], - focal_length * sphere.radius / sphere.center[2], - 0); - } - template - typename Eigen::DenseBase::template FixedSegmentReturnType<2>::Type::PlainObject project(const Eigen::DenseBase& point, typename Eigen::DenseBase::Scalar focal_length) - { - static_assert(Derived::IsVectorAtCompileTime && Derived::SizeAtCompileTime == 3, "Point must be 3 element vector"); - return focal_length * point.template head<2>() / point(2); - } - - - template - std::pair, Circle3D> unproject(const Ellipse2D& ellipse, Scalar circle_radius, Scalar focal_length) - { - using std::sqrt; - using std::abs; - using math::sign; - using math::sq; - typedef Conic Conic; - typedef Conicoid Conicoid; - typedef Circle3D Circle; - typedef Eigen::Matrix Matrix3; - typedef Eigen::Matrix Vector3; - typedef Eigen::Array RowArray3; - typedef Eigen::Translation Translation3; - // Get cone with base of ellipse and vertex at [0 0 -f] - // Safaee-Rad 1992 eq (3) - Conic conic(ellipse); - Vector3 cam_center_in_ellipse(0, 0, -focal_length); - Conicoid pupil_cone(conic, cam_center_in_ellipse); - auto a = pupil_cone.A; - auto b = pupil_cone.B; - auto c = pupil_cone.C; - auto f = pupil_cone.F; - auto g = pupil_cone.G; - auto h = pupil_cone.H; - auto u = pupil_cone.U; - auto v = pupil_cone.V; - auto w = pupil_cone.W; - //auto d = pupil_cone.D; - // Get canonical conic form: - // lambda(1) X^2 + lambda(2) Y^2 + lambda(3) Z^2 = mu - // Safaee-Rad 1992 eq (6) - // Done by solving the discriminating cubic (10) - // Lambdas are sorted descending because order of roots doesn't - // matter, and it later eliminates the case of eq (30), where - // lambda(2) > lambda(1) - RowArray3 lambda; - std::tie(lambda(0), lambda(1), lambda(2)) = solve(1., -(a + b + c), (b * c + c * a + a * b - f * f - g * g - h * h), -(a * b * c + 2 * f * g * h - a * f * f - b * g * g - c * h * h)); - assert(lambda(0) >= lambda(1)); - assert(lambda(1) > 0); - assert(lambda(2) < 0); - // Now want to calculate l,m,n of the plane - // lX + mY + nZ = p - // which intersects the cone to create a circle. - // Safaee-Rad 1992 eq (31) - // [Safaee-Rad 1992 eq (33) comes out of this as a result of lambda(1) == lambda(2)] - auto n = sqrt((lambda(1) - lambda(2)) / (lambda(0) - lambda(2))); - auto m = 0.0; - auto l = sqrt((lambda(0) - lambda(1)) / (lambda(0) - lambda(2))); - // There are two solutions for l, positive and negative, we handle these later - // Want to calculate T1, the rotation transformation from image - // space in the canonical conic frame back to image space in the - // real world - Matrix3 T1; - // Safaee-Rad 1992 eq (8) - auto li = T1.row(0); - auto mi = T1.row(1); - auto ni = T1.row(2); - // Safaee-Rad 1992 eq (12) - RowArray3 t1 = (b - lambda) * g - f * h; - RowArray3 t2 = (a - lambda) * f - g * h; - RowArray3 t3 = -(a - lambda) * (t1 / t2) / g - h / g; - mi = 1 / sqrt(1 + (t1 / t2).square() + t3.square()); - li = (t1 / t2) * mi.array(); - ni = t3 * mi.array(); - - // If li,mi,ni follow the left hand rule, flip their signs - if ((li.cross(mi)).dot(ni) < 0) { - li = -li; - mi = -mi; - ni = -ni; - } - - // Calculate T2, a translation transformation from the canonical - // conic frame to the image space in the canonical conic frame - // Safaee-Rad 1992 eq (14) - Translation3 T2; - T2.translation() = -(u * li + v * mi + w * ni).array() / lambda; - Circle solutions[2]; - Scalar ls[2] = { l, -l }; - - for (int i = 0; i < 2; i++) { - auto l = ls[i]; - // Circle normal in image space (i.e. gaze vector) - Vector3 gaze = T1 * Vector3(l, m, n); - // Calculate T3, a rotation from a frame where Z is the circle normal - // to the canonical conic frame - // Safaee-Rad 1992 eq (19) - // Want T3 = / -m/sqrt(l*l+m*m) -l*n/sqrt(l*l+m*m) l \ - // | l/sqrt(l*l+m*m) -m*n/sqrt(l*l+m*m) m | - // \ 0 sqrt(l*l+m*m) n / - // But m = 0, so this simplifies to - // T3 = / 0 -n*l/sqrt(l*l) l \ - // | l/sqrt(l*l) 0 0 | - // \ 0 sqrt(l*l) n / - // = / 0 -n*sgn(l) l \ - // | sgn(l) 0 0 | - // \ 0 |l| n / - Matrix3 T3; - - if (l == 0) { - // Discontinuity of sgn(l), have to handle explicitly - assert(n == 1); - std::cout << "Warning: l == 0" << std::endl; - T3 << 0, -1, 0, - 1, 0, 0, - 0, 0, 1; - - } else { - //auto sgnl = sign(l); - T3 << 0, -n* sign(l), l, - sign(l), 0, 0, - 0, abs(l), n; - } - - // Calculate the circle center - // Safaee-Rad 1992 eq (38), using T3 as defined in (36) - auto A = lambda.matrix().dot(T3.col(0).cwiseAbs2()); - auto B = lambda.matrix().dot(T3.col(0).cwiseProduct(T3.col(2))); - auto C = lambda.matrix().dot(T3.col(1).cwiseProduct(T3.col(2))); - auto D = lambda.matrix().dot(T3.col(2).cwiseAbs2()); - // Safaee-Rad 1992 eq (41) - Vector3 center_in_Xprime; - center_in_Xprime(2) = A * circle_radius / sqrt(sq(B) + sq(C) - A * D); - center_in_Xprime(0) = -B / A * center_in_Xprime(2); - center_in_Xprime(1) = -C / A * center_in_Xprime(2); - // Safaee-Rad 1992 eq (34) - Translation3 T0; - T0.translation() << 0, 0, focal_length; - // Safaee-Rad 1992 eq (42) using (35) - Vector3 center = T0 * T1 * T2 * T3 * center_in_Xprime; - - // If z is negative (behind the camera), choose the other - // solution of eq (41) [maybe there's a way of calculating which - // solution should be chosen first] - - if (center(2) < 0) { - center_in_Xprime = -center_in_Xprime; - center = T0 * T1 * T2 * T3 * center_in_Xprime; - } - - // Make sure that the gaze vector is toward the camera and is normalised - if (gaze.dot(center) > 0) { - gaze = -gaze; - } - - gaze.normalize(); - // Save the results - solutions[i] = Circle(center, gaze, circle_radius); - } - - return std::make_pair(solutions[0], solutions[1]); - } - -} - -#endif // singleeyefitter_project_h__ diff --git a/pupil_src/shared_cpp/include/solve.h b/pupil_src/shared_cpp/include/solve.h deleted file mode 100644 index 8ac4619c6c..0000000000 --- a/pupil_src/shared_cpp/include/solve.h +++ /dev/null @@ -1,114 +0,0 @@ -#ifndef _SOLVE_H_ -#define _SOLVE_H_ - -#include -#include -#include - -namespace singleeyefitter { - - // a = 0 - template - T solve(T a) - { - if (a == 0) return 0; - else throw std::runtime_error("No solution"); - } - // ax + b = 0 - template - T solve(T a, T b) - { - if (a == 0) return solve(b); - - return -b / a; - } - // ax^2 + bx + c = 0 - template - std::tuple solve(T a, T b, T c) - { - using math::sq; - using std::sqrt; - - if (a == 0) { - auto root = solve(b, c); - return std::tuple(root, root); - } - - // http://www.it.uom.gr/teaching/linearalgebra/NumericalRecipiesInC/c5-6.pdf - // Pg 184 - auto det = sq(b) - 4 * a * c; - - if (det < 0) - throw std::runtime_error("No solution"); - - //auto sqrtdet = sqrt(det); - auto q = -0.5 * (b + (b >= 0 ? 1 : -1) * sqrt(det)); - return std::tuple(q / a, c / q); - } - // ax^2 + bx + c = 0 - template - std::tuple solve(T a, T b, T c, T d) - { - using std::sqrt; - using std::abs; - using math::sq; - using std::cbrt; - - if (a == 0) { - auto roots = solve(b, c, d); - return std::tuple(std::get<0>(roots), std::get<1>(roots), std::get<1>(roots)); - } - - // http://www.it.uom.gr/teaching/linearalgebra/NumericalRecipiesInC/c5-6.pdf - // http://web.archive.org/web/20120321013251/http://linus.it.uts.edu.au/~don/pubs/solving.html - auto p = b / a; - auto q = c / a; - auto r = d / a; - //auto Q = (p*p - 3*q) / 9; - //auto R = (2*p*p*p - 9*p*q + 27*r)/54; - auto u = q - sq(p) / 3; - auto v = r - p * q / 3 + 2 * p * p * p / 27; - auto j = 4 * u * u * u / 27 + v * v; - const auto M = std::numeric_limits::max(); - const auto sqrtM = sqrt(M); - const auto cbrtM = cbrt(M); - - if (b == 0 && c == 0) - return std::tuple(cbrt(-d), cbrt(-d), cbrt(-d)); - - if (abs(p) > 27 * cbrtM) - return std::tuple(-p, -p, -p); - - if (abs(q) > sqrtM) - return std::tuple(-cbrt(v), -cbrt(v), -cbrt(v)); - - if (abs(u) > 3 * cbrtM / 4) - return std::tuple(cbrt(4) * u / 3, cbrt(4) * u / 3, cbrt(4) * u / 3); - - if (j > 0) { - // One real root - auto w = sqrt(j); - T y; - - if (v > 0) - y = (u / 3) * cbrt(2 / (w + v)) - cbrt((w + v) / 2) - p / 3; - else - y = cbrt((w - v) / 2) - (u / 3) * cbrt(2 / (w - v)) - p / 3; - - return std::tuple(y, y, y); - - } else { - // Three real roots - auto s = sqrt(-u / 3); - auto t = -v / (2 * s * s * s); - auto k = acos(t) / 3; - auto y1 = 2 * s * cos(k) - p / 3; - auto y2 = s * (-cos(k) + sqrt(3.) * sin(k)) - p / 3; - auto y3 = s * (-cos(k) - sqrt(3.) * sin(k)) - p / 3; - return std::tuple(y1, y2, y3); - } - } - -} - -#endif//_SOLVE_H_ diff --git a/pupil_src/shared_modules/cython_methods/.gitignore b/pupil_src/shared_modules/cython_methods/.gitignore deleted file mode 100644 index cd92b3c4a6..0000000000 --- a/pupil_src/shared_modules/cython_methods/.gitignore +++ /dev/null @@ -1 +0,0 @@ -*.cpp \ No newline at end of file diff --git a/pupil_src/shared_modules/cython_methods/__init__.py b/pupil_src/shared_modules/cython_methods/__init__.py deleted file mode 100644 index 53dcb12205..0000000000 --- a/pupil_src/shared_modules/cython_methods/__init__.py +++ /dev/null @@ -1,18 +0,0 @@ -""" -(*)~--------------------------------------------------------------------------- -Pupil - eye tracking platform -Copyright (C) 2012-2020 Pupil Labs - -Distributed under the terms of the GNU -Lesser General Public License (LGPL v3.0). -See COPYING and COPYING.LESSER for license details. ----------------------------------------------------------------------------~(*) -""" - -try: - from .methods import * -except ModuleNotFoundError: - # when running from source compile cpp extension if necessary. - from .build import build_cpp_extension - build_cpp_extension() - from .methods import * \ No newline at end of file diff --git a/pupil_src/shared_modules/cython_methods/build.py b/pupil_src/shared_modules/cython_methods/build.py deleted file mode 100644 index 9c1522eb6d..0000000000 --- a/pupil_src/shared_modules/cython_methods/build.py +++ /dev/null @@ -1,39 +0,0 @@ -""" -(*)~--------------------------------------------------------------------------- -Pupil - eye tracking platform -Copyright (C) 2012-2020 Pupil Labs - -Distributed under the terms of the GNU -Lesser General Public License (LGPL v3.0). -See COPYING and COPYING.LESSER for license details. ----------------------------------------------------------------------------~(*) -""" - -import logging - -logger = logging.getLogger(__name__) - - -def build_cpp_extension(): - import subprocess as sp - import os, sys - - src_loc = os.path.dirname(os.path.realpath(__file__)) - install_loc = os.path.split(src_loc)[0] - cwd = os.getcwd() - os.chdir(src_loc) - logger.info("Building extension modules...") - build_cmd = [ - sys.executable, - "setup.py", - "install", - "--install-lib={}".format(install_loc), - ] - ret = sp.check_output(build_cmd).decode(sys.stdout.encoding) - logger.debug("Build log:\n{}".format(ret)) - os.chdir(cwd) - - -if __name__ == "__main__": - logging.basicConfig(level=logging.DEBUG) - build_cpp_extension() diff --git a/pupil_src/shared_modules/cython_methods/methods.pyx b/pupil_src/shared_modules/cython_methods/methods.pyx deleted file mode 100644 index 043e3d733d..0000000000 --- a/pupil_src/shared_modules/cython_methods/methods.pyx +++ /dev/null @@ -1,46 +0,0 @@ -""" -(*)~--------------------------------------------------------------------------- -Pupil - eye tracking platform -Copyright (C) 2012-2020 Pupil Labs - -Distributed under the terms of the GNU -Lesser General Public License (LGPL v3.0). -See COPYING and COPYING.LESSER for license details. ----------------------------------------------------------------------------~(*) -""" - -import numpy as np -cimport numpy as np - - -def cumhist_color_map16(np.ndarray[np.uint16_t, ndim=2] depth_buffer): - cdef int r, c, i, f, width, height - cdef np.uint16_t d - height = depth_buffer.shape[0] - width = depth_buffer.shape[1] - - cdef np.ndarray[np.uint32_t, ndim=1] cumhist = np.zeros(0x10000, dtype=np.uint32) - cdef np.ndarray[np.uint8_t, ndim=3] rgb_image = np.empty((height, width, 3), dtype=np.uint8) - - for r in range(height): - for c in range(width): - cumhist[depth_buffer[r, c]] += 1 - - for i in range(2, 0x10000): - cumhist[i] += cumhist[i-1] - - for r in range(height): - for c in range(width): - d = depth_buffer[r, c] - if d != 0: - # 0-255 based on histogram location - f = cumhist[d] * 255 / cumhist[0xFFFF] - rgb_image[r, c, 0] = f - rgb_image[r, c, 1] = 0 - rgb_image[r, c, 2] = 255 - f - else: - rgb_image[r, c, 0] = 0 - rgb_image[r, c, 1] = 5 - rgb_image[r, c, 2] = 20 - - return rgb_image diff --git a/pupil_src/shared_modules/cython_methods/setup.py b/pupil_src/shared_modules/cython_methods/setup.py deleted file mode 100644 index a1a9e0acf7..0000000000 --- a/pupil_src/shared_modules/cython_methods/setup.py +++ /dev/null @@ -1,42 +0,0 @@ -""" -(*)~--------------------------------------------------------------------------- -Pupil - eye tracking platform -Copyright (C) 2012-2020 Pupil Labs - -Distributed under the terms of the GNU -Lesser General Public License (LGPL v3.0). -See COPYING and COPYING.LESSER for license details. ----------------------------------------------------------------------------~(*) -""" - -from distutils.core import setup -from distutils.extension import Extension -from Cython.Build import cythonize -import numpy as np - - -extensions = [ - Extension( - name="cython_methods.methods", - sources=["methods.pyx"], - include_dirs=[np.get_include()], - libraries=[], - library_dirs=[], - extra_link_args=[], - extra_compile_args=["-D_USE_MATH_DEFINES", "-std=c++11", "-w", "-O2"], - extra_objects=[], - depends=[], - language="c++", - ) -] - -if __name__ == "__main__": - setup( - name="cython_methods", - version="0.1", - url="https://github.com/pupil-labs/pupil", - author="Pupil Labs", - author_email="info@pupil-labs.com", - license="GNU", - ext_modules=cythonize(extensions, quiet=True, nthreads=8), - ) diff --git a/update_license_header.py b/update_license_header.py index dba116b7c0..8b11e34e55 100644 --- a/update_license_header.py +++ b/update_license_header.py @@ -45,7 +45,6 @@ "glfw.py", "version_utils.py", "update_license_header.py", - "shared_cpp*", ] # transform glob patterns to regular expressions