SourceXtractorPlusPlus  0.16
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 
36 
38 
40 
41 namespace SourceXtractor {
42 
43 class MultiThresholdNode : public std::enable_shared_from_this<MultiThresholdNode> {
44 public:
45 
47  : m_pixel_list(pixel_list), m_is_split(false), m_threshold(threshold) {
48  }
49 
51  m_children.push_back(child);
52  child->m_parent = shared_from_this();
53  }
54 
55  bool contains(const Lutz::PixelGroup& pixel_group) const {
56  for (auto pixel : m_pixel_list) {
57  if (pixel_group.pixel_list[0] == pixel) {
58  return true;
59  }
60  }
61  return false;
62  }
63 
65  return m_children;
66  }
67 
69  return m_parent.lock();
70  }
71 
73  DetectionImage::PixelType total_intensity = 0;
74  for (const auto& pixel_coord : m_pixel_list) {
75  total_intensity += (image.getValue(pixel_coord - offset) - m_threshold);
76  }
77 
78  return total_intensity;
79  }
80 
81  bool isSplit() const {
82  return m_is_split;
83  }
84 
85  void flagAsSplit() {
86  m_is_split = true;
87  auto parent = m_parent.lock();
88  if (parent != nullptr) {
89  parent->flagAsSplit();
90  }
91  }
92 
94  return m_pixel_list;
95  }
96 
97  void debugPrint() const {
98  std::cout << "(" << m_pixel_list.size();
99 
100  for (auto& child : m_children) {
101  std::cout << ", ";
102  child->debugPrint();
103  }
104 
105  std::cout << ")";
106  }
107 
108  void addPixel(PixelCoordinate pixel) {
109  m_pixel_list.emplace_back(pixel);
110  }
111 
113  return m_threshold;
114  }
115 
116 private:
118 
121 
123 
125 };
126 
128  std::shared_ptr<SourceInterface> original_source) const {
129 
130  auto parent_source_id = original_source->getProperty<SourceId>().getSourceId();
131 
132  auto& detection_frame = original_source->getProperty<DetectionFrame>();
133 
134  auto& detection_frame_images = original_source->getProperty<DetectionFrameImages>();
135  const auto labelling_image = detection_frame_images.getLockedImage(LayerFilteredImage);
136 
137  auto& pixel_boundaries = original_source->getProperty<PixelBoundaries>();
138 
139  auto& pixel_coords = original_source->getProperty<PixelCoordinateList>().getCoordinateList();
140 
141  auto offset = pixel_boundaries.getMin();
142  auto thumbnail_image = VectorImage<DetectionImage::PixelType>::create(
143  pixel_boundaries.getWidth(), pixel_boundaries.getHeight());
144  thumbnail_image->fillValue(0);
145 
146  auto min_value = original_source->getProperty<PeakValue>().getMinValue() * .8;
147  auto peak_value = original_source->getProperty<PeakValue>().getMaxValue();
148 
149  for (auto pixel_coord : pixel_coords) {
150  auto value = labelling_image->getValue(pixel_coord);
151  thumbnail_image->setValue(pixel_coord - offset, value);
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.getEncapsulatedFrame());
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.getEncapsulatedFrame());
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 
std::shared_ptr< EngineParameter > dx
std::shared_ptr< EngineParameter > dy
T back(T... args)
const std::vector< PixelGroup > & getGroups() const
Definition: Lutz.h:71
void labelImage(const DetectionImage &image, PixelCoordinate offset=PixelCoordinate(0, 0))
Definition: Lutz.h:75
std::vector< PixelCoordinate > pixel_list
Definition: Lutz.h:44
bool contains(const Lutz::PixelGroup &pixel_group) const
std::shared_ptr< MultiThresholdNode > getParent() const
std::weak_ptr< MultiThresholdNode > m_parent
void addChild(std::shared_ptr< MultiThresholdNode > child)
std::vector< std::shared_ptr< MultiThresholdNode > > m_children
MultiThresholdNode(const std::vector< PixelCoordinate > &pixel_list, SeFloat threshold)
const std::vector< PixelCoordinate > & getPixels() const
double getTotalIntensity(VectorImage< DetectionImage::PixelType > &image, const PixelCoordinate &offset) const
const std::vector< std::shared_ptr< MultiThresholdNode > > & getChildren() const
virtual std::vector< std::shared_ptr< SourceInterface > > partition(std::shared_ptr< SourceInterface > source) const
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
The bounding box of all the pixels in the source. Both min and max coordinate are inclusive.
The centroid of all the pixels in the source, weighted by their DetectionImage pixel values.
Definition: PixelCentroid.h:37
static std::shared_ptr< ProcessedImage< T, P > > create(std::shared_ptr< const Image< T >> image_a, std::shared_ptr< const Image< T >> image_b)
Image implementation which keeps the pixel values in memory.
Definition: VectorImage.h:52
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:100
T getValue(int x, int y) const
Definition: VectorImage.h:116
T empty(T... args)
T endl(T... args)
T max(T... args)
SeFloat32 SeFloat
Definition: Types.h:32
@ LayerFilteredImage
Definition: Frame.h:39
T pop_back(T... args)
T pow(T... args)
T push_back(T... args)
T rand(T... args)
T size(T... args)
A pixel coordinate made of two integers m_x and m_y.