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