ShapeParametersTask.cpp

Go to the documentation of this file.

/*
 * ShapeParametersTask.cpp
 *
 *  Created on: Jan 27, 2017
 *      Author: mschefer
 */

#include <iostream>

#include "SEImplementation/Property/PixelCoordinateList.h"
#include "SEImplementation/Plugin/DetectionFramePixelValues/DetectionFramePixelValues.h"
#include "SEImplementation/Plugin/ShapeParameters/ShapeParameters.h"
#include "SEImplementation/Plugin/PixelCentroid/PixelCentroid.h"
#include "SEImplementation/Plugin/PeakValue/PeakValue.h"

#include "SEImplementation/Plugin/ShapeParameters/ShapeParametersTask.h"

namespace SourceXtractor {

void ShapeParametersTask::computeProperties(SourceInterface& source) const {
  const auto& pixel_values = source.getProperty<DetectionFramePixelValues>().getFilteredValues();
  const auto& centroid_x = source.getProperty<PixelCentroid>().getCentroidX();
  const auto& centroid_y = source.getProperty<PixelCentroid>().getCentroidY();
  auto min_value = source.getProperty<PeakValue>().getMinValue();
  auto peak_value = source.getProperty<PeakValue>().getMaxValue();
  auto& coordinates = source.getProperty<PixelCoordinateList>().getCoordinateList();

  SeFloat x_2 = 0.0;
  SeFloat y_2 = 0.0;
  SeFloat x_y = 0.0;

  SeFloat total_intensity = 0;

  DetectionImage::PixelType half_peak_threshold = (peak_value + min_value) / 2.0;
  int nb_of_pixels_above_half = 0;
  int nb_of_pixels = coordinates.size();

  auto i = pixel_values.begin();
  for (auto pixel_coord : coordinates) {
        SeFloat value = *i++;
        SeFloat x_pos = pixel_coord.m_x - centroid_x;
        SeFloat y_pos = pixel_coord.m_y - centroid_y;

    if (value > half_peak_threshold) {
      nb_of_pixels_above_half++;
    }

    x_2 += x_pos * x_pos * value;
    y_2 += y_pos * y_pos * value;
    x_y += x_pos * y_pos * value;

    total_intensity += value;
  }

  x_2 /= total_intensity;
  y_2 /= total_intensity;
  x_y /= total_intensity;

  SeFloat theta;
  if (fabs(x_2 - y_2) > 0.0) {
    theta = atan2(2.0 * x_y, x_2 - y_2) / 2.0;
  } else {
    theta = M_PI/4.0;
  }

  float temp = x_2 - y_2;
  temp = sqrt(0.25 * temp * temp + x_y * x_y);
  SeFloat a = sqrt((x_2 + y_2) / 2 + temp);
  SeFloat b = ((x_2 + y_2) / 2) > temp ? sqrt((x_2 + y_2) / 2 - temp) : 0;

  // From original SExtractor: Handle fully correlated x/y (which cause a singularity...)
  SeFloat tmp = x_2 * y_2 - x_y * x_y;
  bool singu=false;
  if (tmp < 0.00694) {
    x_2 += 0.0833333;
    y_2 += 0.0833333;
    tmp = x_2 * y_2 - x_y * x_y;
    singu=true;
  }

  SeFloat cxx = y_2 / tmp;
  SeFloat cyy = x_2 / tmp;
  SeFloat cxy = -2.0 * x_y / tmp;

  SeFloat area_under_half = (double) pixel_values.size() - nb_of_pixels_above_half;
  SeFloat t1t2 = min_value / half_peak_threshold;

  SeFloat abcor = 1.0;
  if (t1t2 > 0.0) { // && !prefs.dweight_flag
    abcor = (area_under_half > 0.0 ? area_under_half : 1.0) / (2 * M_PI * -log(t1t2 < 1.0 ? t1t2 : 0.99) * a * b);
    if (abcor > 1.0) {
      abcor = 1.0;
    }
  }

  // set the object property
  source.setProperty<ShapeParameters>(a, b, theta, abcor, cxx, cyy, cxy, nb_of_pixels, total_intensity, singu);
}


}