← Back to branch summary

~amcg-stokes/fluidity/block-velocity-nns

« back to all changes in this revision

Viewing changes to spatialindex-1.8.0/test/tprtree/Exhaustive.cc

  • Committer: Rhodri Davies
  • Date: 2014-05-29 00:14:40 UTC
  • mfrom: (3854.1.473 fluidity)
  • Revision ID: rhodri.davies@imperial.ac.uk-20140529001440-0lx74o0djriitwbn
Merge in trunk.

Show diffs side-by-side

added added

removed removed

Lines of Context:
spatialindex-1.8.0/test/tprtree/Exhaustive.cc
 
1
/******************************************************************************
 
2
 * Project:  libspatialindex - A C++ library for spatial indexing
 
3
 * Author:   Marios Hadjieleftheriou, mhadji@gmail.com
 
4
 ******************************************************************************
 
5
 * Copyright (c) 2002, Marios Hadjieleftheriou
 
6
 *
 
7
 * All rights reserved.
 
8
 * 
 
9
 * Permission is hereby granted, free of charge, to any person obtaining a
 
10
 * copy of this software and associated documentation files (the "Software"),
 
11
 * to deal in the Software without restriction, including without limitation
 
12
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 
13
 * and/or sell copies of the Software, and to permit persons to whom the
 
14
 * Software is furnished to do so, subject to the following conditions:
 
15
 *
 
16
 * The above copyright notice and this permission notice shall be included
 
17
 * in all copies or substantial portions of the Software.
 
18
 *
 
19
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 
20
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 
21
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 
22
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 
23
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 
24
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 
25
 * DEALINGS IN THE SOFTWARE.
 
26
******************************************************************************/
 
27
 
 
28
#include <assert.h>
 
29
#include <iostream>
 
30
#include <fstream>
 
31
#include <map>
 
32
#include <queue>
 
33
#include <cmath>
 
34
#include <stdint.h>
 
35
 
 
36
using namespace std;
 
37
 
 
38
#define INSERT 1
 
39
#define DELETE 0
 
40
#define QUERY 2
 
41
 
 
42
#if defined _WIN32 || defined _WIN64 || defined WIN32 || defined WIN64
 
43
#if _MSC_VER <= 1500
 
44
  typedef __int8 int8_t;
 
45
  typedef __int16 int16_t;
 
46
  typedef __int32 int32_t;
 
47
  typedef __int64 int64_t;
 
48
  typedef unsigned __int8 uint8_t;
 
49
  typedef unsigned __int16 uint16_t;
 
50
  typedef unsigned __int32 uint32_t;
 
51
  typedef unsigned __int64 uint64_t;
 
52
#endif
 
53
 
 
54
// Nuke this annoying warning.  See http://www.unknownroad.com/rtfm/VisualStudio/warningC4251.html
 
55
#pragma warning( disable: 4251 )
 
56
 
 
57
#else
 
58
  #include <stdint.h>
 
59
#endif
 
60
 
 
61
class Rectangle
 
