package jpbrt.shapes;
import jpbrt.core.BBox;
import jpbrt.core.DifferentialGeometry;
import jpbrt.core.Global;
import jpbrt.core.Intersection;
import jpbrt.core.Normal;
import jpbrt.core.Point;
import jpbrt.core.Ray;
import jpbrt.core.Shape;
import jpbrt.core.Transform;
import jpbrt.core.Vector;
public class Sphere extends Shape
{
private double radius;
private double phiMax;
private double zmin, zmax;
private double thetaMin, thetaMax;
public Sphere(Transform o2w, Transform w2o, boolean ro,
double radius, double zmin, double zmax, double phiMax)
{
super(o2w, w2o, ro);
assert(radius > 0.0);
this.radius = radius;
this.zmin = Global.clamp(Math.min(zmin, zmax), -radius, radius);
this.zmax = Global.clamp(Math.max(zmin, zmax), -radius, radius);
this.thetaMin = Math.acos( Global.clamp(this.zmin / radius, -1.0, 1.0) );
this.thetaMax = Math.acos( Global.clamp(this.zmax / radius, -1.0, 1.0) );
this.phiMax = Global.clamp(phiMax, 0, 2 * Math.PI);
}
@Override
public BBox objectBound()
{
return new BBox( new Point(-radius, -radius, zmin),
new Point(radius, radius, zmax) );
}
@Override
public Intersection intersect(Ray r)
{
Intersection result = new Intersection();
Ray ray = worldToObject.transform(r);
// compute quadratic sphere coefficients
double a = ray.d.x*ray.d.x +ray.d.y*ray.d.y + ray.d.z*ray.d.z;
double b = 2 * (ray.d.x*ray.o.x + ray.d.y*ray.o.y + ray.d.z* ray.o.z);
double c = ray.o.x*ray.o.x + ray.o.y*ray.o.y + ray.o.z*ray.o.z - radius*radius;
// solve quadratic equation for t values
double[] t = Global.quadratic(a, b, c);
if (t == null)
return result;
// compute intersection distance along ray
if (t[0] > ray.maxt || t[1] < ray.mint)
return result;
double thit = t[0];
if (t[0] < ray.mint)
{
thit = t[1];
if (thit > ray.maxt)
return result;
}
// compute sphere hit position and phi
Point phit = ray.eval(thit);
if (phit.x == 0.0 && phit.y == 0.0)
phit.x = 1e-5 * radius;
double phi = Math.atan2(phit.y, phit.x);
if (phi < 0.0)
phi += 2 * Math.PI;
// test sphere intersection against clipping parameters
if ( (zmin > -radius && phit.z < zmin) ||
(zmax > radius && phit.z > zmax) || phi > phiMax )
{
if (thit == t[1] || t[1] > ray.maxt)
return result;
thit = t[1];
// compute sphere hit position and phi
phit = ray.eval(thit);
if (phit.x == 0.0 && phit.y == 0.0)
phit.x = 1e-5 * radius;
phi = Math.atan2(phit.y, phit.x);
if (phi < 0.0)
phi += 2 * Math.PI;
if ( (zmin > -radius && phit.z < zmin) ||
(zmax < radius && phit.z > zmax) || phi > phiMax )
return result;
}
// find parameteric representation of sphere hit
double u = phi / phiMax;
double theta = Math.acos( Global.clamp(phit.z / radius, -1.0, 1.0) );
double v = (theta - thetaMax) / (thetaMax - thetaMin);
// compute sphere dpdu and dpdv
double zradius = Math.sqrt(phit.x * phit.x + phit.y * phit.y);
double invzradius = 1.0 / zradius;
double cosphi = phit.x / invzradius;
double sinphi = phit.y / invzradius;
Vector dpdu = new Vector(-phiMax * phit.y, phiMax * phit.x, 0.0);
Vector dpdv = new Vector(phit.z * cosphi, phit.z * sinphi, -radius * Math.sin(theta)).mulLocal(thetaMax - thetaMin);
// compute sphere dndu and dndv
Vector d2Pduu = new Vector(phit.x, phit.y, 0).mulLocal(-phiMax * phiMax);
Vector d2Pduv = new Vector(-sinphi, cosphi, 0).mulLocal( (thetaMax - thetaMin) * phit.z * phiMax );
Vector d2Pdvv = new Vector(phit.x, phit.y, phit.z).mulLocal( -(thetaMax - thetaMin) * (thetaMax - thetaMin) );
// compute coefficients for fundamental forms
double E = dpdu.dot(dpdu);
double F = dpdu.dot(dpdv);
double G = dpdv.dot(dpdv);
Vector N = dpdu.cross(dpdv).normalizeLocal();
double e = N.dot(d2Pduu);
double f = N.dot(d2Pduv);
double g = N.dot(d2Pdvv);
// compute dndu and dndv from fundamental form coefficients
double invEGF2 = 1.0 / (E*G - F*F);
Normal dndu = new Normal(dpdu.mul( (f*F - e*G) * invEGF2 ).addLocal( dpdv.mul((e*F - f*E) * invEGF2) ));
Normal dndv = new Normal(dpdu.mul( (g*F - f*G) * invEGF2 ).addLocal( dpdv.mul((f*F - g*E) * invEGF2) ));
// initialize DifferentialGeometry from parametric information
result.dg = new DifferentialGeometry(
objectToWorld.transform(phit),
objectToWorld.transform(dpdu),
objectToWorld.transform(dpdv),
objectToWorld.transform(dndu),
objectToWorld.transform(dndv),
u, v, this);
// update tHit for quadric intersection
result.tHit = thit;
// compute rayEpsilon for quadric intersection
result.rayEpsilon = 5e-4 * thit;
result.hit = true;
return result;
}
@Override
public boolean intersectP(Ray r)
{
// transform ray to object space
Ray ray = worldToObject.transform(r);
// compute quadratic sphere coefficients
double a = ray.d.x*ray.d.x +ray.d.y*ray.d.y + ray.d.z*ray.d.z;
double b = 2 * (ray.d.x*ray.o.x + ray.d.y*ray.o.y + ray.d.z* ray.o.z);
double c = ray.o.x*ray.o.x + ray.o.y*ray.o.y + ray.o.z*ray.o.z - radius*radius;
// solve quadratic equation for t values
double[] t = Global.quadratic(a, b, c);
if (t == null)
return false;
// compute intersection distance along ray
if (t[0] > ray.maxt || t[1] < ray.mint)
return false;
double thit = t[0];
if (t[0] < ray.mint)
{
thit = t[1];
if (thit > ray.maxt)
return false;
}
// compute sphere hit position and phi
Point phit = ray.eval(thit);
if (phit.x == 0.0 && phit.y == 0.0)
phit.x = 1e-5 * radius;
double phi = Math.atan2(phit.y, phit.x);
if (phi < 0.0)
phi += 2 * Math.PI;
// test sphere intersection against clipping parameters
if ( (zmin > -radius && phit.z < zmin) ||
(zmax > radius && phit.z > zmax) || phi > phiMax )
{
if (thit == t[1] || t[1] > ray.maxt)
return false;
thit = t[1];
// compute sphere hit position and phi
phit = ray.eval(thit);
if (phit.x == 0.0 && phit.y == 0.0)
phit.x = 1e-5 * radius;
phi = Math.atan2(phit.y, phit.x);
if (phi < 0.0)
phi += 2 * Math.PI;
if ( (zmin > -radius && phit.z < zmin) ||
(zmax < radius && phit.z > zmax) || phi > phiMax )
return false;
}
return true;
}
@Override
public double area()
{
return phiMax * radius * (zmax - zmin);
}
}