SourceXtractorPlusPlus  0.11
Please provide a description of the project.
MultiThresholdPartitionStep.cpp
Go to the documentation of this file.
1 
17 /*
18  * MultiThresholdPartitionStep.cpp
19  *
20  * Created on: Jan 17, 2017
21  * Author: mschefer
22  */
23 
26 
29 
35 
37 
39 
40 namespace SourceXtractor {
41 
42 class MultiThresholdNode : public std::enable_shared_from_this<MultiThresholdNode> {
43 public:
44 
46  : m_pixel_list(pixel_list), m_is_split(false), m_threshold(threshold) {
47  }
48 
50  m_children.push_back(child);
51  child->m_parent = shared_from_this();
52  }
53 
54  bool contains(const Lutz::PixelGroup& pixel_group) const {
55  for (auto pixel : m_pixel_list) {
56  if (pixel_group.pixel_list[0] == pixel) {
57  return true;
58  }
59  }
60  return false;
61  }
62 
64  return m_children;
65  }
66 
68  return m_parent.lock();
69  }
70 
71  double getTotalIntensity(DetectionImage& image, const PixelCoordinate& offset) const {
72  DetectionImage::PixelType total_intensity = 0;
73  for (const auto& pixel_coord : m_pixel_list) {
74  total_intensity += (image.getValue(pixel_coord - offset) - m_threshold);
75  }
76 
77  return total_intensity;
78  }
79 
80  bool isSplit() const {
81  return m_is_split;
82  }
83 
84  void flagAsSplit() {
85  m_is_split = true;
86  auto parent = m_parent.lock();
87  if (parent != nullptr) {
88  parent->flagAsSplit();
89  }
90  }
91 
93  return m_pixel_list;
94  }
95 
96  void debugPrint() const {
97  std::cout << "(" << m_pixel_list.size();
98 
99  for (auto& child : m_children) {
100  std::cout << ", ";
101  child->debugPrint();
102  }
103 
104  std::cout << ")";
105  }
106 
107  void addPixel(PixelCoordinate pixel) {
108  m_pixel_list.emplace_back(pixel);
109  }
110 
112  return m_threshold;
113  }
114 
115 private:
117 
120 
122 
124 };
125 
127  std::shared_ptr<SourceInterface> original_source) const {
128 
129  auto parent_source_id = original_source->getProperty<SourceId>().getSourceId();
130 
131  auto& detection_frame = original_source->getProperty<DetectionFrame>();
132  const auto labelling_image = detection_frame.getFrame()->getFilteredImage();
133 
134  auto& pixel_boundaries = original_source->getProperty<PixelBoundaries>();
135 
136  auto& pixel_coords = original_source->getProperty<PixelCoordinateList>().getCoordinateList();
137 
138  auto offset = pixel_boundaries.getMin();
139  auto thumbnail_image = VectorImage<DetectionImage::PixelType>::create(
140  pixel_boundaries.getWidth(), pixel_boundaries.getHeight());
141  thumbnail_image->fillValue(0);
142 
143  auto min_value = original_source->getProperty<PeakValue>().getMinValue() * .8;
144  auto peak_value = original_source->getProperty<PeakValue>().getMaxValue();
145 
146  {
148  for (auto pixel_coord : pixel_coords) {
149  auto value = labelling_image->getValue(pixel_coord);
150  thumbnail_image->setValue(pixel_coord - offset, value);
151  }
152  }
153 
154  auto root = std::make_shared<MultiThresholdNode>(pixel_coords, 0);
155 
156  std::list<std::shared_ptr<MultiThresholdNode>> active_nodes { root };
158 
159  // Build the tree
160  for (unsigned int i = 1; i < m_thresholds_nb; i++) {
161 
162  auto threshold = min_value * pow(peak_value / min_value, (double) i / m_thresholds_nb);
163  auto subtracted_image = SubtractImage<DetectionImage::PixelType>::create(thumbnail_image, threshold);
164 
165  LutzList lutz;
166  lutz.labelImage(*subtracted_image, offset);
167 
168  std::list<std::shared_ptr<MultiThresholdNode>> active_nodes_copy(active_nodes);
169  for (auto& node : active_nodes_copy) {
170  int nb_of_groups_inside = 0;
171  for (auto& pixel_group : lutz.getGroups()) {
172  if (pixel_group.pixel_list.size() >= m_min_deblend_area && node->contains(pixel_group)) {
173  nb_of_groups_inside++;
174  }
175  }
176 
177  if (nb_of_groups_inside == 0) {
178  active_nodes.remove(node);
179  }
180 
181  if (nb_of_groups_inside > 1) {
182  active_nodes.remove(node);
183  junction_nodes.push_back(node);
184  for (auto& pixel_group : lutz.getGroups()) {
185  if (pixel_group.pixel_list.size() >= m_min_deblend_area && node->contains(pixel_group)) {
186  auto new_node = std::make_shared<MultiThresholdNode>(pixel_group.pixel_list, threshold);
187  node->addChild(new_node);
188  active_nodes.push_back(new_node);
189  }
190  }
191  }
192  }
193  }
194 
195  // Identify the sources
196  double intensity_threshold = root->getTotalIntensity(*thumbnail_image, offset) * m_contrast;
197 
199  while (!junction_nodes.empty()) {
200  auto node = junction_nodes.back();
201  junction_nodes.pop_back();
202 
203  int nb_of_children_above_threshold = 0;
204 
205  for (auto child : node->getChildren()) {
206  if (child->getTotalIntensity(*thumbnail_image, offset) > intensity_threshold) {
207  nb_of_children_above_threshold++;
208  }
209  }
210 
211  if (nb_of_children_above_threshold >= 2) {
212  node->flagAsSplit();
213  for (auto child : node->getChildren()) {
214  if (child->getTotalIntensity(*thumbnail_image, offset) > intensity_threshold && !child->isSplit()) {
215  source_nodes.push_back(child);
216  }
217  }
218  }
219  }
220 
222  if (source_nodes.empty()) {
223  return { original_source }; // no split, just forward the source unchanged
224  }
225 
226  for (auto source_node : source_nodes) {
227  // remove pixels in the new sources from the image
228  for (auto& pixel : source_node->getPixels()) {
229  thumbnail_image->setValue(pixel - offset, 0);
230  }
231 
232  auto new_source = m_source_factory->createSource();
233 
234  new_source->setProperty<PixelCoordinateList>(source_node->getPixels());
235  new_source->setProperty<DetectionFrame>(detection_frame.getFrame());
236 
237  sources.push_back(new_source);
238  }
239 
240  auto new_sources = reassignPixels(sources, pixel_coords, thumbnail_image, source_nodes, offset);
241 
242  for (auto& new_source : new_sources) {
243  new_source->setProperty<DetectionFrame>(detection_frame.getFrame());
244  new_source->setProperty<SourceId>(parent_source_id);
245  }
246 
247  return new_sources;
248 }
249 
252  const std::vector<PixelCoordinate>& pixel_coords,
255  const PixelCoordinate& offset
256  ) const {
257 
258  std::vector<SeFloat> amplitudes;
259  for (auto& source : sources) {
260  auto& pixel_list = source->getProperty<PixelCoordinateList>().getCoordinateList();
261  auto& shape_parameters = source->getProperty<ShapeParameters>();
262 
263  auto thresh = source->getProperty<PeakValue>().getMinValue();
264  auto peak = source->getProperty<PeakValue>().getMaxValue();
265 
266  auto dist = pixel_list.size() / (2.0 * M_PI * shape_parameters.getAbcor() * shape_parameters.getEllipseA() * shape_parameters.getEllipseB());
267  auto amp = dist < 70.0 ? thresh * expf(dist) : 4.0 * peak;
268 
269  // limit expansion ??
270  if (amp > 4.0 * peak) {
271  amp = 4.0 * peak;
272  }
273 
274  amplitudes.push_back(amp);
275  }
276 
277  for (auto pixel : pixel_coords) {
278  if (image->getValue(pixel - offset) > 0) {
279  SeFloat cumulated_probability = 0;
280  std::vector<SeFloat> probabilities;
281 
283  std::shared_ptr<MultiThresholdNode> closest_source_node;
284 
285  int i = 0;
286  for (auto& source : sources) {
287  auto& shape_parameters = source->getProperty<ShapeParameters>();
288  auto& pixel_centroid = source->getProperty<PixelCentroid>();
289 
290  auto dx = pixel.m_x - pixel_centroid.getCentroidX();
291  auto dy = pixel.m_y - pixel_centroid.getCentroidY();
292 
293  auto dist = 0.5 * (shape_parameters.getEllipseCxx()*dx*dx +
294  shape_parameters.getEllipseCyy()*dy*dy + shape_parameters.getEllipseCxy()*dx*dy) /
295  shape_parameters.getAbcor();
296 
297  if (dist < min_dist) {
298  min_dist = dist;
299  closest_source_node = source_nodes[i];
300  }
301 
302  cumulated_probability += dist < 70.0 ? amplitudes[i] * expf(-dist) : 0.0;
303 
304  probabilities.push_back(cumulated_probability);
305  i++;
306  }
307 
308  if (probabilities.back() > 1.0e-31) {
309  // TODO probably should use a better RNG
310  auto drand = double(probabilities.back()) * double(rand()) / RAND_MAX;
311 
312  unsigned int i=0;
313  for (; i<probabilities.size() && drand >= probabilities[i]; i++);
314  if (i < source_nodes.size()) {
315  source_nodes[i]->addPixel(pixel);
316  } else {
317  std::cout << i << " oops " << drand << " " << probabilities.back() << std::endl;
318  }
319 
320  } else {
321  // select closest source
322  closest_source_node->addPixel(pixel);
323  }
324  }
325  }
326 
327  int total_pixels = 0;
328 
330  for (auto source_node : source_nodes) {
331  // remove pixels in the new sources from the image
332  for (auto& pixel : source_node->getPixels()) {
333  image->setValue(pixel - offset, 0);
334  }
335 
336  auto new_source = m_source_factory->createSource();
337 
338  auto pixels = source_node->getPixels();
339  total_pixels += pixels.size();
340 
341  new_source->setProperty<PixelCoordinateList>(pixels);
342 
343  new_sources.push_back(new_source);
344  }
345 
346  return new_sources;
347 }
348 
349 
350 } // namespace
351 
352 
SourceXtractor::PixelBoundaries
The bounding box of all the pixels in the source. Both min and max coordinate are inclusive.
Definition: PixelBoundaries.h:37
SourceXtractor::PixelCoordinateList
Definition: PixelCoordinateList.h:31
SourceXtractor::MultiThresholdNode::debugPrint
void debugPrint() const
Definition: MultiThresholdPartitionStep.cpp:96
std::lock
T lock(T... args)
SourceXtractor::PixelCoordinate
A pixel coordinate made of two integers m_x and m_y.
Definition: PixelCoordinate.h:37
SourceXtractor::MultiThresholdNode::getPixels
const std::vector< PixelCoordinate > & getPixels() const
Definition: MultiThresholdPartitionStep.cpp:92
PixelBoundaries.h
std::shared_ptr
STL class.
std::list
STL class.
SourceXtractor::Image< SeFloat >::PixelType
SeFloat PixelType
Definition: Image.h:47
SourceId.h
SourceXtractor::PixelCentroid
The centroid of all the pixels in the source, weighted by their DetectionImage pixel values.
Definition: PixelCentroid.h:37
PeakValue.h
SourceXtractor::SeFloat
SeFloat32 SeFloat
Definition: Types.h:32
std::vector
STL class.
std::vector::size
T size(T... args)
SourceXtractor::LutzList::labelImage
void labelImage(const DetectionImage &image, PixelCoordinate offset=PixelCoordinate(0, 0))
Definition: Lutz.h:75
SourceXtractor::MultiThresholdNode::getParent
std::shared_ptr< MultiThresholdNode > getParent() const
Definition: MultiThresholdPartitionStep.cpp:67
std::lock_guard
STL class.
std::list::back
T back(T... args)
SourceXtractor::MultiThresholdNode::addPixel
void addPixel(PixelCoordinate pixel)
Definition: MultiThresholdPartitionStep.cpp:107
SourceXtractor::MultithreadedMeasurement::g_global_mutex
static std::recursive_mutex g_global_mutex
Definition: MultithreadedMeasurement.h:54
SourceXtractor::MultiThresholdPartitionStep::m_contrast
SeFloat m_contrast
Definition: MultiThresholdPartitionStep.h:68
SourceXtractor::Image< SeFloat >
VectorImage.h
DetectionFramePixelValues.h
MultiThresholdPartitionStep.h
std::list::push_back
T push_back(T... args)
SourceXtractor::SourceId
Definition: SourceId.h:31
SourceXtractor
Definition: Aperture.h:30
SourceXtractor::DetectionFrame::getFrame
std::shared_ptr< DetectionImageFrame > getFrame() const
Definition: DetectionFrame.h:38
std::enable_shared_from_this< MultiThresholdNode >::shared_from_this
T shared_from_this(T... args)
SourceXtractor::Image::getValue
virtual T getValue(int x, int y) const =0
Returns the value of the pixel with the coordinates (x,y)
SourceXtractor::MultiThresholdPartitionStep::m_thresholds_nb
unsigned int m_thresholds_nb
Definition: MultiThresholdPartitionStep.h:69
SourceXtractor::MultiThresholdNode
Definition: MultiThresholdPartitionStep.cpp:42
dy
std::shared_ptr< EngineParameter > dy
Definition: MoffatModelFittingTask.cpp:92
std::cout
SourceXtractor::MultiThresholdPartitionStep::m_source_factory
std::shared_ptr< SourceFactory > m_source_factory
Definition: MultiThresholdPartitionStep.h:67
SourceXtractor::LutzList
Definition: Lutz.h:65
SourceXtractor::MultiThresholdPartitionStep::partition
virtual std::vector< std::shared_ptr< SourceInterface > > partition(std::shared_ptr< SourceInterface > source) const
Definition: MultiThresholdPartitionStep.cpp:126
std::enable_shared_from_this
SourceXtractor::DetectionFrame
Definition: DetectionFrame.h:33
SourceXtractor::VectorImage::create
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:89
SourceXtractor::MultiThresholdNode::m_parent
std::weak_ptr< MultiThresholdNode > m_parent
Definition: MultiThresholdPartitionStep.cpp:118
SourceXtractor::MultiThresholdNode::isSplit
bool isSplit() const
Definition: MultiThresholdPartitionStep.cpp:80
ProcessedImage.h
SourceXtractor::MultiThresholdNode::getChildren
const std::vector< std::shared_ptr< MultiThresholdNode > > & getChildren() const
Definition: MultiThresholdPartitionStep.cpp:63
SourceXtractor::PeakValue
Definition: PeakValue.h:32
std::list::pop_back
T pop_back(T... args)
Lutz.h
SourceXtractor::MultiThresholdNode::m_is_split
bool m_is_split
Definition: MultiThresholdPartitionStep.cpp:121
SourceXtractor::MultiThresholdNode::m_children
std::vector< std::shared_ptr< MultiThresholdNode > > m_children
Definition: MultiThresholdPartitionStep.cpp:119
SourceXtractor::Lutz::PixelGroup
Definition: Lutz.h:39
std::rand
T rand(T... args)
std::weak_ptr
STL class.
SourceXtractor::ProcessedImage::create
static std::shared_ptr< ProcessedImage< T, P > > create(std::shared_ptr< const Image< T >> image_a, std::shared_ptr< const Image< T >> image_b)
Definition: ProcessedImage.h:54
SourceXtractor::LutzList::getGroups
const std::vector< PixelGroup > & getGroups() const
Definition: Lutz.h:71
std::endl
T endl(T... args)
SourceXtractor::MultiThresholdNode::MultiThresholdNode
MultiThresholdNode(const std::vector< PixelCoordinate > &pixel_list, SeFloat threshold)
Definition: MultiThresholdPartitionStep.cpp:45
SourceXtractor::VectorImage< DetectionImage::PixelType >
MultithreadedMeasurement.h
SourceXtractor::MultiThresholdNode::getTotalIntensity
double getTotalIntensity(DetectionImage &image, const PixelCoordinate &offset) const
Definition: MultiThresholdPartitionStep.cpp:71
std::list::empty
T empty(T... args)
SourceXtractor::MultiThresholdNode::addChild
void addChild(std::shared_ptr< MultiThresholdNode > child)
Definition: MultiThresholdPartitionStep.cpp:49
SourceXtractor::MultiThresholdNode::contains
bool contains(const Lutz::PixelGroup &pixel_group) const
Definition: MultiThresholdPartitionStep.cpp:54
ShapeParameters.h
std::numeric_limits::max
T max(T... args)
dx
std::shared_ptr< EngineParameter > dx
Definition: MoffatModelFittingTask.cpp:92
SourceXtractor::MultiThresholdPartitionStep::reassignPixels
std::vector< std::shared_ptr< SourceInterface > > reassignPixels(const std::vector< std::shared_ptr< SourceInterface >> &sources, const std::vector< PixelCoordinate > &pixel_coords, std::shared_ptr< VectorImage< DetectionImage::PixelType >> image, const std::vector< std::shared_ptr< MultiThresholdNode >> &source_nodes, const PixelCoordinate &offset) const
Definition: MultiThresholdPartitionStep.cpp:250
SourceXtractor::MultiThresholdNode::m_pixel_list
std::vector< PixelCoordinate > m_pixel_list
Definition: MultiThresholdPartitionStep.cpp:116
SourceXtractor::ShapeParameters::getAbcor
SeFloat getAbcor() const
Definition: ShapeParameters.h:91
SourceXtractor::Lutz::PixelGroup::pixel_list
std::vector< PixelCoordinate > pixel_list
Definition: Lutz.h:44
SourceXtractor::ShapeParameters
Definition: ShapeParameters.h:32
PixelCentroid.h
SourceXtractor::MultiThresholdNode::getThreshold
SeFloat getThreshold() const
Definition: MultiThresholdPartitionStep.cpp:111
SourceXtractor::MultiThresholdNode::flagAsSplit
void flagAsSplit()
Definition: MultiThresholdPartitionStep.cpp:84
SourceXtractor::MultiThresholdNode::m_threshold
SeFloat m_threshold
Definition: MultiThresholdPartitionStep.cpp:123
SourceXtractor::MultiThresholdPartitionStep::m_min_deblend_area
unsigned int m_min_deblend_area
Definition: MultiThresholdPartitionStep.h:70
std::pow
T pow(T... args)