|
|
Post by likeaglove on Apr 15, 2020 1:52:36 GMT
I have implemented the Cone Shape from Chapter 13 All tests and examples pass except one. In the test scenario: Intersecting a cone with a ray. With inputs point(0, 0, -5) | vector(1, 1, 1) the expected output is t0 = 8.66025 & t1 = 8.66025;
However, the LocalIntersect function is returning an empty Array, when debugging I see that the "disc" value is a rather large negative value.
Here is the set up code:
Cone con = new Cone();
Tuple Origin = new Tuple(0f,0f,-5f,1);
Tuple Direction = new Tuple(1f,1f,1f,0).Normalise();
Ray r = new Ray(Origin,Direction);
Intersection[] xs = con.LocalIntersect(r);
Console.WriteLine(xs.Length);
Console.WriteLine(xs[0].T_Value);
Console.WriteLine(xs[1].T_Value);
This results in a out of bounds exception, due to the Intersection Array being empty, where it should be filled with duplicated 8.66025 values.
I was wondering if anybody could point out an error in my Cone implementation. All other tests in the chapter pass, including the ones beyond this example.
using System;
using System.Collections.Generic;
namespace The_Ray_Tracer_Challenge
{
public class Cone : Shape
{
public float minimum { get; set; } = float.NegativeInfinity;
public float maximum { get; set; } = float.PositiveInfinity;
public bool closed { get; set; } = false;
public override Intersection[] LocalIntersect(Ray localRay)
{
List<Intersection> xsList = new List<Intersection>();
Tuple o = localRay.Origin;
Tuple d = localRay.Direction;
float a = localRay.Direction.x * localRay.Direction.x - localRay.Direction.y * localRay.Direction.y + localRay.Direction.z * localRay.Direction.z;
float b = 2.0f * localRay.Origin.x * localRay.Direction.x - 2.0f * localRay.Origin.y * localRay.Direction.y + 2.0f * localRay.Origin.z * localRay.Direction.z;
float c = localRay.Origin.x * localRay.Origin.x - localRay.Origin.y * localRay.Origin.y + localRay.Origin.z * localRay.Origin.z;
if (a <= 0 && MathF.Abs(b) > 0)
{
float t = -c / (2 * b);
xsList.Add(new Intersection(t, this));
}
if (a != 0)
{
float disc = (b * b) - 4.0f * a * c;
if (disc >= 0)
{
float t0 = (-b - MathF.Sqrt(disc)) / (2 * a);
float t1 = (-b + MathF.Sqrt(disc)) / (2 * a);
if (t0 > t1)
{
float temp = t0;
t0 = t1;
t1 = temp;
}
float y0 = localRay.Origin.y + t0 * localRay.Direction.y;
if (this.minimum < y0 && y0 < this.maximum)
{
xsList.Add(new Intersection(t0, this));
}
float y1 = localRay.Origin.y + t1 * localRay.Direction.y;
if (this.minimum < y1 && y1 < this.maximum)
{
xsList.Add(new Intersection(t1, this));
}
}
}
IntersectCaps(this, localRay, xsList);
return xsList.ToArray();
}
public override Tuple LocalNormalAt(Tuple localPoint)
{
float dist = localPoint.x * localPoint.x + localPoint.z * localPoint.z;
if (dist < 1 && localPoint.y >= this.maximum - Arithmetic.EPSILON)
{
return new Tuple(0, 1, 0, 0);
}
if (dist < 1 && localPoint.y <= this.minimum + Arithmetic.EPSILON)
{
return new Tuple(0, -1, 0, 0);
}
float y = MathF.Sqrt(dist);
if (localPoint.y > 0)
{
y *= -1;
}
return new Tuple(localPoint.x, y, localPoint.z);
}
private static bool CheckCap(Ray r, float t, float y)
{
float x = r.Origin.x + t * r.Direction.x;
float z = r.Origin.z + t * r.Direction.z;
return (x * x) + (z * z) <= y * y;
}
private static void IntersectCaps(Cone cyl, Ray r, List<Intersection> xsList)
{
if (!cyl.closed || (r.Direction.y) <= Arithmetic.EPSILON)
{
return;
}
float t = (cyl.minimum - r.Origin.y) / r.Direction.y;
if (CheckCap(r, t, cyl.minimum))
{
xsList.Add(new Intersection(t, cyl));
}
t = (cyl.maximum - r.Origin.y) / r.Direction.y;
if (CheckCap(r, t, cyl.minimum))
{
xsList.Add(new Intersection(t, cyl));
}
return;
}
}
}
|
|
|
|
Post by likeaglove on Apr 15, 2020 2:10:16 GMT
I now have all test passing by changing the line float disc = (b * b) - 4.0f * a * c;to float disc = MathF.Abs((b * b) - 4.0f * a * c);It seems to be working fine, but im unsure this could cause some errors down the road. Edit: Yup, cause some errors in rendering, however, when I take away this MathF.Abs() call everything seems to be rendering just fine. However, the test has gone back to failing. ![]() |
|
fremag
Junior Member
Posts: 73
|
Post by fremag on Apr 15, 2020 8:30:50 GMT
Hi, In LocalIntersect, I see these lines: if (a <= 0 && MathF.Abs(b) > 0)
if (a != 0)You should check:
if( MathF.Abs(a) <= Arithmetic.EPSILON ...
and:
if( MathF.Abs(a) > Arithmetic.EPSILON)
The rest looks good to me. |
|
|
|
Post by likeaglove on Apr 17, 2020 4:21:18 GMT
Hello,
I have also tried those solutions, but the result is the same.
When debugging the value for my solution the values I receive are:
a: 0.3333333
b: -5.7735023
c: 25
This then follows on to a disc value of: -3.8146973E-06
causing it to only run the IntersectCaps function, and as the cone is capless, this returns immediately.
|
|
fremag
Junior Member
Posts: 73
|
Post by fremag on Apr 19, 2020 14:51:27 GMT
Here is a screenshot of my debugger running the test, you have all the values displayed so you can compare with yours.
The C# code is here, if you want to have a look, it's very close to yours.
|
|
fremag
Junior Member
Posts: 73
|
Post by fremag on Apr 19, 2020 14:59:53 GMT
There is a bug IntersectCaps:
t = (cyl.maximum - r.Origin.y) / r.Direction.y;
if (CheckCap(r, t, cyl.minimum))
The second call to CheckCap should be done with cyl.maximum, not cyl.minimum
|
|
|
|
Post by jaredp on Sept 5, 2020 18:58:08 GMT
I too had the exact same issue as likeaglove. My values for a,b,c and the discriminant were as follows:
a: 0.333333 b: -5.7735 c: 25
discr: -3.08061e-15
It's merely a floating point precision issue. I can actually get the test to pass (i.e, discr == 0) by re-compiling without the "-march=native" GCC flag. I'm on Linux with an AMD CPU and my raytracer is in C++.
Hopefully that saves some time for anyone who might run into the same problem.
|
|
|
|
Post by ulfben on Mar 28, 2023 5:38:06 GMT
I now have all test passing by changing the line float disc = (b * b) - 4.0f * a * c;to float disc = MathF.Abs((b * b) - 4.0f * a * c);It seems to be working fine, but im unsure this could cause some errors down the road. Edit: Yup, cause some errors in rendering, however, when I take away this MathF.Abs() call everything seems to be rendering just fine. However, the test has gone back to failing. ![]() If changing the calculation of the discriminant to always return a positive value makes the unit test pass, it indicates that there were some precision errors in the original calculation of the discriminant. In general, the discriminant of the quadratic equation should be non-negative if there are real solutions, and negative if there are complex solutions. However, due to floating-point roundoff errors, the discriminant may sometimes become slightly negative even when there are real solutions, leading to incorrect results. I had the exact same issues, and solved it like so; constexpr auto local_intersect([[maybe_unused]] const Cone& cone, const Ray& local_ray) noexcept {
using math::square, math::is_zero, math::sqrt, math::is_between, math::max, math::abs;
std::vector<Real> result;
const auto a = square(local_ray.dx()) - square(local_ray.dy()) + square(local_ray.dz());
const auto b = 2 * ((local_ray.x() * local_ray.dx()) - (local_ray.y() * local_ray.dy()) + (local_ray.z() * local_ray.dz()));
const auto c = square(local_ray.x()) - square(local_ray.y()) + square(local_ray.z());
//When a is zero, it means the ray is parallel to one of the cone’s halves
//this still means the ray might intersect the other half of the cone.
if (is_zero(a)){ //In this case the ray will miss when both a and b are zero.
if (!is_zero(b)) { //If a is zero but b isn’t, calc the single point of intersection:
const auto t = -c / (2.0f * b);
result.push_back(t);
}
}
else {// a is non-zero
//In general, the discriminant of the quadratic equation should be non-negative if there are real solutions,
// and negative if there are complex solutions. However, due to floating-point roundoff errors, the discriminant
// may sometimes become slightly negative even when there are real solutions, leading to incorrect results.
//Here, rel_eps is a small relative epsilon based on the magnitudes of the coefficients, making it more robust to different combinations of coefficients
const auto rel_eps = math::EPSILON * max(abs(a), abs(b), abs(c));
const auto discriminant = square(b) - 4.0f * a * c;
if (discriminant >= -rel_eps) { //The discriminant is considered to be non-negative if it is greater than or equal to -rel_eps
const auto sqrt_discriminant = sqrt(max(discriminant, 0.0f)); //round to 0 if need be.
const auto t1 = (-b - sqrt_discriminant) / (2 * a);
const auto t2 = (-b + sqrt_discriminant) / (2 * a);
const auto y1 = local_ray.y() + t1 * local_ray.dy();
if (is_between(y1, cone.minimum, cone.maximum)) {
result.push_back(t1);
}
const auto y2 = local_ray.y() + t2 * local_ray.dy();
if (is_between(y2, cone.minimum, cone.maximum)) {
result.push_back(t2);
}
}
}
intersect_caps(cone, local_ray, result);
return result;
};In my case an epsilon of 0.000001f / 1.0e-6f worked fine and passes all the tests from the book. To help future readers find this thread: |
|
|
|
Post by loop23 on Jul 18, 2023 1:10:48 GMT
Hi, I'm having the same issue where the test "Intersecting a cone with a ray" will not pass. I implemented the code for getting comparing the discriminant against a "-rel_eps" value and calculating a square root discriminant to try fix the rounding errors which got me a little closer. Currently my program outputs these values for for each of the 3 examples in that test.
A: 1
B: -10
C: 25
Discriminant: 0
-rel_eps: -2.5e-05
Square Root Discriminant: 0
XS Count: 2 Expected: 2
xs[0] time: 5 Expected t0: 5
xs[1] time: 5 Expected t1: 5
A: 0.333333
B: -5.7735
C: 25
Discriminant: -3.88971e-06
-rel_eps: -2.5e-05
Square Root Discriminant: 0
XS Count: 2 Expected: 2
xs[0] time: 8.66025 Expected t0: 8.66025
xs[1] time: 8.66025 Expected t1: 8.66025
A: 0.111111
B: -6
C: 25
Discriminant: 24.8889
-rel_eps: -2.5e-05
Square Root Discriminant: 4.98888
XS Count: 2 Expected: 2
xs[0] time: 16.9489 Expected t0: 4.55006
xs[1] time: 37.0511 Expected t1: 49.4499
My code implementation is below.
#include "cone.h"
#include <limits>
#include <algorithm>
#include <cmath>
Cone::Cone(float minimum, float maximum, bool closed)
: Shape("cone", ObjectType::kCone),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
Cone::Cone(const std::string& name, float minimum, float maximum, bool closed)
: Shape(name, ObjectType::kCone),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
Cone::Cone(const Material& material, float minimum, float maximum, bool closed)
: Shape("cone", ObjectType::kCone, material),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
Cone::Cone(const Matrix4& transform, float minimum, float maximum, bool closed)
: Shape("cone", ObjectType::kCone, transform),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
Cone::Cone(const std::string& name, const Material& material, float minimum, float maximum, bool closed)
: Shape(name, ObjectType::kCone, material),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
Cone::Cone(const std::string& name, const Matrix4& transform, float minimum, float maximum, bool closed)
: Shape(name, ObjectType::kCone, transform),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
Cone::Cone(const Material& material, const Matrix4& transform, float minimum, float maximum, bool closed)
: Shape("cone", ObjectType::kCone, material, transform),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
Cone::Cone(const std::string& name, const Material& material, const Matrix4& transform, float minimum, float maximum, bool closed)
: Shape(name, ObjectType::kCone, material, transform),
minimum_(minimum), maximum_(maximum), closed_(closed) {}
bool Cone::operator==(const Object& object) {
Cone* other = (Cone*)&object;
return(this->GetName() == other->GetName() &&
this->GetTransform() == other->GetTransform() &&
this->GetObjectType() == other->GetObjectType() &&
this->GetMaterial() == other->GetMaterial());
}
Cone& Cone::operator=(const Object& object) {
Cone* other = (Cone*)&object;
this->SetName(other->GetName());
this->SetObjectType(other->GetObjectType());
this->SetTransform(other->GetTransform());
this->SetMaterial(other->GetMaterial());
return *this;
}
std::vector<Intersection> Cone::local_intersect(const utils::RayStruct& local_ray) {
std::vector<Intersection> xs;
const auto origin_x = local_ray.origin[0];
const auto direction_x = local_ray.direction[0];
const auto origin_y = local_ray.origin[1];
const auto direction_y = local_ray.direction[1];
const auto origin_z = local_ray.origin[2];
const auto direction_z = local_ray.direction[2];
const auto a = pow(direction_x, 2) - pow(direction_y, 2) + pow(direction_z, 2);
const auto b = 2 * origin_x * direction_x - 2 * origin_y * direction_y + 2 * origin_z * direction_z;
const auto c = pow(origin_x, 2) - pow(origin_y, 2) + pow(origin_z, 2);
const auto discriminant = (b * b) - 4.0f * a * c;
std::cout << "A: " << a << std::endl;
std::cout << "B: " << b << std::endl;
std::cout << "C: " << c << std::endl;
std::cout << "Discriminant: " << discriminant << std::endl;
// When a is zero, it means the ray is parallel to one of the cone’s halves
// this still means the ray might intersect the other half of the cone.
if (utils::equal(a, 0)) {
if (!utils::equal(b, 0)) {
const auto t = -c / (2 * b);
xs.push_back(Intersection(t, this));
}
}
// In general, the discriminant of the quadratic equation should be non-negative if there are real solutions,
// and negative if there are complex solutions. However, due to floating-point roundoff errors, the discriminant
// may sometimes become slightly negative even when there are real solutions, leading to incorrect results.
// Here, rel_eps is a small relative epsilon based on the magnitudes of the coefficients, making it more robust to0
// different combinations of coefficients
const float eps = 0.000001f;
const auto rel_eps = eps * std::max<double>({ abs(a), abs(b), abs(c) });
std::cout << "-rel_eps: " << -rel_eps << std::endl;
// ray does not intersect the cone
if (discriminant < -rel_eps) { //The discriminant is considered to be negative if it is less than -rel_ep
return {};
}
const auto sqrt_discriminant = sqrt(std::max(discriminant, 0.0)); //round to 0 if need be.
std::cout << "Square Root Discriminant: " << sqrt_discriminant << std::endl;
auto t0 = (-b - sqrt(sqrt_discriminant)) / (2 * a);
auto t1 = (-b + sqrt(sqrt_discriminant)) / (2 * a);
//if (t0 > t1) {
// auto temp = t0;
// t0 = t1;
// t1 = temp;
//}
//xs = {};
const auto y0 = origin_y + t0 * direction_y;
const auto y1 = origin_y + t1 * direction_y;
if (this->minimum_ < y0 && y0 < this->maximum_) {
xs.push_back(Intersection(t0, this));
}
if (this->minimum_ < y1 && y1 < this->maximum_) {
xs.push_back(Intersection(t1, this));
}
this->intersect_caps(local_ray, xs);
return xs;
}
Vector Cone::local_normal_at(const Point& local_point) const {
// compute the square of the distance from the y axis
const auto dist = pow(local_point[0], 2) + pow(local_point[2], 2);
if (dist < 1 && local_point[1] >= this->maximum_ - utils::kEPSILON) {
return Vector(0, 1, 0);
}
else if (dist < 1 && local_point[1] <= this->minimum_ + utils::kEPSILON) {
return Vector(0, -1, 0);
}
return Vector(local_point[0], 0, local_point[2]);
}
float Cone::GetMinimum() const { return minimum_; }
float Cone::GetMaximum() const { return maximum_; }
bool Cone::GetClosed() const { return closed_; }
void Cone::SetMinimum(float minimum) {
minimum_ = minimum;
}
void Cone::SetMaximum(float maximum) {
maximum_ = maximum;
}
void Cone::SetClosed(bool closed) {
closed_ = closed;
}
void Cone::intersect_caps(const utils::RayStruct& local_ray, std::vector<Intersection>& xs) {
// caps only matter if the cone is closed, and might possibly be
// intersected by the ray.
if (this->closed_ == false || utils::is_close_to_zero(local_ray.direction[1])) return;
// check for an intersection with the lower end cap by intersecting
// the ray with the plane at y = cyl.minimum
const auto lower_t = (this->minimum_ - local_ray.origin[1]) / local_ray.direction[1];
if (utils::check_cap(local_ray, lower_t)) {
xs.push_back(Intersection(lower_t, this));
}
// check for an intersection with the upper end cap by intersecting
// the ray with the plane at y = cyl.maximum
const auto upper_t = (this->maximum_ - local_ray.origin[1]) / local_ray.direction[1];
if (utils::check_cap(local_ray, upper_t)) {
xs.push_back(Intersection(upper_t, this));
}
}
Any idea why the values for the last example are off?
Thanks
|
|
|
|
Post by icolomby on Jul 31, 2023 19:19:15 GMT
Hi,
You are taking the square root of the discriminate twice when calculating t0 and t1.
auto t0 = (-b - sqrt(sqrt_discriminant)) / (2 * a);
auto t1 = (-b + sqrt(sqrt_discriminant)) / (2 * a);
sqrt_discriminant is already the square root of the discriminant no need to call sqrt again.
Ian.
|
|
|
|
Post by loop23 on Aug 3, 2023 1:02:22 GMT
Thanks, that worked and all tests are passing now.
|
|
