tools/HTMesh-0.01/HTMesh.xs

0001 #ifdef __cplusplus
0002 extern "C" {
0003 #endif
0004 #include "EXTERN.h"
0005 #include "perl.h"
0006 #include "XSUB.h"
0007 #ifdef __cplusplus
0008 }
0009 #endif
0010
0011 #include "ppport.h"
0012
0013 #include <iostream>
0014
0015 #include "SpatialVector.h"
0016 #include "SpatialIndex.h"
0017 #include "RangeConvex.h"
0018 #include "HtmRange.h"
0019 #include "HtmRangeIterator.h"
0020
0021
0022 class HTMesh {
0023 int level, build_level;
0024 SpatialIndex *htm;
0025 char name_buffer[128];
0026 HtmRange* range;
0027 HtmRangeIterator* iterator;
0028 long numTriangles;
0029
0030 public:
0031         HTMesh(int level_in, int build_level_in) {
0032             level = level_in;
0033             build_level = build_level_in;
0034             if (build_level > 0) {
0035                 if (build_level > level) {
0036                     build_level = level;
0037                 }
0038                 htm = new SpatialIndex(level, build_level);
0039             }
0040             else {
0041                 htm = new SpatialIndex(level);
0042             }
0043
0044             numTriangles = 8;
0045             for (int i = level; i--;) {
0046                 numTriangles *= 4;
0047             }
0048         }
0049
0050         ~HTMesh() {
0051             delete htm;
0052         }
0053
0054         /* long lookupId(double ra, double dec) { */
0055         /*      return htm->idByPoint(ra * 15.0, dec); */
0056         /* } */
0057
0058         const char* idToName(long id) {
0059             htm->nameById(id, name_buffer);
0060             return name_buffer;
0061         }
0062
0063         /* const char* lookupName(double ra, double dec) { */
0064         /*     long id = lookupId(ra, dec); */
0065         /*     return idToName(id); */
0066         /* } */
0067
0068         void intersectCircle(double ra, double dec, double rad) {
0069             ra = ra * 15.0;
0070             RangeConvex convex;
0071             double d = cos( 3.1415926535897932385E0 * rad/180.0);
0072             SpatialConstraint c(SpatialVector(ra, dec), d);
0073             convex.add(c); // [ed:RangeConvex::add]
0074             convex.setOlevel(level);
0075             range = new HtmRange();
0076             convex.intersect(htm, range); // commit fec61ac2c400e7138425961e779fed3bbf5687a5 removed unmanipulated bool parameter from RangeConvex::intersect -- asimha
0077             iterator = new HtmRangeIterator(range);
0078         }
0079
0080         void intersectRectangle(double ra, double dec, double dra, double ddec) {
0081             dra = dra / 15.0;
0082             double ra1 = ra - dra;
0083             double ra2 = ra + dra;
0084             double dec1 = dec - ddec;
0085             double dec2 = dec + ddec;
0086             intersectPolygon4(ra1, dec1, ra1, dec2, ra2, dec2, ra2, dec1);
0087         }
0088
0089         void intersectPolygon4(double ra1, double dec1, double ra2, double dec2,
0090                                double ra3, double dec3, double ra4, double dec4) {
0091             SpatialVector* corner[4];
0092             corner[0] = new SpatialVector(ra1 * 15.0, dec1);
0093             corner[1] = new SpatialVector(ra2 * 15.0, dec2);
0094             corner[2] = new SpatialVector(ra3 * 15.0, dec3);
0095             corner[3] = new SpatialVector(ra4 * 15.0, dec4);
0096             RangeConvex* convex = new RangeConvex(corner[0], corner[1], corner[2], corner[3]);
0097             convex->setOlevel(level);
0098             range = new HtmRange();
0099             convex->intersect(htm, range); // commit fec61ac2c400e7138425961e779fed3bbf5687a5 removed unmanipulated bool parameter from RangeConvex::intersect -- asimha
0100             iterator = new HtmRangeIterator(range);
0101         }
0102
0103         double timeRectangle(int iter, double ra, double dec, double dra, double ddec) {
0104             int i = iter;
0105             time_t t0 = clock();
0106             while(i--) {
0107                 intersectRectangle(ra, dec, dra, ddec);
0108             }
0109             time_t t1 = clock();
0110             double nsec = (double)(t1-t0)/(double)CLOCKS_PER_SEC ;
0111             return nsec;
0112         }
0113
0114         double timeCircle(int iter, double ra, double dec, double rad) {
0115             int i = iter;
0116             time_t t0 = clock();
0117             while(i--) {
0118                 intersectCircle(ra, dec, rad);
0119             }
0120             time_t t1 = clock();
0121             double nsec = (double)(t1-t0)/(double)CLOCKS_PER_SEC ;
0122             return nsec;
0123         }
0124
0125         HtmRangeIterator* getIterator() {
0126             return iterator;
0127         }
0128
0129         SpatialIndex* getHtm() {
0130             return htm;
0131         }
0132
0133         bool hasNext() {
0134             return iterator->hasNext();
0135         }
0136
0137         long nextId() {
0138             return iterator->next();
0139         }
0140         const char * nextName() {
0141             return iterator->nextSymbolic(name_buffer);
0142         }
0143
0144         long totalTriangles() {
0145             return numTriangles;
0146         }
0147
0148         int resultSize() {
0149             HtmRangeIterator* iter;
0150             int size = 0;
0151             iter = new HtmRangeIterator(range);
0152             while (iter->hasNext()) {
0153                 size++;
0154                 iter->next();
0155             }
0156             return size;
0157         }
0158 };
0159
0160 MODULE = HTMesh         PACKAGE = HTMesh
0161
0162 HTMesh *
0163 HTMesh::new(int level, int build_level=0)
0164
0165 void
0166 HTMesh::DESTROY()
0167
0168 void
0169 HTMesh::intersect_rect(double ra, double dec, double dra, double ddec)
0170 CODE:
0171     THIS->intersectRectangle(ra, dec, dra, ddec);
0172
0173 void
0174 HTMesh::intersect_poly4(double ra1, double dec1, double ra2, double dec2, double ra3, double dec3, double ra4, double dec4)
0175 CODE:
0176     THIS->intersectPolygon4(ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4);
0177
0178 void
0179 HTMesh::intersect_circle(double ra, double dec, double rad)
0180 CODE:
0181     THIS->intersectCircle(ra, dec, rad);
0182
0183 double
0184 HTMesh::time_rect(int iterations, double ra, double dec, double dra, double ddec)
0185 CODE:
0186     RETVAL = THIS->timeRectangle(iterations, ra, dec, dra, ddec);
0187 OUTPUT:
0188     RETVAL
0189
0190 double
0191 HTMesh::time_circle(int iter, double ra, double dec, double rad)
0192 CODE:
0193     RETVAL = THIS->timeCircle(iter, ra, dec, rad);
0194 OUTPUT:
0195     RETVAL
0196
0197 void
0198 HTMesh::results()
0199 PPCODE:
0200         HtmRangeIterator* iterator = THIS->getIterator();
0201         while ( iterator->hasNext() ) {
0202                 XPUSHs(sv_2mortal(newSVnv(iterator->next())));
0203         }
0204
0205 void
0206 HTMesh::next_triangle()
0207 PPCODE:
0208         long id = THIS->nextId();
0209         SpatialIndex* htm = THIS->getHtm();
0210         SpatialVector v1, v2, v3;
0211         htm->nodeVertex(id, v1, v2, v3);
0212
0213         XPUSHs(sv_2mortal(newSVnv( v1.ra() / 15.0 )));
0214         XPUSHs(sv_2mortal(newSVnv( v1.dec() )));
0215
0216         XPUSHs(sv_2mortal(newSVnv( v2.ra() / 15.0 )));
0217         XPUSHs(sv_2mortal(newSVnv( v2.dec() )));
0218
0219         XPUSHs(sv_2mortal(newSVnv( v3.ra() / 15.0 )));
0220         XPUSHs(sv_2mortal(newSVnv( v3.dec() )));
0221
0222         XPUSHs(sv_2mortal(newSVnv( id )));
0223
0224 bool
0225 HTMesh::has_next()
0226 CODE:
0227     RETVAL = THIS->hasNext();
0228 OUTPUT:
0229     RETVAL
0230
0231 const char*
0232 HTMesh::id_to_name(long id)
0233 CODE:
0234     RETVAL = THIS->idToName(id);
0235 OUTPUT:
0236     RETVAL
0237
0238 long
0239 HTMesh::next_id()
0240 CODE:
0241     RETVAL = THIS->nextId();
0242 OUTPUT:
0243     RETVAL
0244
0245 const char *
0246 HTMesh::next_name()
0247 CODE:
0248     RETVAL = THIS->nextName();
0249 OUTPUT:
0250     RETVAL
0251
0252 long
0253 HTMesh::total_triangles()
0254 CODE:
0255     RETVAL = THIS->totalTriangles();
0256 OUTPUT:
0257     RETVAL
0258
0259 int
0260 HTMesh::result_size()
0261 CODE:
0262     RETVAL = THIS->resultSize();
0263 OUTPUT:
0264     RETVAL
0265