SourceXtractorPlusPlus  0.11
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GrowthCurveTask.cpp
Go to the documentation of this file.
1 
27 #include "SEUtils/Mat22.h"
28 
29 
30 namespace SourceXtractor {
31 
32 using SExtractor::Mat22;
33 
34 static const SeFloat GROWTH_NSIG = 6.;
35 static const size_t GROWTH_NSAMPLES = 64;
36 
37 GrowthCurveTask::GrowthCurveTask(unsigned instance, bool use_symmetry)
38  : m_instance{instance}, m_use_symmetry{use_symmetry} {}
39 
41  // Measurement frame and properties defined on it
42  auto measurement_frame = source.getProperty<MeasurementFrame>(m_instance).getFrame();
43  auto image = measurement_frame->getSubtractedImage();
44  auto variance_map = measurement_frame->getVarianceMap();
45  auto variance_threshold = measurement_frame->getVarianceThreshold();
46  auto centroid_x = source.getProperty<MeasurementFramePixelCentroid>(m_instance).getCentroidX();
47  auto centroid_y = source.getProperty<MeasurementFramePixelCentroid>(m_instance).getCentroidY();
48  Mat22 jacobian{source.getProperty<JacobianSource>(m_instance).asTuple()};
49 
50  // These are from the detection frame
51  auto& shape = source.getProperty<ShapeParameters>();
52  auto& kron = source.getProperty<KronRadius>();
53 
54  // Radius for computing the growth curve and step size *on the detection frame*
55  double detection_rlim = std::max(GROWTH_NSIG * shape.getEllipseA(), kron.getKronRadius());
56 
57  // Now we need to compute the rlim for the measurement frame
58  // We take two vectors defined by the radius on the detection frame along the X and Y,
59  // transform them, and we use as a limit now the longest of the two after the transformation
60  Mat22 radius_22{detection_rlim, 0, 0, detection_rlim};
61  radius_22 = radius_22 * jacobian;
62  double r1 = radius_22[0] * radius_22[0] + radius_22[1] * radius_22[1];
63  double r2 = radius_22[2] * radius_22[2] + radius_22[3] * radius_22[3];
64  double rlim = std::sqrt(std::max(r1, r2));
65 
66  double step_size = rlim / GROWTH_NSAMPLES;
67 
68  // List of apertures
71  apertures.reserve(GROWTH_NSAMPLES);
72  for (size_t step = 1; step <= GROWTH_NSAMPLES; ++step) {
73  apertures.emplace_back(step_size * step);
74  }
75 
76  // Boundaries for the computation
77  // We know the last aperture is the widest, so we take the limits from it
78  auto min_coord = apertures.back().getMinPixel(centroid_x, centroid_y);
79  auto max_coord = apertures.back().getMaxPixel(centroid_x, centroid_y);
80 
81  // Compute fluxes for each ring
83  for (auto y = min_coord.m_y; y <= max_coord.m_y; ++y) {
84  for (auto x = min_coord.m_x; x <= max_coord.m_x; ++x) {
85  if (x >= 0 && x < image->getWidth() && y >= 0 && y < image->getHeight()) {
86  // Get the pixel value
87  WeightImage::PixelType variance_tmp = variance_map ? variance_map->getValue(x, y) : 1;
88  DetectionImage::PixelType pixel_value = 0;
89  // Masked out
90  if (variance_tmp > variance_threshold) {
91  if (m_use_symmetry) {
92  auto mirror_x = 2 * centroid_x - x + 0.49999;
93  auto mirror_y = 2 * centroid_y - y + 0.49999;
94  if (mirror_x >= 0 && mirror_y >= 0 && mirror_x < image->getWidth() && mirror_y < image->getHeight()) {
95  variance_tmp = variance_map ? variance_map->getValue(mirror_x, mirror_y) : 1;
96  if (variance_tmp < variance_threshold) {
97  // mirror pixel is OK: take the value
98  pixel_value = image->getValue(mirror_x, mirror_y);
99  }
100  }
101  }
102  }
103  // Not masked
104  else {
105  pixel_value = image->getValue(x, y);
106  }
107 
108  // Assign the pixel value according to the affected area
109  auto dx = x - centroid_x;
110  auto dy = y - centroid_y;
111  double r = std::sqrt(dx * dx + dy * dy);
112  // The pixel may be affected by multiple rings, so we look for those
113  // that overlap the start and the end of the pixels (which has size sqrt(2) on the diagonal)
114  size_t idx = 0;
115  if (r > 1.42 / 2.) {
116  idx = static_cast<size_t>((r - 1.42 / 2.) / step_size);
117  }
118  size_t outer_idx = static_cast<size_t>(std::ceil((r + 1.42 / 2.) / step_size));
119 
120  double inner = 0;
121  for (; idx <= outer_idx; ++idx) {
122  if (idx < apertures.size()) {
123  auto& aperture = apertures[idx];
124  auto area = aperture.getArea(centroid_x, centroid_y, x, y);
125 
126  fluxes[idx] += area * pixel_value - inner;
127  inner = area * pixel_value;
128  }
129  }
130  }
131  }
132  }
133 
134  // Accumulate
135  for (size_t i = 1; i < fluxes.size(); ++i) {
136  fluxes[i] += fluxes[i - 1];
137  }
138 
139  // Set property
140  source.setIndexedProperty<GrowthCurve>(m_instance, std::move(fluxes), rlim);
141 }
142 
143 } // end of namespace SourceXtractor
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
T ceil(T...args)
void computeProperties(SourceInterface &source) const override
Computes one or more properties for the Source.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
void setIndexedProperty(std::size_t index, Args...args)
Convenience template method to call setProperty() with a more user-friendly syntax.
SeFloat32 SeFloat
Definition: Types.h:32
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
GrowthCurveTask(unsigned instance, bool use_symmetry)
T lock(T...args)
static const std::string GROWTH_NSAMPLES
T max(T...args)
T move(T...args)
T size(T...args)
STL class.
static const SeFloat GROWTH_NSIG
T back(T...args)
T sqrt(T...args)
std::shared_ptr< EngineParameter > dx
The SourceInterface is an abstract &quot;source&quot; that has properties attached to it.
T reserve(T...args)
std::shared_ptr< EngineParameter > dy
T emplace_back(T...args)