Fréchet View  1.6.0
A Tool for Exploring Fréchet Distance Algorithms
freespace.cpp
Go to the documentation of this file.
1 
2 #include <freespace.h>
3 #include <numeric.h>
4 #include <tbb/tbb.h>
5 
6 using namespace frechet;
7 using namespace data;
8 using namespace numeric;
9 using namespace fs;
10 
11 const QRectF FreeSpace::UNIT_RECT(0,0,1,1);
12 const double FreeSpace::DET_MINIMUM = 1e-12;
13 
14 FreeSpace::FreeSpace(const Curve &ap, const Curve &aq)
15  : P(ap), Q(aq),
16  n(P.length()), m(Q.length()),
17  cells(n,m),
18  _components(n,m)
19 { }
20 
21 
23  : P(that.P), Q(that.Q),
24  n(that.P.length()), m(that.Q.length()),
25  cells(that.cells),
26  _components(that._components)
27 { }
28 
35 void FreeSpace::calculateFreeSpace(double epsilon)
36 {
37  // inner cells
38  tbb::static_partitioner sp;
39  tbb::parallel_for(0, n-1, [&] (int i)
40  {
41  const Point& p1 = P[i];
42  const Point& p2 = P[i+1];
43 
44  for(int j=0; j < (m-1); ++j)
45  {
46  // Calculate L,B
47  const Point& q1 = Q[j];
48  const Point& q2 = Q[j+1];
49 
50  Cell& cl = cell(i,j);
51  calculateFreeSpaceSegment<accurate_float>(
52  p1, q1,q2,
53  epsilon,
54  cl.L);
55 
56  calculateFreeSpaceSegment<accurate_float>(
57  q1, p1,p2,
58  epsilon,
59  cl.B);
60  }
61  }, sp);
62 
63  // rightmost column
64  if (n > 0)
65  {
66  int i = n-1;
67  const Point& p1 = P[i];
68  for(int j=0; j < (m-1); ++j)
69  {
70  // Calculate L
71  Cell& cl = cell(i,j);
72  calculateFreeSpaceSegment<accurate_float>(
73  p1, Q[j],Q[j+1],
74  epsilon,
75  cl.L);
76  }
77  }
78 
79  // topmost row
80  if (m > 0)
81  {
82  int j = m-1;
83  const Point& q1 = Q[j];
84  for(int i=0; i < (n-1); ++i)
85  {
86  // Calculate B
87  Cell& cl = cell(i,j);
88  calculateFreeSpaceSegment<accurate_float>(
89  q1, P[i],P[i+1],
90  epsilon,
91  cl.B);
92  }
93  }
94 
95  //
96  // Fix border cases
97  //
98 
99  tbb::parallel_for (0, n - 1, [&] (int i)
100  {
101  for (int j = 0; j <= (m - 1); ++j) {
102  Cell& cl = cell(i, j);
104  (j < (m-1)) ? &cl.L : nullptr,
105  (i < (n-1)) ? &cl.B : nullptr,
106  (j>0) ? &cell(i, j - 1).L : nullptr,
107  (i>0) ? &cell(i - 1, j).B : nullptr);
108  }
109  }, sp);
110 }
111 
112 void FreeSpace::calculateBounds(double epsilon)
113 {
114  int i,j;
115  for(i=0; i < (n-1); ++i)
116  {
117  for(j=0; j < (m-1); ++j)
118  {
119  Cell& cl = cell(i,j);
120  const Interval& L = cl.L;
121  const Interval& B = cell(i,j).B;
122  const Interval& R = cell(i+1,j).L; // neighbor to the right
123  const Interval& T = cell(i,j+1).B; // neighbor at top
124 
125  Interval& H = cl.bounds.H;
126  Interval& V = cl.bounds.V;
127 
128  if (L)
129  H.setLower(0.0); // easy: left boundary = 0.0
130  else
131  H.setLower(min(B.lower(),T.lower())); // check minimum of intersections (possibly undefined)
132 
133  if (R)
134  H.setUpper(1.0); // easy: right boundary = 1.0
135  else
136  H.setUpper(max(B.upper(),T.upper())); // check maximum of intersections (possibly undefined)
137 
138  if (B)
139  V.setLower(0.0); // easy: bottom boundary = 0.0
140  else
141  V.setLower(min(L.lower(),R.lower()));
142 
143  if (T)
144  V.setUpper(1.0);
145  else
146  V.setUpper(max(L.upper(),R.upper()));
147 
148  // do we have to calculate the ellipse boundaries ??
149  if (H != Interval::UNIT)
150  {
151  double lambda[2];
152  calculateEllipseBoundaries<accurate_float>(
153  P[i],P[i+1], Q[j],Q[j+1],
154  epsilon, lambda);
155 
156  H.setLower(min(H.lower(),lambda[0]));
157  H.setUpper(max(H.upper(),lambda[1]));
158  }
159 
160  if (V != Interval::UNIT)
161  {
162  double lambda[2];
163  calculateEllipseBoundaries<accurate_float>(
164  Q[j],Q[j+1], P[i],P[i+1],
165  epsilon, lambda);
166 
167  V.setLower(min(V.lower(),lambda[0]));
168  V.setUpper(max(V.upper(),lambda[1]));
169  }
170  }
171  }
172 }
173 
174 bool FreeSpace::cellEmptyBounds(int i, int j) const {
175  const Cell& cl = cell(i,j);
176  // easy to calculate but assumes that bounds have been calculate !!
177  // (only k-Frechet does it!)
178  return !cl.bounds.H.valid()
179  && !cl.bounds.V.valid();
180 }
181 
182 bool FreeSpace::cellEmptyIntervals(int i, int j) const {
183  const Cell& cl = cell(i,j);
184  return ! cell(i,j).L
185  && ! cell(i,j).B
186  && ! cell(i+1,j).L
187  && ! cell(i,j+1).B;
188 }
189 
190 bool FreeSpace::cellFull(int i, int j) const {
191  // is this cell completely free?
192  const Cell& cl = cell(i,j);
193  return cl.L.contains(Interval::UNIT)
194  && cl.B.contains(Interval::UNIT)
195  && cell(i+1,j).L.contains(Interval::UNIT)
196  && cell(i,j+1).B.contains(Interval::UNIT);
197 }
198 
199 
200 double FreeSpace::determinant(int i, int j) const {
201  Point s1 = P[i];
202  Point s2 = P[i+1];
203 
204  Point t1 = Q[j];
205  Point t2 = Q[j+1];
206 
207  return (s1.x()-s2.x())*(t1.y()-t2.y()) - (s1.y()-s2.y())*(t1.x()-t2.x());
208 }
209 
210 
211 QRectF FreeSpace::segmentBounds(int i, int j)
212 {
213  const Cell& cl = cell(i,j);
214  return cl.bounds.boundingRect();
215 }
216 
217 std::list<Interval> FreeSpace::collectFreeVerticalIntervals(int i) const
218 {
219  std::list<Interval> result;
220  const Interval& L = cell(i,0).L;
221  bool inside = (L.lower()==0.0);
222  if (inside) result.push_back(L);
223 
224  for(int j=0; j < m-1; ++j)
225  {
226  const Interval& L = cell(i,j).L;
227  const Interval& L2 = cell(i,j+1).L;
228 
229  if (!L) {
230  inside=false;
231  }
232  else if (L.lower() > 0.0) {
233  result.push_back(L + (double)j);
234  inside = (L.upper() >= 1.0) && (L2.lower()==0.0);
235  }
236  else if (inside) {
237  Q_ASSERT(L.lower()==0.0);
238  result.back().upper() = L.upper() + j;
239  inside = (L.upper() >= 1.0) && (L2.lower()==0.0);
240  }
241  }
242  return result;
243 }
244 
253 {
254 
255  bool fix = (N && (N->lower()==0.0))
256  || (E && (E->lower()==0.0))
257  || (S && (S->upper()==1.0))
258  || (W && (W->upper()==1.0));
259 
260  if (fix) {
261  if (N) {
262  N->lower()=0.0;
263  if (std::isnan(N->upper())) N->upper()=0.0;
264  }
265  if (E) {
266  E->lower()=0.0;
267  if (std::isnan(E->upper())) E->upper()=0.0;
268  }
269  if (S) {
270  S->upper()=1.0;
271  if (std::isnan(S->lower())) S->lower()=1.0;
272  }
273  if (W && (std::isnan(W->upper()) || (W->upper()!=1.0))) {
274  W->upper()=1.0;
275  if (std::isnan(W->lower())) W->lower()=1.0;
276  }
277  }
278 }
279 
280 template<typename Float>
282  const Point& p1, const Point& p2,
283  const Point& q1, const Point& q2,
284  double epsilon,
285  double lambda[2])
286 {
287  Float A,B,D;
288 
289  Float p1x=p1.x();
290  Float p1y=p1.y();
291  Float p2x=p2.x();
292  Float p2y=p2.y();
293  Float q1x=q1.x();
294  Float q1y=q1.y();
295  Float q2x=q2.x();
296  Float q2y=q2.y();
297 
298  lambda[0]=lambda[1] = NAN;
299 
300  A = q2x*(p1y-q1y) + q1x*(q2y-p1y) + p1x*(q1y-q2y);
301  B = euclidean_distance_sq<Float>(q1,q2);
302  D = (q2y-q1y)*(p2x-p1x) - (p2y-p1y)*(q2x-q1x);
303 
304  if (D==0.0) return;
305 
306  // lambda_1/2
307  Float lambda0 = (A - epsilon*sqrt(B)) / D;
308  Float lambda1 = (A + epsilon*sqrt(B)) / D;
309 
310  if (lambda0 > lambda1)
311  std::swap(lambda0,lambda1);
312 
313  // clip to unit rect
314  if ((lambda0<0.0) || (lambda0>1.0))
315  lambda0 = NAN;
316  if ((lambda1<0.0) || (lambda1>1.0))
317  lambda1 = NAN;
318 
319  if (std::isnan(lambda0) && std::isnan(lambda1)) return;
320 
321  // mu_1/2
322  Float E8 = (p2x-p1x)*(q2x-q1x) + (p2y-p1y)*(q2y-q1y);
323  Float E9 = (q2x-q1x)*(p1x-q1x) + (q2y-q1y)*(p1y-q1y);
324 
325  Float mu0 = (lambda0*E8+E9)/B;
326  Float mu1 = (lambda1*E8+E9)/B;
327 
328  // clip to unit rect
329  if ((mu0>=0.0) && (mu0<=1.0))
330  lambda[0] = (double)lambda0;
331 
332  if ((mu1>=0.0) && (mu1<=1.0))
333  lambda[1] = (double)lambda1;
334 }
335 
void swap(gpuword **A, gpuword **B)
static const QRectF UNIT_RECT
the unit rectangle [0,1]x[0,1]
Definition: freespace.h:95
bool cellFull(int i, int j) const
test if a cell is completely covered by the free-space area
Definition: freespace.cpp:190
static void fixBorderCase(data::Interval *N, data::Interval *E, data::Interval *S, data::Interval *W)
fix border case at the joint of four intervals
Definition: freespace.cpp:252
QRectF segmentBounds(int i, int j)
Definition: freespace.cpp:211
Interval & setLower(double value)
update lower bound
Definition: interval.h:111
Interval H
horizontal interval
Definition: interval.h:459
global definitions for all algorithms.
double determinant(int i, int j) const
Definition: freespace.cpp:200
void calculateFreeSpace(double epsilon)
calculate all free space intervals
Definition: freespace.cpp:35
static const double DET_MINIMUM
minimum determinant to avoid roundoff errors
Definition: freespace.h:97
data::Interval B
Definition: freespace.h:31
Float sqrt(const Float &x)
square-root function template for floating point types
data::IntervalPair bounds
Definition: freespace.h:40
bool cellEmptyIntervals(int i, int j) const
test if a cell is empty
Definition: freespace.cpp:182
double upper() const
Definition: interval.h:96
data::Interval L
Definition: freespace.h:31
Interval V
vertical interval
Definition: interval.h:460
QPointF Point
a point in the plane; with double floating point precision. This type is heavily used throughout all ...
Definition: types.h:14
QPolygonF Curve
a polygonal curve in the plane; with double floating point precision. This type is heavily used throu...
Definition: types.h:20
double lower() const
Definition: interval.h:92
std::list< data::Interval > collectFreeVerticalIntervals(int i) const
traverse one column of the free-space and find all connected intervals. This method is needed by the ...
Definition: freespace.cpp:217
double min(double a, double b)
minimum function with checks for NAN
Definition: numeric.h:222
Interval & setUpper(double value)
update upper bound
Definition: interval.h:117
const int m
number of vertexes in Q
Definition: freespace.h:89
bool cellEmptyBounds(int i, int j) const
test if the bounding box of a cell is empty
Definition: freespace.cpp:174
static const Interval UNIT
the unit interval [0,1]
Definition: interval.h:37
void calculateBounds(double epsilon)
calculate the bounding boxes of the free-space areas (k-Frechet only)
Definition: freespace.cpp:112
The FreeSpace class models the Free-Space Diagram calculated from two curves.
Definition: freespace.h:83
holds data for one cell of the free-space diagram
Definition: freespace.h:27
bool contains(double x) const
containment test (assumes closed interval, bounds inclusive)
Definition: interval.h:206
const int n
number of vertexes in P
Definition: freespace.h:88
an interval of two double values.
Definition: interval.h:31
FreeSpace(const Curve &ap, const Curve &aq)
default constructor with two input curves
Definition: freespace.cpp:14
bool valid() const
validity test
Definition: interval.h:74
const Curve & Q
second input curve
Definition: freespace.h:87
void calculateEllipseBoundaries(const Point &p1, const Point &p2, const Point &q1, const Point &q2, double epsilon, double lambda[2])
compute the bounding box of one cell in the free-space are. As parameters we expect two segments of P...
Cell & cell(int x, int y)
Definition: freespace.h:160
const Curve & P
first input curve
Definition: freespace.h:86
double max(double a, double b)
maximum function with checks for NAN
Definition: numeric.h:233
QRectF boundingRect() const
calculate the bounding rectangle, i.e. a rectangle that is defined by the horizontal and vertical int...
Definition: interval.h:480