001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math3.geometry.euclidean.threed; 018 019 020import java.io.Serializable; 021 022import org.apache.commons.math3.util.FastMath; 023 024/** This class provides conversions related to <a 025 * href="http://mathworld.wolfram.com/SphericalCoordinates.html">spherical coordinates</a>. 026 * <p> 027 * The conventions used here are the mathematical ones, i.e. spherical coordinates are 028 * related to Cartesian coordinates as follows: 029 * </p> 030 * <ul> 031 * <li>x = r cos(θ) sin(Φ)</li> 032 * <li>y = r sin(θ) sin(Φ)</li> 033 * <li>z = r cos(Φ)</li> 034 * </ul> 035 * <ul> 036 * <li>r = √(x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>)</li> 037 * <li>θ = atan2(y, x)</li> 038 * <li>Φ = acos(z/r)</li> 039 * </ul> 040 * <p> 041 * r is the radius, θ is the azimuthal angle in the x-y plane and Φ is the polar 042 * (co-latitude) angle. These conventions are <em>different</em> from the conventions used 043 * in physics (and in particular in spherical harmonics) where the meanings of θ and 044 * Φ are reversed. 045 * </p> 046 * <p> 047 * This class provides conversion of coordinates and also of gradient and Hessian 048 * between spherical and Cartesian coordinates. 049 * </p> 050 * @since 3.2 051 */ 052public class SphericalCoordinates implements Serializable { 053 054 /** Serializable UID. */ 055 private static final long serialVersionUID = 20130206L; 056 057 /** Cartesian coordinates. */ 058 private final Vector3D v; 059 060 /** Radius. */ 061 private final double r; 062 063 /** Azimuthal angle in the x-y plane θ. */ 064 private final double theta; 065 066 /** Polar angle (co-latitude) Φ. */ 067 private final double phi; 068 069 /** Jacobian of (r, θ &Phi). */ 070 private double[][] jacobian; 071 072 /** Hessian of radius. */ 073 private double[][] rHessian; 074 075 /** Hessian of azimuthal angle in the x-y plane θ. */ 076 private double[][] thetaHessian; 077 078 /** Hessian of polar (co-latitude) angle Φ. */ 079 private double[][] phiHessian; 080 081 /** Build a spherical coordinates transformer from Cartesian coordinates. 082 * @param v Cartesian coordinates 083 */ 084 public SphericalCoordinates(final Vector3D v) { 085 086 // Cartesian coordinates 087 this.v = v; 088 089 // remaining spherical coordinates 090 this.r = v.getNorm(); 091 this.theta = v.getAlpha(); 092 this.phi = FastMath.acos(v.getZ() / r); 093 094 } 095 096 /** Build a spherical coordinates transformer from spherical coordinates. 097 * @param r radius 098 * @param theta azimuthal angle in x-y plane 099 * @param phi polar (co-latitude) angle 100 */ 101 public SphericalCoordinates(final double r, final double theta, final double phi) { 102 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); 107 108 // spherical coordinates 109 this.r = r; 110 this.theta = theta; 111 this.phi = phi; 112 113 // Cartesian coordinates 114 this.v = new Vector3D(r * cosTheta * sinPhi, 115 r * sinTheta * sinPhi, 116 r * cosPhi); 117 118 } 119 120 /** Get the Cartesian coordinates. 121 * @return Cartesian coordinates 122 */ 123 public Vector3D getCartesian() { 124 return v; 125 } 126 127 /** Get the radius. 128 * @return radius r 129 * @see #getTheta() 130 * @see #getPhi() 131 */ 132 public double getR() { 133 return r; 134 } 135 136 /** Get the azimuthal angle in x-y plane. 137 * @return azimuthal angle in x-y plane θ 138 * @see #getR() 139 * @see #getPhi() 140 */ 141 public double getTheta() { 142 return theta; 143 } 144 145 /** Get the polar (co-latitude) angle. 146 * @return polar (co-latitude) angle Φ 147 * @see #getR() 148 * @see #getTheta() 149 */ 150 public double getPhi() { 151 return phi; 152 } 153 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θ, df/dΦ} 158 * @return gradient with respect to Cartesian coordinates 159 * {df/dx, df/dy, df/dz} 160 */ 161 public double[] toCartesianGradient(final double[] sGradient) { 162 163 // lazy evaluation of Jacobian 164 computeJacobian(); 165 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 }; 173 174 } 175 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θ, d<sup>2</sup>f/drdΦ}, 185 * {d<sup>2</sup>f/drdθ, d<sup>2</sup>f/dθ<sup>2</sup>, d<sup>2</sup>f/dθdΦ}, 186 * {d<sup>2</sup>f/drdΦ, d<sup>2</sup>f/dθdΦ, d<sup>2</sup>f/dΦ<sup>2</sup>} 187 * @param sGradient gradient with respect to spherical coordinates 188 * {df/dr, df/dθ, df/dΦ} 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) { 195 196 computeJacobian(); 197 computeHessians(); 198 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]; 204 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]; 216 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]; 224 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]; 232 233 // ensure symmetry 234 cHessian[0][1] = cHessian[1][0]; 235 cHessian[0][2] = cHessian[2][0]; 236 cHessian[1][2] = cHessian[2][1]; 237 238 return cHessian; 239 240 } 241 242 /** Lazy evaluation of (r, θ, φ) Jacobian. 243 */ 244 private void computeJacobian() { 245 if (jacobian == null) { 246 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; 254 255 jacobian = new double[3][3]; 256 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; 261 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 266 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; 271 272 } 273 } 274 275 /** Lazy evaluation of Hessians. 276 */ 277 private void computeHessians() { 278 279 if (rHessian == null) { 280 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; 299 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; 308 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]; 313 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; 319 320 // upper-right part is symmetric 321 thetaHessian[0][1] = thetaHessian[1][0]; 322 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; 336 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]; 341 342 } 343 344 } 345 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 } 353 354 /** Internal class used only for serialization. */ 355 private static class DataTransferObject implements Serializable { 356 357 /** Serializable UID. */ 358 private static final long serialVersionUID = 20130206L; 359 360 /** Abscissa. 361 * @serial 362 */ 363 private final double x; 364 365 /** Ordinate. 366 * @serial 367 */ 368 private final double y; 369 370 /** Height. 371 * @serial 372 */ 373 private final double z; 374 375 /** Simple constructor. 376 * @param x abscissa 377 * @param y ordinate 378 * @param z height 379 */ 380 DataTransferObject(final double x, final double y, final double z) { 381 this.x = x; 382 this.y = y; 383 this.z = z; 384 } 385 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 } 392 393 } 394 395}