/* 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 . */ #ifndef HELPERS_H #define HELPERS_H #include #include #include #include #ifndef M_PI #define M_PI 3.14159265358979323 #endif typedef float Float; typedef Eigen::Vector2f Vector2; typedef Eigen::Vector3f Vector3; typedef std::complex 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 inline T randUniform() { return rand() / (T) RAND_MAX; } template inline T normrnd(T mu, T sigma) { T U = randUniform(); T V = randUniform(); 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(mu(0), sigma), normrnd(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