ProlateHyperspheroid.cpp
1 /*********************************************************************
2 * Software License Agreement (BSD License)
3 *
4 * Copyright (c) 2014, University of Toronto
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 *
11 * * Redistributions of source code must retain the above copyright
12 * notice, this list of conditions and the following disclaimer.
13 * * Redistributions in binary form must reproduce the above
14 * copyright notice, this list of conditions and the following
15 * disclaimer in the documentation and/or other materials provided
16 * with the distribution.
17 * * Neither the name of the University of Toronto nor the names of its
18 * contributors may be used to endorse or promote products derived
19 * from this software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
24 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
25 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
27 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
28 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
29 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
31 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
32 * POSSIBILITY OF SUCH DAMAGE.
33 *********************************************************************/
34 
35 /* Author: Jonathan Gammell*/
36 
37 // The class's header
38 #include "ompl/util/ProlateHyperspheroid.h"
39 // For OMPL exceptions
40 #include "ompl/util/Exception.h"
41 // For OMPL information
42 #include "ompl/util/Console.h"
43 // For geometric equations like prolateHyperspheroidMeasure
44 #include "ompl/util/GeometricEquations.h"
45 
46 // For std::make_shared
47 #include <memory>
48 
49 // Eigen core:
50 #include <Eigen/Core>
51 // Inversion and determinants
52 #include <Eigen/LU>
53 // SVD decomposition
54 #include <Eigen/SVD>
55 
56 struct ompl::ProlateHyperspheroid::PhsData
57 {
59  unsigned int dim_;
61  bool isTransformUpToDate_;
63  double minTransverseDiameter_;
65  double transverseDiameter_;
67  double phsMeasure_;
70  Eigen::VectorXd xFocus1_;
73  Eigen::VectorXd xFocus2_;
76  Eigen::VectorXd xCentre_;
79  Eigen::MatrixXd rotationWorldFromEllipse_;
82  Eigen::MatrixXd transformationWorldFromEllipse_;
83 };
84 
85 ompl::ProlateHyperspheroid::ProlateHyperspheroid(unsigned int n, const double focus1[], const double focus2[])
86  : dataPtr_(std::make_shared<PhsData>())
87 {
88  // Initialize the data:
89  dataPtr_->dim_ = n;
90  dataPtr_->transverseDiameter_ = 0.0; // Initialize to something.
91  dataPtr_->isTransformUpToDate_ = false;
92 
93  // Copy the arrays into their Eigen containers via the Eigen::Map "view"
94  dataPtr_->xFocus1_ = Eigen::Map<const Eigen::VectorXd>(focus1, dataPtr_->dim_);
95  dataPtr_->xFocus2_ = Eigen::Map<const Eigen::VectorXd>(focus2, dataPtr_->dim_);
96 
97  // Calculate the minimum transverse diameter
98  dataPtr_->minTransverseDiameter_ = (dataPtr_->xFocus1_ - dataPtr_->xFocus2_).norm();
99 
100  // Calculate the centre:
101  dataPtr_->xCentre_ = 0.5 * (dataPtr_->xFocus1_ + dataPtr_->xFocus2_);
102 
103  // Calculate the rotation
104  updateRotation();
105 }
106 
108 {
109  if (transverseDiameter < dataPtr_->minTransverseDiameter_)
110  {
111  OMPL_ERROR("%g < %g", transverseDiameter, dataPtr_->minTransverseDiameter_);
112  throw Exception("Transverse diameter cannot be less than the distance between the foci.");
113  }
114 
115  // Store and update if changed
116  if (dataPtr_->transverseDiameter_ != transverseDiameter)
117  {
118  // Mark as out of date
119  dataPtr_->isTransformUpToDate_ = false;
120 
121  // Store
122  dataPtr_->transverseDiameter_ = transverseDiameter;
123 
124  // Update the transform
125  updateTransformation();
126  }
127  // No else, the diameter didn't change
128 }
129 
130 void ompl::ProlateHyperspheroid::transform(const double sphere[], double phs[]) const
131 {
132  if (!dataPtr_->isTransformUpToDate_)
133  {
134  throw Exception("The transformation is not up to date in the PHS class. Has the transverse diameter been set?");
135  }
136 
137  // Calculate the tranformation and offset, using Eigen::Map views of the data
138  Eigen::Map<Eigen::VectorXd>(phs, dataPtr_->dim_) =
139  dataPtr_->transformationWorldFromEllipse_ * Eigen::Map<const Eigen::VectorXd>(sphere, dataPtr_->dim_);
140  Eigen::Map<Eigen::VectorXd>(phs, dataPtr_->dim_) += dataPtr_->xCentre_;
141 }
142 
143 bool ompl::ProlateHyperspheroid::isInPhs(const double point[]) const
144 {
145  if (!dataPtr_->isTransformUpToDate_)
146  {
147  // The transform is not up to date until the transverse diameter has been set
148  throw Exception("The transverse diameter has not been set");
149  }
150 
151  return (getPathLength(point) < dataPtr_->transverseDiameter_);
152 }
153 
154 bool ompl::ProlateHyperspheroid::isOnPhs(const double point[]) const
155 {
156  if (!dataPtr_->isTransformUpToDate_)
157  {
158  // The transform is not up to date until the transverse diameter has been set
159  throw Exception("The transverse diameter has not been set");
160  }
161 
162  return (getPathLength(point) == dataPtr_->transverseDiameter_);
163 }
164 
166 {
167  return dataPtr_->dim_;
168 }
169 
171 {
172  if (!dataPtr_->isTransformUpToDate_)
173  {
174  // The transform is not up to date until the transverse diameter has been set, therefore we have no transverse
175  // diameter and we have infinite measure
176  return std::numeric_limits<double>::infinity();
177  }
178 
179  // Calculate and return:
180  return dataPtr_->phsMeasure_;
181 }
182 
183 double ompl::ProlateHyperspheroid::getPhsMeasure(double tranDiam) const
184 {
185  return prolateHyperspheroidMeasure(dataPtr_->dim_, dataPtr_->minTransverseDiameter_, tranDiam);
186 }
187 
189 {
190  return dataPtr_->minTransverseDiameter_;
191 }
192 
193 double ompl::ProlateHyperspheroid::getPathLength(const double point[]) const
194 {
195  return (dataPtr_->xFocus1_ - Eigen::Map<const Eigen::VectorXd>(point, dataPtr_->dim_)).norm() +
196  (Eigen::Map<const Eigen::VectorXd>(point, dataPtr_->dim_) - dataPtr_->xFocus2_).norm();
197 }
198 
200 {
201  return dataPtr_->dim_;
202 }
203 
204 void ompl::ProlateHyperspheroid::updateRotation()
205 {
206  // Mark the transform as out of date
207  dataPtr_->isTransformUpToDate_ = false;
208 
209  // If the minTransverseDiameter_ is too close to 0, we treat this as a circle.
210  double circleTol = 1E-9;
211  if (dataPtr_->minTransverseDiameter_ < circleTol)
212  {
213  dataPtr_->rotationWorldFromEllipse_.setIdentity(dataPtr_->dim_, dataPtr_->dim_);
214  }
215  else
216  {
217  // Variables
218  // The transverse axis of the PHS expressed in the world frame.
219  Eigen::VectorXd transverseAxis(dataPtr_->dim_);
220  // The matrix representation of the Wahba problem
221  Eigen::MatrixXd wahbaProb(dataPtr_->dim_, dataPtr_->dim_);
222  // The middle diagonal matrix in the SVD solution to the Wahba problem
223  Eigen::VectorXd middleM(dataPtr_->dim_);
224 
225  // Calculate the major axis, storing as the first eigenvector
226  transverseAxis = (dataPtr_->xFocus2_ - dataPtr_->xFocus1_) / dataPtr_->minTransverseDiameter_;
227 
228  // Calculate the rotation that will allow us to generate the remaining eigenvectors
229  // Formulate as a Wahba problem, first forming the matrix a_j*a_i' where a_j is the transverse axis if the
230  // ellipse in the world frame, and a_i is the first basis vector of the world frame (i.e., [1 0 .... 0])
231  wahbaProb = transverseAxis * Eigen::MatrixXd::Identity(dataPtr_->dim_, dataPtr_->dim_).col(0).transpose();
232 
233  // Then run it through the SVD solver
234  Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::NoQRPreconditioner> svd(wahbaProb,
235  Eigen::ComputeFullV | Eigen::ComputeFullU);
236 
237  // Then calculate the rotation matrix from the U and V components of SVD
238  // Calculate the middle diagonal matrix
239  middleM = Eigen::VectorXd::Ones(dataPtr_->dim_);
240  // Make the last value equal to det(U)*det(V) (zero-based indexing remember)
241  middleM(dataPtr_->dim_ - 1) = svd.matrixU().determinant() * svd.matrixV().determinant();
242 
243  // Calculate the rotation
244  dataPtr_->rotationWorldFromEllipse_ = svd.matrixU() * middleM.asDiagonal() * svd.matrixV().transpose();
245  }
246 }
247 
248 void ompl::ProlateHyperspheroid::updateTransformation()
249 {
250  // Variables
251  // The radii of the ellipse
252  Eigen::VectorXd diagAsVector(dataPtr_->dim_);
253  // The conjugate diameters:
254  double conjugateDiamater;
255 
256  // Calculate the conjugate radius
257  conjugateDiamater = std::sqrt(dataPtr_->transverseDiameter_ * dataPtr_->transverseDiameter_ -
258  dataPtr_->minTransverseDiameter_ * dataPtr_->minTransverseDiameter_);
259 
260  // Store into the diagonal matrix
261  // All the elements but one are the conjugate radius
262  diagAsVector.fill(conjugateDiamater / 2.0);
263 
264  // The first element in diagonal is the transverse radius
265  diagAsVector(0) = 0.5 * dataPtr_->transverseDiameter_;
266 
267  // Calculate the transformation matrix
268  dataPtr_->transformationWorldFromEllipse_ = dataPtr_->rotationWorldFromEllipse_ * diagAsVector.asDiagonal();
269 
270  // Calculate the measure:
271  dataPtr_->phsMeasure_ =
272  prolateHyperspheroidMeasure(dataPtr_->dim_, dataPtr_->minTransverseDiameter_, dataPtr_->transverseDiameter_);
273 
274  // Mark as up to date
275  dataPtr_->isTransformUpToDate_ = true;
276 }
bool isOnPhs(const double point[]) const
Check if the given point lies on the PHS.
double getPhsMeasure() const
The measure of the PHS.
double getPathLength(const double point[]) const
Calculate length of a line that originates from one focus, passes through the given point...
unsigned int getPhsDimension() const
The dimension of the PHS.
STL namespace.
void setTransverseDiameter(double transverseDiameter)
Set the transverse diameter of the PHS.
double getMinTransverseDiameter() const
The minimum transverse diameter of the PHS, i.e., the distance between the foci.
#define OMPL_ERROR(fmt,...)
Log a formatted error string.
Definition: Console.h:64
double prolateHyperspheroidMeasure(unsigned int N, double dFoci, double dTransverse)
The Lebesgue measure (i.e., "volume") of an n-dimensional prolate hyperspheroid (a symmetric hyperell...
The exception type for ompl.
Definition: Exception.h:46
bool isInPhs(const double point[]) const
Check if the given point lies in the PHS.
ProlateHyperspheroid(unsigned int n, const double focus1[], const double focus2[])
The description of an n-dimensional prolate hyperspheroid.
void transform(const double sphere[], double phs[]) const
Transform a point from a sphere to PHS. The return variable phs is expected to already exist...
unsigned int getDimension() const
The state dimension of the PHS.