ConvHull, C++ (gcc)

ConvHull

//
//  main.cpp
//  ConvHull
//
//  Created by Nik on 07/11/2018.
//
//

#include <algorithm>
#include <cmath>
#include <iostream>
#include <queue>
#include <string>
#include <functional>
#include <set>
#include <unordered_map>
#include <vector>

struct point;

struct vector3d {
public:
    friend class point;
    vector3d (double x, double y, double z): x(x), y(y), z(z) {}
    double operator * (const vector3d& otherVector) {
        return x * otherVector.x + y * otherVector.y + z * otherVector.z;
    }
    vector3d cross(vector3d secondVector) {
        double newX = y * secondVector.z - z * secondVector.y;
        double newY = z * secondVector.x - x * secondVector.z;
        double newZ = x * secondVector.y - y * secondVector.x;
        return vector3d(newX, newY, newZ);
    }
    inline double norm() const {
        return x * x + y * y + z * z;
    }
    inline double length() const {
        return sqrt(norm());
    }
    inline vector3d operator * (double num) {
        return vector3d(x * num, y * num, z * num);
    }
    inline vector3d operator / (double num) const {
        return vector3d(x / num, y / num, z / num);
    }
    inline vector3d operator - (const vector3d& secondVector) const {
        return vector3d(secondVector.x - x, secondVector.y - y, secondVector.z -z);
    }
    inline double findAngleWith(const vector3d& secondVector) const {
        vector3d tmpVector1 = *this / length();
        vector3d tmpVector2 = secondVector / secondVector.length();
        return -(tmpVector1 * tmpVector2);
    }
    vector3d operator = (const vector3d& other) {
        x = other.x;
        y = other.y;
        z = other.z;
        return *this;
    }
    inline double triple (const vector3d& fVec, const vector3d& sVec) {
        return x * (fVec.y * sVec.z - fVec.z * sVec.y) + y * (fVec.z * sVec.x - fVec.x * sVec.z) + z * (fVec.x * sVec.y - fVec.y * sVec.x);
    }
    
    double x, y, z;
    
};

struct point {
    double x;
    double y;
    double z;
    point (double x, double y, double z): x(x), y(y) , z(z){}
    inline double distanceTo(point otherPoint) {
        return sqrt((x - otherPoint.x) * (x - otherPoint.x) + (y - otherPoint.y) * (y - otherPoint.y) + (z - otherPoint.z) * (z - otherPoint.z));
    }
    inline vector3d operator - (const point& otherPoint) const {
        return vector3d(x - otherPoint.x, y - otherPoint.y, z - otherPoint.z);
    }
    inline point operator - (const vector3d& vec) {
        return point(x - vec.x, y - vec.y, z - vec.z);
    }
    point operator = (const point& point) {
        x = point.x;
        y = point.y;
        z = point.z;
        return *this;
    }
    inline point operator + (const vector3d& vec) const {
        return point(x+vec.x, y+vec.y, z+vec.z);
    }
    inline bool operator == (const point otherPoint) const {
        return (x == otherPoint.x && y == otherPoint.y && z == otherPoint.z);
    }
    inline bool operator < (const point& otherPoint) const {
        return z < otherPoint.z;
    }
    inline vector3d findNorm (const point& secondPoint, const point& thirdPoint) {
        vector3d firstVector = secondPoint - *this;
        vector3d secondVector = thirdPoint - *this;
        return firstVector.cross(secondVector);
    }
    inline bool operator != (const point& otherPoint){
        return !(*this == otherPoint);
    }
    
    inline point operator + (const point& otherPoint) const {
        double tmpX = x + otherPoint.x;
        double tmpY = y + otherPoint.y;
        double tmpZ = z + otherPoint.z;
        return point(tmpX, tmpY, tmpZ);
    }
};

