SphericalCoordinates.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.geometry.euclidean.threed;


  20. import java.io.Serializable;

  22. import org.apache.commons.math3.util.FastMath;

  24. /** This class provides conversions related to <a
  25.  * href="http://mathworld.wolfram.com/SphericalCoordinates.html">spherical coordinates</a>.
  26.  * <p>
  27.  * The conventions used here are the mathematical ones, i.e. spherical coordinates are
  28.  * related to Cartesian coordinates as follows:
  29.  * </p>
  30.  * <ul>
  31.  *   <li>x = r cos(&theta;) sin(&Phi;)</li>
  32.  *   <li>y = r sin(&theta;) sin(&Phi;)</li>
  33.  *   <li>z = r cos(&Phi;)</li>
  34.  * </ul>
  35.  * <ul>
  36.  *   <li>r       = &radic;(x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>)</li>
  37.  *   <li>&theta; = atan2(y, x)</li>
  38.  *   <li>&Phi;   = acos(z/r)</li>
  39.  * </ul>
  40.  * <p>
  41.  * r is the radius, &theta; is the azimuthal angle in the x-y plane and &Phi; is the polar
  42.  * (co-latitude) angle. These conventions are <em>different</em> from the conventions used
  43.  * in physics (and in particular in spherical harmonics) where the meanings of &theta; and
  44.  * &Phi; are reversed.
  45.  * </p>
  46.  * <p>
  47.  * This class provides conversion of coordinates and also of gradient and Hessian
  48.  * between spherical and Cartesian coordinates.
  49.  * </p>
  50.  * @since 3.2
  51.  */
  52. public class SphericalCoordinates implements Serializable {

  54.     /** Serializable UID. */
  55.     private static final long serialVersionUID = 20130206L;

  57.     /** Cartesian coordinates. */
  58.     private final Vector3D v;

  60.     /** Radius. */
  61.     private final double r;

  63.     /** Azimuthal angle in the x-y plane &theta;. */
  64.     private final double theta;

  66.     /** Polar angle (co-latitude) &Phi;. */
  67.     private final double phi;

  69.     /** Jacobian of (r, &theta; &Phi). */
  70.     private double[][] jacobian;

  72.     /** Hessian of radius. */
  73.     private double[][] rHessian;

  75.     /** Hessian of azimuthal angle in the x-y plane &theta;. */
  76.     private double[][] thetaHessian;

  78.     /** Hessian of polar (co-latitude) angle &Phi;. */
  79.     private double[][] phiHessian;

  81.     /** Build a spherical coordinates transformer from Cartesian coordinates.
  82.      * @param v Cartesian coordinates
  83.      */
  84.     public SphericalCoordinates(final Vector3D v) {

  86.         // Cartesian coordinates
  87.         this.v = v;

  89.         // remaining spherical coordinates
  90.         this.r     = v.getNorm();
  91.         this.theta = v.getAlpha();
  92.         this.phi   = FastMath.acos(v.getZ() / r);

  94.     }

  96.     /** Build a spherical coordinates transformer from spherical coordinates.
  97.      * @param r radius
  98.      * @param theta azimuthal angle in x-y plane
  99.      * @param phi polar (co-latitude) angle
  100.      */
  101.     public SphericalCoordinates(final double r, final double theta, final double phi) {

  103.         final double cosTheta = FastMath.cos(theta);
  104.         final double sinTheta = FastMath.sin(theta);
  105.         final double cosPhi   = FastMath.cos(phi);
  106.         final double sinPhi   = FastMath.sin(phi);

  108.         // spherical coordinates
  109.         this.r     = r;
  110.         this.theta = theta;
  111.         this.phi   = phi;

  113.         // Cartesian coordinates
  114.         this.v  = new Vector3D(r * cosTheta * sinPhi,
  115.                                r * sinTheta * sinPhi,
  116.                                r * cosPhi);

  118.     }

  120.     /** Get the Cartesian coordinates.
  121.      * @return Cartesian coordinates
  122.      */
  123.     public Vector3D getCartesian() {
  124.         return v;
  125.     }

  127.     /** Get the radius.
  128.      * @return radius r
  129.      * @see #getTheta()
  130.      * @see #getPhi()
  131.      */
  132.     public double getR() {
  133.         return r;
  134.     }

  136.     /** Get the azimuthal angle in x-y plane.
  137.      * @return azimuthal angle in x-y plane &theta;
  138.      * @see #getR()
  139.      * @see #getPhi()
  140.      */
  141.     public double getTheta() {
  142.         return theta;
  143.     }

  145.     /** Get the polar (co-latitude) angle.
  146.      * @return polar (co-latitude) angle &Phi;
  147.      * @see #getR()
  148.      * @see #getTheta()
  149.      */
  150.     public double getPhi() {
  151.         return phi;
  152.     }

  154.     /** Convert a gradient with respect to spherical coordinates into a gradient
  155.      * with respect to Cartesian coordinates.
  156.      * @param sGradient gradient with respect to spherical coordinates
  157.      * {df/dr, df/d&theta;, df/d&Phi;}
  158.      * @return gradient with respect to Cartesian coordinates
  159.      * {df/dx, df/dy, df/dz}
  160.      */
  161.     public double[] toCartesianGradient(final double[] sGradient) {

  163.         // lazy evaluation of Jacobian
  164.         computeJacobian();

  166.         // compose derivatives as gradient^T . J
  167.         // the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0
  168.         return new double[] {
  169.             sGradient[0] * jacobian[0][0] + sGradient[1] * jacobian[1][0] + sGradient[2] * jacobian[2][0],
  170.             sGradient[0] * jacobian[0][1] + sGradient[1] * jacobian[1][1] + sGradient[2] * jacobian[2][1],
  171.             sGradient[0] * jacobian[0][2]                                 + sGradient[2] * jacobian[2][2]
  172.         };

  174.     }

  176.     /** Convert a Hessian with respect to spherical coordinates into a Hessian
  177.      * with respect to Cartesian coordinates.
  178.      * <p>
  179.      * As Hessian are always symmetric, we use only the lower left part of the provided
  180.      * spherical Hessian, so the upper part may not be initialized. However, we still
  181.      * do fill up the complete array we create, with guaranteed symmetry.
  182.      * </p>
  183.      * @param sHessian Hessian with respect to spherical coordinates
  184.      * {{d<sup>2</sup>f/dr<sup>2</sup>, d<sup>2</sup>f/drd&theta;, d<sup>2</sup>f/drd&Phi;},
  185.      *  {d<sup>2</sup>f/drd&theta;, d<sup>2</sup>f/d&theta;<sup>2</sup>, d<sup>2</sup>f/d&theta;d&Phi;},
  186.      *  {d<sup>2</sup>f/drd&Phi;, d<sup>2</sup>f/d&theta;d&Phi;, d<sup>2</sup>f/d&Phi;<sup>2</sup>}
  187.      * @param sGradient gradient with respect to spherical coordinates
  188.      * {df/dr, df/d&theta;, df/d&Phi;}
  189.      * @return Hessian with respect to Cartesian coordinates
  190.      * {{d<sup>2</sup>f/dx<sup>2</sup>, d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dxdz},
  191.      *  {d<sup>2</sup>f/dxdy, d<sup>2</sup>f/dy<sup>2</sup>, d<sup>2</sup>f/dydz},
  192.      *  {d<sup>2</sup>f/dxdz, d<sup>2</sup>f/dydz, d<sup>2</sup>f/dz<sup>2</sup>}}
  193.      */
  194.     public double[][] toCartesianHessian(final double[][] sHessian, final double[] sGradient) {

  196.         computeJacobian();
  197.         computeHessians();

  199.         // compose derivative as J^T . H_f . J + df/dr H_r + df/dtheta H_theta + df/dphi H_phi
  200.         // the expressions have been simplified since we know jacobian[1][2] = dTheta/dZ = 0
  201.         // and H_theta is only a 2x2 matrix as it does not depend on z
  202.         final double[][] hj = new double[3][3];
  203.         final double[][] cHessian = new double[3][3];

  205.         // compute H_f . J
  206.         // beware we use ONLY the lower-left part of sHessian
  207.         hj[0][0] = sHessian[0][0] * jacobian[0][0] + sHessian[1][0] * jacobian[1][0] + sHessian[2][0] * jacobian[2][0];
  208.         hj[0][1] = sHessian[0][0] * jacobian[0][1] + sHessian[1][0] * jacobian[1][1] + sHessian[2][0] * jacobian[2][1];
  209.         hj[0][2] = sHessian[0][0] * jacobian[0][2]                                   + sHessian[2][0] * jacobian[2][2];
  210.         hj[1][0] = sHessian[1][0] * jacobian[0][0] + sHessian[1][1] * jacobian[1][0] + sHessian[2][1] * jacobian[2][0];
  211.         hj[1][1] = sHessian[1][0] * jacobian[0][1] + sHessian[1][1] * jacobian[1][1] + sHessian[2][1] * jacobian[2][1];
  212.         // don't compute hj[1][2] as it is not used below
  213.         hj[2][0] = sHessian[2][0] * jacobian[0][0] + sHessian[2][1] * jacobian[1][0] + sHessian[2][2] * jacobian[2][0];
  214.         hj[2][1] = sHessian[2][0] * jacobian[0][1] + sHessian[2][1] * jacobian[1][1] + sHessian[2][2] * jacobian[2][1];
  215.         hj[2][2] = sHessian[2][0] * jacobian[0][2]                                   + sHessian[2][2] * jacobian[2][2];

  217.         // compute lower-left part of J^T . H_f . J
  218.         cHessian[0][0] = jacobian[0][0] * hj[0][0] + jacobian[1][0] * hj[1][0] + jacobian[2][0] * hj[2][0];
  219.         cHessian[1][0] = jacobian[0][1] * hj[0][0] + jacobian[1][1] * hj[1][0] + jacobian[2][1] * hj[2][0];
  220.         cHessian[2][0] = jacobian[0][2] * hj[0][0]                             + jacobian[2][2] * hj[2][0];
  221.         cHessian[1][1] = jacobian[0][1] * hj[0][1] + jacobian[1][1] * hj[1][1] + jacobian[2][1] * hj[2][1];
  222.         cHessian[2][1] = jacobian[0][2] * hj[0][1]                             + jacobian[2][2] * hj[2][1];
  223.         cHessian[2][2] = jacobian[0][2] * hj[0][2]                             + jacobian[2][2] * hj[2][2];

  225.         // add gradient contribution
  226.         cHessian[0][0] += sGradient[0] * rHessian[0][0] + sGradient[1] * thetaHessian[0][0] + sGradient[2] * phiHessian[0][0];
  227.         cHessian[1][0] += sGradient[0] * rHessian[1][0] + sGradient[1] * thetaHessian[1][0] + sGradient[2] * phiHessian[1][0];
  228.         cHessian[2][0] += sGradient[0] * rHessian[2][0]                                     + sGradient[2] * phiHessian[2][0];
  229.         cHessian[1][1] += sGradient[0] * rHessian[1][1] + sGradient[1] * thetaHessian[1][1] + sGradient[2] * phiHessian[1][1];
  230.         cHessian[2][1] += sGradient[0] * rHessian[2][1]                                     + sGradient[2] * phiHessian[2][1];
  231.         cHessian[2][2] += sGradient[0] * rHessian[2][2]                                     + sGradient[2] * phiHessian[2][2];

  233.         // ensure symmetry
  234.         cHessian[0][1] = cHessian[1][0];
  235.         cHessian[0][2] = cHessian[2][0];
  236.         cHessian[1][2] = cHessian[2][1];

  238.         return cHessian;

  240.     }

  242.     /** Lazy evaluation of (r, &theta;, &phi;) Jacobian.
  243.      */
  244.     private void computeJacobian() {
  245.         if (jacobian == null) {

  247.             // intermediate variables
  248.             final double x    = v.getX();
  249.             final double y    = v.getY();
  250.             final double z    = v.getZ();
  251.             final double rho2 = x * x + y * y;
  252.             final double rho  = FastMath.sqrt(rho2);
  253.             final double r2   = rho2 + z * z;

  255.             jacobian = new double[3][3];

  257.             // row representing the gradient of r
  258.             jacobian[0][0] = x / r;
  259.             jacobian[0][1] = y / r;
  260.             jacobian[0][2] = z / r;

  262.             // row representing the gradient of theta
  263.             jacobian[1][0] = -y / rho2;
  264.             jacobian[1][1] =  x / rho2;
  265.             // jacobian[1][2] is already set to 0 at allocation time

  267.             // row representing the gradient of phi
  268.             jacobian[2][0] = x * z / (rho * r2);
  269.             jacobian[2][1] = y * z / (rho * r2);
  270.             jacobian[2][2] = -rho / r2;

  272.         }
  273.     }

  275.     /** Lazy evaluation of Hessians.
  276.      */
  277.     private void computeHessians() {

  279.         if (rHessian == null) {

  281.             // intermediate variables
  282.             final double x      = v.getX();
  283.             final double y      = v.getY();
  284.             final double z      = v.getZ();
  285.             final double x2     = x * x;
  286.             final double y2     = y * y;
  287.             final double z2     = z * z;
  288.             final double rho2   = x2 + y2;
  289.             final double rho    = FastMath.sqrt(rho2);
  290.             final double r2     = rho2 + z2;
  291.             final double xOr    = x / r;
  292.             final double yOr    = y / r;
  293.             final double zOr    = z / r;
  294.             final double xOrho2 = x / rho2;
  295.             final double yOrho2 = y / rho2;
  296.             final double xOr3   = xOr / r2;
  297.             final double yOr3   = yOr / r2;
  298.             final double zOr3   = zOr / r2;

  300.             // lower-left part of Hessian of r
  301.             rHessian = new double[3][3];
  302.             rHessian[0][0] = y * yOr3 + z * zOr3;
  303.             rHessian[1][0] = -x * yOr3;
  304.             rHessian[2][0] = -z * xOr3;
  305.             rHessian[1][1] = x * xOr3 + z * zOr3;
  306.             rHessian[2][1] = -y * zOr3;
  307.             rHessian[2][2] = x * xOr3 + y * yOr3;

  309.             // upper-right part is symmetric
  310.             rHessian[0][1] = rHessian[1][0];
  311.             rHessian[0][2] = rHessian[2][0];
  312.             rHessian[1][2] = rHessian[2][1];

  314.             // lower-left part of Hessian of azimuthal angle theta
  315.             thetaHessian = new double[2][2];
  316.             thetaHessian[0][0] = 2 * xOrho2 * yOrho2;
  317.             thetaHessian[1][0] = yOrho2 * yOrho2 - xOrho2 * xOrho2;
  318.             thetaHessian[1][1] = -2 * xOrho2 * yOrho2;

  320.             // upper-right part is symmetric
  321.             thetaHessian[0][1] = thetaHessian[1][0];

  323.             // lower-left part of Hessian of polar (co-latitude) angle phi
  324.             final double rhor2       = rho * r2;
  325.             final double rho2r2      = rho * rhor2;
  326.             final double rhor4       = rhor2 * r2;
  327.             final double rho3r4      = rhor4 * rho2;
  328.             final double r2P2rho2    = 3 * rho2 + z2;
  329.             phiHessian = new double[3][3];
  330.             phiHessian[0][0] = z * (rho2r2 - x2 * r2P2rho2) / rho3r4;
  331.             phiHessian[1][0] = -x * y * z * r2P2rho2 / rho3r4;
  332.             phiHessian[2][0] = x * (rho2 - z2) / rhor4;
  333.             phiHessian[1][1] = z * (rho2r2 - y2 * r2P2rho2) / rho3r4;
  334.             phiHessian[2][1] = y * (rho2 - z2) / rhor4;
  335.             phiHessian[2][2] = 2 * rho * zOr3 / r;

  337.             // upper-right part is symmetric
  338.             phiHessian[0][1] = phiHessian[1][0];
  339.             phiHessian[0][2] = phiHessian[2][0];
  340.             phiHessian[1][2] = phiHessian[2][1];

  342.         }

  344.     }

  346.     /**
  347.      * Replace the instance with a data transfer object for serialization.
  348.      * @return data transfer object that will be serialized
  349.      */
  350.     private Object writeReplace() {
  351.         return new DataTransferObject(v.getX(), v.getY(), v.getZ());
  352.     }

  354.     /** Internal class used only for serialization. */
  355.     private static class DataTransferObject implements Serializable {

  357.         /** Serializable UID. */
  358.         private static final long serialVersionUID = 20130206L;

  360.         /** Abscissa.
  361.          * @serial
  362.          */
  363.         private final double x;

  365.         /** Ordinate.
  366.          * @serial
  367.          */
  368.         private final double y;

  370.         /** Height.
  371.          * @serial
  372.          */
  373.         private final double z;

  375.         /** Simple constructor.
  376.          * @param x abscissa
  377.          * @param y ordinate
  378.          * @param z height
  379.          */
  380.         public DataTransferObject(final double x, final double y, final double z) {
  381.             this.x = x;
  382.             this.y = y;
  383.             this.z = z;
  384.         }

  386.         /** Replace the deserialized data transfer object with a {@link SphericalCoordinates}.
  387.          * @return replacement {@link SphericalCoordinates}
  388.          */
  389.         private Object readResolve() {
  390.             return new SphericalCoordinates(new Vector3D(x, y, z));
  391.         }

  393.     }

  395. }
Created with JaCoCo 0.7.2.201409121644