geometry.h

Go to the documentation of this file.

/***************************************************************************
*
* Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
*
* Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
* 02111-1307  USA
*
*  All comments concerning this program package may be sent to the
*  e-mail address 'xmipp@cnb.csic.es'
***************************************************************************/

#ifndef CORE_GEOMETRY_H
#define CORE_GEOMETRY_H

#include "multidim_array.h"
#include "multidim_array_generic.h"

#ifndef FLT_EPSILON
#define FLT_EPSILON 1.19209e-07
#endif

#include <vector>
#include <iostream>



void Uproject_to_plane(const Matrix1D< double >& point,
                       const Matrix1D< double >& direction,
                       double distance,
                       Matrix1D< double >& result);

void Uproject_to_plane(const Matrix1D< double >& r,
                       double rot,
                       double tilt,
                       double psi,
                       Matrix1D< double >& result);

void Uproject_to_plane(const Matrix1D< double >& r,
                       const Matrix2D< double >& euler,
                       Matrix1D< double >& result);

double spherical_distance(const Matrix1D< double >& r1,
                          const Matrix1D< double >& r2);

double point_line_distance_3D(const Matrix1D< double >& p,
                              const Matrix1D< double >& a,
                              const Matrix1D< double >& v);

inline double point_plane_distance_3D(const Matrix1D< double >& p,
                                      const Matrix1D< double >& a,
                                      const Matrix1D< double >& v)
{
    static Matrix1D< double > p_a(3);

    V3_MINUS_V3(p_a, p, a);
    return (dotProduct(p_a, v) / v.module());
}

struct FitPoint
{
    double x;
    double y;
    double z;
    double w;
};

void least_squares_plane_fit(FitPoint *IN_points,
                             int Npoints,
                             double& plane_A,
                             double& plane_B,
                             double& plane_C);

void least_squares_plane_fit_All_Points(const MultidimArray<double> &Image,
                             double& plane_A,
                             double& plane_B,
                             double& plane_C);

struct fit_point2D
{
    double x;
    double y;
    double w;
};

void least_squares_line_fit(const std::vector< fit_point2D >& IN_points,
                            double& line_A,
                            double& line_B);

class Bspline_model
{
public:
    int l0;
    int lF;
    int m0;
    int mF;

    double x0;
    double y0;

    int SplineDegree;

    double h_x;
    double h_y;

    MultidimArray< double > c_ml;

    double evaluate(double x, double y) const;
};

void Bspline_model_fitting(const std::vector< FitPoint >& IN_points,
                           int SplineDegree,
                           int l0,
                           int lF,
                           int m0,
                           int mF,
                           double h_x,
                           double h_y,
                           double x0,
                           double y0,
                           Bspline_model& result);

void rectangle_enclosing(const Matrix1D< double >& v0,
                         const Matrix1D< double >& vF,
                         const Matrix2D< double >& V,
                         Matrix1D< double >& corner1,
                         Matrix1D< double >& corner2);

void box_enclosing(const Matrix1D< double >& v0,
                   const Matrix1D< double >& vF,
                   const Matrix2D< double >& V,
                   Matrix1D< double >& corner1,
                   Matrix1D< double >& corner2);

bool point_inside_polygon(const std::vector< Matrix1D< double > > & polygon,
                          const Matrix1D< double >& point);

void def_affinity(double u1x, double u1y, double u2x, double u2y, double u3x, double u3y, double t1x,
        double t1y, double t2x, double t2y, double t3x, double t3y, Matrix2D<double> &A, Matrix1D<double> &T, Matrix2D<double> &invW);



double triangle_area(double x1, double y1, double x2, double y2, double x3, double y3);

int line_plane_intersection(const Matrix1D< double > normal_plane,
                            const Matrix1D< double > vector_line,
                            Matrix1D< double >& intersection_point,
                            const Matrix1D< double > point_line,
                            double point_plane_at_x_y_zero = 0.);


#define EULER_CLIPPING(rot,tilt,psi) \
    rot = realWRAP(rot, 0, 360); \
    tilt = realWRAP(tilt, 0, 360); \
    psi = realWRAP(psi, 0, 360);

#define EULER_CLIPPING_RAD(rot, tilt, psi) \
    rot = realWRAP(rot, 0, 2.0*PI); \
    tilt = realWRAP(tilt, 0, 2.0*PI); \
    psi = realWRAP(psi, 0, 2.0*PI);

template<typename T>
void Euler_angles2matrix(T a, T b, T g, Matrix2D<T>& A,
                         bool homogeneous=false);

void Euler_anglesZXZ2matrix(double a, double b, double g, Matrix2D< double >& A,
                         bool homogeneous=false);

double Euler_distanceBetweenMatrices(const Matrix2D<double> &E1,
                                     const Matrix2D<double> &E2);

template<typename T>
T Euler_distanceBetweenAngleSets(T rot1, T tilt1, T psi1,
                                      T rot2, T tilt2, T psi2,
                                      bool only_projdir);

template<typename T>
T Euler_distanceBetweenAngleSets_fast(const Matrix2D<T> &E1,
                                           T rot2, T tilt2, T psi2,
                                           bool only_projdir, Matrix2D<T> &E2);

void Euler_Angles_after_compresion(const double rot,
                                   double tilt,
                                   double psi,
                                   double& new_rot,
                                   double& new_tilt,
                                   double& new_psi,
                                   Matrix2D< double >& D);

void Euler_direction(double alpha,
                     double beta,
                     double gamma,
                     Matrix1D< double >& v);

void Euler_direction2angles(Matrix1D< double >& v,
                            double& alpha,
                            double& beta,
                            double& gamma);

void Euler_matrix2angles(const Matrix2D< double >& A,
                         double& alpha,
                         double& beta,
                         double& gamma,bool homogeneous=false);

void Euler_up_down(double rot,
                   double tilt,
                   double psi,
                   double& newrot,
                   double& newtilt,
                   double& newpsi);

void Euler_another_set(double rot,
                       double tilt,
                       double psi,
                       double& newrot,
                       double& newtilt,
                       double& newpsi);

void Euler_mirrorY(double rot,
                   double tilt,
                   double psi,
                   double& newrot,
                   double& newtilt,
                   double& newpsi);

void Euler_mirrorX(double rot,
                   double tilt,
                   double psi,
                   double& newrot,
                   double& newtilt,
                   double& newpsi);

void Euler_mirrorXY(double rot,
                    double tilt,
                    double psi,
                    double& newrot,
                    double& newtilt,
                    double& newpsi);

void Euler_apply_transf(const Matrix2D< double >& L,
                        const Matrix2D< double >& R,
                        double rot,
                        double tilt,
                        double psi,
                        double& newrot,
                        double& newtilt,
                        double& newpsi);

void Euler_rotate(const MultidimArray< double >& V,
                  double rot,
                  double tilt,
                  double psi,
                  MultidimArray< double >& result);

void Euler_rotate(const MultidimArrayGeneric &V,
                  double rot,
                  double tilt,
                  double psi,
                  MultidimArray<double> &result);
void computeCircleAroundE(const Matrix2D<double> &E,
                          double angCircle, double angStep, std::vector<double> &outputEulerAngles);


double intersection_unit_sphere(const Matrix1D< double >& u,
                                const Matrix1D< double >& r);

double intersection_unit_cylinder(const Matrix1D< double >& u,
                                  const Matrix1D< double >& r);

double intersection_unit_cube(const Matrix1D< double >& u,
                              const Matrix1D< double >& r);

#endif