62
{
 
63
public:
 
64
        Rectangle(double xlow, double xhigh, double ylow, double yhigh)
 
65
                : m_xlow(xlow), m_xhigh(xhigh), m_ylow(ylow), m_yhigh(yhigh) {}
 
66
 
 
67
        size_t computeCode(double x, double y)
 
68
        {
 
69
                size_t c = 0;
 
70
                if(y > m_yhigh) c |= TOP;
 
71
                else if(y < m_ylow) c |= BOTTOM;
 
72
                if( x > m_xhigh) c |= RIGHT;
 
73
                else if(x < m_xlow) c |= LEFT;
 
74
                return c;
 
75
        }
 
76
 
 
77
        bool intersectsPoint(double x, double y)
 
78
        {
 
79
                if (m_xlow <= x && x <= m_xhigh &&
 
80
                        m_ylow <= y && y <= m_yhigh) return true;
 
81
                return false;
 
82
        }
 
83
 
 
84
        bool intersectsSegment(double x0, double x1, double y0, double y1)
 
85
        {
 
86
                size_t C0, C1, C;
 
87
                double x,y;
 
88
 
 
89
                C0 = computeCode(x0, y0);
 
90
                C1 = computeCode(x1, y1);
 
91
 
 
92
                for(;;)
 
93
                {
 
94
                        /*
 
95
                         * trivial accept: both ends inside rectangle
 
96
                         */
 
97
                        if((C0 | C1) == 0)
 
98
                        {
 
99
                                return true;
 
100
                        }
 
101
 
 
102
                        /*
 
103
                         * trivial reject: both ends on the external side
 
104
                         * of the rectanlge
 
105
                         */
 
106
                        if((C0 & C1) != 0)
 
107
                        {
 
108
                                return false;
 
109
                        }
 
110
 
 
111
                        /*
 
112
                         * normal case: clip end outside rectangle
 
113
                         */
 
114
                        C = C0 ? C0 : C1;
 
115
                        if(C & TOP)
 
116
                        {
 
117
                                x = x0 + (x1 - x0) * (m_yhigh - y0) / (y1 - y0);
 
118
                                y = m_yhigh;
 
119
                        }
 
120
                        else if(C & BOTTOM)
 
121
                        {
 
122
                                x = x0 + (x1 - x0) * (m_ylow  - y0) / (y1 - y0);
 
123
                                y = m_ylow;
 
124
                        }
 
125
                        else if(C & RIGHT)
 
126
                        {
 
127
                                x = m_xhigh;
 
128
                                y = y0 + (y1 - y0) * (m_xhigh - x0) / (x1 - x0);
 
129
                        }
 
130
                        else
 
131
                        {
 
132
                                x = m_xlow;
 
133
                                y = y0 + (y1 - y0) * (m_xlow - x0) / (x1 - x0);
 
134
                        }
 
135
 
 
136
                        /*
 
137
                         * set new end point and iterate
 
138
                         */
 
139
                        if(C == C0)
 
140
                        {
 
141
                                x0 = x; y0 = y;
 
142
                                C0 = computeCode(x0, y0);
 
143
                        }
 
144
                        else
 
145
                        {
 
146
                                x1 =x; y1 = y;
 
147
                                C1 = computeCode(x1, y1);
 
148
                        }
 
149
                }
 
150
        }
 
151
 
 
152
        static const size_t TOP = 0x1;
 
153
        static const size_t BOTTOM = 0x2;
 
154
        static const size_t RIGHT = 0x4;
 
155
        static const size_t LEFT = 0x8;
 
156
 
 
157
        double m_xlow, m_xhigh, m_ylow, m_yhigh;
 
158
};
 
159
 
 
160
class MovingPoint
 
161
{
 
162
public:
 
163
        MovingPoint(double ax, double vx, double ay, double vy, double rt)
 
164
                : m_ax(ax), m_vx(vx), m_ay(ay), m_vy(vy), m_rt(rt) {}
 
165
 
 
166
        double getX(double t) { return m_ax + m_vx * (t - m_rt); }
 
167
        double getY(double t) { return m_ay + m_vy * (t - m_rt); }
 
168
 
 
169
        double m_ax, m_vx, m_ay, m_vy;
 
170
        double m_rt;
 
171
};
 
172
 
 
173
class TimeRectangle
 
