Loading...
Searching...
No Matches
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
56struct 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
87ompl::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
132void 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
145bool 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
156bool 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
185double ompl::ProlateHyperspheroid::getPhsMeasure(double tranDiam) const
186{
187 return prolateHyperspheroidMeasure(dataPtr_->dim_, dataPtr_->minTransverseDiameter_, tranDiam);
188}
189
191{
192 return dataPtr_->minTransverseDiameter_;
193}
194
195double 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
206void 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
250void 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}
The exception type for ompl.
Definition Exception.h:47
bool isInPhs(const double point[]) const
Check if the given point lies in 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.
void setTransverseDiameter(double transverseDiameter)
Set the transverse diameter of the PHS.
unsigned int getPhsDimension() const
The dimension of the PHS.
double getPathLength(const double point[]) const
Calculate length of a line that originates from one focus, passes through the given point,...
double getPhsMeasure() const
The measure of the PHS.
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.
ProlateHyperspheroid(unsigned int n, const double focus1[], const double focus2[])
The description of an n-dimensional prolate hyperspheroid.
unsigned int getDimension() const
The state dimension of the PHS.
#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...
STL namespace.