Merge pull request #82 from liuyun-0002/glint

feature: Add glint material

src/CoreLayer/Geometry/BoundingBox.h

@@ -116,3 +116,84 @@ static BoundingBox3<BaseType> BoundingBoxPointUnion(const BoundingBox3<BaseType>
 }
 
 using BoundingBox3f = BoundingBox3<double>;
+
+//Axis-aligned bounding box 2D base type
+template <typename BaseType>
+class BoundingBox2
+{
+public:
+	TPoint2<BaseType> pMin;			
+	TPoint2<BaseType> pMax;
+
+	BoundingBox2(){
+		pMin[0] = pMin[1] = std::numeric_limits<BaseType>::max();
+		pMax[0] = pMax[1] = std::numeric_limits<BaseType>::lowest();
+	}
+
+	BoundingBox2(TPoint2<BaseType> _p)
+	{
+		pMin = pMax = _p;
+	}
+
+	BoundingBox2(TPoint2<BaseType> _pMin, TPoint2<BaseType> _pMax)
+	{
+		pMin = _pMin;
+		pMax = _pMax;
+
+		// corner case: bounding box not exist
+		if(pMin.x > pMax.x || pMin.y > pMax.y){
+			pMin[0] = pMin[1] = std::numeric_limits<BaseType>::max();
+			pMax[0] = pMax[1] = std::numeric_limits<BaseType>::lowest();
+		}
+	}
+
+	// caculate area
+	inline double getSurfaceArea(){
+		double x = pMax[0] - pMin[0];
+		double y = pMax[1] - pMin[1];
+
+		return x * y;
+	}
+
+	// judge if overlap
+	inline bool overlap(BoundingBox2 &box){
+		if(pMin[0] > box.pMax[0] || pMax[0] < box.pMin[0])
+			return false;
+		if(pMin[1] > box.pMax[1] || pMax[1] < box.pMin[1])
+			return false;
+		return true;
+	}
+
+	// judge if contains some box
+	inline bool contains(BoundingBox2 &box){
+		if(pMin[0] > box.pMin[0] || pMax[0] < box.pMax[0])
+			return false;
+		if(pMin[1] > box.pMin[1] || pMax[1] < box.pMax[1])
+			return false;
+		return true;
+	}
+
+	inline TPoint2<BaseType> getCenter(){
+		return TPoint2<BaseType>((pMin[0] + pMax[0]) / 2, (pMin[1] + pMax[1]) / 2);
+	}
+};
+
+template <typename BaseType>
+static BoundingBox2<BaseType> BoundingBoxUnionIntersect(const BoundingBox2<BaseType>& b1, 
+														const BoundingBox2<BaseType>& b2)
+{
+	TPoint2<BaseType> _pMin(std::max(b1.pMin[0], b2.pMin[0]), std::max(b1.pMin[1], b2.pMin[1]));
+	TPoint2<BaseType> _pMax(std::min(b1.pMax[0], b2.pMax[0]), std::min(b1.pMax[1], b2.pMax[1]));
+	return BoundingBox2(_pMin, _pMax);
+}
+
+template <typename BaseType>
+static BoundingBox2<BaseType> BoundingBoxPointUnion(const BoundingBox2<BaseType>& b1,
+                                                    const TPoint2<BaseType>& p)
+{
+    TPoint2<BaseType> _pMin(std::min(b1.pMin[0], p[0]), std::min(b1.pMin[1], p[1]));
+    TPoint2<BaseType> _pMax(std::max(b1.pMax[0], p[0]), std::max(b1.pMax[1], p[1]));
+    return BoundingBox2(_pMin, _pMax);
+}
+
+using BoundingBox2d = BoundingBox2<double>;

src/CoreLayer/Ray/Ray.cpp

