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_ == false)
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_ == false)
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_ == false)
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_ == false)
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  else
179  {
180  // Calculate and return:
181  return dataPtr_->phsMeasure_;
182  }
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 }
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.