Eirik Hatlen / gammalib

Commit cc862b5b376e131a8b99287b31f3f2bb525767f2

Authored by Andreas Specovius
Committed by Jürgen Knödlseder
1 parent de78ac38

Implemented first approach to build a GSkyRegionMap from a GSkyRegionRect (no rotation angle yet).

Inline Side-by-side
Showing 2 changed files with 97 additions and 20 deletions
include/GSkyRegionMap.hpp
  View file @cc862b5
... ... @@ -37,6 +37,7 @@
37 37 /* __ Forward declarations _______________________________________________ */
38 38 class GSkyDir;
39 39 class GSkyRegionMap;
  40 +class GSkyRegionRect;
40 41  
41 42  
42 43 /***********************************************************************//**
... ... @@ -87,6 +88,7 @@ protected:
87 88 void compute_solid_angle(void);
88 89 void set_nonzero_indices(void);
89 90 void set_region_circle(const GSkyRegionCircle* circle);
  91 + void set_region_rect(const GSkyRegionRect* rect);
90 92  
91 93 // Protected members
92 94 GSkyMap m_map; //!< The region map
... ... @@ -108,7 +110,7 @@ std::string GSkyRegionMap::classname(void) const
108 110  
109 111 /***********************************************************************//**
110 112 * @brief Return sky map
111   - *
  113 + *
112 114 * @return region sky map.
113 115 ***************************************************************************/
114 116 inline
... ... @@ -120,7 +122,7 @@ const GSkyMap&amp; GSkyRegionMap::map(void) const
120 122  
121 123 /***********************************************************************//**
122 124 * @brief Get non-zero index vector
123   - *
  125 + *
124 126 * @return Reference to non-zero pixel indices vector.
125 127 ***************************************************************************/
126 128 inline
... ...
src/sky/GSkyRegionMap.cpp
  View file @cc862b5
... ... @@ -34,6 +34,7 @@
34 34 #include "GSkyMap.hpp"
35 35 #include "GSkyRegion.hpp"
36 36 #include "GSkyRegionCircle.hpp"
  37 +#include "GSkyRegionRect.hpp"
37 38 #include "GSkyRegionMap.hpp"
38 39  
39 40 /* __ Method name definitions ____________________________________________ */
... ... @@ -120,6 +121,7 @@ GSkyRegionMap::GSkyRegionMap(const GSkyRegion* region) : GSkyRegion(*region)
120 121 // Get dynamic casts
121 122 const GSkyRegionMap* map = dynamic_cast<const GSkyRegionMap*>(region);
122 123 const GSkyRegionCircle* circle = dynamic_cast<const GSkyRegionCircle*>(region);
  124 + const GSkyRegionRect* rect = dynamic_cast<const GSkyRegionRect*>(region);
123 125  
124 126 // If region is of type GSkyRegionMap then copy class members
125 127 if (map != NULL) {
... ... @@ -132,6 +134,11 @@ GSkyRegionMap::GSkyRegionMap(const GSkyRegion* region) : GSkyRegion(*region)
132 134 set_region_circle(circle);
133 135 }
134 136  
  137 + // If region is of type GSkyRegionRect than set region map fron region rect
  138 + else if (rect != NULL) {
  139 + set_region_rect(rect);
  140 + }
  141 +
