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