SemiDiscreteOT 1.0
Semi-Discrete Optimal Transport Library
Loading...
Searching...
No Matches
Distance.h
Go to the documentation of this file.
1#ifndef DISTANCE_H
2#define DISTANCE_H
3
4#include <deal.II/base/point.h>
5
11#include <deal.II/base/point.h>
12
13using namespace dealii;
14
15// TODO implement factory
16
23template <int spacedim>
25 const Point<spacedim> a, const Point<spacedim> b) {
26 return (a-b).norm();
27}
28
35template <int spacedim>
37 const Point<spacedim> a, const Point<spacedim> b) {
38 double dot_product = 0.0;
39 for (unsigned int i = 0; i < spacedim; ++i) {
40 dot_product += a[i] * b[i];
41 }
42
43 // Normalize by the product of magnitudes
44 double norm_a = a.norm();
45 double norm_b = b.norm();
46
47 // Ensure dot product is between -1 and 1 to avoid numerical issues
48 double cosine = std::min(1.0, std::max(-1.0, dot_product / (norm_a * norm_b)));
49
50 return std::acos(cosine);
51}
52
53// distance gradients: grad of d^2
54
61template <int spacedim>
63 const Point<spacedim> a, const Point<spacedim> b) {
64 Vector<double> gradient(spacedim);
65 for (unsigned int i = 0; i < spacedim; ++i) {
66 gradient[i] = (a[i] - b[i]);
67 }
68 return gradient;
69}
70
71// it is actually the log map and grad of d^2, mind the order, x is the evaluation point
78template <int spacedim>
80 const Point<spacedim> a, const Point<spacedim> b) {
81 Vector<double> gradient(spacedim);
82 double dist = spherical_distance(a, b);
83
84 double dot_product = 0.0;
85 for (unsigned int i = 0; i < spacedim; ++i) {
86 dot_product += a[i] * b[i];
87 }
88
89 for (unsigned int i = 0; i < spacedim; ++i) {
90 gradient[i] = -2 * (dist/std::sin(dist)) * (b[i]-dot_product*a[i]);
91 }
92 return gradient;
93}
94
95// exponential maps
102template <int spacedim>
104 const Point<spacedim> a, const Vector<double> v) {
105 Point<spacedim> b;
106 for (unsigned int i = 0; i < spacedim; ++i) {
107 b[i] = a[i] + v[i];
108 }
109 return b;
110}
111
112// v must be different form 0, a perpendicular to v
119template <int spacedim>
121 const Point<spacedim> a, const Vector<double> v) {
122 Point<spacedim> b;
123 for (unsigned int i = 0; i < spacedim; ++i) {
124 b[i] = std::cos(v.l2_norm())*a[i] + std::sin(v.l2_norm()) * (v[i]/v.l2_norm());
125 }
126 return b;
127}
128
129#endif // DISTANCE_H
Point< spacedim > spherical_distance_exp_map(const Point< spacedim > a, const Vector< double > v)
Computes the exponential map for the spherical distance.
Definition Distance.h:120
Vector< double > euclidean_distance_gradient(const Point< spacedim > a, const Point< spacedim > b)
Computes the gradient of the squared Euclidean distance.
Definition Distance.h:62
double euclidean_distance(const Point< spacedim > a, const Point< spacedim > b)
Computes the Euclidean distance between two points.
Definition Distance.h:24
double spherical_distance(const Point< spacedim > a, const Point< spacedim > b)
Computes the spherical distance between two points.
Definition Distance.h:36
Point< spacedim > euclidean_distance_exp_map(const Point< spacedim > a, const Vector< double > v)
Computes the exponential map for the Euclidean distance.
Definition Distance.h:103
Vector< double > spherical_distance_gradient(const Point< spacedim > a, const Point< spacedim > b)
Computes the gradient of the squared spherical distance.
Definition Distance.h:79