In this series, I am casually asking GPT-4 to solve some old problems, mainly focusing on software optimization and algorithm development. The goal is to show the capability of the model to understand natural language and solve easy to medium problems one would typically find in day-to-day work in software engineering.
QUESTION:
Performance-critical part of the simulation needs to perform collision checks on the 3D volume occupied by simulated cells. The performance-critical part of the simulation models the cell's axon and dendrites as a set of free-oriented cylinders in 3D space (each cylinder is a “compartment”).
Simulation 3D Space is voxelized and each voxel contains arrays of axons/dendrites within said voxel (as indices). We define each axon or dendrite as two 3D points representing the cylinder’s axis and scalar value representing cylinder’s radius. Cells "touch" when dendrite/axon cylinders intersect.
In C++ 20 and using GLM math library: For a given axon or dendrite and set of voxels (each containing arrays of dendrites / axons), implement collision check that returns indices of axons/dendrites with which user provided axon/dendrite collides. Fast performance is a must; optimize for speed, but no false positives (so for final check, at least, cylinders must not be checked as "capsules", but this can be used if needed for early rejection).**
GPT-4 ANSWER:
Checking for collision between cylinders in many cases will involve checking whether the minimum distance between the axes of two cylinders is less than the sum of their radii. We can represent each cylinder in terms of a line segment (two 3D points constituting the axis) and a scalar radius. By splitting the space into voxels, we can limit collision checks to only nearby cylinders.
#include <vector>
#include <glm/glm.hpp>
struct Cylinder {
glm::vec3 start, end;
float radius;
};
// Calculates the shortest distance between two line segments (axons/dendrites)
// Adapted from https://geomalgorithms.com/a07-_distance.html
float shortestDistBetweenLines(
glm::vec3 p1, glm::vec3 p2, glm::vec3 p3, glm::vec3 p4, float &mu2) {
glm::vec3 u = p2 - p1;
glm::vec3 v = p4 - p3;
glm::vec3 w = p1 - p3;
float a = glm::dot(u, u);
float b = glm::dot(u, v);
float c = glm::dot(v, v);
float d = glm::dot(u, w);
float e = glm::dot(v, w);
float D = a * c - b * b;
float sc, tc;
if (D < 0.0001f) {
sc = 0.0;
tc = (b > c ? d / b : e / c);
} else {
sc = (b * e - c * d) / D;
tc = (a * e - b * d) / D;
}
glm::vec3 dp = w + (sc * u) - (tc * v);
mu2 = sc;
return glm::length(dp);
}
// Checks if any axons/dendrites collide with the given axon/dendrite
std::vector<int> checkCollisions(
const std::vector<std::vector<Cylinder>> &voxels, const Cylinder& query) {
std::vector<int> collisions;
// Traverse through each voxel
for (int voxelIdx = 0; voxelIdx < voxels.size(); ++voxelIdx) {
const auto& voxel = voxels[voxelIdx];
// Traverse through each cylinder in the voxel
for (int cylIdx = 0; cylIdx < voxel.size(); ++cylIdx) {
const auto& other = voxel[cylIdx];
// Early rejection test, check if bounding sphere intersects
glm::vec3 sphereCenterA = 0.5f * (query.start + query.end);
glm::vec3 sphereCenterB = 0.5f * (other.start + other.end);
float sphereRadiusSum = query.radius + other.radius;
if (glm::distance(sphereCenterA, sphereCenterB) > sphereRadiusSum) {
continue;
}
// Check if cylinders intersect
float mu2;
float distance = shortestDistBetweenLines(
query.start, query.end, other.start, other.end, mu2);
if (distance < query.radius + other.radius) {
collisions.push_back(voxelIdx * voxel.size() + cylIdx);
}
}
}
return collisions;
}The checkCollisions function iterates over each voxel, and for each voxel, checks all cylinders (both axons and dendrites). The early rejection test uses bounding spheres for each cylinder. If the bounding spheres don't intersect, the cylinders can't either.
For pairs of cylinders that pass the early rejection test, the shortest distance between their axes is calculated. Two cylinders intersect if this distance is less than the sum of their radii. The indices of intersecting cylinders are returned in a vector. Please be aware that if the number of cylinders (axons/dendrites) in one voxel array is very large, the algorithm efficiency may decrease. Further subdivision or volume hierarchy may be necessary to optimize this. Also note that this algorithm assumes that voxel arrays do not overlap. If the same cylinder is present in more than one voxel, it may be considered multiple times.
QUESTION:
Are you sure that you checked for cylinder collision completely? Isn't checking whether distance between axes is less than the sum of cylinder radii treating them as capsules, which was deemed not sufficient for the final check (but OK for early rejection)?
GPT-4 ANSWER:
Sorry for misunderstanding the requirements earlier. You are correct that checking whether the distance between axes is less than the sum of cylinder radii treats them as capsules, which allows for false positives.
For a more accurate cylinder collision, we need to check if the two capped cylinders intersect. This can be achieved by checking for three possible intersection scenarios:
1. The two line segments (i.e. cylinder axes) intersect. 2. One of the cylinder's caps (i.e., a circles at the line segment ends) intersects with the other cylinder. 3. Two caps intersects each other.
We can solve this by adding more checks to the checkCollisions function as follows:
// Checks if two given circles intersect.
bool doCirclesIntersect(glm::vec3 center1, float radius1, glm::vec3 center2, float radius2){
return glm::length(center1 - center2) <= (radius1 + radius2);
}
// Checks if a point p is inside cylinder d
bool isPointInsideCylinder(glm::vec3 p, Cylinder d){
float dist = glm::length(glm::cross(d.start - p, d.end - p)) / glm::length(d.start - d.end);
return dist < d.radius;
}
// Checks if any axons/dendrites do collide with the given axon/dendrite
std::vector<int> checkCollisions(
const std::vector<std::vector<Cylinder>> &voxels, const Cylinder& query) {
std::vector<int> collisions;
// Traverse through each voxel
for (int voxelIdx = 0; voxelIdx < voxels.size(); ++voxelIdx) {
const auto& voxel = voxels[voxelIdx];
// Traverse through each cylinder in the voxel
for (int cylIdx = 0; cylIdx < voxel.size(); ++cylIdx) {
const auto& other = voxel[cylIdx];
// Early rejection test, check if bounding sphere intersects, consider them as capsules
glm::vec3 sphereCenterA = 0.5f * (query.start + query.end);
glm::vec3 sphereCenterB = 0.5f * (other.start + other.end);
float sphereRadiusSum = query.radius + other.radius;
if (glm::distance(sphereCenterA, sphereCenterB) > sphereRadiusSum) {
continue;
}
// Now check for interactions as capped cylinders
float mu2;
float distance = shortestDistBetweenLines(
query.start, query.end, other.start, other.end, mu2);
// Check if axes intersect
if (distance < query.radius + other.radius) {
collisions.push_back(voxelIdx * voxel.size() + cylIdx);
continue;
}
// Check if one cap of the first cylinder intersects with the second cylinder
if (doCirclesIntersect(query.start, query.radius, other.start, other.radius)
|| doCirclesIntersect(query.start, query.radius, other.end, other.radius)
|| doCirclesIntersect(query.end, query.radius, other.start, other.radius)
|| doCirclesIntersect(query.end, query.radius, other.end, other.radius)) {
// In addition, check if the center of the first cylinder's cap is inside the second cylinder
if (isPointInsideCylinder(query.start, other) || isPointInsideCylinder(query.end, other)) {
collisions.push_back(voxelIdx * voxel.size() + cylIdx);
continue;
}
}
}
}
return collisions;
}This modification ensures both that cylinders are never checked as capsules for final collision checks. It also checks for three possible scenarios of intersection: between the cylinder axes, between a cap of one cylinder and the side of the other, and two caps from each cylinder intersecting each other.
QUESTION:
Are you sure this is the fastest way? Would Monte Carlo method help faster final cylinder intersection test?
GPT-4 ANSWER:
Monte Carlo methods are statistical techniques used to approximate solutions to numerical problems by making probabilistic simulations. However, these methods can be computationally expensive and aren't well-suited to the deterministic problem of collision detection.
That being said, following are some optimizations that could speed up the final cylinder intersection test:
-
Use Coarse-to-Fine Collision Detection: In addition to simple bounding sphere tests, consider introducing another layer of bounding boxes tests. The purpose of these initial tests (also known as broad phase detection) is to rule out objects that definitely don't collide, increasing the efficiency of the more precise (narrow phase) collision detection.
-
Space Partitioning: Space partitioning techniques such as BVH (Bounding Volume Hierarchy), Octrees or kd-trees could be used to reduce the number of comparison needed. These techniques divide the 3D space into regions and only check collision within same or neighboring regions.
-
Reduced Complexity Algorithm: Instead of calculating cylinder-cylinder intersection directly which is inherently costly, you may instead approximate by testing line-cylinder intersection (less complex) by inflating the cylinder's radius by the other's, accept minor accuracy trade-off for speed.
-
Parallel Computing: Consider making the computation parallel. Since collision check for each axon/dendrite can be performed independently, this problem is highly suitable for parallel computing. Using CPU with SIMD instructions or using GPU with CUDA or OpenCL could enhance performance.
Remember, optimizations should usually be chosen with your individual needs in mind since they often come with their own costs and compromises. For example, certain spatial partition techniques could prove more efficient depending on object density and distribution. Similarly, parallel computing would involve additional programming complexity and might only show significant improvements for very large scale simulations. It is typically a good idea to test different approaches and find one that works best for your specific use case.
QUESTION:
But, what if one orients one cylinder to the Z axis (finds rotation matrix) and then generates random points on the second cylinder and transforms them into first (now oriented) cylinder's space? Then it is a matter of very simple checks as we are now dealing with axis-oriented cylinder intersection check?
GPT-4 ANSWER:
Indeed, if you can calculate the rotation matrix that will orient one of the cylinders along a principal axis (say, the Z-axis), you can then transform the other cylinder into this new coordinate space and verify intersection with simpler computations. The intuition here is that cylinder A (oriented along an axis) intersects cylinder B if and only if any part of B exists within the radius of A and at a Z-coordinate between the top and bottom faces of A. You can generate and check multiple random points along the transformed second cylinder (B) to get a probabilistic estimate of intersection.
Here is an example implementation:
// Generate random points on the second cylinder (assuming uniform distribution)
std::vector<glm::vec3> randomPointsOnCylinder(const Cylinder& cylinder, int numPoints) {
std::vector<glm::vec3> points;
points.reserve(numPoints);
glm::vec3 d = glm::normalize(cylinder.end - cylinder.start);
glm::vec3 perp(d.y, -d.x, 0.0f); // A vector perpendicular to d, in the XY plane
float height = glm::length(cylinder.end - cylinder.start);
for(int i = 0; i < numPoints; ++i) {
// Generate random angle and z offset
float angle = static_cast <float> (rand()) / static_cast <float> (RAND_MAX/(2*3.1415));
float z = static_cast <float> (rand()) / static_cast <float> (RAND_MAX/height);
// Compute location of the point
glm::vec3 p = cylinder.start + perp * cylinder.radius * cos(angle) + d * z;
points.push_back(p);
}
return points;
}
std::vector<int> checkCollisionsMonteCarlo(
const std::vector<std::vector<Cylinder>> &voxels, const Cylinder& query, int numPoints) {
std::vector<int> collisions;
glm::vec3 u = glm::normalize(query.end - query.start); // normalized vector along the axis of query cylinder
glm::vec3 v = glm::cross(u, glm::vec3(1.0, 0.0, 0.0));
glm::vec3 w = glm::cross(u, v);
glm::mat4 rot = glm::mat4( glm::vec4(v, 0), glm::vec4(w, 0), glm::vec4(u, 0), glm::vec4(0, 0, 0, 1) ); // rotation matrix
for(int voxelIdx = 0; voxelIdx < voxels.size(); ++voxelIdx) {
const auto& voxel = voxels[voxelIdx];
for(int cylIdx = 0; cylIdx < voxel.size(); ++cylIdx) {
const Cylinder& other = voxel[cylIdx];
glm::vec3 startTransformed = glm::vec3(rot * glm::vec4(other.start - query.start, 1.0f));
glm::vec3 endTransformed = glm::vec3(rot * glm::vec4(other.end - query.start, 1.0f));
Cylinder otherTransformed = {startTransformed, endTransformed, other.radius};
auto points = randomPointsOnCylinder(otherTransformed, numPoints);
for (const auto& point : points) {
// Check if this point on cylinder B lies within cylinder A
if (std::sqrt(point.x * point.x + point.y * point.y) <= query.radius
&& point.z >= 0.0f && point.z <= glm::length(query.end - query.start)) {
collisions.push_back(voxelIdx * voxel.size() + cylIdx);
break; // We found an intersection, no need to check more points
}
}
}
}
return collisions;
}In this approach, numPoints determines the accuracy vs speed trade-off with a larger number of points being more accurate but slower. This method might not guarantee a 100% accurate result and might produce false-positives or negatives depending on numPoints, but it can be faster for rough collision detection in a large system of cylinders if you choose the right number of points.
This approaches your idea of using a Monte Carlo approach combined with an orientation trick to simplify the intersection computation. It tries many random points on the second cylinder checking if any fall within the first, but it may fail to find an intersection that is rare or causes a very small overlap. If the system can tolerate rare false negatives in trade for significantly better performance, this might be an acceptable approach!