struct facet {
    using it = std::vector<point>::iterator;
    it fPoint, sPoint, tPoint;
    bool operator < (const facet& other) const {
        if (fPoint != other.fPoint) {
            return fPoint < other.fPoint;
        }
        else if (sPoint != other.sPoint) {
            return sPoint < other.sPoint;
        }
        else {
            return tPoint < other.tPoint;
        }
    }
    facet(it fPoint, it sPoint, it tPoint): fPoint(fPoint), sPoint(sPoint), tPoint(tPoint) {}
    bool operator == (const facet& other) const {
        return (fPoint == other.fPoint && sPoint == other.sPoint && tPoint == other.tPoint);
    }
};

namespace std {
    template <>
    class hash<point> {
    public:
        size_t operator()(const point& pointToHash) const
        {
            auto tmpString = std::to_string(pointToHash.x) + '$' + std::to_string(pointToHash.y)
            + '$' + std::to_string(pointToHash.z);
            std::hash<string> ourHash;
            return ourHash(tmpString);
            
        }
    };
}

struct triangle {
    point fPoint, sPoint, tPoint;
    triangle (point fPoint, point sPoint, point tPoint): fPoint(fPoint), sPoint(sPoint), tPoint(tPoint) {}
    
};

class ConvHull {
public:
    
