/*

Copyright 2018 Lingqi Yan

This file is part of WaveOpticsBrdf.

WaveOpticsBrdf 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 3 of the License, or
(at your option) any later version.

WaveOpticsBrdf 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 WaveOpticsBrdf.  If not, see <https://www.gnu.org/licenses/>.

*/

#ifndef HELPERS_H
#define HELPERS_H

#define M_PIf Float(M_PI)
#define M_2PIf 6.28318530718f
#define M_2PIf_inv 0.15915494309f

#include <cmath>
#include <complex>
#include <cstdlib>
#include <Eigen/Dense>

typedef float Float;
typedef Eigen::Vector2f Vector2;
typedef Eigen::Vector3f Vector3;
typedef std::complex<Float> comp;

const Float eps = 1e-6f;

inline Float dist(Float x1, Float y1, Float x2, Float y2) {
    return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}

inline Float G(Float x, Float mu, Float sigma) {
    return 1.0 / (sqrt(2.0 * M_PI) * sigma) * exp(-0.5 * pow((x - mu) / sigma, 2.0));
}

inline Float G(Vector2 x, Vector2 mu, Float sigma) {
    return G(x(0), mu(0), sigma) * G(x(1), mu(1), sigma);
}

inline Float GRot(Float x, Float mu, Float sigma, Float period) {
    return G(x - period, mu, sigma) + G(x, mu, sigma) + G(x + period, mu, sigma);
}

inline Float GRot(Vector2 x, Vector2 mu, Float sigma, Float period) {
    return GRot(x(0), mu(0), sigma, period) * GRot(x(1), mu(1), sigma, period);
}

// e^(ix) = cos x + i sin x
inline comp cis(Float x) {
    Float sinx, cosx;
    // sincosf(x, &sinx, &cosx);
    // return comp(cosx, sinx);
   return comp(std::cos(x), std::sin(x));
}

// e^(-ix) = cos x - i sin x
inline comp cnis(Float x) {
    Float sinx, cosx;
    // sincosf(x, &sinx, &cosx);
    // return comp(cosx, -sinx);
   return comp(std::cos(x), -std::sin(x));
}

template <class T>
inline T randUniform() {
    return rand() / (T) RAND_MAX;
}

template <class T>
inline T normrnd(T mu, T sigma) {
    T U = randUniform<T>();
    T V = randUniform<T>();
    T X = sqrt(-2.0 * log(U)) * cos(2.0 * M_PI * V);
    return X * sigma + mu;
}

inline Vector2 normrnd2D(Vector2 mu, Float sigma) {
    return Vector2(normrnd<Float>(mu(0), sigma), normrnd<Float>(mu(1), sigma));
}


struct AABB {
    Float xMin, xMax, yMin, yMax;
    AABB(Float xMin_ = Float(0.0f), Float xMax_ = Float(0.0f), Float yMin_ = Float(0.0f), Float yMax_ = Float(0.0f)) :
        xMin(xMin_), xMax(xMax_), yMin(yMin_), yMax(yMax_) {}
};

inline bool intersectAABB(const AABB &aabb1, const AABB &aabb2) {
    if (aabb1.xMax < aabb2.xMin)
        return false;
    if (aabb1.xMin > aabb2.xMax)
        return false;
    if (aabb1.yMax < aabb2.yMin)
        return false;
    if (aabb1.yMin > aabb2.yMax)
        return false;
    return true;
}

inline bool intersectAABBRot(const AABB &aabb1, const AABB &aabb2, Float period) {
    for (int i = -1; i <= 1; i++) {
        for (int j = -1; j <= 1; j++) {
            AABB aabb0(aabb1.xMin + i * period, aabb1.xMax + i * period,
                       aabb1.yMin + j * period, aabb1.yMax + j * period);
            if (intersectAABB(aabb0, aabb2))
                return true;
        }
    }
    return false;
}

inline bool insideAABB(const AABB &aabb, Float x, Float y) {
    if (x < aabb.xMin || x > aabb.xMax)
        return false;
    if (y < aabb.yMin || y > aabb.yMax)
        return false;
    return true;
}

inline bool insideAABB(const AABB &aabb, Vector2 xy) {
    return insideAABB(aabb, xy(0), xy(1));
}

inline AABB combineAABB(const AABB &aabb1, const AABB &aabb2) {
    AABB aabb;
    aabb.xMin = std::min(aabb1.xMin, aabb2.xMin);
    aabb.xMax = std::max(aabb1.xMax, aabb2.xMax);
    aabb.yMin = std::min(aabb1.yMin, aabb2.yMin);
    aabb.yMax = std::max(aabb1.yMax, aabb2.yMax);
    return aabb;
}

inline AABB combineAABB(const AABB &aabb1, const AABB &aabb2, const AABB &aabb3, const AABB &aabb4) {
    AABB aabbA = combineAABB(aabb1, aabb2);
    AABB aabbB = combineAABB(aabb3, aabb4);
    return combineAABB(aabbA, aabbB);
}

inline Float distSqrPeriod(Vector2 p1, Vector2 p2, Float period) {
    Float minDistX = std::abs(p1(0) - p2(0));
    minDistX = std::min(minDistX, std::abs(p1(0) + period - p2(0)));
    minDistX = std::min(minDistX, std::abs(p1(0) - period - p2(0)));
    Float minDistY = std::abs(p1(1) - p2(1));
    minDistY = std::min(minDistY, std::abs(p1(1) + period - p2(1)));
    minDistY = std::min(minDistY, std::abs(p1(1) - period - p2(1)));
    return minDistX * minDistX + minDistY * minDistY;
}

#endif