1 #ifdef __cplusplus
2 extern "C" {
3 #endif
4 #include "EXTERN.h"
5 #include "perl.h"
6 #include "XSUB.h"
7 #ifdef __cplusplus
8 }
9 #endif
10 
11 #include "ppport.h"
12 
13 #include <iostream>
14 
15 #include "SpatialVector.h"
16 #include "SpatialIndex.h"
17 #include "RangeConvex.h"
18 #include "HtmRange.h"
19 #include "HtmRangeIterator.h"
20 
21 
22 class HTMesh {
23 int level, build_level;
24 SpatialIndex *htm;
25 char name_buffer[128];
26 HtmRange* range;
27 HtmRangeIterator* iterator;
28 long numTriangles;
29 
30 public:
HTMesh(int level_in,int build_level_in)31         HTMesh(int level_in, int build_level_in) {
32             level = level_in;
33             build_level = build_level_in;
34             if (build_level > 0) {
35                 if (build_level > level) {
36                     build_level = level;
37                 }
38                 htm = new SpatialIndex(level, build_level);
39             }
40             else {
41                 htm = new SpatialIndex(level);
42             }
43 
44             numTriangles = 8;
45             for (int i = level; i--;) {
46                 numTriangles *= 4;
47             }
48         }
49 
~HTMesh()50         ~HTMesh() {
51             delete htm;
52         }
53 
54         /* long lookupId(double ra, double dec) { */
55         /*      return htm->idByPoint(ra * 15.0, dec); */
56         /* } */
57 
idToName(long id)58         const char* idToName(long id) {
59             htm->nameById(id, name_buffer);
60             return name_buffer;
61         }
62 
63         /* const char* lookupName(double ra, double dec) { */
64         /*     long id = lookupId(ra, dec); */
65         /*     return idToName(id); */
66         /* } */
67 
intersectCircle(double ra,double dec,double rad)68         void intersectCircle(double ra, double dec, double rad) {
69             ra = ra * 15.0;
70             RangeConvex convex;
71             double d = cos( 3.1415926535897932385E0 * rad/180.0);
72             SpatialConstraint c(SpatialVector(ra, dec), d);
73             convex.add(c); // [ed:RangeConvex::add]
74             convex.setOlevel(level);
75             range = new HtmRange();
76             convex.intersect(htm, range); // commit fec61ac2c400e7138425961e779fed3bbf5687a5 removed unmanipulated bool parameter from RangeConvex::intersect -- asimha
77             iterator = new HtmRangeIterator(range);
78         }
79 
intersectRectangle(double ra,double dec,double dra,double ddec)80         void intersectRectangle(double ra, double dec, double dra, double ddec) {
81             dra = dra / 15.0;
82             double ra1 = ra - dra;
83             double ra2 = ra + dra;
84             double dec1 = dec - ddec;
85             double dec2 = dec + ddec;
86             intersectPolygon4(ra1, dec1, ra1, dec2, ra2, dec2, ra2, dec1);
87         }
88 
intersectPolygon4(double ra1,double dec1,double ra2,double dec2,double ra3,double dec3,double ra4,double dec4)89         void intersectPolygon4(double ra1, double dec1, double ra2, double dec2,
90                                double ra3, double dec3, double ra4, double dec4) {
91             SpatialVector* corner[4];
92             corner[0] = new SpatialVector(ra1 * 15.0, dec1);
93             corner[1] = new SpatialVector(ra2 * 15.0, dec2);
94             corner[2] = new SpatialVector(ra3 * 15.0, dec3);
95             corner[3] = new SpatialVector(ra4 * 15.0, dec4);
96             RangeConvex* convex = new RangeConvex(corner[0], corner[1], corner[2], corner[3]);
97             convex->setOlevel(level);
98             range = new HtmRange();
99             convex->intersect(htm, range); // commit fec61ac2c400e7138425961e779fed3bbf5687a5 removed unmanipulated bool parameter from RangeConvex::intersect -- asimha
100             iterator = new HtmRangeIterator(range);
101         }
102 
timeRectangle(int iter,double ra,double dec,double dra,double ddec)103         double timeRectangle(int iter, double ra, double dec, double dra, double ddec) {
104             int i = iter;
105             time_t t0 = clock();
106             while(i--) {
107                 intersectRectangle(ra, dec, dra, ddec);
108             }
109             time_t t1 = clock();
110             double nsec = (double)(t1-t0)/(double)CLOCKS_PER_SEC ;
111             return nsec;
112         }
113 
timeCircle(int iter,double ra,double dec,double rad)114         double timeCircle(int iter, double ra, double dec, double rad) {
115             int i = iter;
116             time_t t0 = clock();
117             while(i--) {
118                 intersectCircle(ra, dec, rad);
119             }
120             time_t t1 = clock();
121             double nsec = (double)(t1-t0)/(double)CLOCKS_PER_SEC ;
122             return nsec;
123         }
124 
getIterator()125         HtmRangeIterator* getIterator() {
126             return iterator;
127         }
128 
getHtm()129         SpatialIndex* getHtm() {
130             return htm;
131         }
132 
hasNext()133         bool hasNext() {
134             return iterator->hasNext();
135         }
136 
nextId()137         long nextId() {
138             return iterator->next();
139         }
nextName()140         const char * nextName() {
141             return iterator->nextSymbolic(name_buffer);
142         }
143 
totalTriangles()144         long totalTriangles() {
145             return numTriangles;
146         }
147 
resultSize()148         int resultSize() {
149             HtmRangeIterator* iter;
150             int size = 0;
151             iter = new HtmRangeIterator(range);
152             while (iter->hasNext()) {
153                 size++;
154                 iter->next();
155             }
156             return size;
157         }
158 };
159 
160 MODULE = HTMesh		PACKAGE = HTMesh
161 
162 HTMesh *
163 HTMesh::new(int level, int build_level=0)
164 
165 void
166 HTMesh::DESTROY()
167 
168 void
169 HTMesh::intersect_rect(double ra, double dec, double dra, double ddec)
170 CODE:
171     THIS->intersectRectangle(ra, dec, dra, ddec);
172 
173 void
174 HTMesh::intersect_poly4(double ra1, double dec1, double ra2, double dec2, double ra3, double dec3, double ra4, double dec4)
175 CODE:
176     THIS->intersectPolygon4(ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4);
177 
178 void
179 HTMesh::intersect_circle(double ra, double dec, double rad)
180 CODE:
181     THIS->intersectCircle(ra, dec, rad);
182 
183 double
184 HTMesh::time_rect(int iterations, double ra, double dec, double dra, double ddec)
185 CODE:
186     RETVAL = THIS->timeRectangle(iterations, ra, dec, dra, ddec);
187 OUTPUT:
188     RETVAL
189 
190 double
191 HTMesh::time_circle(int iter, double ra, double dec, double rad)
192 CODE:
193     RETVAL = THIS->timeCircle(iter, ra, dec, rad);
194 OUTPUT:
195     RETVAL
196 
197 void
198 HTMesh::results()
199 PPCODE:
200         HtmRangeIterator* iterator = THIS->getIterator();
201         while ( iterator->hasNext() ) {
202                 XPUSHs(sv_2mortal(newSVnv(iterator->next())));
203         }
204 
205 void
206 HTMesh::next_triangle()
207 PPCODE:
208         long id = THIS->nextId();
209         SpatialIndex* htm = THIS->getHtm();
210         SpatialVector v1, v2, v3;
211         htm->nodeVertex(id, v1, v2, v3);
212 
213         XPUSHs(sv_2mortal(newSVnv( v1.ra() / 15.0 )));
214         XPUSHs(sv_2mortal(newSVnv( v1.dec() )));
215 
216         XPUSHs(sv_2mortal(newSVnv( v2.ra() / 15.0 )));
217         XPUSHs(sv_2mortal(newSVnv( v2.dec() )));
218 
219         XPUSHs(sv_2mortal(newSVnv( v3.ra() / 15.0 )));
220         XPUSHs(sv_2mortal(newSVnv( v3.dec() )));
221 
222         XPUSHs(sv_2mortal(newSVnv( id )));
223 
224 bool
225 HTMesh::has_next()
226 CODE:
227     RETVAL = THIS->hasNext();
228 OUTPUT:
229     RETVAL
230 
231 const char*
232 HTMesh::id_to_name(long id)
233 CODE:
234     RETVAL = THIS->idToName(id);
235 OUTPUT:
236     RETVAL
237 
238 long
239 HTMesh::next_id()
240 CODE:
241     RETVAL = THIS->nextId();
242 OUTPUT:
243     RETVAL
244 
245 const char *
246 HTMesh::next_name()
247 CODE:
248     RETVAL = THIS->nextName();
249 OUTPUT:
250     RETVAL
251 
252 long
253 HTMesh::total_triangles()
254 CODE:
255     RETVAL = THIS->totalTriangles();
256 OUTPUT:
257     RETVAL
258 
259 int
260 HTMesh::result_size()
261 CODE:
262     RETVAL = THIS->resultSize();
263 OUTPUT:
264     RETVAL
265 
266