Indoor/geo/Sphere3.h

/*
 * © Copyright 2014 – Urheberrechtshinweis
 * Alle Rechte vorbehalten / All Rights Reserved
 *
 * Programmcode ist urheberrechtlich geschuetzt.
 * Das Urheberrecht liegt, soweit nicht ausdruecklich anders gekennzeichnet, bei Frank Ebner.
 * Keine Verwendung ohne explizite Genehmigung.
 * (vgl. § 106 ff UrhG / § 97 UrhG)
 */

#ifndef GEO_SPHERE3_H
#define GEO_SPHERE3_H

#include "Point3.h"
#include "Ray3.h"

#include <vector>


struct Sphere3 {

	Point3 center;

	float radius;

	/** empty ctor */
	Sphere3() : center(), radius(0) {;}

	/** ctor */
	Sphere3(const Point3 center, const float radius) : center(center), radius(radius) {;}



	/** does this sphere contain the given sphere? */
	bool contains(const Sphere3& o) const {
		return (o.center.getDistance(this->center) + o.radius) <= this->radius;
	}

	/** does the sphere contain the given point? */
	bool contains(const Point3& p) const {
		return center.getDistance(p) <= radius;
	}

	/** does the sphere intersect with the given sphere? */
	bool intersects(const Sphere3& other) const {
		return center.getDistance(other.center) < (radius + other.radius);
	}

	/** does the sphere intersect with the given ray? */
	bool intersects(const Ray3& ray) const {

		// if the sphere contains the ray's start -> done
		//if (contains(ray.start)) {return true;}


		// https://stackoverflow.com/questions/6533856/ray-sphere-intersection

		/*
		const float xA = ray.start.x;
		const float yA = ray.start.y;
		const float zA = ray.start.z;

		const float xB = (ray.start + ray.dir).x;
		const float yB = (ray.start + ray.dir).y;
		const float zB = (ray.start + ray.dir).z;

		const float xC = center.x;
		const float yC = center.y;
		const float zC = center.z;

		const float a = (xB-xA)*(xB-xA)+(yB-yA)*(yB-yA)+(zB-zA)*(zB-zA);
		const float b = 2*((xB-xA)*(xA-xC)+(yB-yA)*(yA-yC)+(zB-zA)*(zA-zC));
		const float c = (xA-xC)*(xA-xC)+(yA-yC)*(yA-yC)+(zA-zC)*(zA-zC)-(radius*radius);

		const float delta = (b*b) - (4*a*c);

		// intersection?
		return delta >= 0;
	*/


		// http://www.ccs.neu.edu/home/fell/CSU540/programs/RayTracingFormulas.htm
		/*
		const float x0 = ray.start.x;
		const float y0 = ray.start.y;
		const float z0 = ray.start.z;

		const float cx = center.x;
		const float cy = center.y;
		const float cz = center.z;

		const float dx = ray.dir.x;
		const float dy = ray.dir.y;
		const float dz = ray.dir.z;

		const float a = dx*dx + dy*dy + dz*dz;
		const float b = 2*dx*(x0-cx) +  2*dy*(y0-cy) +  2*dz*(z0-cz);
		const float c = cx*cx + cy*cy + cz*cz + x0*x0 + y0*y0 + z0*z0 + -2*(cx*x0 + cy*y0 + cz*z0) - radius*radius;

		const float discriminant = (b*b) - (4*a*c);
		return discriminant >= 0;
		*/


		//https://gamedev.stackexchange.com/questions/96459/fast-ray-sphere-collision-code

		const Point3 m = ray.start - center;
		const float b = dot(m, ray.dir);
		const float c = dot(m, m) - radius*radius;

		// Exit if r’s origin outside s (c > 0) and r pointing away from s (b > 0)
		if (c > 0.0f && b > 0.0f) {return false;}
		const float discr = b*b - c;

		// A negative discriminant corresponds to ray missing sphere
		if (discr < 0.0f) {return false;}

		return true;

//		// Ray now found to intersect sphere, compute smallest t value of intersection
//		t = -b - Sqrt(discr);

//		// If t is negative, ray started inside sphere so clamp t to zero
//		if (t < 0.0f) t = 0.0f;
//		q = p + t * d;

//		return 1;


	}

	/** configure this sphere to contain the given point-set */
	void adjustToPointSet(const std::vector<Point3>& lst) {

		// NOT OPTIMAL but fast

		// calculate the point set's center
		Point3 sum(0,0,0);
		for (const Point3 p : lst) {sum += p;}
		const Point3 center = sum / lst.size();

		// calculate the sphere's radius (maximum distance from center
		float radius = 0;
		for (const Point3 p : lst) {
			const float dist = center.getDistance(p);
			if (dist > radius) {radius = dist;}
		}

		this->radius = radius;
		this->center = center;

	}

public:

	/** create a sphere that fully contains the given point-set */
	static Sphere3 around(const std::vector<Point3>& lst) {
		Sphere3 sphere;
		sphere.adjustToPointSet(lst);
		return sphere;
	}

	/** combine two spheres into a new one containing both */
	static Sphere3 join(const Sphere3& a, const Sphere3& b) {

		// if one already contains the other, just return it as-is
		if (a.contains(b)) {return a;}
		if (b.contains(a)) {return b;}

		// calculate the new center and radius
//		Point3 newCen = (a.center + b.center) / 2;
//		float newRad = (a.center.getDistance(b.center) + std::max(a.radius, b.radius) * 2) / 2;
//		return Sphere3(newCen, newRad);

//		Point3 newCen = (a.center*a.radius + b.center*b.radius) / (a.radius+b.radius);
//		float newRad = (a.center.getDistance(b.center) + a.radius + b.radius) / 2 * 1.02f;	// slightly larger to prevent rounding issues
//		Sphere3 res(newCen, newRad);

		// create both maximum ends
		const Point3 out1 = a.center + (a.center-b.center).normalized() * a.radius;
		const Point3 out2 = b.center + (b.center-a.center).normalized() * b.radius;

		// center is within both ends, so is the radius
		Point3 newCen = (out1 + out2) / 2;
		float newRad = out1.getDistance(out2) / 2 * 1.0001;	// slightly larger to prevent rounding issues
		Sphere3 res(newCen, newRad);

		// check
		Assert::isTrue(res.contains(a), "sphere joining error. base-spheres are not contained.");
		Assert::isTrue(res.contains(b), "sphere joining error. base-spheres are not contained.");

		return res;

	}

};

#endif // GEO_SPHERE3_H