174
{
 
175
public:
 
176
        TimeRectangle(double xlow, double xhigh, double ylow, double yhigh, double tlow, double thigh)
 
177
                : m_xlow(xlow), m_xhigh(xhigh), m_ylow(ylow), m_yhigh(yhigh), m_tlow(tlow), m_thigh(thigh) {}
 
178
 
 
179
        bool intersects(MovingPoint& mp)
 
180
        {
 
181
                double x0 = mp.getX(m_tlow);
 
182
                double x1 = mp.getX(m_thigh);
 
183
                double y0 = mp.getY(m_tlow);
 
184
                double y1 = mp.getY(m_thigh);
 
185
                //double t0 = m_tlow;
 
186
                //double t1 = m_thigh;
 
187
 
 
188
                Rectangle rxy(m_xlow, m_xhigh, m_ylow, m_yhigh);
 
189
                return rxy.intersectsSegment(x0, x1, y0, y1);
 
190
 
 
191
/*
 
192
                // not needed to check all planes since it is
 
193
                // guaranteed that on the time dimension
 
194
                // the line segment and the query cube have
 
195
                // exactly the same length (thus, if they intersect
 
196
                // they should intersect on the X-Y projection for sure).
 
197
 
 
198
                Rectangle rxy(m_xlow, m_xhigh, m_ylow, m_yhigh);
 
199
                if (rxy.intersectsSegment(x0, x1, y0, y1))
 
200
                {
 
201
                        Rectangle rxt(m_xlow, m_xhigh, m_tlow, m_thigh);
 
202
                        if (rxt.intersectsSegment(x0, x1, t0, t1))
 
203
                        {
 
204
                                Rectangle ryt(m_ylow, m_yhigh, m_tlow, m_thigh);
 
205
                                if (ryt.intersectsSegment(y0, y1, t0, t1)) return true;
 
206
                        }
 
207
                }
 
208
*/
 
209
 
 
210
                return false;
 
211
        }
 
212
 
 
213
        bool intersectsStupid(MovingPoint& mp)
 
214
        {
 
215
                size_t t0 = static_cast<size_t>(std::floor(m_tlow));
 
216
                size_t t1 = static_cast<size_t>(std::floor(m_thigh));
 
217
 
 
218
                Rectangle rxy(m_xlow, m_xhigh, m_ylow, m_yhigh);
 
219
 
 
220
                for (size_t T = t0; T <= t1; T++)
 
221
                {
 
222
                        if (rxy.intersectsPoint(mp.getX(T), mp.getY(T))) return true;
 
223
                }
 
224
                return false;
 
225
        }
 
226
 
 
227
        double m_xlow, m_xhigh, m_ylow, m_yhigh;
 
228
        double m_tlow, m_thigh;
 
229
};
 
230
 
 
231
int main(int argc, char** argv)
 
232
{
 
233
        if (argc != 2)
 
234
        {
 
235
                cerr << "Usage: " << argv[0] << " data_file." << endl;
 
236
                return -1;
 
237
        }
 
238
 
 
239
        ifstream fin(argv[1]);
 
240
        if (! fin)
 
241
        {
 
242
                cerr << "Cannot open data file" << argv[1] << "." << endl;
 
243
                return -1;
 
244
        }
 
245
 
 
246
        map<size_t, MovingPoint> data;
 
247
        size_t id, op;
 
248
        double ax, vx, ay, vy, ct, rt, unused;
 
249
 
 
250
        while (fin)
 
251
        {
 
252
                fin >> id >> op >> ct >> rt >> unused >> ax >> vx >> unused >> ay >> vy;
 
253
                if (! fin.good()) continue;
 
254
 
 
255
                if (op == INSERT)
 
256
                {
 
257
                        data.insert(pair<size_t, MovingPoint>(id, MovingPoint(ax, vx, ay, vy, ct)));
 
258
                }
 
259
                else if (op == DELETE)
 
260
                {
 
261
                        data.erase(id);
 
262
                }
 
263
                else if (op == QUERY)
 
264
                {
 
265
                        TimeRectangle query = TimeRectangle(ax, vx, ay, vy, ct, rt);
 
266
                        std::map<size_t, MovingPoint>::iterator it;
 
267
                        for (it = data.begin(); it != data.end(); it++)
 
268
                        {
 
269
                                //assert(query.intersects((*it).second) == query.intersectsStupid((*it).second));
 
270
                                if (query.intersects((*it).second) == false && query.intersectsStupid((*it).second) == true)
 
271
                                {
 
272
                                        cerr << "Something is wrong: " << ct << " " << (*it).first << endl;
 
273
                                        return -1;
 
274
                                }
 
275
                                if (query.intersects((*it).second)) cout << (*it).first << endl;
 
276
                        }
 
277
                }
 
278
        }
 
279
 
 
280
        return 0;
 
281
}