Hard work on Curved 3D tracer.

Notes/formulae.md

@@ -87,4 +87,47 @@ $$
 \vec{r} = \vec{v} - 2\vec{n}\langle \vec{n}, \vec{v}\rangle
 $$
 
-If we wish to compute the reflection *off* the surface, then we simply mirror the incoming vector pointing towards the normal. So this formula is a *mirror* reflection when $\vec{v}$ points in the same direction as the normal, and a *surface* reflection if $\vec{v}$ points in the opposite direction.
+If we wish to compute the reflection *off* the surface, then we simply mirror the incoming vector pointing towards the normal. So this formula is a *mirror* reflection when $\vec{v}$ points in the same direction as the normal, and a *surface* reflection if $\vec{v}$ points in the opposite direction.
+
+## 4. Raytracer3D overview
+
+The rays are cast by the fragment shader, which computes the ray direction using the projection and view matrices. NDC Coordinates are inverted to local space by multiplying with the inverse of the model-view-projection matrix.
+
+```mermaid
+graph LR;
+3dm["3D Model"]-->vs["Vertex Shader"]--NDC coords [-1, 1]-->fs["Fragment Shader"]
+```
+
+The fragment shader is where the raytracing happens. In general for raytracing, the essence is to
+
+1. Compute the local space ray origin $\vec{o}$ and direction $\vec{d}$ from NDC coordinates.
+2. For each reflection (numBounces) do the following:
+2.1 Determine which mirror has been collided with. This can be optimized but we simply iterate every mirror.
+2.2. Given the collided mirror and collision point, determine if there is an edge intersection. If so, stop and give color. This edge intersection test requires a distance check for each edge (straight or curved) which adds up to the cost of the algorithm. Edge visualising is not mandatory and can be substituted for corner visualisation or even single prop visualiation.
+2.3. Reflect the ray in the normal and set new origin.
+3. Set the fragment color.
+
+```mermaid
+graph TD
+Ray-->rayloop["For each reflection"]-->mloop["For each mirror"]-->det["Determine collision step"]
+det--Set new closest-->mloop
+mloop--Closest hit found-->edgeloop["For each edge"]-->detdi["Determine hit distance to edge"]
+detdi--Edge collision-->break["Break loops and set color"]
+edgeloop--No edge collision-->compref["Set new ray origin and direction"]
+compref--Next reflection-->rayloop
+rayloop--End of loop-->black["Set background color (Black)"]
+```
+
+## 5. Determine If Point In Polygon
+
+For the 3D curved mirrors we use a region of a circle as the surface. We must be able to determine of a collision occured inside that surface, which is determine by a polygon on a sphere. So given a point and a polygon, we have to determine if the point is inside. 
+
+There are two ways to do this
+
+### Convex Polygons
+
+Here we determine the winding number as a sum of internal angles of the vertices to the point. The sum of internal angles must be 360.
+
+### Any Polygon
+
+Here we must trace a ray from the point to any constant direction. The number of intersections with the polygon must be an odd number. This is more costly.

Shaders/ray3.frag

@@ -3,15 +3,23 @@
 // Output fragment color
 out vec4 finalColor;
 
+struct MirrorInfo {
+    vec4 normal;
+    vec4 center;
+    int offset;
+    int vertexCount;
+    float p2, p3;
+};
+
 // First buffer for vertices.
 // Second buffer for polygon counts.
 layout(std430, binding=0) buffer ssbo0 { vec4 vertices[]; };
-layout(std430, binding=1) buffer ssbo1 { vec4 normals[]; };
-layout(std430, binding=2) buffer ssbo2 { int sizes[]; };
+layout(std430, binding=1) buffer ssbo1 { MirrorInfo mirrors[]; };
 
 layout(location=0) uniform int numMirrors;
 layout(location=1) uniform float edgeThickness;
 layout(location=2) uniform mat4 invMvp;
+layout(location=3) uniform float sphereFocus;
 uniform mat4 mvp;
 
 const float STEP_MIN = 0.001f;
