From 415d1c57d1eaad75f0636980862b33efe47713cb Mon Sep 17 00:00:00 2001 From: Don Cross Date: Tue, 14 Dec 2021 20:35:42 -0500 Subject: [PATCH] Fixed raytracer problems due to numeric scale. I had lots of problems with using AU as my scale units. By changing to 10000*km units, the vector equation solver works correctly again, and I actually get an image of Jupiter and its moons. However, it does not match the test photo: Ganymede does not appear close to the planet, nor is there a shadow of it on the planet. I will have to debug that separately. --- demo/c/raytrace/debug.cpp | 88 ++++++++++++++++++++++++++++++++++++ demo/c/raytrace/main.cpp | 35 ++++++++------ demo/c/raytrace/run | 2 + demo/c/raytrace/spheroid.cpp | 12 ++--- 4 files changed, 117 insertions(+), 20 deletions(-) create mode 100644 demo/c/raytrace/debug.cpp create mode 100755 demo/c/raytrace/run diff --git a/demo/c/raytrace/debug.cpp b/demo/c/raytrace/debug.cpp new file mode 100644 index 00000000..c2c9dc3d --- /dev/null +++ b/demo/c/raytrace/debug.cpp @@ -0,0 +1,88 @@ +/* + debug.cpp + + Copyright (C) 2013 by Don Cross - http://cosinekitty.com/raytrace + + This software is provided 'as-is', without any express or implied + warranty. In no event will the author be held liable for any damages + arising from the use of this software. + + Permission is granted to anyone to use this software for any purpose, + including commercial applications, and to alter it and redistribute it + freely, subject to the following restrictions: + + 1. The origin of this software must not be misrepresented; you must not + claim that you wrote the original software. If you use this software + in a product, an acknowledgment in the product documentation would be + appreciated but is not required. + + 2. Altered source versions must be plainly marked as such, and must not be + misrepresented as being the original software. + + 3. This notice may not be removed or altered from any source + distribution. + + ------------------------------------------------------------------------- + + Implements output operators and other stuff + useful for debugging the ray tracer. +*/ + +#include +#include +#include +#include "imager.h" + +namespace Imager +{ + void Indent(std::ostream& output, int depth) + { + const int numSpaces = 4 * depth; + for (int i=0; i < numSpaces; ++i) + { + output << ' '; + } + } + + std::ostream& operator<< (std::ostream &output, const Color& color) + { + using namespace std; + + output.precision(3); + output << scientific; + output << "(R " << color.red << ", G " << color.green << ", B " << color.blue << ")"; + return output; + } + + std::ostream& operator<< (std::ostream& output, const Intersection& intersection) + { + using namespace std; + + output << "{ "; + if (intersection.solid) + { + output << intersection.solid->GetTag(); + } + if (intersection.tag) + { + output << ":" << intersection.tag; + } + output.precision(3); + output << fixed; + output << " d=" << intersection.distanceSquared; + output << ", p=" << intersection.point; + output << ", n=" << intersection.surfaceNormal; + output << " }"; + return output; + } + + std::ostream& operator<< (std::ostream& output, const Vector& vector) + { + using namespace std; + + output.precision(3); + output << fixed; + output << "(" << vector.x << ", " << vector.y << ", " << vector.z << ")"; + return output; + } +} diff --git a/demo/c/raytrace/main.cpp b/demo/c/raytrace/main.cpp index f6fa4cda..66c8ecf9 100644 --- a/demo/c/raytrace/main.cpp +++ b/demo/c/raytrace/main.cpp @@ -88,12 +88,16 @@ int JupiterImage(const char *filename, int width, astro_time_t time) Scene scene(Color(0.0, 0.0, 0.0)); - const double equ_radius_au = JUPITER_EQUATORIAL_RADIUS_KM / KM_PER_AU; - const double pol_radius_au = JUPITER_POLAR_RADIUS_KM / KM_PER_AU; - Spheroid *planet = new Spheroid(equ_radius_au, equ_radius_au, pol_radius_au); + const double km_scale = 10000.0; // makes Jupiter's radius about 7 units. + const double au_scale = km_scale / KM_PER_AU; + + const double equ_radius = JUPITER_EQUATORIAL_RADIUS_KM / km_scale; + const double pol_radius = JUPITER_POLAR_RADIUS_KM / km_scale; + Spheroid *planet = new Spheroid(equ_radius, equ_radius, pol_radius); scene.AddSolidObject(planet); planet->SetFullMatte(Color(1.0, 1.0, 0.7)); - planet->Move(Vector(jupiter.x, jupiter.y, jupiter.z)); + planet->Move(Vector(jupiter.x, jupiter.y, jupiter.z) / au_scale); + planet->SetTag("Jupiter"); // Reorient Jupiter's rotation axis to match the calculated orientation. planet->RotateY(90.0 - axis.dec); // ?? not sure about direction @@ -102,51 +106,54 @@ int JupiterImage(const char *filename, int width, astro_time_t time) // Add Jupiter's moons to the scene. // Io - const double io_radius = IO_RADIUS_KM / KM_PER_AU; + const double io_radius = IO_RADIUS_KM / km_scale; Spheroid *io = new Spheroid(io_radius, io_radius, io_radius); scene.AddSolidObject(io); io->SetFullMatte(Color(1.0, 1.0, 1.0)); // ??? Actual moon colors ??? io->Move(Vector( jm.moon[JM_IO].x + jupiter.x, jm.moon[JM_IO].y + jupiter.y, - jm.moon[JM_IO].z + jupiter.z) + jm.moon[JM_IO].z + jupiter.z) / au_scale ); + io->SetTag("Io"); // Europa - const double eu_radius = EUROPA_RADIUS_KM / KM_PER_AU; + const double eu_radius = EUROPA_RADIUS_KM / km_scale; Spheroid *europa = new Spheroid(eu_radius, eu_radius, eu_radius); scene.AddSolidObject(europa); europa->SetFullMatte(Color(1.0, 1.0, 1.0)); // ??? Actual moon colors ??? europa->Move(Vector( jm.moon[JM_EUROPA].x + jupiter.x, jm.moon[JM_EUROPA].y + jupiter.y, - jm.moon[JM_EUROPA].z + jupiter.z) + jm.moon[JM_EUROPA].z + jupiter.z) / au_scale ); + europa->SetTag("Europa"); // Ganymede - const double gan_radius = GANYMEDE_RADIUS_KM / KM_PER_AU; + const double gan_radius = GANYMEDE_RADIUS_KM / km_scale; Spheroid *ganymede = new Spheroid(gan_radius, gan_radius, gan_radius); scene.AddSolidObject(ganymede); ganymede->SetFullMatte(Color(1.0, 1.0, 1.0)); // ??? Actual moon colors ??? ganymede->Move(Vector( jm.moon[JM_GANYMEDE].x + jupiter.x, jm.moon[JM_GANYMEDE].y + jupiter.y, - jm.moon[JM_GANYMEDE].z + jupiter.z) + jm.moon[JM_GANYMEDE].z + jupiter.z) / au_scale ); // Callisto - const double cal_radius = CALLISTO_RADIUS_KM / KM_PER_AU; + const double cal_radius = CALLISTO_RADIUS_KM / km_scale; Spheroid *callisto = new Spheroid(cal_radius, cal_radius, cal_radius); scene.AddSolidObject(callisto); callisto->SetFullMatte(Color(1.0, 1.0, 1.0)); // ??? Actual moon colors ??? callisto->Move(Vector( jm.moon[JM_CALLISTO].x + jupiter.x, jm.moon[JM_CALLISTO].y + jupiter.y, - jm.moon[JM_CALLISTO].z + jupiter.z) + jm.moon[JM_CALLISTO].z + jupiter.z) / au_scale ); + callisto->SetTag("Callisto"); // Add the Sun as the point light source. - scene.AddLightSource(LightSource(Vector(sun.x, sun.y, sun.z), Color(1.0, 1.0, 1.0))); + scene.AddLightSource(LightSource(Vector(sun.x, sun.y, sun.z) / au_scale, Color(1.0, 1.0, 1.0))); // Aim the camera at the center of Jupiter. // Start with an identity matrix, which leaves the camera pointing in the -z direction, @@ -178,7 +185,7 @@ int JupiterImage(const char *filename, int width, astro_time_t time) // Magnify the image as appropriate for Jupiter's closest approach (opposition). // Render the image. - const double zoom = 200.0; + const double zoom = 300.0; scene.SaveImage(filename, (size_t)width, (size_t)width, zoom, 4); return 0; diff --git a/demo/c/raytrace/run b/demo/c/raytrace/run new file mode 100755 index 00000000..75d2238e --- /dev/null +++ b/demo/c/raytrace/run @@ -0,0 +1,2 @@ +#!/bin/bash +./raytrace jupiter.png 1000 Jupiter 2021-12-11T23:12:00Z diff --git a/demo/c/raytrace/spheroid.cpp b/demo/c/raytrace/spheroid.cpp index 8dc4fac9..cfdd1cf0 100644 --- a/demo/c/raytrace/spheroid.cpp +++ b/demo/c/raytrace/spheroid.cpp @@ -24,7 +24,7 @@ ------------------------------------------------------------------------- Implements class Spheroid. - A spheroid is like a sphere, only it may have three different + A spheroid is like a sphere, only it may have three different diameters in the x-, y-, and z-directions. */ @@ -34,8 +34,8 @@ namespace Imager { void Spheroid::ObjectSpace_AppendAllIntersections( - const Vector& vantage, - const Vector& direction, + const Vector& vantage, + const Vector& direction, IntersectionList& intersectionList) const { double u[2]; @@ -55,9 +55,9 @@ namespace Imager intersection.distanceSquared = displacement.MagnitudeSquared(); intersection.point = vantage + displacement; - // The surface normal vector was calculated by expressing the spheroid as a + // The surface normal vector was calculated by expressing the spheroid as a // function z(x,y) = sqrt(1 - (x/a)^2 - (y/b)^2), - // taking partial derivatives dz/dx = (c*c*x)/(a*a*z), dz/dy = (c*c*y)/(b*b*z), + // taking partial derivatives dz/dx = (c*c*x)/(a*a*z), dz/dy = (c*c*y)/(b*b*z), // and using these to calculate the vectors <1, 0, dz/dx> and <0, 1, dy,dz>. // The normalized cross product of these two vectors yields the surface normal vector. const double x = intersection.point.x; @@ -74,7 +74,7 @@ namespace Imager } else { - // The equation devolves to an ellipse on the xy plane : + // The equation devolves to an ellipse on the xy plane : // (x^2)/(a^2) + (y^2)/(b^2) = 1. intersection.surfaceNormal = Vector(-1.0, -(a2*y)/(b2*x), 0.0).UnitVector(); }