Skip to content
Snippets Groups Projects
facet.cc 4.67 KiB
Newer Older
Maik Nijhuis's avatar
Maik Nijhuis committed
// Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later

#include "facet.h"

#include <aocommon/imagecoordinates.h>

#include <cassert>
#include <cmath>
#include <stdexcept>

namespace schaapcommon {
namespace facets {

BoundingBox::BoundingBox(const std::vector<Pixel>& pixels, size_t align) {
  if (pixels.empty()) {
    throw std::invalid_argument("Cannot create boundingbox for 0 pixels");
  }

  min_ = max_ = pixels.front();
  for (auto i = pixels.begin() + 1; i != pixels.end(); ++i) {
    min_.x = std::min(min_.x, i->x);
    max_.x = std::max(max_.x, i->x);
    min_.y = std::min(min_.y, i->y);
    max_.y = std::max(max_.y, i->y);
  }

  if (align > 1) {
    min_.x -= min_.x % align;
    min_.y -= min_.y % align;
    max_.x += align - 1 - (max_.x + align - 1) % align;
    max_.y += align - 1 - (max_.y + align - 1) % align;
  }
}

std::pair<int, int> Facet::HorizontalIntersections(
    const int y_intersect) const {
  if (y_intersect < min_y_ || y_intersect >= max_y_) {
    return std::make_pair<int, int>(0, 0);
  }

  bool found_x1 = false;
  int x1 = 0;
  int x2 = 0;
  size_t i;
  for (i = 0; i != pixels_.size(); ++i) {
    const Pixel& p1 = pixels_[i];
    const Pixel& p2 = pixels_[(i + 1) % pixels_.size()];
    if ((p1.y <= y_intersect && p2.y > y_intersect) ||
        (p2.y <= y_intersect && p1.y > y_intersect)) {
      int x;
      if (p1.y == y_intersect) {
        x = p1.x;
      } else if (p2.y == y_intersect) {
        x = p2.x;
      } else {
        const double beta = double(p2.x - p1.x) / double(p2.y - p1.y);
        const double xfl = p1.x + beta * (y_intersect - p1.y);
        x = round(xfl);
      }
      if (!found_x1) {
        x1 = x;
        found_x1 = true;
      } else {
        x2 = x;
        break;
      }
    }
  }

  // The loop should have found x1 and x2, and then stopped using 'break'.
  assert(i != pixels_.size());
  return std::minmax(x1, x2);
}

void Facet::CalculatePixels(double phase_centre_ra, double phase_centre_dec,
                            double px_scale_x, double px_scale_y,
                            size_t image_width, size_t image_height,
                            double shift_l, double shift_m,
                            bool origin_at_centre, double padding, size_t align,
                            bool make_square) {
  if (padding < 1.0) {
    throw std::invalid_argument("Padding factor should be >= 1.0");
  }
  if ((align > 1) && ((image_width % align) || (image_height % align))) {
    throw std::invalid_argument("Image is not aligned");
  }
  if (coords_.size() < 3) {
    throw std::runtime_error("Number of coordinates < 3, facet incomplete!");
  }

  pixels_.clear();
  pixels_.reserve(coords_.size());
  for (const Coord& coord : coords_) {
    double l, m;
    int x, y;
    aocommon::ImageCoordinates::RaDecToLM(coord.ra, coord.dec, phase_centre_ra,
                                          phase_centre_dec, l, m);
    l -= shift_l;
    m -= shift_m;
    aocommon::ImageCoordinates::LMToXY(l, m, px_scale_x, px_scale_y,
                                       image_width, image_height, x, y);

    // Clip x and y to the image size.
    x = std::max(x, 0);
    y = std::max(y, 0);
    x = std::min(x, static_cast<int>(image_width));
    y = std::min(y, static_cast<int>(image_height));

    if (origin_at_centre) {
      x = static_cast<int>(image_width / 2) - x;
      y = y - static_cast<int>(image_height / 2);
    }

    pixels_.emplace_back(x, y);
  }

  // Calculate bounding box for the pixels only, and set min_y_ and max_y_.
  BoundingBox pixel_box(pixels_);
  min_y_ = pixel_box.Min().y;
  max_y_ = pixel_box.Max().y;

  // Calculate the trimmed_box_.
  trimmed_box_ = BoundingBox({pixel_box.Min(), pixel_box.Max()}, align);
  if (make_square) {
    size_t trimmed_size = std::max(trimmed_box_.Width(), trimmed_box_.Height());
    const Pixel square_pad((trimmed_size - trimmed_box_.Width()) / 2,
                           (trimmed_size - trimmed_box_.Height()) / 2);
    trimmed_box_ = BoundingBox(
        {trimmed_box_.Min() - square_pad, trimmed_box_.Max() + square_pad});
  }

  // Calculate the untrimmed_box_.
  size_t width = static_cast<size_t>(std::ceil(trimmed_box_.Width() * padding));
  size_t height =
      static_cast<size_t>(std::ceil(trimmed_box_.Height() * padding));

  // Calculate padding. Divide by two since the padding occurs at both sides.
  const Pixel pad((width - trimmed_box_.Width()) / 2,
                  (height - trimmed_box_.Height()) / 2);

  // Create the padded and aligned bounding box for the facet.
  untrimmed_box_ =
      BoundingBox({trimmed_box_.Min() - pad, trimmed_box_.Max() + pad}, align);
}

}  // namespace facets
}  // namespace schaapcommon