Thierry Louge / gammalib

Commit 08ec3fe0aaa6a831d566dd0137cbbee1622712fa

Authored by Maria Krause
1 parent 37bbe9c2

corrected contains and overlaps function

Inline Side-by-side
Showing 3 changed files with 79 additions and 96 deletions
include/GSkyRegionMap.hpp
  View file @08ec3fe
... ... @@ -33,6 +33,7 @@
33 33 #include "GMath.hpp"
34 34 #include "GSkyDir.hpp"
35 35 #include "GSkyRegion.hpp"
  36 +#include "GSkyRegions.hpp"
36 37 #include "GFits.hpp"
37 38 #include "GSkymap.hpp"
38 39 #include "GSkyRegionCircle.hpp"
... ... @@ -71,9 +72,7 @@ public:
71 72 const GSkymap& map(void) const;
72 73 bool contains(const GSkyDir& dir) const;
73 74 bool contains(const GSkyRegion& reg) const;
74   - bool contains(const GSkyPixel& pixel) const;
75 75 bool overlaps(const GSkyRegion& reg) const;
76   - bool overlaps(const GSkyDir& dir) const;
77 76 virtual std::string print(const GChatter& chatter = NORMAL) const;
78 77 std::string write(void) const;
79 78 virtual void read(const std::string& regstring);
... ... @@ -86,30 +85,14 @@ protected:
86 85 void copy_members(const GSkyRegionMap& region);
87 86 void free_members(void);
88 87 void compute_solid_angle(void);
89   - void pixels_sky_region_map(void);
  88 + void get_array_nonzero(std::vector<int> pix_array) const;
90 89  
91 90 // Protected members
92 91 GSkymap m_map; //!< The map
93 92 GSkyPixel dir2pix(const GSkyDir& dir) const;
94   - GSkyDir m_centre; //!< Centre or reference point of the region
95   - double m_radius; //!< Radius of circular the region [deg]
96 93  
97 94 };
98 95  
99   -
100   -/***********************************************************************//**
101   - * @brief Return region sky map
102   - *
103   - * @return region sky map.
104   - *
105   - * Returns the region sky map.
106   - ***************************************************************************/
107   -/*inline
108   -const double& GSkyRegionMap::map(void) const
109   -{
110   - return (m_map);
111   -}
112   -*/
113 96 /***********************************************************************//**
114 97 * @brief Set sky map region
115 98 *
... ...
src/sky/GSkyRegionMap.cpp
  View file @08ec3fe
... ... @@ -86,6 +86,10 @@ GSkyRegionMap::GSkyRegionMap(const GSkyRegionMap&amp; region)
86 86  
87 87 // Compute solid angle
88 88 compute_solid_angle();
  89 +
  90 + // Get non-zero indices
  91 + std::vector<int> pix_array;
  92 + get_array_nonzero(pix_array);
89 93  
90 94 // Return
91 95 return;
... ... @@ -386,58 +390,70 @@ void GSkyRegionMap::read(const GFitsHDU&amp; hdu)
386 390 }
387 391 */
388 392 /***********************************************************************//**
389   - * @brief Verifies if sky direction falls in map
  393 + * @brief Tells if direction is contained in one of the region maps
  394 + *
  395 + * @param[in] dir A Sky direction.
390 396 *
391   - * @param[in] dir Sky direction.
  397 + * @return True or False
392 398 *
393   - * This method checks if the specified sky direction falls within the pixels
394   - * covered by the skymap region. The method uses the dir2xy method to convert
395   - * the sky direction into 2D pixel indices, and then checks whether the pixel
396   - * indices fall in the skymap.
  399 + * Tells if direction is contained in one of the regions
397 400 ***************************************************************************/
398 401 bool GSkyRegionMap::contains(const GSkyDir& dir) const
399 402 {
  403 + bool fully_inside = false;
  404 +
  405 + // Get non zeros indices of sky region map
  406 + std::vector<int> pix_array;
  407 + get_array_nonzero(pix_array);
  408 +
400 409 // Convert sky direction into sky pixel
401 410 GSkyPixel pixel = dir2pix(dir);
402 411  
403   - // Return location flag
404   - return (contains(pixel));
  412 + for (int i = 0; i < pix_array.size(); i++ ) {
  413 + GSkyDir direc = m_map.pix2dir(pix_array[i]);
  414 + if (i < pix_array.size() && dir==direc) {
  415 + fully_inside = true;
  416 + }
  417 + else {
  418 + throw GException::feature_not_implemented(G_CONTAINS,
  419 + "Cannot compare two different region types yet");
  420 + }
  421 + }
  422 +
  423 + // Return value
  424 + return fully_inside;
405 425 }
406 426  
407 427 /***********************************************************************//**
408   - * @brief Checks if sky map pixel falls in map
  428 + * @brief Checks if region is fully contained within this region
409 429 *
410   - * @param[in] pixel Sky map pixel.
411   - * @return Trus if pixels is within map, false otherwise.
  430 + * @param[in] reg Sky region map.
412 431 *
413   - * Checks whether the specified sky map @p pixel falls within the skymap
414   - * or not.
  432 + * @exception GException::feature_not_implemented
  433 + * Regions differ in type.
415 434 ***************************************************************************/
416   -bool GSkymap::contains(const GSkyPixel& pixel) const
  435 +bool GSkyRegionMap::contains(const GSkyRegion& reg) const
