Upgraded dot product and cross product-related functions to doubles, to improve accuracy of various calculations

bsp.h

@@ -238,7 +238,7 @@ typedef struct World
             const plane_t* plane = &planes[face->plane_id];
 
             // Check if the point lies on the face's plane (it's not strictly necessary to do this, but this check makes the whole function a lot faster)
-            double pointPlaneDist = ((double)point.x * plane->normal.x + (double)point.y * plane->normal.y + (double)point.z * plane->normal.z) - plane->dist;
+            double pointPlaneDist = point.dotProduct(plane->normal) - plane->dist;
             if (fabs(pointPlaneDist) > 0.001)
                 continue;
 

common.h

@@ -19,6 +19,7 @@ typedef struct Vec3                 // Vector or Position
 
     Vec3() : x(0), y(0), z(0) { }
     Vec3(float x, float y, float z) : x(x), y(y), z(z) { }
+    Vec3(double x, double y, double z) : x((float)x), y((float)y), z((float)z) { }
 
     Vec3 operator+(const Vec3& other) const
     {
@@ -45,25 +46,29 @@ typedef struct Vec3                 // Vector or Position
         return Vec3(-x, -y, -z);
     }
 
-    float magnitude() const
+    double magnitude() const
     {
-        return sqrtf(x * x + y * y + z * z);
+        return sqrt((double)x * x + (double)y * y + (double)z * z);
     }
 
     Vec3 normalized() const
     {
-        float invMag = 1.0f / magnitude();
-        return Vec3(x * invMag, y * invMag, z * invMag);
+        double invMag = 1.0 / magnitude();
+        return Vec3(invMag * x, invMag * y, invMag * z);
     }
 
-    float dotProduct(const Vec3 &other) const
+    double dotProduct(const Vec3 &other) const
     {
-        return x * other.x + y * other.y + z * other.z;
+        return (double)x * other.x + (double)y * other.y + (double)z * other.z;
     }
 
     Vec3 crossProduct(const Vec3 &other) const
     {
-        return Vec3(y * other.z - z * other.y, z * other.x - x * other.z, x * other.y - y * other.x);
+        return Vec3(
+            (double)y * other.z - (double)z * other.y,
+            (double)z * other.x - (double)x * other.z,
+            (double)x * other.y - (double)y * other.x
+        );
     }
 
     Vec3 floor() const
@@ -76,3 +81,26 @@ typedef struct Vec3                 // Vector or Position
         return Vec3(ceilf(x), ceilf(y), ceilf(z));
     }
 } vec3_t;
+
+static float halton(int32_t index, int32_t base)
+{
+    float f = 1.0f, result = 0.0f;
+
+    for (int32_t currentIndex = index; currentIndex > 0;) {
+
+        f /= (float)base;
+        result = result + f * (float)(currentIndex % base);
+        currentIndex = (uint32_t)(floorf((float)(currentIndex) / (float)(base)));
+    }
+
+    return result;
+}
+
+static void getJitterOffset(float* outX, float* outY, int32_t index, int32_t phaseCount)
+{
+    const float x = halton((index % phaseCount) + 1, 2) - 0.5f;
+    const float y = halton((index % phaseCount) + 1, 3) - 0.5f;
+
+    *outX = x;
+    *outY = y;
+}

lighting.cpp

@@ -160,8 +160,8 @@ std::unordered_map<const edge_t*, EdgeData> analyze_edges(const world_t* world)
 			const plane_t* planeB = &world->planes[faceB->plane_id];
 			vec3_t normalA = faceA->side ? -planeA->normal : planeA->normal;
 			vec3_t normalB = faceB->side ? -planeB->normal : planeB->normal;
-			float dot = normalA.dotProduct(normalB);
-			bool isSmooth = dot >= 0.5f;//&& dot <= 1;
+			double dot = normalA.dotProduct(normalB);
+			bool isSmooth = dot >= 0.5;//&& dot <= 1;
 			iter->second.isSharpEdge = !isSmooth;
 			break;
 		}
