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  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
85 };
86 
87 ompl::ProlateHyperspheroid::ProlateHyperspheroid(unsigned int n, const double focus1[], const double focus2[])
88  : dataPtr_(std::make_shared<PhsData>())
89 {
90  // Initialize the data:
91  dataPtr_->dim_ = n;
92  dataPtr_->transverseDiameter_ = 0.0; // Initialize to something.
93  dataPtr_->isTransformUpToDate_ = false;
94 
95  // Copy the arrays into their Eigen containers via the Eigen::Map "view"
96  dataPtr_->xFocus1_ = Eigen::Map<const Eigen::VectorXd>(focus1, dataPtr_->dim_);
97  dataPtr_->xFocus2_ = Eigen::Map<const Eigen::VectorXd>(focus2, dataPtr_->dim_);
98 
99  // Calculate the minimum transverse diameter
100  dataPtr_->minTransverseDiameter_ = (dataPtr_->xFocus1_ - dataPtr_->xFocus2_).norm();
101 
102  // Calculate the centre:
103  dataPtr_->xCentre_ = 0.5 * (dataPtr_->xFocus1_ + dataPtr_->xFocus2_);
104 
105  // Calculate the rotation
106  updateRotation();
107 }
108 
110 {
111  if (transverseDiameter < dataPtr_->minTransverseDiameter_)
112  {
113  OMPL_ERROR("%g < %g", transverseDiameter, dataPtr_->minTransverseDiameter_);
114  throw Exception("Transverse diameter cannot be less than the distance between the foci.");
115  }
116 
117  // Store and update if changed
118  if (dataPtr_->transverseDiameter_ != transverseDiameter)
119  {
120  // Mark as out of date
121  dataPtr_->isTransformUpToDate_ = false;
122 
123  // Store
124  dataPtr_->transverseDiameter_ = transverseDiameter;
125 
126  // Update the transform
127  updateTransformation();
128  }
129  // No else, the diameter didn't change
130 }
131 
132 void ompl::ProlateHyperspheroid::transform(const double sphere[], double phs[]) const
133 {
134  if (!dataPtr_->isTransformUpToDate_)
135  {
136  throw Exception("The transformation is not up to date in the PHS class. Has the transverse diameter been set?");
137  }
138 
139  // Calculate the tranformation and offset, using Eigen::Map views of the data
140  Eigen::Map<Eigen::VectorXd>(phs, dataPtr_->dim_) =
141  dataPtr_->transformationWorldFromEllipse_ * Eigen::Map<const Eigen::VectorXd>(sphere, dataPtr_->dim_);
142  Eigen::Map<Eigen::VectorXd>(phs, dataPtr_->dim_) += dataPtr_->xCentre_;
143 }
144 
145 bool ompl::ProlateHyperspheroid::isInPhs(const double point[]) const
146 {
147  if (!dataPtr_->isTransformUpToDate_)
148  {
149  // The transform is not up to date until the transverse diameter has been set
150  throw Exception("The transverse diameter has not been set");
151  }
152 
153  return (getPathLength(point) < dataPtr_->transverseDiameter_);
154 }
155 
156 bool ompl::ProlateHyperspheroid::isOnPhs(const double point[]) const
157 {
158  if (!dataPtr_->isTransformUpToDate_)
159  {
160  // The transform is not up to date until the transverse diameter has been set
161  throw Exception("The transverse diameter has not been set");
162  }
163 
164  return (getPathLength(point) == dataPtr_->transverseDiameter_);
165 }
166 
168 {
169  return dataPtr_->dim_;
170 }
171 
173 {
174  if (!dataPtr_->isTransformUpToDate_)
175  {
176  // The transform is not up to date until the transverse diameter has been set, therefore we have no transverse
177  // diameter and we have infinite measure
178  return std::numeric_limits<double>::infinity();
179  }
180 
181  // Calculate and return:
182  return dataPtr_->phsMeasure_;
183 }
184 
185 double ompl::ProlateHyperspheroid::getPhsMeasure(double tranDiam) const
186 {
187  return prolateHyperspheroidMeasure(dataPtr_->dim_, dataPtr_->minTransverseDiameter_, tranDiam);
188 }
189 
191 {
192  return dataPtr_->minTransverseDiameter_;
193 }
194 
195 double ompl::ProlateHyperspheroid::getPathLength(const double point[]) const
196 {
197  return (dataPtr_->xFocus1_ - Eigen::Map<const Eigen::VectorXd>(point, dataPtr_->dim_)).norm() +
198  (Eigen::Map<const Eigen::VectorXd>(point, dataPtr_->dim_) - dataPtr_->xFocus2_).norm();
199 }
200 
202 {
203  return dataPtr_->dim_;
204 }
205 
206 void ompl::ProlateHyperspheroid::updateRotation()
207 {
208  // Mark the transform as out of date
209  dataPtr_->isTransformUpToDate_ = false;
210 
211  // If the minTransverseDiameter_ is too close to 0, we treat this as a circle.
212  double circleTol = 1E-9;
213  if (dataPtr_->minTransverseDiameter_ < circleTol)
214  {
215  dataPtr_->rotationWorldFromEllipse_.setIdentity(dataPtr_->dim_, dataPtr_->dim_);
216  }
217  else
218  {
219  // Variables
220  // The transverse axis of the PHS expressed in the world frame.
221  Eigen::VectorXd transverseAxis(dataPtr_->dim_);
222  // The matrix representation of the Wahba problem
223  Eigen::MatrixXd wahbaProb(dataPtr_->dim_, dataPtr_->dim_);
224  // The middle diagonal matrix in the SVD solution to the Wahba problem
225  Eigen::VectorXd middleM(dataPtr_->dim_);
226 
227  // Calculate the major axis, storing as the first eigenvector
228  transverseAxis = (dataPtr_->xFocus2_ - dataPtr_->xFocus1_) / dataPtr_->minTransverseDiameter_;
229 
230  // Calculate the rotation that will allow us to generate the remaining eigenvectors
231  // Formulate as a Wahba problem, first forming the matrix a_j*a_i' where a_j is the transverse axis if the
232  // ellipse in the world frame, and a_i is the first basis vector of the world frame (i.e., [1 0 .... 0])
233  wahbaProb = transverseAxis * Eigen::MatrixXd::Identity(dataPtr_->dim_, dataPtr_->dim_).col(0).transpose();
234 
235  // Then run it through the SVD solver
236  Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::NoQRPreconditioner> svd(wahbaProb,
237  Eigen::ComputeFullV | Eigen::ComputeFullU);
238 
239  // Then calculate the rotation matrix from the U and V components of SVD
240  // Calculate the middle diagonal matrix
241  middleM = Eigen::VectorXd::Ones(dataPtr_->dim_);
242  // Make the last value equal to det(U)*det(V) (zero-based indexing remember)
243  middleM(dataPtr_->dim_ - 1) = svd.matrixU().determinant() * svd.matrixV().determinant();
244 
245  // Calculate the rotation
246  dataPtr_->rotationWorldFromEllipse_ = svd.matrixU() * middleM.asDiagonal() * svd.matrixV().transpose();
247  }
248 }
249 
250 void ompl::ProlateHyperspheroid::updateTransformation()
251 {
252  // Variables
253  // The radii of the ellipse
254  Eigen::VectorXd diagAsVector(dataPtr_->dim_);
255  // The conjugate diameters:
256  double conjugateDiamater;
257 
258  // Calculate the conjugate radius
259  conjugateDiamater = std::sqrt(dataPtr_->transverseDiameter_ * dataPtr_->transverseDiameter_ -
260  dataPtr_->minTransverseDiameter_ * dataPtr_->minTransverseDiameter_);
261 
262  // Store into the diagonal matrix
263  // All the elements but one are the conjugate radius
264  diagAsVector.fill(conjugateDiamater / 2.0);
265 
266  // The first element in diagonal is the transverse radius
267  diagAsVector(0) = 0.5 * dataPtr_->transverseDiameter_;
268 
269  // Calculate the transformation matrix
270  dataPtr_->transformationWorldFromEllipse_ = dataPtr_->rotationWorldFromEllipse_ * diagAsVector.asDiagonal();
271 
272  // Calculate the measure:
273  dataPtr_->phsMeasure_ =
274  prolateHyperspheroidMeasure(dataPtr_->dim_, dataPtr_->minTransverseDiameter_, dataPtr_->transverseDiameter_);
275 
276  // Mark as up to date
277  dataPtr_->isTransformUpToDate_ = true;
278 }
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.
bool isInPhs(const double point[]) const
Check if the given point lies in the PHS.
double getPathLength(const double point[]) const
Calculate length of a line that originates from one focus, passes through the given point,...
void setTransverseDiameter(double transverseDiameter)
Set the transverse diameter of the PHS.
unsigned int getDimension() const
The state dimension of the PHS.
double prolateHyperspheroidMeasure(unsigned int N, double dFoci, double dTransverse)
The Lebesgue measure (i.e., "volume") of an n-dimensional prolate hyperspheroid (a symmetric hyperell...
#define OMPL_ERROR(fmt,...)
Log a formatted error string.
Definition: Console.h:64
double getPhsMeasure() const
The measure of the PHS.
The exception type for ompl.
Definition: Exception.h:78
unsigned int getPhsDimension() const
The dimension of the PHS.
bool isOnPhs(const double point[]) const
Check if the given point lies on the PHS.
double getMinTransverseDiameter() const
The minimum transverse diameter of the PHS, i.e., the distance between the foci.