Torus.h

Go to the documentation of this file.

#ifndef TORUS_HEADER
#define TORUS_HEADER
#include <stdio.h>
 
#include <istream>
#include <ostream>
 
#include "PointCloud.h"
#include "basic.h"
#include <GfxTL/HyperplaneCoordinateSystem.h>
#include <GfxTL/MathHelper.h>
#include <GfxTL/VectorXD.h>
#include <MiscLib/NoShrinkVector.h>
#include <MiscLib/Vector.h>
 
#ifndef DLL_LINKAGE
#define DLL_LINKAGE
#endif
 
class DLL_LINKAGE Torus
{
public:
    enum
    {
        RequiredSamples = 4
    };
 
    bool Init(const MiscLib::Vector<Vec3f>& samples);
    bool InitAverage(const MiscLib::Vector<Vec3f>& samples);
    bool Init(bool binary, std::istream* i);
    void Init(FILE* i);
    void Init(float* array);
    inline float Distance(const Vec3f& p) const;
    inline void Normal(const Vec3f& p, Vec3f* n) const;
    inline float DistanceAndNormal(const Vec3f& p, Vec3f* n) const;
    inline float SignedDistance(const Vec3f& p) const;
    inline float SignedDistanceAndNormal(const Vec3f& p, Vec3f* n) const;
    inline void Project(const Vec3f& p, Vec3f* pp) const;
    void Transform(float scale, const Vec3f& translate);
 
    const Vec3f&
    Center() const
    {
        return m_center;
    }
 
    const Vec3f&
    AxisDirection() const
    {
        return m_normal;
    }
 
    const float
    MinorRadius() const
    {
        return m_rminor;
    }
 
    const float
    MajorRadius() const
    {
        return m_rmajor;
    }
 
    bool LeastSquaresFit(const PointCloud& pc,
                         MiscLib::Vector<size_t>::const_iterator begin,
                         MiscLib::Vector<size_t>::const_iterator end);
 
    bool
    Fit(const PointCloud& pc,
        MiscLib::Vector<size_t>::const_iterator begin,
        MiscLib::Vector<size_t>::const_iterator end)
    {
        return LeastSquaresFit(pc, begin, end);
    }
 
    bool
    IsAppleShaped() const
    {
        return m_appleShaped;
    }
 
    float
    AppleCutOffAngle() const
    {
        return m_cutOffAngle;
    }
 
    void Serialize(bool binary, std::ostream* o) const;
    static size_t SerializedSize();
    void Serialize(FILE* o) const;
    void Serialize(float* array) const;
    static size_t SerializedFloatSize();
 
private:
    void ComputeAppleParams();
 
private:
    Vec3f m_normal;
    Vec3f m_center;
    float m_rminor;
    float m_rmajor;
    bool m_appleShaped; // an apple shaped torus has rminor > rmajor
    float m_cutOffAngle; // for an apple shaped torus
    // the minor circle is cut off
    float m_appleHeight; // height of the "apple" point
};
 
float
Torus::Distance(const Vec3f& p) const
{
    Vec3f s = p - m_center;
    float spin1 = m_normal.dot(s);
    float spin0 = (s - spin1 * m_normal).length();
    spin0 -= m_rmajor;
    if (!m_appleShaped)
    {
        return abs(std::sqrt(spin0 * spin0 + spin1 * spin1) - m_rminor);
    }
    // apple shaped torus distance
    float minorAngle = std::atan2(spin1, spin0); // minor angle
    if (abs(minorAngle) < m_cutOffAngle)
    {
        return abs(std::sqrt(spin0 * spin0 + spin1 * spin1) - m_rminor);
    }
    spin0 += 2 * m_rmajor - m_rminor;
    if (minorAngle < 0)
    {
        spin1 += m_appleHeight;
    }
    else
    {
        spin1 -= m_appleHeight;
    }
    return std::sqrt(spin0 * spin0 + spin1 * spin1);
}
 
void
Torus::Normal(const Vec3f& p, Vec3f* n) const
{
    Vec3f s = p - m_center, tmp;
    float spin1 = m_normal.dot(s);
    float spin0 = (s - (tmp = spin1 * m_normal)).length();
    spin0 -= m_rmajor;
    if (m_appleShaped)
    {
        float minorAngle = std::atan2(spin1, spin0); // minor angle
        if (abs(minorAngle) > m_cutOffAngle)
        {
            *n = m_normal;
            if (minorAngle < 0)
            {
                *n *= -1;
            }
            return;
        }
    }
    Vec3f pln = s.cross(m_normal);
    Vec3f plx = m_normal.cross(pln);
    plx.normalize();
    *n = spin0 * plx + tmp;
    *n /= std::sqrt(spin0 * spin0 + spin1 * spin1);
}
 