@@ -252,8 +252,8 @@ unsigned char compute_faceVertex_light2(const world_t* world, const face_t* face
 		const plane_t* otherPlane = &world->planes[otherFace->plane_id];
 		vec3_t otherNormal = otherFace->side ? -otherPlane->normal : otherPlane->normal;
 
-		float dot = thisNormal.dotProduct(otherNormal);
-		if (dot < 0.5f)
+		double dot = thisNormal.dotProduct(otherNormal);
+		if (dot < 0.5)
 			continue;	// Sharp edge, we don't want light contribution from this face
 
 		light += sample_lightmap(world, otherFace, faceBounds.find(otherFace)->second, point) + (0xFF - otherFace->baselight);
@@ -318,8 +318,8 @@ SurfaceList group_surfaces(const world_t* world, const VertexFaces& vertexFaces)
 
 						const plane_t* otherPlane = &world->planes[otherFace->plane_id];
 						vec3_t otherNormal = otherFace->side ? -otherPlane->normal : otherPlane->normal;
-						float dot = thisNormal.dotProduct(otherNormal);
-						if (dot < 0.5f)
+						double dot = thisNormal.dotProduct(otherNormal);
+						if (dot < 0.5)
 							continue;	// Sharp edge, face belongs to a different surface
 
 						// Add face to this surface and make sure it won't be reconsidered for any other surfaces
@@ -358,8 +358,8 @@ unsigned char compute_faceVertex_light4(const world_t* world, const face_t* refF
 		vec3_t normal = face->side ? -plane->normal : plane->normal;
 
 		// Check if the face is at a shallow angle with the reference face
-		float dot = normal.dotProduct(refNormal);
-		if (dot < 0.5f)
+		double dot = normal.dotProduct(refNormal);
+		if (dot < 0.5)
 			continue;
 
 		unsigned char sample;

main.cpp

@@ -61,8 +61,8 @@ static float computeFaceArea(const world_t* world, const face_t* face)
 		vertex_t* vertex = &world->vertices[vertIndex];
 		Vec3 vec = vertex->toVec();
 
-		float x = tangent.dotProduct(vec);
-		float y = bitangent.dotProduct(vec);
+		double x = tangent.dotProduct(vec);
+		double y = bitangent.dotProduct(vec);
 		bounds.includePoint(Vec3(x, y, 0));
 	}
 
@@ -125,8 +125,8 @@ int process_faces(const world_t* world, const std::vector<ps1bsp_texture_t>& tex
 		outFace.firstFaceVertex = (unsigned short)outFaceVertices.size();
 		outFace.textureId = (unsigned char)texinfo->texture_id;
 
-		float minS = FLT_MAX, minT = FLT_MAX;
-		float maxS = FLT_MIN, maxT = FLT_MIN;
+		double minS = DBL_MAX, minT = DBL_MAX;
+		double maxS = DBL_MIN, maxT = DBL_MIN;
 
 		// Traverse the list of face edges to collect all of the face's vertices
 		Vec3 vertexSum;
@@ -143,8 +143,8 @@ int process_faces(const world_t* world, const std::vector<ps1bsp_texture_t>& tex
 			Vec3 vertexPoint = vertex->toVec();
 
 			// Calculate texture UV bounds
-			float s = (vertexPoint.dotProduct(texinfo->vectorS) + texinfo->distS) / miptex->width;
-			float t = (vertexPoint.dotProduct(texinfo->vectorT) + texinfo->distT) / miptex->height;
+			double s = (vertexPoint.dotProduct(texinfo->vectorS) + texinfo->distS) / miptex->width;
+			double t = (vertexPoint.dotProduct(texinfo->vectorT) + texinfo->distT) / miptex->height;
 			if (s > maxS) maxS = s; if (s < minS) minS = s;
 			if (t > maxT) maxT = t; if (t < minT) minT = t;
 
@@ -159,8 +159,8 @@ int process_faces(const world_t* world, const std::vector<ps1bsp_texture_t>& tex
 		}
 
 		// If the texture doesn't tile, we don't need to correct the UVs as much
-		float sRange = maxS - minS;
-		float tRange = maxT - minT;
+		double sRange = maxS - minS;
+		double tRange = maxT - minT;
 		if (sRange < 1) sRange = 1;
 		if (tRange < 1) tRange = 1;
 
@@ -181,8 +181,8 @@ int process_faces(const world_t* world, const std::vector<ps1bsp_texture_t>& tex
 			faceVertex.light = 0;
 
 			// Calculate texture UVs
-			float s = (vertexPoint.dotProduct(texinfo->vectorS) + texinfo->distS) / miptex->width;
-			float t = (vertexPoint.dotProduct(texinfo->vectorT) + texinfo->distT) / miptex->height;
+			double s = (vertexPoint.dotProduct(texinfo->vectorS) + texinfo->distS) / miptex->width;
+			double t = (vertexPoint.dotProduct(texinfo->vectorT) + texinfo->distT) / miptex->height;
 			if (minS < 0 || maxS > 1) s = (s - minS) / sRange;
 			if (minT < 0 || maxT > 1) t = (t - minT) / tRange;
 
@@ -218,7 +218,7 @@ int process_faces(const world_t* world, const std::vector<ps1bsp_texture_t>& tex
 			ps1bsp_facevertex_t& faceVertex = outFaceVertices[outFace.firstFaceVertex + faceVertIdx];
 			const vertex_t* vertex = &world->vertices[faceVertex.index];
 			faceVertex.light = compute_faceVertex_light4(world, face, faceBounds, vertex->toVec());
-			faceVertex.light <<= 1;
+			faceVertex.light = (short)((float)faceVertex.light * 1.5f);
 			if (faceVertex.light > 255)
 				faceVertex.light = 255;
 		}