/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.geometry.euclidean.threed;
import java.io.Serializable;
import org.apache.commons.math3.util.FastMath;
/** This class provides conversions related to <a
* href="http://mathworld.wolfram.com/SphericalCoordinates.html">spherical coordinates</a>.
* <p>
* The conventions used here are the mathematical ones, i.e. spherical coordinates are
* related to Cartesian coordinates as follows:
* </p>
* <ul>
* <li>x = r cos(θ) sin(Φ)</li>
* <li>y = r sin(θ) sin(Φ)</li>
* <li>z = r cos(Φ)</li>
* </ul>
* <ul>
* <li>r = √(x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>)</li>
* <li>θ = atan2(y, x)</li>
* <li>Φ = acos(z/r)</li>
* </ul>
* <p>
* r is the radius, θ is the azimuthal angle in the x-y plane and Φ is the polar
* (co-latitude) angle. These conventions are <em>different</em> from the conventions used
* in physics (and in particular in spherical harmonics) where the meanings of θ and
* Φ are reversed.
* </p>
* <p>
* This class provides conversion of coordinates and also of gradient and Hessian
* between spherical and Cartesian coordinates.
* </p>
* @since 3.2
*/
public class SphericalCoordinates implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20130206L;
/** Cartesian coordinates. */
private final Vector3D v;
/** Radius. */
private final double r;
/** Azimuthal angle in the x-y plane θ. */
private final double theta;
/** Polar angle (co-latitude) Φ. */
private final double phi;
/** Jacobian of (r, θ &Phi). */
private double[][] jacobian;
/** Hessian of radius. */
private double[][] rHessian;
/** Hessian of azimuthal angle in the x-y plane θ. */
private double[][] thetaHessian;
/** Hessian of polar (co-latitude) angle Φ. */
private double[][] phiHessian;
/** Build a spherical coordinates transformer from Cartesian coordinates.
* @param v Cartesian coordinates
*/
public SphericalCoordinates(final Vector3D v) {
// Cartesian coordinates
this.v = v;
// remaining spherical coordinates
this.r = v.getNorm();
this.theta = v.getAlpha();
this.phi = FastMath.acos(v.getZ() / r);
}
/** Build a spherical coordinates transformer from spherical coordinates.
* @param r radius
* @param theta azimuthal angle in x-y plane
* @param phi polar (co-latitude) angle
*/
public SphericalCoordinates(final double r, final double theta, final double phi) {
final double cosTheta = FastMath.cos(theta);
final double sinTheta = FastMath.sin(theta);
final double cosPhi = FastMath.cos(phi);
final double sinPhi = FastMath.sin(phi);
// spherical coordinates
this.r = r;
this.theta = theta;
this.phi = phi;
// Cartesian coordinates
this.v = new Vector3D(r * cosTheta * sinPhi,
r * sinTheta * sinPhi,
r * cosPhi);
}
/** Get the Cartesian coordinates.
* @return Cartesian coordinates
*/
public Vector3D getCartesian() {
return v;
}
/** Get the radius.
* @return radius r
* @see #getTheta()
* @see #getPhi()
*/
public double getR() {
return r;
}
/** Get the azimuthal angle in x-y plane.
* @return azimuthal angle in x-y plane θ
* @see #getR()
* @see #getPhi()
*/
public double getTheta() {
return theta;
}
/** Get the polar (co-latitude) angle.
* @return polar (co-latitude) angle Φ
* @see #getR()
* @see #getTheta()
*/
public double getPhi() {
return phi;
}
/** Convert a gradient with respect to spherical coordinates into a gradient
* with respect to Cartesian coordinates.
* @param sGradient gradient with respect to spherical coordinates
* {df/dr, df/dθ, df/dΦ}
* @return gradient with respect to Cartesian coordinates
* {df/dx, df/dy, df/dz}
*/
public double[] toCartesianGradient(final double[] sGradient) {
// lazy evaluation of Jacobian
computeJacobian();
// compose derivatives as gradient^T . J
// the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0
return new double[] {
sGradient[0] * jacobian[0][0] + sGradient[1] * jacobian[1][0] + sGradient[2] * jacobian[2][0],
sGradient[0] * jacobian[0][1] + sGradient[1] * jacobian[1][1] + sGradient[2] * jacobian[2][1],
sGradient[0] * jacobian[0][2] + sGradient[2] * jacobian[2][2]
};
}
/** Convert a Hessian with respect to spherical coordinates into a Hessian
* with respect to Cartesian coordinates.
* <p>
* As Hessian are always symmetric, we use only the lower left part of the provided
* spherical Hessian, so the upper part may not be initialized. However, we still
* do fill up the complete array we create, with guaranteed symmetry.
* </p>
* @param sHessian Hessian with respect to spherical coordinates
* {{d<sup>2</sup>f/dr<sup>2</sup>, d<sup>2</sup>f/drdθ, d<sup>2</sup>f/drdΦ},
* {d<sup>2</sup>f/drdθ, d<sup>2</sup>f/dθ<sup>2</sup>, d<sup>2</sup>f/dθdΦ},
* {d<sup>2</sup>f/drdΦ, d<sup>2</sup>f/dθdΦ, d<sup>2</sup>f/dΦ<sup>2</sup>}
* @param sGradient gradient with respect to spherical coordinates
* {df/dr, df/dθ, df/dΦ}
* @return Hessian with respect to Cartesian coordinates
* {{d<sup>2</sup>f/dx<sup>2</sup>, d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dxdz},
* {d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dy<sup>2</sup>, d<sup>2</sup>f/dydz},
* {d<sup>2</sup>f/dxdz, d<sup>2</sup>f/dydz, d<sup>2</sup>f/dz<sup>2</sup>}}
*/
public double[][] toCartesianHessian(final double[][] sHessian, final double[] sGradient) {
computeJacobian();
computeHessians();
// compose derivative as J^T . H_f . J + df/dr H_r + df/dtheta H_theta + df/dphi H_phi
// the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0
// and H_theta is only a 2x2 matrix as it does not depend on z
final double[][] hj = new double[3][3];
final double[][] cHessian = new double[3][3];
// compute H_f . J
// beware we use ONLY the lower-left part of sHessian
hj[0][0] = sHessian[0][0] * jacobian[0][0] + sHessian[1][0] * jacobian[1][0] + sHessian[2][0] * jacobian[2][0];
hj[0][1] = sHessian[0][0] * jacobian[0][1] + sHessian[1][0] * jacobian[1][1] + sHessian[2][0] * jacobian[2][1];
hj[0][2] = sHessian[0][0] * jacobian[0][2] + sHessian[2][0] * jacobian[2][2];
hj[1][0] = sHessian[1][0] * jacobian[0][0] + sHessian[1][1] * jacobian[1][0] + sHessian[2][1] * jacobian[2][0];
hj[1][1] = sHessian[1][0] * jacobian[0][1] + sHessian[1][1] * jacobian[1][1] + sHessian[2][1] * jacobian[2][1];
// don't compute hj[1][2] as it is not used below
hj[2][0] = sHessian[2][0] * jacobian[0][0] + sHessian[2][1] * jacobian[1][0] + sHessian[2][2] * jacobian[2][0];
hj[2][1] = sHessian[2][0] * jacobian[0][1] + sHessian[2][1] * jacobian[1][1] + sHessian[2][2] * jacobian[2][1];
hj[2][2] = sHessian[2][0] * jacobian[0][2] + sHessian[2][2] * jacobian[2][2];
// compute lower-left part of J^T . H_f . J
cHessian[0][0] = jacobian[0][0] * hj[0][0] + jacobian[1][0] * hj[1][0] + jacobian[2][0] * hj[2][0];
cHessian[1][0] = jacobian[0][1] * hj[0][0] + jacobian[1][1] * hj[1][0] + jacobian[2][1] * hj[2][0];
cHessian[2][0] = jacobian[0][2] * hj[0][0] + jacobian[2][2] * hj[2][0];
cHessian[1][1] = jacobian[0][1] * hj[0][1] + jacobian[1][1] * hj[1][1] + jacobian[2][1] * hj[2][1];
cHessian[2][1] = jacobian[0][2] * hj[0][1] + jacobian[2][2] * hj[2][1];
cHessian[2][2] = jacobian[0][2] * hj[0][2] + jacobian[2][2] * hj[2][2];
// add gradient contribution
cHessian[0][0] += sGradient[0] * rHessian[0][0] + sGradient[1] * thetaHessian[0][0] + sGradient[2] * phiHessian[0][0];
cHessian[1][0] += sGradient[0] * rHessian[1][0] + sGradient[1] * thetaHessian[1][0] + sGradient[2] * phiHessian[1][0];
cHessian[2][0] += sGradient[0] * rHessian[2][0] + sGradient[2] * phiHessian[2][0];
cHessian[1][1] += sGradient[0] * rHessian[1][1] + sGradient[1] * thetaHessian[1][1] + sGradient[2] * phiHessian[1][1];
cHessian[2][1] += sGradient[0] * rHessian[2][1] + sGradient[2] * phiHessian[2][1];
cHessian[2][2] += sGradient[0] * rHessian[2][2] + sGradient[2] * phiHessian[2][2];
// ensure symmetry
cHessian[0][1] = cHessian[1][0];
cHessian[0][2] = cHessian[2][0];
cHessian[1][2] = cHessian[2][1];
return cHessian;
}
/** Lazy evaluation of (r, θ, φ) Jacobian.
*/
private void computeJacobian() {
if (jacobian == null) {
// intermediate variables
final double x = v.getX();
final double y = v.getY();
final double z = v.getZ();
final double rho2 = x * x + y * y;
final double rho = FastMath.sqrt(rho2);
final double r2 = rho2 + z * z;
jacobian = new double[3][3];
// row representing the gradient of r
jacobian[0][0] = x / r;
jacobian[0][1] = y / r;
jacobian[0][2] = z / r;
// row representing the gradient of theta
jacobian[1][0] = -y / rho2;
jacobian[1][1] = x / rho2;
// jacobian[1][2] is already set to 0 at allocation time
// row representing the gradient of phi
jacobian[2][0] = x * z / (rho * r2);
jacobian[2][1] = y * z / (rho * r2);
jacobian[2][2] = -rho / r2;
}
}
/** Lazy evaluation of Hessians.
*/
private void computeHessians() {
if (rHessian == null) {
// intermediate variables
final double x = v.getX();
final double y = v.getY();
final double z = v.getZ();
final double x2 = x * x;
final double y2 = y * y;
final double z2 = z * z;
final double rho2 = x2 + y2;
final double rho = FastMath.sqrt(rho2);
final double r2 = rho2 + z2;
final double xOr = x / r;
final double yOr = y / r;
final double zOr = z / r;
final double xOrho2 = x / rho2;
final double yOrho2 = y / rho2;
final double xOr3 = xOr / r2;
final double yOr3 = yOr / r2;
final double zOr3 = zOr / r2;
// lower-left part of Hessian of r
rHessian = new double[3][3];
rHessian[0][0] = y * yOr3 + z * zOr3;
rHessian[1][0] = -x * yOr3;
rHessian[2][0] = -z * xOr3;
rHessian[1][1] = x * xOr3 + z * zOr3;
rHessian[2][1] = -y * zOr3;
rHessian[2][2] = x * xOr3 + y * yOr3;
// upper-right part is symmetric
rHessian[0][1] = rHessian[1][0];
rHessian[0][2] = rHessian[2][0];
rHessian[1][2] = rHessian[2][1];
// lower-left part of Hessian of azimuthal angle theta
thetaHessian = new double[2][2];
thetaHessian[0][0] = 2 * xOrho2 * yOrho2;
thetaHessian[1][0] = yOrho2 * yOrho2 - xOrho2 * xOrho2;
thetaHessian[1][1] = -2 * xOrho2 * yOrho2;
// upper-right part is symmetric
thetaHessian[0][1] = thetaHessian[1][0];
// lower-left part of Hessian of polar (co-latitude) angle phi
final double rhor2 = rho * r2;
final double rho2r2 = rho * rhor2;
final double rhor4 = rhor2 * r2;
final double rho3r4 = rhor4 * rho2;
final double r2P2rho2 = 3 * rho2 + z2;
phiHessian = new double[3][3];
phiHessian[0][0] = z * (rho2r2 - x2 * r2P2rho2) / rho3r4;
phiHessian[1][0] = -x * y * z * r2P2rho2 / rho3r4;
phiHessian[2][0] = x * (rho2 - z2) / rhor4;
phiHessian[1][1] = z * (rho2r2 - y2 * r2P2rho2) / rho3r4;
phiHessian[2][1] = y * (rho2 - z2) / rhor4;
phiHessian[2][2] = 2 * rho * zOr3 / r;
// upper-right part is symmetric
phiHessian[0][1] = phiHessian[1][0];
phiHessian[0][2] = phiHessian[2][0];
phiHessian[1][2] = phiHessian[2][1];
}
}
/**
* Replace the instance with a data transfer object for serialization.
* @return data transfer object that will be serialized
*/
private Object writeReplace() {
return new DataTransferObject(v.getX(), v.getY(), v.getZ());
}
/** Internal class used only for serialization. */
private static class DataTransferObject implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20130206L;
/** Abscissa.
* @serial
*/
private final double x;
/** Ordinate.
* @serial
*/
private final double y;
/** Height.
* @serial
*/
private final double z;
/** Simple constructor.
* @param x abscissa
* @param y ordinate
* @param z height
*/
public DataTransferObject(final double x, final double y, final double z) {
this.x = x;
this.y = y;
this.z = z;
}
/** Replace the deserialized data transfer object with a {@link SphericalCoordinates}.
* @return replacement {@link SphericalCoordinates}
*/
private Object readResolve() {
return new SphericalCoordinates(new Vector3D(x, y, z));
}
}
}