@@ -45,59 +53,107 @@ void main()
 
     // Iterate all mirrors and keep track of individual
     // mirror offset indices in the main vertex buffer.
-    for (int b = 0; b < 16; b++) {
-        int offset = 0;
-        int mirrorIndex = 0;
-        int mirrorOffset = 0;
-        int mirrorSize = 0;
-        HitInfo mirrorHit;
-        mirrorHit.t = -1;
+    for (int b = 0; b < 1; b++) {
+        // We keep track of the closest mirror hit.
+        HitInfo closestMirrorHit;
+        closestMirrorHit.t = -1;
+        MirrorInfo closestMirror;
+        float closestEdgeProj = -1;
 
         // Find the nearest plane to intersect with valid normal.
-        for (int i = 0; i < numMirrors; i++) {
-            int size = sizes[i];
-            vec3 normal = normals[i].xyz;
-            vec3 surfacePoint = vertices[offset].xyz;
-            HitInfo hit = RayPlaneHit(rayOrigin, rayDir, normal, surfacePoint);
-
-            bool validity = dot(normal, rayDir) < 0 && hit.t > STEP_MIN;
-            bool optimality = mirrorHit.t < 0 || hit.t < mirrorHit.t;
+        for (int i = 0; i < 1; i++) {
+            // For spherical collision we must test three things.
+            // 1. The ray must collide with the sphere.
+            //    We take the second collision and test if it is on the correct side.
+            // 2. We test if the collision point is within the bounds of the edge planes.
+            MirrorInfo mirror = mirrors[i];
+            vec3 normal = mirror.normal.xyz;
+            vec3 center = mirror.center.xyz;
+            vec3 surfacePoint = vertices[mirror.offset].xyz;
+
+            // First we collide with the sphere.
+            vec3 sphereCenter = center + sphereFocus*normal;
+            float rsq = DistanceSq(sphereCenter, surfacePoint);
+            HitInfo hit = RaySphereHit(rayOrigin, rayDir, sphereCenter, rsq, false);
+
+            // Test if the collision occured at the correct side.
+            // For focus > 0, we have spherical space and we collide with the outside.
+            // For focus < 0, we have hyperbolic space and we collide with the inside.
+            // Hyperbolic space does not require any special tests.
+            float planeProj = dot(hit.point-center, normal);
+
+            // Now we compute the closest edge projection.
+            // We iterate each edge and project on the normal given by the plane spanned by
+            // the outward normal at a point and its edge vector.
+            // We face the normal of the plane toward our side, then a projection > 0 is correct. 
+            float edgeProj = -1;
+            for (int j = 0; j < mirror.vertexCount; j++) {
+                vec3 u = vertices[mirror.offset + j].xyz;
+                vec3 v = vertices[mirror.offset + (j+1)%mirror.vertexCount].xyz;
+
+                // Now the normal is given by the (v-u) x u, since our vertices lie on the sphere.
+                // We flip the normal to face the center of the face.
+                vec3 edgeNormal = normalize(cross(v-u, u));
+                edgeNormal *= sign(dot(edgeNormal, center-u));
+
+                // Now we project with the vector from u to the hit point.
+                float proj = dot(edgeNormal, hit.point-u);
+                if (proj >= 0 && (proj < edgeProj || edgeProj < 0)) {
+                    edgeProj = proj;
+                }
+            }
+            edgeProj = 0.01;
+
+            bool validity = planeProj < 0 && hit.t > STEP_MIN;
+            bool optimality = closestMirrorHit.t < 0 || hit.t < closestMirrorHit.t;
+
+            // HitInfo hit = RayPlaneHit(rayOrigin, rayDir, normal, surfacePoint);
+            // bool validity = hit.t > STEP_MIN && dot(rayDir, normal) < 0;
+            // bool optimality = closestMirrorHit.t < 0 || hit.t < closestMirrorHit.t;
+
             if (optimality && validity) {
-                mirrorHit = hit;
-                mirrorIndex = i;
-                mirrorOffset = offset;
-                mirrorSize = size;
+                closestMirrorHit = hit;
+                closestMirror = mirror;
+                closestEdgeProj = edgeProj;
             }
-            offset += size;
         }
-
-        if (mirrorHit.t < 0)
+
+        // if (closestMirrorHit.t < 0) {
+        //     break;
+        // }
+
+        if (closestEdgeProj > 0 && closestEdgeProj < 0.05) {
+            //col = vec3(0, (1-0.5*closestEdgeProj/edgeThickness)*pow(0.85, float(b)), 0);
+            col = vec3(0, 0, 1);
             break;
-
-        // Test collision with an edge.
-        bool hasHit = false;
-        for (int j = 0; j < mirrorSize; j++) {
-            vec3 u = vertices[mirrorOffset+j].xyz;
-            vec3 v = vertices[mirrorOffset+(j+1)%mirrorSize].xyz;
-            LineProjection lp = PointLineProject(mirrorHit.point, u, v);
-            bool bounds = 0 <= lp.utov && lp.utov <= 1;
-            float ds = DistanceSq(lp.projection, mirrorHit.point);
-            if (ds < edgeThickness && bounds) {
-
-                col = vec3(0, (1-0.5*ds/edgeThickness)*pow(0.85, float(b+1)), 0);
-                hasHit = true;
-            }
         }
-        if (hasHit) break;
-
-        rayDir = reflect(rayDir, mirrorHit.normal);
-        rayOrigin = mirrorHit.point;
+        // bool hasHit = false;
+        // for (int j = 0; j < closestMirror.vertexCount; j++) {
+        //     vec3 u = vertices[closestMirror.offset+j].xyz;
+        //     vec3 v = vertices[closestMirror.offset+(j+1)%closestMirror.vertexCount].xyz;
+
+        //     // Project point on the line and use distance as color.
+        //     LineProjection lp = PointLineProject(closestMirrorHit.point, u, v);
+        //     bool bounds = 0 <= lp.utov && lp.utov <= 1;
+        //     float ds = DistanceSq(lp.projection, closestMirrorHit.point);
+        //     if (ds < edgeThickness && bounds) {
+        //         col = vec3(0, (1-0.5*ds/edgeThickness)*pow(0.85, float(b)), 0);
+        //         hasHit = true;
+        //         break;
+        //     }
+        // }
+        // if (hasHit) break;
+
+        // Dont reflect on backward facing normals, but do allow collision
+        // with the edges of that mirror.
+        rayOrigin = closestMirrorHit.point;
+        rayDir = reflect(rayDir, closestMirrorHit.normal);
     }
 
     finalColor = vec4(col, 1);
 }
 
-HitInfo RaySphereHit(vec3 o, vec3 d, vec3 cc, float rsq, bool convex)
+HitInfo RaySphereHit(vec3 o, vec3 d, vec3 cc, float rsq, bool outside)
 {
     // Make the point relative to the circle.
     vec3 rel = o-cc;
@@ -113,10 +169,8 @@ HitInfo RaySphereHit(vec3 o, vec3 d, vec3 cc, float rsq, bool convex)
     float t1 = (-b - sqrtD) / (2*a);
     float t2 = (-b + sqrtD) / (2*a);
 
-    // Convex means that we take the last collision point, which relative
-    // to the ray is a convex curve. This is default spherical.
-    // Hyperbolic is not convex but arched toward the ray so concave.
-    float t = float(convex)*max(t1, t2) + (1-float(convex))*min(t1, t2);
+    // Outside collides only from outside. Inside only after entering.
+    float t = float(outside)*min(t1, t2) + (1-float(outside))*max(t1, t2);
 
     HitInfo hit;
     hit.point = o + t*d;

Shaders/ray3.vert

@@ -11,5 +11,5 @@ uniform mat4 mvp;
 
 void main()
 {
-    gl_Position = mvp*vec4(vertexPosition, 1.0);
+    gl_Position = mvp*vec4(vertexPosition*1.3, 1.0);
 }