SourceXtractorPlusPlus  0.11
Please provide a description of the project.
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
MeasurementFrame.h
std::lock
T lock(T... args)
Jacobian.h
std::move
T move(T... args)
SourceXtractor::Image< SeFloat >::PixelType
SeFloat PixelType
Definition: Image.h:47
MeasurementFramePixelCentroid.h
SourceXtractor::MeasurementFramePixelCentroid
Definition: MeasurementFramePixelCentroid.h:31
std::vector::reserve
T reserve(T... args)
SourceXtractor::SeFloat
SeFloat32 SeFloat
Definition: Types.h:32
std::vector
STL class.
std::vector::size
T size(T... args)
SourceXtractor::KronRadius
Kron radius.
Definition: KronRadius.h:36
std::lock_guard
STL class.
SourceXtractor::GrowthCurveTask::computeProperties
void computeProperties(SourceInterface &source) const override
Computes one or more properties for the Source.
Definition: GrowthCurveTask.cpp:40
Mat22.h
std::vector::back
T back(T... args)
SourceXtractor::MultithreadedMeasurement::g_global_mutex
static std::recursive_mutex g_global_mutex
Definition: MultithreadedMeasurement.h:54
std::sqrt
T sqrt(T... args)
CircularAperture.h
SourceXtractor::GrowthCurveTask::m_instance
unsigned m_instance
Definition: GrowthCurveTask.h:34
SourceXtractor
Definition: Aperture.h:30
dy
std::shared_ptr< EngineParameter > dy
Definition: MoffatModelFittingTask.cpp:92
SourceXtractor::GROWTH_NSIG
static const SeFloat GROWTH_NSIG
Definition: GrowthCurveTask.cpp:34
SourceXtractor::SourceInterface::setIndexedProperty
void setIndexedProperty(std::size_t index, Args... args)
Convenience template method to call setProperty() with a more user-friendly syntax.
Definition: SourceInterface.h:64
SourceXtractor::GROWTH_NSAMPLES
static const std::string GROWTH_NSAMPLES
Definition: GrowthCurveConfig.cpp:26
SourceXtractor::MeasurementFrame
Definition: MeasurementFrame.h:36
std::ceil
T ceil(T... args)
SourceXtractor::JacobianSource
Definition: Jacobian.h:51
SourceXtractor::GrowthCurveTask::GrowthCurveTask
GrowthCurveTask(unsigned instance, bool use_symmetry)
Definition: GrowthCurveTask.cpp:37
std::vector::emplace_back
T emplace_back(T... args)
GrowthCurve.h
SourceXtractor::GrowthCurve
Definition: GrowthCurve.h:30
x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
Definition: MoffatModelFittingTask.cpp:93
MultithreadedMeasurement.h
GrowthCurveTask.h
SourceXtractor::SourceInterface::getProperty
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
Definition: SourceInterface.h:57
ShapeParameters.h
SourceXtractor::SourceInterface
The SourceInterface is an abstract "source" that has properties attached to it.
Definition: SourceInterface.h:46
std::max
T max(T... args)
dx
std::shared_ptr< EngineParameter > dx
Definition: MoffatModelFittingTask.cpp:92
SourceXtractor::GrowthCurveTask::m_use_symmetry
bool m_use_symmetry
Definition: GrowthCurveTask.h:35
y
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
Definition: MoffatModelFittingTask.cpp:93
SourceXtractor::ShapeParameters
Definition: ShapeParameters.h:32
SExtractor::Mat22
Definition: Mat22.h:19
KronRadius.h