135 142 // Return
136 143 return;
137 144 }
... ... @@ -239,7 +246,7 @@ GSkyRegionMap* GSkyRegionMap::clone(void) const
239 246  
240 247 /***********************************************************************//**
241 248 * @brief Set sky map
242   - *
  249 + *
243 250 * @param[in] map Sky map
244 251 *
245 252 * Sets the sky map of a region map.
... ... @@ -251,7 +258,7 @@ void GSkyRegionMap::map(const GSkyMap&amp; map)
251 258  
252 259 // Set non-zero pixel indices
253 260 set_nonzero_indices();
254   -
  261 +
255 262 // Compute solid angle
256 263 compute_solid_angle();
257 264  
... ... @@ -275,10 +282,10 @@ void GSkyRegionMap::load(const GFilename&amp; filename)
275 282  
276 283 // Load map
277 284 m_map.load(filename);
278   -
  285 +
279 286 // Set non-zero indices
280 287 set_nonzero_indices();
281   -
  288 +
282 289 // Compute solid angle
283 290 compute_solid_angle();
284 291  
... ... @@ -365,13 +372,13 @@ bool GSkyRegionMap::contains(const GSkyDir&amp; dir) const
365 372  
366 373 // If direction is within map boundaries
367 374 if (m_map.contains(dir)) {
368   -
  375 +
369 376 // Convert sky direction into pixel index
370 377 int i = m_map.dir2inx(dir);
371 378  
372 379 // Set containment flag
373 380 contain = (m_map(i,0) != 0.0);
374   -
  381 +
375 382 } // endif: map contains direction
376 383  
377 384 // Return value
... ... @@ -393,17 +400,17 @@ bool GSkyRegionMap::contains(const GSkyRegion&amp; reg) const
393 400 {
394 401 // Initialise return value
395 402 bool inside = true;
396   -
  403 +
397 404 // Case of a region map
398 405 if (reg.type() == "Map") {
399   -
  406 +
400 407 // Cast to map type
401 408 const GSkyRegionMap* inreg = dynamic_cast<const GSkyRegionMap*>(&reg);
402   -
  409 +
403 410 // Retrieve map data
404 411 const GSkyMap& regmap = inreg->map();
405 412 const std::vector<int>& regnzvec = inreg->nonzero_indices();
406   -
  413 +
407 414 // Loop over input map non-zero pixels and check if they are all in
408 415 for (int i = 0; i < regnzvec.size(); ++i) {
409 416 if (!contains(regmap.inx2dir(regnzvec[i]))) {
... ... @@ -416,14 +423,14 @@ bool GSkyRegionMap::contains(const GSkyRegion&amp; reg) const
416 423  
417 424 // Case of a circular region
418 425 else if (reg.type() == "Circle") {
419   -
  426 +
420 427 // Create circular region from reg
421 428 const GSkyRegionCircle* inreg = dynamic_cast<const GSkyRegionCircle*>(&reg);
422 429  
423 430 // Retrieve circle data
424 431 const GSkyDir& centre = inreg->centre();
425 432 double rad = inreg->radius();
426   -
  433 +
427 434 // Test a number of points along the circle
428 435 const int numtest = 360;
429 436 double rotangle = 360.0 / numtest;
... ... @@ -444,7 +451,7 @@ bool GSkyRegionMap::contains(const GSkyRegion&amp; reg) const
444 451 throw GException::feature_not_implemented(G_CONTAINS,
445 452 "Method only implemented for cicular and map regions.");
446 453 }
447   -
  454 +
448 455 // Return result
449 456 return inside;
450 457 }
... ... @@ -460,8 +467,8 @@ bool GSkyRegionMap::overlaps(const GSkyRegion&amp; reg) const
460 467 {
461 468 // Initialise return value
462 469 bool overlap = false;
463   -
464   - // Loop over non-zero pixels until direction is contained into region
  470 +
  471 + // Loop over non-zero pixels until direction is contained into region
465 472 for (int i = 0; i < m_nonzero_indices.size(); ++i) {
466 473 GSkyDir dir = m_map.inx2dir(m_nonzero_indices[i]);
467 474 if (reg.contains(dir)) {
... ... @@ -508,7 +515,7 @@ void GSkyRegionMap::copy_members(const GSkyRegionMap&amp; map)
508 515 // Copy members
509 516 m_map = map.m_map;
510 517 m_nonzero_indices = map.m_nonzero_indices;
511   -
  518 +
512 519 // Compute solid angle
513 520 compute_solid_angle();
514 521  
... ... @@ -539,7 +546,7 @@ void GSkyRegionMap::compute_solid_angle(void)
539 546 for (int i = 0; i < m_nonzero_indices.size(); ++i) {
540 547 m_solid += m_map.solidangle(m_nonzero_indices[i]);
541 548 }
542   -
  549 +
543 550 // Return
544 551 return;
545 552 }
... ... @@ -559,7 +566,7 @@ void GSkyRegionMap::set_nonzero_indices(void)
559 566 m_nonzero_indices.push_back(i);
560 567 }
561 568 }
562   -
  569 +
563 570 // Return
564 571 return;
565 572 }
... ... @@ -620,3 +627,71 @@ void GSkyRegionMap::set_region_circle(const GSkyRegionCircle* circle)
620 627 // Return
621 628 return;
622 629 }
  630 +
  631 +
  632 +/***********************************************************************//**
  633 + * @brief Set region map from region rectangle
  634 + *
  635 + * @param[in] rect Region rectangle
  636 + *
  637 + * Sets a region map from a region rectangle.
  638 + * A skymap is allocated with a resolution goal of 0.01 degrees (same for both
  639 + * axes).
  640 + * The particular resolution of either axes depends on the w/h of the rectangle.
  641 + * The sky map will be in celestial coordinates and TAN projection.
  642 + ***************************************************************************/
  643 +void GSkyRegionMap::set_region_rect(const GSkyRegionRect* rect)
  644 +{
  645 + // Define desired resolution (in degrees)
  646 + const double resolution_goal = 0.01;
  647 +
  648 + // Compute num of bins in x direction, round
  649 + int nbins_x = int(rect->width() / resolution_goal +0.5);
  650 + int nbins_y = int(rect->height() / resolution_goal +0.5);
  651 +
  652 + // Ensure a minimum num of bins
  653 + if (nbins_x < 10) {
  654 + nbins_x = 10;
  655 + }
  656 + if (nbins_y < 10) {
  657 + nbins_y = 10;
  658 + }
  659 +
  660 + // Compute actual resolution using nbins
  661 + double resolution_x = rect->width() / nbins_x;
  662 + double resolution_y = rect->height() / nbins_y;
  663 +
  664 + // Add leakage pixels
  665 + nbins_x += 6;
  666 + nbins_y += 6;
  667 +
  668 + // Create coarse 9-pixel map to later avoid "bending" of fine box-map
  669 + GSkyMap map_coarse("TAN", "CEL", rect->ra(), rect->dec(),
  670 + rect->width(), rect->height(),
  671 + 3, 3
  672 + );
  673 + map_coarse(4) = 1.0;
  674 + GSkyRegionMap regionmap_coarse = GSkyRegionMap(map_coarse);
  675 +
  676 + // Create fine sky map
  677 + GSkyMap map_fine("TAN", "CEL", rect->ra(), rect->dec(),
  678 + resolution_x, resolution_y,
  679 + nbins_x, nbins_y
  680 + );
  681 +
  682 + // Set sky map pixels
  683 + for (int i = 0; i < map_fine.npix(); ++i) {
  684 +
  685 + // If pixel is contained in rectangle then set map pixel to 1
  686 + if (regionmap_coarse.contains(map_fine.inx2dir(i))) {
  687 + map_fine(i) = 1.0;
  688 + }
  689 +
  690 + } // endfor: looped over sky map pixels
  691 +
  692 + // Set sky map as region map
  693 + this->map(map_fine);
  694 +
  695 + // Return
  696 + return;
  697 +}
623 698 \ No newline at end of file
... ...