@@ -15,6 +15,7 @@
 Ray::Ray(const Point3d &_origin, const Vec3d &_direction, double _timeMin, double _timeMax) : origin(_origin),
                                                                                               direction(_direction)
 {
+    hasDifferential = false;    
 	timeMin = _timeMin;
 	timeMax = _timeMax;
 }

src/CoreLayer/Ray/Ray.h

@@ -22,6 +22,11 @@ struct Ray
 	Point3d origin;		///< The origin of the ray
 	Vec3d direction;	///< The direction of the ray
 
+	// ray differential
+	bool hasDifferential; 
+	Point3d origin_x, origin_y;
+	Vec3d direction_x, direction_y;
+
 	double timeMin;		///< The time when the shutter open
 	double timeMax;		///< The time when the shutter close
 	Point3d at(double t) const ;

src/FunctionLayer/Camera/CameraFactory.cpp

@@ -1,6 +1,7 @@
 #include "CameraFactory.h"
 #include "Camera.h"
 #include "Pinhole.h"
+#include "DifferentialPinhole.h"
 
 namespace CameraFactory {
 
@@ -10,7 +11,9 @@ std::shared_ptr<T> cameraCreatorHelper(const Json &json) {
 }
 
 static std::unordered_map<std::string, std::function<std::shared_ptr<Camera>(const Json &)>> cameraMap{
-    {"pinhole", cameraCreatorHelper<PinholeCamera>}};
+    {"pinhole", cameraCreatorHelper<PinholeCamera>},
+    {"differential_pinhole", cameraCreatorHelper<DifferentialPinholeCamera>}
+};
 
 std::shared_ptr<Camera> LoadCameraFromJson(const Json &json) {
     return cameraMap[json["type"]](json);

src/FunctionLayer/Camera/DifferentialPinhole.cpp

@@ -0,0 +1,27 @@
+#include"DifferentialPinhole.h"
+
+Ray DifferentialPinholeCamera::generateRay(const Point2i &filmResolution, const Point2i &pixelPosition, const CameraSample &sample) const {
+    double x = (double)(pixelPosition.x + sample.xy.x) / filmResolution.x,
+           y = (double)(pixelPosition.y + sample.xy.y) / filmResolution.y;
+
+    Point3d pointOnFilm = sampleToFilm * Point3d {x, y, 0},
+            origin = Point3d {0, 0, 0};
+
+    Vec3d dir = cameraToWorld * normalize(pointOnFilm - origin);
+    Vec3d dirX = cameraToWorld * normalize((sampleToFilm * Point3d {x + 1.f, y, 0}) - origin);
+    Vec3d dirY = cameraToWorld * normalize((sampleToFilm * Point3d {x, y + 1.f, 0}) - origin);
+
+    origin = cameraToWorld * origin;
+
+    Ray ret = Ray(origin, dir);
+
+    ret.hasDifferential = true;
+    ret.direction_x = dirX;
+    ret.direction_y = dirY;
+    ret.origin_x = ret.origin_y = origin;
+
+    return ret;
+}
+
+DifferentialPinholeCamera::DifferentialPinholeCamera(const Json & json) : PerspectiveCamera(json) {
+}

src/FunctionLayer/Camera/DifferentialPinhole.h

@@ -0,0 +1,24 @@
+#pragma once
+
+#include "Perspective.h"
+#include "CoreLayer/Adapter/JsonUtil.h"
+
+// pinhole with differential
+class DifferentialPinholeCamera : public PerspectiveCamera {
+public:
+    DifferentialPinholeCamera() = default;
+    DifferentialPinholeCamera(
+        const Point3d &lookFrom, 
+        const Point3d &lookAt, 
+        const Vec3d &up,
+        double xFov,
+        double aspectRatio,
+        double distToFilm
+    ) : PerspectiveCamera(lookFrom, lookAt, up, xFov, aspectRatio, distToFilm) { }
+
+    DifferentialPinholeCamera(const Json & json);
+
+    virtual Ray generateRay(const Point2i &filmResolution, 
+                            const Point2i &pixelPosition, 
+                            const CameraSample &sample) const override;
+};

src/FunctionLayer/Integrator/AbstractPathIntegrator.cpp

@@ -66,6 +66,9 @@ Spectrum AbstractPathIntegrator::Li(const Ray &initialRay, std::shared_ptr<Scene
         }
 
         Intersection its = itsOpt.value();
+        if(ray.hasDifferential){
+            its.computeRayDifferential(ray);
+        }
 
         // RUSSIAN ROULETTE
         nBounce++;

src/FunctionLayer/Intersection.h

@@ -36,6 +36,10 @@ struct Intersection
 	Vec3d dpdu, dpdv;
 	Normal3d dndu, dndv;
 
+	// raydifferential
+	float dudx, dvdx, dudy, dvdy;
+	Vec3d dpdx, dpdy;
+
 	// std::shared_ptr<Entity> object;
 	const Entity *object;
 
@@ -54,4 +58,68 @@ struct Intersection
 	{
 		return shFrame.toWorld(d);
 	}
+
+	// compute differential, reference implementation in lite
+	void computeRayDifferential(const Ray &ray){
+
+		dudx = 0.f;
+		dudy = 0.f;
+		dvdx = 0.f;
+		dvdy = 1.f;
+		dpdx = Vec3d(0.f);
+		dpdy = Vec3d(0.f);
+
+		Point3d p = position;
+		Vec3d n = geometryNormal;
+		Vec3d ox = Vec3d{ray.origin_x[0], ray.origin_x[1], ray.origin_x[2]};
+		Vec3d oy = Vec3d{ray.origin_y[0], ray.origin_y[1], ray.origin_y[2]};
+
+		double distance = dot(n, Vec3d(p[0], p[1], p[2]));
+		double tx = (distance - dot(n, ox)) / (dot(n, ray.direction_x));
+		if(std::isinf(tx) || std::isnan(tx))
+			return;
+		double ty = (distance - dot(n, oy)) / (dot(n, ray.direction_y));
+		if(std::isinf(ty) || std::isnan(ty))
+			return;
+
+		dpdx = normalize((ray.origin_x + tx * ray.direction_x) - p);
+		dpdy = normalize((ray.origin_y + ty * ray.direction_y) - p);
+
+		int dim[2];
+		if (std::abs(n[0]) > std::abs(n[1]) && std::abs(n[0]) > std::abs(n[2])) {
+			dim[0] = 1;
+			dim[1] = 2;
+		} else if (std::abs(n[1]) > std::abs(n[2])) {
+			dim[0] = 0;
+			dim[1] = 2;
+		} else {
+			dim[0] = 0;
+			dim[1] = 1;
+		}
+
+		float A[2][2] = {{dpdu[dim[0]], dpdv[dim[0]]},
+						{dpdu[dim[1]], dpdv[dim[1]]}};
+		float Bx[2] = {dpdx[dim[0]], dpdx[dim[1]]};
+		float By[2] = {dpdy[dim[0]], dpdy[dim[1]]};
+
+		auto solveLinearSystem2x2 = [](const float A[2][2], const float B[2],
+										float *x0, float *x1) {
+			float det = A[0][0] * A[1][1] - A[0][1] * A[1][0];
+			if (std::abs(det) < 1e-10f)
+				return false;
+			*x0 = (A[1][1] * B[0] - A[0][1] * B[1]) / det;
+			*x1 = (A[0][0] * B[1] - A[1][0] * B[0]) / det;
+			if (std::isnan(*x0) || std::isnan(*x1))
+				return false;
+			return true;
+		};
+
+		if (!solveLinearSystem2x2(A, Bx, &dudx, &dvdx))
+			dudx = dvdx = .0f;
+		if (!solveLinearSystem2x2(A, Bx, &dudy, &dvdy))
+			dudy = dvdy = .0f;
+
+		return;
+
+	}
 };