417 436 {
418   - // Initialise containment flag
419   - bool inmap = false;
  437 + // Initialise return value
  438 + bool fully_inside = false;
420 439  
421   - // Test 1D pixels
422   - if (pixel.is_1D()) {
423   - if (pixel.index()+0.5 >= 0.0 && pixel.index()-0.5 < m_num_pixels) {
424   - inmap = true;
  440 + // Get non zeros indices of sky region map
  441 + std::vector<int> pix_array;
  442 + get_array_nonzero(pix_array);
  443 +
  444 + for (int i = 0; i < pix_array.size(); i++ ) {
  445 + GSkyDir dir = m_map.pix2dir(pix_array[i]);
  446 + if (i < pix_array.size() && reg.contains(dir)) {
  447 + fully_inside = true;
425 448 }
426   - }
427   -
428   - // Test 2D pixels
429   - else if (pixel.is_2D()) {
430   -
431   - // If pixel is in range then set containment flag to true
432   - if ((pixel.x()+0.5 >= 0.0 && pixel.x()-0.5 < m_num_x) &&
433   - (pixel.y()+0.5 >= 0.0 && pixel.y()-0.5 < m_num_y)) {
434   - inmap = true;
  449 + else {
  450 + throw GException::feature_not_implemented(G_CONTAINS,
  451 + "Cannot compare two different region types yet");
435 452 }
436   -
437 453 }
438   -
439   - // Return containment flag
440   - return inmap;
  454 +
  455 + // Return value
  456 + return fully_inside;
441 457 }
442 458  
443 459 /***********************************************************************//**
... ... @@ -445,43 +461,27 @@ bool GSkymap::contains(const GSkyPixel&amp; pixel) const
445 461 *
446 462 * @param[in] reg Sky region.
447 463 *
448   - * @TODO: region sky map overlap with ring or circle region
449   - *
450 464 * @exception GException::feature_not_implemented
451 465 * Regions differ in type.
452 466 ***************************************************************************/
453   -bool GSkyRegionMap::overlaps(const GSkyDir& dir) const
  467 +bool GSkyRegionMap::overlaps(const GSkyRegion& reg) const
454 468 {
455 469 // Initialise return value
456 470 bool overlap = false;
457 471  
458   - // Convert sky direction into sky pixel
459   - GSkyPixel pixel = dir2pix(dir);
460   -
461   - // Return location flag
462   - return (contains(pixel));
  472 + // Get non zeros indices of sky region map
  473 + std::vector<int> pix_array;
  474 + get_array_nonzero(pix_array);
463 475  
464   -/* GSkyDir dir;
465   - for (int i=0; i<360; i++){
466   - for (int j=0; j<90; j++){
467   -
468   - const GSkyDir* regdir =
469   - dynamic_cast<const GSkyDir*>(&dir);
470   - double direction = dir->radec_deg(i,j);
471   -
472   - if (reg->contains(dir)) {
473   - overlap = true;
474   - }
475   - else {
476   - overlap = false;
477   - }
478   - }
  476 + for (int i = 0; i < pix_array.size(); i++ ) {
  477 + GSkyDir dir = m_map.pix2dir(pix_array[i]);
  478 + reg.contains(dir);
  479 + if (overlap) break;
479 480 }
480   - // Return value
481   - return overlap;
482   - */
483   -}
484 481  
  482 + // Return if region is in map or not
  483 + return overlap;
  484 +}
485 485  
486 486  
487 487 /*==========================================================================
... ... @@ -539,7 +539,6 @@ void GSkyRegionMap::compute_solid_angle(void)
539 539 for (int i =0; i < m_map.npix(); ++i) {
540 540 GSkyPixel pix=m_map.inx2pix(i);
541 541 std::cout << pix.size() << std::endl;
542   - //value += m_map.solidangle(pix);
543 542 if (pix.size() > 0){
544 543 value += m_map.solidangle(pix);
545 544 }
... ... @@ -557,20 +556,21 @@ void GSkyRegionMap::compute_solid_angle(void)
557 556  
558 557  
559 558 /***********************************************************************//**
560   - * @brief Get an array of pixels of the sky region map
  559 + * @brief Create an array of pixels
561 560 ***************************************************************************/
562   - void GSkyRegionMap::pixels_sky_region_map(void)
563   - {
564   - double num_pix_map[m_map.npix()];
565   -
566   - for (int i =0; i < m_map.npix(); ++i) {
567   - GSkyPixel pix=m_map.inx2pix(i);
568   - std::cout << pix.size() << std::endl;
  561 +void GSkyRegionMap::get_array_nonzero(std::vector<int> pix_array) const
  562 +{
  563 + //std::vector<int> pix_array;
  564 +
  565 + for (int i = 0; i < m_map.npix(); i++)
  566 + {
  567 + GSkyPixel pix=m_map.inx2pix(i);
  568 + std::cout << pix.size() << std::endl;
569 569 if (pix.size() > 0){
570   - num_pix_map[i] = i;
571   - }
572   - }
573   -
574   - // Return
575   - return;
576   -}
577 570 \ No newline at end of file
  571 + pix_array.push_back(pix);
  572 + }
  573 + }
  574 +
  575 + std::cout << pix_array.size() << std::endl;
  576 +
  577 +}
... ...
src/sky/GSkyRegions.cpp
  View file @08ec3fe
... ... @@ -842,7 +842,7 @@ std::string GSkyRegions::print(const GChatter&amp; chatter) const
842 842  
843 843  
844 844 /***********************************************************************//**
845   - * @brief Tells if region overlaps one of the regions
  845 + * @brief Tells if region one of the regions
846 846 *
847 847 * @param[in] reg A sky region
848 848 * @return True or False
... ...