    ConvHull(std::istream& inputStream) {
        double tmpX, tmpY, tmpZ;
        size_t s;
        
        inputStream >> s;
        for (size_t j = 0; j < s; ++j) {
            inputStream >> tmpX >> tmpY >> tmpZ;
            InputPoints.emplace_back(tmpX, tmpY, tmpZ);
        }
        for (size_t i = 0; i < InputPoints.size(); ++i) {
            posInInp[InputPoints[i]] = i;
        }
        
        preprocessing();
        iterations();
        
    }
    void printTo(std::ostream& outputStream) {
        std::vector<size_t> bufArray(3);
        outputStream << HullsFacets.size() << std::endl;
        for (auto&& triang : HullsFacets) {
            bufArray[0] = posInInp[*triang.fPoint];
            bufArray[1] = posInInp[*triang.sPoint];
            bufArray[2] = posInInp[*triang.tPoint];
            outputStream <<"3 " << bufArray[0] << " " << bufArray[1] << " " << bufArray[2] << std::endl;
        }
    }
private:
    std::unordered_map<point, size_t> posInInp;
    std::vector<point> InputPoints;
    std::vector<point> shiftedPoints;
    std::vector<point> hull;
    std::set<facet> HullsFacets;
    point center = point(0,0,0);
    
    
    void preprocessing() {
        auto lowestPoint = *std::min_element(InputPoints.begin(), InputPoints.end());
        hull.push_back(point(lowestPoint));
        std::vector<double> tmpAngles;
        point secondPoint = lowestPoint + point(0,1,0);
        bool tmpFlag = false;
        for (auto&& curPoint :  InputPoints) {
            if (curPoint != lowestPoint) {
                if (curPoint.z == lowestPoint.z && lowestPoint.y < curPoint.y && curPoint.x == lowestPoint.x) {
                    secondPoint = curPoint;
                    tmpFlag = true;
                }
            }
        }
        point thirdPoint(0,0,0);
        if (tmpFlag) {
            thirdPoint = findThirdPoint(lowestPoint, secondPoint);
        }
        else {
            secondPoint = findThirdPoint(lowestPoint, secondPoint);
            thirdPoint = nextPoint(lowestPoint, secondPoint, lowestPoint + point(0,1,0));
        }
        hull.push_back(secondPoint);
        hull.push_back(thirdPoint);
    }
    
    
    void permutateFacet(facet& inpFacet) {
        if (inpFacet.fPoint > inpFacet.sPoint) {
            std::swap(inpFacet.fPoint, inpFacet.sPoint);
        }
        if (inpFacet.fPoint > inpFacet.tPoint) {
            std::swap(inpFacet.fPoint, inpFacet.tPoint);
        }
        
        vector3d dir = center - *(inpFacet.fPoint);
        if (dir.triple(*inpFacet.sPoint - *inpFacet.fPoint, *inpFacet.tPoint - *inpFacet.fPoint) > 0) {
            std::swap(inpFacet.sPoint, inpFacet.tPoint);
        }
    }
    
    
    void iterations(){
        std::queue<facet> triangles;
        auto firstNextIt = InputPoints.begin() + posInInp[hull[0]];
        auto secondNextIt = InputPoints.begin() + posInInp[hull[1]];
        auto thirdNextIt = InputPoints.begin() + posInInp[hull[2]];
        triangles.emplace(secondNextIt, firstNextIt, thirdNextIt);
        triangles.emplace(thirdNextIt, secondNextIt, firstNextIt);
        triangles.emplace(firstNextIt, thirdNextIt, secondNextIt);
        for (auto&& curPoint : InputPoints) {
            center = curPoint + center;
        }
        center.x /= InputPoints.size();
        center.y /= InputPoints.size();
        center.z /= InputPoints.size();
        size_t tmpPos1 = posInInp[hull[0]];
        size_t tmpPos2 = posInInp[hull[1]];
        size_t tmpPos3 = posInInp[hull[2]];
        facet firstFacet(InputPoints.begin() + tmpPos1, InputPoints.begin() + tmpPos2, InputPoints.begin() + tmpPos3);
        permutateFacet(firstFacet);
        HullsFacets.insert(firstFacet);
        while (!triangles.empty()) {
            auto curF = triangles.front();
            triangles.pop();
            point nextP = nextPoint(*curF.fPoint, *curF.sPoint, *curF.tPoint);
            firstNextIt = InputPoints.begin() + posInInp[*curF.fPoint];
            secondNextIt = InputPoints.begin() + posInInp[*curF.sPoint];
            thirdNextIt = InputPoints.begin() + posInInp[nextP];
            facet nextF(firstNextIt, secondNextIt, thirdNextIt);
            permutateFacet(nextF);
            bool is_in = HullsFacets.find(nextF) != HullsFacets.end();
            if (!is_in) {
                triangles.emplace(nextF.tPoint, nextF.sPoint, nextF.fPoint);
                triangles.emplace(nextF.fPoint, nextF.tPoint, nextF.sPoint);
                HullsFacets.insert(nextF);
            }
            
        }
    }
    
    
    point nextPoint(const point& fPoint, const point& sPoint, const point& tPoint) {
        double minAngle = 2;
        point candidate(0,0,0);
        auto mainNorm = (tPoint - sPoint).cross(tPoint - fPoint);  // fix me
        for(auto&& curPoint : InputPoints) {
            if(curPoint != fPoint && curPoint != sPoint && curPoint != tPoint) {
                auto normal = (curPoint - fPoint).cross(curPoint - sPoint);
                if (mainNorm.findAngleWith(normal) < minAngle) {
                    minAngle = mainNorm.findAngleWith(normal); // check angle
                    candidate = curPoint;
                }
            }
        }
        return candidate;
    }
    
    
    point findThirdPoint(point& firstPoint, point& secondPoint) {
        double minAngle = 2;
        point candidate = InputPoints[0];
        for (auto&& curPoint : InputPoints) {
            if (curPoint != firstPoint && curPoint != secondPoint) {
                vector3d proj(curPoint.x, curPoint.y, firstPoint.z);
                vector3d firstVec(proj.x - firstPoint.x, proj.y - firstPoint.y, proj.z - firstPoint.z);
                vector3d secondVec(proj.x - secondPoint.x, proj.y - secondPoint.y, proj.z - secondPoint.z);
                auto normal1 = firstVec.cross(secondVec);
                vector3d thirdVec(curPoint.x - firstPoint.x, curPoint.y - firstPoint.y, curPoint.z - firstPoint.z);
                vector3d fourthVec(curPoint.x - secondPoint.x, curPoint.y - secondPoint.y, curPoint.z - secondPoint.z);
                auto normal2 = thirdVec.cross(fourthVec);
                if (normal1.findAngleWith(normal2) <= minAngle) {
                    minAngle = normal1.findAngleWith(normal2);
                    candidate = curPoint;
                }
            }
        }
        return candidate;
    }
    
};

int main(int argc, const char * argv[]) {
    size_t n;
    std:: cin >> n;
    for (int i = 0; i < n; ++i) {
        ConvHull ourHull(std::cin);
        ourHull.printTo(std::cout);
    }
    return 0;
}


run \| edit \| history \| help	0