float
Torus::DistanceAndNormal(const Vec3f& p, Vec3f* n) const
{
    Vec3f s = p - m_center, tmp;
    float spin1 = m_normal.dot(s);
    float spin0 = (s - (tmp = spin1 * m_normal)).length();
    spin0 -= m_rmajor;
    if (m_appleShaped)
    {
        float minorAngle = std::atan2(spin1, spin0); // minor angle
        if (abs(minorAngle) > m_cutOffAngle)
        {
            *n = m_normal;
            if (minorAngle < 0)
            {
                *n *= -1;
            }
            spin0 += 2 * m_rmajor - m_rminor;
            if (minorAngle < 0)
            {
                spin1 += m_appleHeight;
            }
            else
            {
                spin1 -= m_appleHeight;
            }
            return std::sqrt(spin0 * spin0 + spin1 * spin1);
        }
    }
    Vec3f pln = s.cross(m_normal);
    Vec3f plx = m_normal.cross(pln);
    plx.normalize();
    *n = spin0 * plx + tmp;
    float d = std::sqrt(spin0 * spin0 + spin1 * spin1);
    *n /= d;
    return abs(d - m_rminor);
}
 
float
Torus::SignedDistance(const Vec3f& p) const
{
    Vec3f s = p - m_center;
    float spin1 = m_normal.dot(s);
    float spin0 = (s - spin1 * m_normal).length();
    spin0 -= m_rmajor;
    if (!m_appleShaped)
    {
        return std::sqrt(spin0 * spin0 + spin1 * spin1) - m_rminor;
    }
    // apple shaped torus distance
    float minorAngle = std::atan2(spin1, spin0); // minor angle
    if (abs(minorAngle) < m_cutOffAngle)
    {
        return std::sqrt(spin0 * spin0 + spin1 * spin1) - m_rminor;
    }
    spin0 += 2 * m_rmajor - m_rminor;
    if (minorAngle < 0)
    {
        spin1 += m_appleHeight;
    }
    else
    {
        spin1 -= m_appleHeight;
    }
    return -std::sqrt(spin0 * spin0 + spin1 * spin1);
}
 
float
Torus::SignedDistanceAndNormal(const Vec3f& p, Vec3f* n) const
{
    Vec3f s = p - m_center, tmp;
    float spin1 = m_normal.dot(s);
    float spin0 = (s - (tmp = spin1 * m_normal)).length();
    spin0 -= m_rmajor;
    if (m_appleShaped)
    {
        float minorAngle = std::atan2(spin1, spin0); // minor angle
        if (abs(minorAngle) > m_cutOffAngle)
        {
            *n = m_normal;
            if (minorAngle < 0)
            {
                *n *= -1;
            }
            spin0 += 2 * m_rmajor - m_rminor;
            if (minorAngle < 0)
            {
                spin1 += m_appleHeight;
            }
            else
            {
                spin1 -= m_appleHeight;
            }
            return -std::sqrt(spin0 * spin0 + spin1 * spin1);
        }
    }
    Vec3f pln = s.cross(m_normal);
    Vec3f plx = m_normal.cross(pln);
    plx.normalize();
    *n = spin0 * plx + tmp;
    float d = std::sqrt(spin0 * spin0 + spin1 * spin1);
    *n /= d;
    return d - m_rminor;
}
 
void
Torus::Project(const Vec3f& p, Vec3f* pp) const
{
    Vec3f s = p - m_center, tmp;
    float spin1 = m_normal.dot(s);
    float spin0 = (s - (tmp = spin1 * m_normal)).length();
    spin0 -= m_rmajor;
    if (m_appleShaped)
    {
        float minorAngle = std::atan2(spin1, spin0); // minor angle
        if (abs(minorAngle) > m_cutOffAngle)
        {
            *pp = m_center + GfxTL::Math<float>::Sign(minorAngle) * m_normal;
            return;
        }
    }
    Vec3f pln = s.cross(m_normal);
    Vec3f plx = m_normal.cross(pln);
    plx.normalize();
    float d = std::sqrt(spin0 * spin0 + spin1 * spin1);
    *pp = m_center + (m_rminor / d) * (spin0 * plx + tmp) + m_rmajor * plx;
}
 
#endif