ecg-prog-filtered/u02/src/util.cpp

// SPDX-License-Identifier: LGPL-3.0-or-later
#include "util.h"
#include <cmath>

void transform_mut(Transformation transformation, int &x, int &y) {
  if (transformation & TRANSFORM_ROTATE_CW) {
    std::swap(x, y);
    x = -x;
  }
  if (transformation & TRANSFORM_ROTATE_CCW) {
    std::swap(x, y);
    y = -y;
  }
  if (transformation & TRANSFORM_MIRROR_X) {
    y = -y;
  }
  if (transformation & TRANSFORM_MIRROR_Y) {
    x = -x;
  }
}

std::pair<int, int> transform(Transformation transformation, int x, int y) {
  transform_mut(transformation, x, y);
  return std::make_pair(x, y);
}

void transform_inv_mut(Transformation transformation, int &x, int &y) {
  if (transformation & TRANSFORM_MIRROR_Y) {
    x = -x;
  }
  if (transformation & TRANSFORM_MIRROR_X) {
    y = -y;
  }
  if (transformation & TRANSFORM_ROTATE_CCW) {
    // does clockwise rotation
    std::swap(x, y);
    x = -x;
  }
  if (transformation & TRANSFORM_ROTATE_CW) {
    // does counterclockwise rotation
    std::swap(x, y);
    y = -y;
  }
}

std::pair<int, int> transform_inv(Transformation transformation, int x, int y) {
  transform_inv_mut(transformation, x, y);
  return std::make_pair(x, y);
}

/*
 * After it took me many hours to get this right,
 * I at least want to document how I got to it:
 *
 * N.B. In the following,
 * Modulo is *not* the remainder of the euclidean division,
 * but instead the remainder of truncated division
 * (i.e., negative quotients produce negative results).
 *
 * There are two main cases:
 * The simple one is where angle % 90° ≤ 45°.
 * To transform this into the standard case,
 * only mirrors are needed.
 * The more complicated is when angle % 90° ≥ 45°.
 * To transform this into the standard case,
 * a rotation has to be done, followed by a mirror in some cases.
 *
 * The following matrices show what must be done when.
 * The ASCII art arrows show the line as it should be drawn,
 * the the column/row headings show how they can be identified in code,
 * the capital letters in the field show what needs to be done
 * to reach the standard case
 * (X/Y: mirror X/Y; CW/CCW: rotate CW/CCW).
 * Because there is no nice way to draw arrows with an angle < 45°,
 * they are just differentiated by the heading.
 *
 * Let (x_0, y_0) be the starting point and (x_1, y_1) the end point.
 * Let Δx = x_1 - x_0, Δy = y_1 - y_0.
 * Let m = Δy/Δx, α = atan(m).
 *
 * α ≤ 45°:
 *
 *      Δx>0   Δx<0
 *
 *         A | A
 * Δy<0   /  |  \
 *       /   |   \
 *      /  - | Y  \
 *     ------+------
 *      \  X | XY /
 *       \   |   /
 * Δy>0   \  |  /
 *         V | V
 *
 * α ≥ 45°:
 *
 *      Δx>0   Δx<0                       Δx>0   Δx<0
 *
 *         A | A                          \    |    A
 * Δy<0   /  |  \                    Δy<0  \   |   /
 *       /   |   \                          \  |  /
 *      / CW | CW \                       X  V | /
 *    -------+------- → after rotation → ------+-----
 *     \ CCW | CCW /                         A | \
 *      \    |    /                         /  |  \
 * Δy>0  \   |   /                   Δy>0  /   |   \
 *        V  |  V                         /    | X  V
 */
Transformation transformation_to_standard_case(int x0, int y0, int x1, int y1) {
  Transformation transformation = 0;

  int delta_y = y1 - y0;
  int delta_x = x1 - x0;
  // checks if angle ∈ (-90°, 90°) is ≥ 45°
  // this is a simplified version of atan(Δy/Δx) > π/4:
  //   atan(Δy/Δx) > π/4 | tan(…)
  // ⇔ Δx/Δy       > 1 | Δy
  // ⇔ Δx          > Δy
  if (abs(delta_y) > abs(delta_x)) {
    // if-else is needed, because of special case Δy = 0
    if (delta_y < 0) {
      transformation |= TRANSFORM_ROTATE_CW;
    } else if (delta_y > 0) {
      transformation |= TRANSFORM_ROTATE_CCW;
    }
    // the sign of Δx and Δy (pre-rotation!) differ,
    // an additional mirror is needed
    if (delta_x * delta_y < 0) {
      transformation |= TRANSFORM_MIRROR_X;
    }
  } else {
    if (delta_x < 0) {
      transformation |= TRANSFORM_MIRROR_Y;
    }
    if (delta_y > 0) {
      transformation |= TRANSFORM_MIRROR_X;
    }
  }

  return transformation;
}

std::tuple<float, float, float> barycentric_coordinates(int x0, int y0, int x1,
                                                        int y1, int x2, int y2,
                                                        int xp, int yp) {
  // Source:
  // https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Vertex_approach
  float b1 = x1 * y2 - x2 * y1 + xp * (y1 - y2) + yp * (x2 - x1);
  float b2 = x2 * y0 - x0 * y2 + xp * (y2 - y0) + yp * (x0 - x2);
  float b3 = x0 * y1 - x1 * y0 + xp * (y0 - y1) + yp * (x1 - x0);

  // reciprocal computed directly for performance
  float area_factor =
      1 / static_cast<float>(x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1));

  b1 *= area_factor;
  b2 *= area_factor;
  b3 *= area_factor;

  return {b1, b2, b3};
}

bool point_in_triangle(int x0, int y0, int x1, int y1, int x2, int y2, int xp,
                       int yp) {
  float b1, b2, b3;
  std::tie(b1, b2, b3) =
      barycentric_coordinates(x0, y0, x1, y1, x2, y2, xp, yp);
  return b1 >= 0.0 && b1 <= 1.0 && b2 >= 0.0 && b2 <= 1.0 && b3 >= 0.0 &&
         b3 <= 1.0;
}

void sort_triangle_points(int &x0, int &y0, int &x1, int &y1, int &x2,
                          int &y2) {
  // Bubble sort is not really ideal in general.
  // It could be changed to use a more efficient algorithm,
  // but for only 3 values, it should suffice.
  // Moreover, implementing sorting on an array/vector of tuples
  // is probably more overhead.
  if (y0 > y1) {
    std::swap(x0, x1);
    std::swap(y0, y1);
  }
  if (y0 > y2) {
    std::swap(x0, x2);
    std::swap(y0, y2);
  }
  if (y1 > y2) {
    std::swap(x1, x2);
    std::swap(y1, y2);
  }
}

float slope(int x0, int y0, int x1, int y1) {
  float m = ((float)(y1 - y0)) / ((float)(x1 - x0));
  if (std::isinf(m) || std::isnan(m)) {
    // This is a special case for two things:
    //
    // IEEE 754 specifies ∞ × 0 / 0 × ∞ to be an invalid operation,
    // and therefore return NaN.
    // That makes the computation of Δy fail when x0 == x1.
    //
    // In the case that additionally y0 == y1,
    // the expression is 0/0, also defined in IEEE 754 as invalid.
    m = 0;
  }
  return m;
}