7

I'm really close to getting a thick ellipse algorithm working but I'm having a bit of trouble. I've taken the midpoint thick circle algorithm from here, and the midpoint ellipse algorithm from here, and I'm trying to combine them together to get the midpoint thick ellipse algorithm. I'm doing this because Googling "midpoint thick ellipse algorithm" didn't show what I'm looking for. The output from my attempt resembles a thick circle (images are at the bottom of the post).

This is the image code (just a placeholder):

struct Point {
  int x, y;
};

struct Image {};
using Color = int;

void setPixel(Image &, Color, Point) {
  // ...
}

void horiLine(Image &image, Color color, Point first, int last) {
  while (first.x <= last) {
    setPixel(image, color, first);
    first.x++;
  }
}

void vertLine(Image &image, Color color, Point first, int last) {
  while (first.y <= last) {
    setPixel(image, color, first);
    first.y++;
  }
}

Here's the midpoint thick circle algorithm:

void midpointCircleThick(
  Image &image,
  Color color,
  Point center,
  int innerRadius,
  int outerRadius
) {
  int innerX = innerRadius;
  int outerX = outerRadius;
  int posY = 0;
  int innerErr = 1 - innerRadius;
  int outerErr = 1 - outerRadius;

  while (outerX >= posY) {
    horiLine(image, color, {center.x + innerX, center.y + posY},   center.x + outerX);
    vertLine(image, color, {center.x + posY,   center.y + innerX}, center.y + outerX);
    horiLine(image, color, {center.x - outerX, center.y + posY},   center.x - innerX);
    vertLine(image, color, {center.x - posY,   center.y + innerX}, center.y + outerX);

    horiLine(image, color, {center.x - outerX, center.y - posY},   center.x - innerX);
    vertLine(image, color, {center.x - posY,   center.y - outerX}, center.y - innerX);
    horiLine(image, color, {center.x + innerX, center.y - posY},   center.x + outerX);
    vertLine(image, color, {center.x + posY,   center.y - outerX}, center.y - innerX);

    posY++;

    if (outerErr < 0) {
      outerErr += 2 * posY + 1;
    } else {
      outerX--;
      outerErr += 2 * (posY - outerX) + 1;
    }

    if (posY > innerRadius) {
      innerX = posY;
    } else {
      if (innerErr < 0) {
        innerErr += 2 * posY + 1;
      } else {
        innerX--;
        innerErr += 2 * (posY - innerX) + 1;
      }
    }
  }
}

Here's the midpoint ellipse algorithm:

void midpointEllipse(
  Image &image,
  Color color,
  Point center,
  Point radius
) {
  Point pos = {radius.x, 0};
  Point delta = {
    2 * radius.y * radius.y * pos.x,
    2 * radius.x * radius.x * pos.y
  };
  int err = radius.x * radius.x
          - radius.y * radius.y * radius.x
          + (radius.y * radius.y) / 4;

  while (delta.y < delta.x) {
    setPixel(image, color, {center.x + pos.x, center.y + pos.y});
    setPixel(image, color, {center.x + pos.x, center.y - pos.y});
    setPixel(image, color, {center.x - pos.x, center.y + pos.y});
    setPixel(image, color, {center.x - pos.x, center.y - pos.y});

    pos.y++;

    if (err < 0) {
      delta.y += 2 * radius.x * radius.x;
      err += delta.y + radius.x * radius.x;
    } else {
      pos.x--;
      delta.y += 2 * radius.x * radius.x;
      delta.x -= 2 * radius.y * radius.y;
      err += delta.y - delta.x + radius.x * radius.x;
    }
  }

  err = radius.x * radius.x * (pos.y * pos.y + pos.y)
      + radius.y * radius.y * (pos.x - 1) * (pos.x - 1)
      - radius.y * radius.y * radius.x * radius.x;

  while (pos.x >= 0) {
    setPixel(image, color, {center.x + pos.x, center.y + pos.y});
    setPixel(image, color, {center.x + pos.x, center.y - pos.y});
    setPixel(image, color, {center.x - pos.x, center.y + pos.y});
    setPixel(image, color, {center.x - pos.x, center.y - pos.y});

    pos.x--;

    if (err > 0) {
      delta.x -= 2 * radius.y * radius.y;
      err += radius.y * radius.y - delta.x;
    } else {
      pos.y++;
      delta.y += 2 * radius.x * radius.x;
      delta.x -= 2 * radius.y * radius.y;
      err += delta.y - delta.x + radius.y * radius.y;
    }
  }
}

I attempted to combine the two algorithms and this is what I have so far. I left some ? where I'm not sure about the code. I am well aware of the messiness and duplication here. I just want to get it working before I worry about what the code looks like.

void midpointEllipseThick(
  Image &image,
  Color color,
  Point center,
  Point innerRadius,
  Point outerRadius
) {
  int innerX = innerRadius.x;
  int outerX = outerRadius.x;
  int posY = 0;
  Point innerDelta = {
    2 * innerRadius.y * innerRadius.y * innerX,
    2 * innerRadius.x * innerRadius.x * posY
  };
  Point outerDelta = {
    2 * outerRadius.y * outerRadius.y * outerX,
    2 * outerRadius.x * outerRadius.x * posY
  };
  int innerErr = innerRadius.x * innerRadius.x
               - innerRadius.y * innerRadius.y * innerRadius.x
               + (innerRadius.y * innerRadius.y) / 4;
  int outerErr = outerRadius.x * outerRadius.x
               - outerRadius.y * outerRadius.y * outerRadius.x
               + (outerRadius.y * outerRadius.y) / 4;

  while (outerDelta.y < outerDelta.x) { // ?
    horiLine(image, color, {center.x + innerX, center.y + posY},   center.x + outerX);
    vertLine(image, color, {center.x + posY,   center.y + innerX}, center.y + outerX);
    horiLine(image, color, {center.x - outerX, center.y + posY},   center.x - innerX);
    vertLine(image, color, {center.x - posY,   center.y + innerX}, center.y + outerX);

    horiLine(image, color, {center.x - outerX, center.y - posY},   center.x - innerX);
    vertLine(image, color, {center.x - posY,   center.y - outerX}, center.y - innerX);
    horiLine(image, color, {center.x + innerX, center.y - posY},   center.x + outerX);
    vertLine(image, color, {center.x + posY,   center.y - outerX}, center.y - innerX);

    posY++;

    if (outerErr < 0) {
      outerDelta.y += 2 * outerRadius.x * outerRadius.x;
      outerErr += outerDelta.y + outerRadius.x * outerRadius.x;
    } else {
      outerX--;
      outerDelta.y += 2 * outerRadius.x * outerRadius.x;
      outerDelta.x -= 2 * outerRadius.y * outerRadius.y;
      outerErr += outerDelta.y - outerDelta.x + outerRadius.x * outerRadius.x;
    }

    // ?
    // if (posY > innerRadius.y) {
    //   innerX = posY;
    // } else {
      if (innerErr < 0) {
        innerDelta.y += 2 * innerRadius.x * innerRadius.x;
        innerErr += innerDelta.y + innerRadius.x * innerRadius.x;
      } else {
        innerX--;
        innerDelta.y += 2 * innerRadius.x * innerRadius.x;
        innerDelta.x -= 2 * innerRadius.y * innerRadius.y;
        innerErr += innerDelta.y - innerDelta.x + innerRadius.x * innerRadius.x;
      }
    // }
  }

  innerErr = innerRadius.x * innerRadius.x * (posY * posY + posY)
           + innerRadius.y * innerRadius.y * (innerX - 1) * (innerX - 1)
           - innerRadius.y * innerRadius.y * innerRadius.x * innerRadius.x;
  outerErr = outerRadius.x * outerRadius.x * (posY * posY + posY)
           + outerRadius.y * outerRadius.y * (outerX - 1) * (outerX - 1)
           - outerRadius.y * outerRadius.y * outerRadius.x * outerRadius.x;

  while (outerX >= 0) { // ?
    horiLine(image, color, {center.x + innerX, center.y + posY},   center.x + outerX);
    vertLine(image, color, {center.x + posY,   center.y + innerX}, center.y + outerX);
    horiLine(image, color, {center.x - outerX, center.y + posY},   center.x - innerX);
    vertLine(image, color, {center.x - posY,   center.y + innerX}, center.y + outerX);

    horiLine(image, color, {center.x - outerX, center.y - posY},   center.x - innerX);
    vertLine(image, color, {center.x - posY,   center.y - outerX}, center.y - innerX);
    horiLine(image, color, {center.x + innerX, center.y - posY},   center.x + outerX);
    vertLine(image, color, {center.x + posY,   center.y - outerX}, center.y - innerX);

    outerX--; // ?
    innerX--;

    if (outerErr > 0) {
      outerDelta.x -= 2 * outerRadius.y * outerRadius.y;
      outerErr += outerRadius.y * outerRadius.y - outerDelta.x;
    } else {
      posY++;
      outerDelta.y += 2 * outerRadius.x * outerRadius.x;
      outerDelta.x -= 2 * outerRadius.y * outerRadius.y;
      outerErr += outerDelta.y - outerDelta.x + outerRadius.y * outerRadius.y;
    }

    // ?
    // if (innerX < -innerRadius.x) {

    // } else {
      if (outerErr > 0) {
        innerDelta.x -= 2 * innerRadius.y * innerRadius.y;
        innerErr += innerRadius.y * innerRadius.y - innerDelta.x;
      } else {
        posY++;
        innerDelta.y += 2 * innerRadius.x * innerRadius.x;
        innerDelta.x -= 2 * innerRadius.y * innerRadius.y;
        outerErr += innerDelta.y - innerDelta.x + innerRadius.y * innerRadius.y;
      }
    // }
  }
}

Here's a thick circle with innerRadius = 22; outerRadius = 24:

Thick circle

Here's an ellipse with radius = {32, 24}:

Ellipse

Here's (what's supposed to be) a thick ellipse with innerRadius = {30, 22}; outerRadius = {32, 24}:

Thick ellipse

I'm close but not quite there. Could someone who knows more about this stuff than I do get me over the finish line?

Indiana Kernick
  • 5,041
  • 2
  • 20
  • 50
  • @Scheff I’m using Qt for testing this as well actually! All you have to do is create type aliases for `Image` and `Color` and then implement `setPixel` to call `QImage::setPixel` – Indiana Kernick May 04 '19 at 08:58
  • How about drawing an ellipse, followed by similar but slightly smaller ellipses to create the thickness? – גלעד ברקן May 04 '19 at 11:42
  • 1
    @גלעדברקן That would probably leave holes. Also it isn’t as efficient as drawing horizontal and vertical lines. – Indiana Kernick May 04 '19 at 11:45

2 Answers2

11

I must admit that I strongly believe there is more symmetry in a circle than in an ellipse. Where a circle might be mirrored on any axis through center, for an ellipse, this is possible only with x and y axis in general. Hence, I believe that the midPointCircleThick() cannot be adapted for an ellipse.

So, I started my implementation with the midpointEllipse() provided by the OP.

These were my basic thoughts:

  • IMHO, the Bresenham Line algorithm is the origin of the Midpoint Circle algorithm as well as the Midpoint Ellipse algorithm. This might be helpful to understand the error/delta magic which is used. It's much simpler for a line but follows the same idea adapted to x²/a² + y²/b² = 1 (the ellipse equation).

  • With origin in center of ellipse, the midpointEllipse() renders all 4 quadrants simultaneously (exploiting the symmetry). Hence, only the curve in one quadrant has to be computed effectively. The curve is in this area monotonic.

  • The midpointEllipse() has two regions:

    1. Starting at points on x axis, ∆y > ∆x until cross-even.
    2. Afterwards, ∆x > ∆y.

My concept was to adapt the midpointEllipse() that way, that the code is "duplicated" to manage two points (one for inner border, one for outer) with identical y coordinates to draw horizontal lines (span lines).

My first observation was that the new algorithm has to manage a final phase (for innerRadius.y < y ≤ outerRadius.y where only the points on outer border have to be considered.

Remembering that the original algorithm has two regions, there are now two regions for the outer border, two regions for the inner border, and the two phases mentioned above. This allows a variety of combinations. (To get this managed was the main effort in my implementation.)

The sample implementation (based on Qt to have a simple visualization):

#include <functional>

#include <QtWidgets>

class View: public QLabel {

  public:
    View(QWidget *pQParent = nullptr):
      QLabel(pQParent)
    { }
    virtual ~View() = default;

    View(const View&) = delete;
    View& operator=(const View&) = delete;

  protected:

    virtual void paintEvent(QPaintEvent *pQEvent) override;
};

struct Point { int x, y; };

using Color = QColor;

void midpointEllipse(
  Point center,
  Point radius,
  std::function<void(const Color&, const Point&)> setPixel)
{
  Point pos = { radius.x, 0 };
  Point delta = {
    2 * radius.y * radius.y * pos.x,
    2 * radius.x * radius.x * pos.y
  };
  int err = radius.x * radius.x
    - radius.y * radius.y * radius.x
    + (radius.y * radius.y) / 4;

  while (delta.y < delta.x) {
    setPixel(Qt::blue, { center.x + pos.x, center.y + pos.y });
    setPixel(Qt::blue, { center.x + pos.x, center.y - pos.y });
    setPixel(Qt::blue, { center.x - pos.x, center.y + pos.y });
    setPixel(Qt::blue, { center.x - pos.x, center.y - pos.y });

    pos.y++;

    if (err < 0) {
      delta.y += 2 * radius.x * radius.x;
      err += delta.y + radius.x * radius.x;
    } else {
      pos.x--;
      delta.y += 2 * radius.x * radius.x;
      delta.x -= 2 * radius.y * radius.y;
      err += delta.y - delta.x + radius.x * radius.x;
    }
  }

  err = radius.x * radius.x * (pos.y * pos.y + pos.y)
    + radius.y * radius.y * (pos.x - 1) * (pos.x - 1)
    - radius.y * radius.y * radius.x * radius.x;

  while (pos.x >= 0) {
    setPixel(Qt::yellow, { center.x + pos.x, center.y + pos.y });
    setPixel(Qt::yellow, { center.x + pos.x, center.y - pos.y });
    setPixel(Qt::yellow, { center.x - pos.x, center.y + pos.y });
    setPixel(Qt::yellow, { center.x - pos.x, center.y - pos.y });

    pos.x--;

    if (err > 0) {
      delta.x -= 2 * radius.y * radius.y;
      err += radius.y * radius.y - delta.x;
    } else {
      pos.y++;
      delta.y += 2 * radius.x * radius.x;
      delta.x -= 2 * radius.y * radius.y;
      err += delta.y - delta.x + radius.y * radius.y;
    }
  }
}

void midpointEllipseThick(
  Point center,
  Point innerRadius,
  Point outerRadius,
  std::function<void(const Color&, const Point&, int)> horiLine)
{
  /// @todo validate/correct innerRadius and outerRadius
  Point pos = { outerRadius.x, 0 };
  Point deltaOuter = {
    2 * outerRadius.y * outerRadius.y * pos.x,
    2 * outerRadius.x * outerRadius.x * pos.y
  };
  auto errOuterYX
    = [&]() {
      return outerRadius.x * outerRadius.x
        - outerRadius.y * outerRadius.y * outerRadius.x
        + (outerRadius.y * outerRadius.y) / 4;
    };
  auto errOuterXY
    = [&]() {
      return outerRadius.x * outerRadius.x * (pos.y * pos.y + pos.y)
        + outerRadius.y * outerRadius.y * (pos.x - 1) * (pos.x - 1)
        - outerRadius.y * outerRadius.y * outerRadius.x * outerRadius.x;
    };
  int errOuter = errOuterYX();
  int xInner = innerRadius.x;
  Point deltaInner = {
    2 * innerRadius.y * innerRadius.y * xInner,
    2 * innerRadius.x * innerRadius.x * pos.y
  };
  auto errInnerYX
    = [&]() {
      return innerRadius.x * innerRadius.x
        - innerRadius.y * innerRadius.y * innerRadius.x
        + (innerRadius.y * innerRadius.y) / 4;
    };
  auto errInnerXY
    = [&]() {
      return innerRadius.x * innerRadius.x * (pos.y * pos.y + pos.y)
        + innerRadius.y * innerRadius.y * (xInner - 1) * (xInner - 1)
        - innerRadius.y * innerRadius.y * innerRadius.x * innerRadius.x;
    };
  int errInner = errInnerYX();
  // helpers (to reduce code duplication)
  auto stepOuterYX
    = [&]() {
      ++pos.y;
      if (errOuter < 0) {
        deltaOuter.y += 2 * outerRadius.x * outerRadius.x;
        errOuter += deltaOuter.y + outerRadius.x * outerRadius.x;
      } else {
        --pos.x;
        deltaOuter.y += 2 * outerRadius.x * outerRadius.x;
        deltaOuter.x -= 2 * outerRadius.y * outerRadius.y;
        errOuter += deltaOuter.y - deltaOuter.x + outerRadius.x * outerRadius.x;
      }
    };
  auto stepOuterXY
    = [&]() {
      while (--pos.x > 0) {
        if (errOuter > 0) {
          deltaOuter.x -= 2 * outerRadius.y * outerRadius.y;
          errOuter += outerRadius.y * outerRadius.y - deltaOuter.x;
        } else {
          ++pos.y;
          deltaOuter.y += 2 * outerRadius.x * outerRadius.x;
          deltaOuter.x -= 2 * outerRadius.y * outerRadius.y;
          errOuter += deltaOuter.y - deltaOuter.x + outerRadius.y * outerRadius.y;
          break;
        }
      }
    };
  auto stepInnerYX
    = [&]() {
      if (errInner < 0) {
        deltaInner.y += 2 * innerRadius.x * innerRadius.x;
        errInner += deltaInner.y + innerRadius.x * innerRadius.x;
      } else {
        --xInner;
        deltaInner.y += 2 * innerRadius.x * innerRadius.x;
        deltaInner.x -= 2 * innerRadius.y * innerRadius.y;
        errInner += deltaInner.y - deltaInner.x + innerRadius.x * innerRadius.x;
      }
    };
  auto stepInnerXY
    = [&]() {
      while (--xInner >= 0) {
        if (errInner > 0) {
          deltaInner.x -= 2 * innerRadius.y * innerRadius.y;
          errInner += innerRadius.y * innerRadius.y - deltaInner.x;
        } else {
          deltaInner.y += 2 * innerRadius.x * innerRadius.x;
          deltaInner.x -= 2 * innerRadius.y * innerRadius.y;
          errInner += deltaInner.y - deltaInner.x + innerRadius.y * innerRadius.y;
          break;
        }
      }
    };
  // 1st phase
  while (deltaOuter.y < deltaOuter.x && deltaInner.y < deltaInner.x) {
    horiLine(Qt::blue, { center.x - pos.x, center.y + pos.y }, center.x - xInner);
    horiLine(Qt::blue, { center.x + pos.x, center.y + pos.y }, center.x + xInner);
    horiLine(Qt::blue, { center.x - pos.x, center.y - pos.y }, center.x - xInner);
    horiLine(Qt::blue, { center.x + pos.x, center.y - pos.y }, center.x + xInner);
    stepOuterYX();
    stepInnerYX();
  }

  // 2nd phase
  if (deltaOuter.y < deltaOuter.x) { // inner flipped
    //errOuter = errOuterYX();
    errInner = errInnerXY();
    while (deltaOuter.y < deltaOuter.x && xInner >= 0) {
      horiLine(Qt::green, { center.x - pos.x, center.y + pos.y }, center.x - xInner);
      horiLine(Qt::green, { center.x + pos.x, center.y + pos.y }, center.x + xInner);
      horiLine(Qt::green, { center.x - pos.x, center.y - pos.y }, center.x - xInner);
      horiLine(Qt::green, { center.x + pos.x, center.y - pos.y }, center.x + xInner);
      stepOuterYX();
      stepInnerXY();
    }
    //errOuter = errOuterYX();
    while (deltaOuter.y < deltaOuter.x) {
      horiLine(Qt::red, { center.x - pos.x, center.y + pos.y }, center.x + pos.x);
      horiLine(Qt::red, { center.x - pos.x, center.y - pos.y }, center.x + pos.x);
      stepOuterYX();
    }
  } else { // outer flipped
    errOuter = errOuterXY();
    //errInner = errInnerYX();
    while (deltaInner.y < deltaInner.x) {
      horiLine(Qt::cyan, { center.x - pos.x, center.y + pos.y }, center.x - xInner);
      horiLine(Qt::cyan, { center.x + pos.x, center.y + pos.y }, center.x + xInner);
      horiLine(Qt::cyan, { center.x - pos.x, center.y - pos.y }, center.x - xInner);
      horiLine(Qt::cyan, { center.x + pos.x, center.y - pos.y }, center.x + xInner);
      stepOuterXY();
      stepInnerYX();
    }
    //errOuter = errOuterXY();
  }
  // 3rd phase
  errOuter = errOuterXY();
  errInner = errInnerXY();
  while (xInner >= 0) {
    horiLine(Qt::yellow, { center.x - pos.x, center.y + pos.y }, center.x - xInner);
    horiLine(Qt::yellow, { center.x + pos.x, center.y + pos.y }, center.x + xInner);
    horiLine(Qt::yellow, { center.x - pos.x, center.y - pos.y }, center.x - xInner);
    horiLine(Qt::yellow, { center.x + pos.x, center.y - pos.y }, center.x + xInner);
    stepOuterXY();
    stepInnerXY();
  }
  // 4th phase
  //errOuter = errOuterXY();
  while (pos.x >= 0) {
    horiLine(Qt::magenta, { center.x - pos.x, center.y + pos.y }, center.x + pos.x);
    horiLine(Qt::magenta, { center.x - pos.x, center.y - pos.y }, center.x + pos.x);
    stepOuterXY();
  }
}

void View::paintEvent(QPaintEvent*)
{
  QPainter qPainter(this);
#if 0 // warm up
  auto setPixel
    = [&](const Color &color, const Point &point)
    {
      qPainter.setPen(color);
      qPainter.drawPoint(point.x, point.y);
    };
  Point center = { 0.5 * width(), 0.5 * height() };
  midpointEllipse(center, center, setPixel);
#else // my attempt to adapt it to thick ellipses
  auto horiLine
    = [&](const Color &color, const Point &pos0, int x1)
    {
      qPainter.setPen(color);
      qPainter.drawLine(pos0.x, pos0.y, x1, pos0.y);
    };
  Point center = { 0.5 * width(), 0.5 * height() };
  Point innerRadius = { 0.5 * center.x, 0.5 * center.y };
  Point outerRadius = { 0.9 * center.x, 0.9 * center.y };
  midpointEllipseThick(center, innerRadius, outerRadius, horiLine);
#endif // 0
}

int main(int argc, char **argv)
{
  qDebug() << "Qt Version:" << QT_VERSION_STR;
  QApplication app(argc, argv);
  // setup UI
  View qWin;
  qWin.setWindowTitle(QString::fromUtf8("Draw Thick Ellipse"));
  qWin.resize(320, 240);
  qWin.show();
  // runtime loop
  return app.exec();
}

Compiled an tested in VS2017 (Qt 5.11.2):

Snapshot of testQThickEllipse Snapshot of testQThickEllipse (with different width/height) Snapshot of testQThickEllipse (with yet another width/height)

I used colors to visualize the different combinations of regions and phases. This is intended to simply illustrate which part of code was responsible to render which part of ellipse.


I was a bit uncertain about the else case in // 2nd phase. I tested with

  Point center = { 0.5 * width(), 0.5 * height() };
  Point innerRadius = { 0.3 * center.x, 0.8 * center.y };
  Point outerRadius = { 0.9 * center.x, 0.9 * center.y };
  midpointEllipseThick(center, innerRadius, outerRadius, horiLine);

and got this:

Snapshot of testQThickEllipse (with different arguments)

Now, the // 1st phase stops due to failing deltaOuter.y < deltaOuter.x (and cyan areas appear).


OP complained about poor handling of edge cases like e.g. innerRadius = outerRadius;. I checked it with the following test set:

  Point center = { 0.5 * width(), 0.5 * height() };
  // test edge cases
  { Point outerRadius = { 0.9 * center.x, 0.9 * center.y };
    Point innerRadius = { outerRadius.x, outerRadius.y };
    Old::midpointEllipseThick(center, innerRadius, outerRadius, horiLine);
  }
  { Point outerRadius = { 0.8 * center.x, 0.8 * center.y };
    Point innerRadius = { outerRadius.x - 1, outerRadius.y };
    Old::midpointEllipseThick(center, innerRadius, outerRadius, horiLine);
  }
  { Point outerRadius = { 0.7 * center.x, 0.7 * center.y };
    Point innerRadius = { outerRadius.x, outerRadius.y - 1 };
    Old::midpointEllipseThick(center, innerRadius, outerRadius, horiLine);
  }
  { Point outerRadius = { 0.6 * center.x, 0.6 * center.y };
    Point innerRadius = { outerRadius.x - 1, outerRadius.y - 1 };
    Old::midpointEllipseThick(center, innerRadius, outerRadius, horiLine);
  }
  { Point outerRadius = { 0.5 * center.x, 0.5 * center.y };
    Point innerRadius = { outerRadius.x - 2, outerRadius.y - 2 };
    Old::midpointEllipseThick(center, innerRadius, outerRadius, horiLine);
  }

changed Qt::yellow to Qt::darkgray (for a better contrast) and got this:

Snapshot of testQThickEllipse (edge tests)

It becomes obvious that gaps appear when ∆xy →y+1 > xOuter - xInner.

To fix this issue, the ∆xy →y+1 has to be considered as well for the generation of span lines. To achieve this, I modified the iterations for ∆x ≥ ∆y (in the bottom part of function):

void midpointEllipseThick(
  Point center,
  Point innerRadius,
  Point outerRadius,
  std::function<void(const Color&, const Point&, int)> horiLine)
{
  /// @todo validate/correct innerRadius and outerRadius
  Point pos = { outerRadius.x, 0 };
  Point deltaOuter = {
    2 * outerRadius.y * outerRadius.y * pos.x,
    2 * outerRadius.x * outerRadius.x * pos.y
  };
  auto errOuterYX
    = [&]() {
      return outerRadius.x * outerRadius.x
        - outerRadius.y * outerRadius.y * outerRadius.x
        + (outerRadius.y * outerRadius.y) / 4;
    };
  auto errOuterXY
    = [&]() {
      return outerRadius.x * outerRadius.x * (pos.y * pos.y + pos.y)
        + outerRadius.y * outerRadius.y * (pos.x - 1) * (pos.x - 1)
        - outerRadius.y * outerRadius.y * outerRadius.x * outerRadius.x;
    };
  int errOuter;
  int xInner = innerRadius.x;
  Point deltaInner = {
    2 * innerRadius.y * innerRadius.y * xInner,
    2 * innerRadius.x * innerRadius.x * pos.y
  };
  auto errInnerYX
    = [&]() {
      return innerRadius.x * innerRadius.x
        - innerRadius.y * innerRadius.y * innerRadius.x
        + (innerRadius.y * innerRadius.y) / 4;
    };
  auto errInnerXY
    = [&]() {
      return innerRadius.x * innerRadius.x * (pos.y * pos.y + pos.y)
        + innerRadius.y * innerRadius.y * (xInner - 1) * (xInner - 1)
        - innerRadius.y * innerRadius.y * innerRadius.x * innerRadius.x;
    };
  int errInner;
  // helpers (to reduce code duplication)
  auto stepOuterYX
    = [&]() {
      ++pos.y;
      if (errOuter < 0) {
        deltaOuter.y += 2 * outerRadius.x * outerRadius.x;
        errOuter += deltaOuter.y + outerRadius.x * outerRadius.x;
      } else {
        --pos.x;
        deltaOuter.y += 2 * outerRadius.x * outerRadius.x;
        deltaOuter.x -= 2 * outerRadius.y * outerRadius.y;
        errOuter += deltaOuter.y - deltaOuter.x + outerRadius.x * outerRadius.x;
      }
    };
  auto stepInnerYX
    = [&]() {
      if (errInner < 0) {
        deltaInner.y += 2 * innerRadius.x * innerRadius.x;
        errInner += deltaInner.y + innerRadius.x * innerRadius.x;
      } else {
        --xInner;
        deltaInner.y += 2 * innerRadius.x * innerRadius.x;
        deltaInner.x -= 2 * innerRadius.y * innerRadius.y;
        errInner += deltaInner.y - deltaInner.x + innerRadius.x * innerRadius.x;
      }
    };
  auto stepOuterXY
    = [&]() {
      while (--pos.x >= 0) {
        if (errOuter > 0) {
          deltaOuter.x -= 2 * outerRadius.y * outerRadius.y;
          errOuter += outerRadius.y * outerRadius.y - deltaOuter.x;
        } else {
          ++pos.y;
          deltaOuter.y += 2 * outerRadius.x * outerRadius.x;
          deltaOuter.x -= 2 * outerRadius.y * outerRadius.y;
          errOuter += deltaOuter.y - deltaOuter.x + outerRadius.y * outerRadius.y;
          break;
        }
      }
    };
  auto stepInnerXY
    = [&]() {
      while (--xInner >= 0) {
        if (errInner > 0) {
          deltaInner.x -= 2 * innerRadius.y * innerRadius.y;
          errInner += innerRadius.y * innerRadius.y - deltaInner.x;
        } else {
          deltaInner.y += 2 * innerRadius.x * innerRadius.x;
          deltaInner.x -= 2 * innerRadius.y * innerRadius.y;
          errInner += deltaInner.y - deltaInner.x + innerRadius.y * innerRadius.y;
          break;
        }
      }
    };
  auto min
    = [](int x1, int x2, int x3) {
      return std::min(std::min(x1, x2), x3);
    };
  // 1st phase
  errOuter = errOuterYX(); // init error for delta y < delta x
  errInner = errInnerYX(); // init error for delta y < delta x
  while (deltaOuter.y < deltaOuter.x && deltaInner.y < deltaInner.x) {
    horiLine(Qt::blue, { center.x - pos.x, center.y + pos.y }, center.x - xInner);
    horiLine(Qt::blue, { center.x + pos.x, center.y + pos.y }, center.x + xInner);
    horiLine(Qt::blue, { center.x - pos.x, center.y - pos.y }, center.x - xInner);
    horiLine(Qt::blue, { center.x + pos.x, center.y - pos.y }, center.x + xInner);
    stepOuterYX();
    stepInnerYX();
  }

  // 2nd phase
  if (deltaOuter.y < deltaOuter.x) { // inner flipped
    //errOuter = errOuterYX(); // still delta y < delta x
    errInner = errInnerXY(); // init error for delta x < delta y
    while (deltaOuter.y < deltaOuter.x && xInner >= 0) {
      horiLine(Qt::green, { center.x - pos.x, center.y + pos.y }, center.x - xInner);
      horiLine(Qt::green, { center.x + pos.x, center.y + pos.y }, center.x + xInner);
      horiLine(Qt::green, { center.x - pos.x, center.y - pos.y }, center.x - xInner);
      horiLine(Qt::green, { center.x + pos.x, center.y - pos.y }, center.x + xInner);
      stepOuterYX();
      stepInnerXY();
    }
    //errOuter = errOuterYX(); // still delta y < delta x
    while (deltaOuter.y < deltaOuter.x) {
      horiLine(Qt::red, { center.x - pos.x, center.y + pos.y }, center.x + pos.x);
      horiLine(Qt::red, { center.x - pos.x, center.y - pos.y }, center.x + pos.x);
      stepOuterYX();
    }
  } else { // outer flipped
    errOuter = errOuterXY(); // init error for delta x < delta y
    //errInner = errInnerYX(); // still delta y < delta x
    while (deltaInner.y < deltaInner.x) {
      Point pos_ = pos;
      stepOuterXY();
      stepInnerYX();
      int xInner_ = std::min(pos.x, xInner);
      horiLine(Qt::cyan, { center.x - pos_.x, center.y + pos_.y }, center.x - xInner_);
      horiLine(Qt::cyan, { center.x + pos_.x, center.y + pos_.y }, center.x + xInner_);
      horiLine(Qt::cyan, { center.x - pos_.x, center.y - pos_.y }, center.x - xInner_);
      horiLine(Qt::cyan, { center.x + pos_.x, center.y - pos_.y }, center.x + xInner_);
    }
  }
  // 3rd phase
  errOuter = errOuterXY(); // init error for delta x < delta y
  errInner = errInnerXY(); // init error for delta x < delta y
  while (xInner >= 0) {
    Point pos_ = pos;
    stepOuterXY();
    int xInner_ = std::min(pos.x, xInner);
    horiLine(Qt::darkGray, { center.x - pos_.x, center.y + pos_.y }, center.x - xInner_);
    horiLine(Qt::darkGray, { center.x + pos_.x, center.y + pos_.y }, center.x + xInner_);
    horiLine(Qt::darkGray, { center.x - pos_.x, center.y - pos_.y }, center.x - xInner_);
    horiLine(Qt::darkGray, { center.x + pos_.x, center.y - pos_.y }, center.x + xInner_);
    stepInnerXY();
  }
  // 4th phase
  //errOuter = errOuterXY(); // still delta x < delta y
  while (pos.x >= 0) {
    horiLine(Qt::magenta, { center.x - pos.x, center.y + pos.y }, center.x + pos.x + 1);
    horiLine(Qt::magenta, { center.x - pos.x, center.y - pos.y }, center.x + pos.x + 1);
    stepOuterXY();
  }
}

The result looks not that bad:

Snapshot of testQThickEllipse (fixed for edge tests)

The gaps are removed.

I realized that there is still the other complained issue about off-by-one error:

The thickness at the top and bottom part of the ellipse seems to be one pixel too small.

Hmmm… That's a question of definition. Whenever a range has to be given, there has to be said whether start and end are (each) inclusive or exclusive. (Compare e.g. with iterator ranges in standard containers – start → inclusive, end → exclusive.)

The Qt doc. dedicates a whole extra chapter to this topic Coordinate System.

What I have to admit: My current algorithm handles this different for horizontal and vertical direction which I would consider as “ugliness”. IMHO, the easiest fix is to make it consistent horizontally and vertically. Afterwards the doc. might be adjusted respectively.

Employee: “Boss! Our recently produced buckets have a hole and lose water.”
Boss: “Good to know. We should mention that in the manual.”

Thus, I fixed the horizontal border size by tweaking the horiLine helper lambda:

  auto horiLine
    = [&](const Color &color, const Point &pos0, int x1)
    {
      qPainter.setPen(color);
      if (x1 != pos0.x) x1 += x1 < pos0.x ? +1 : -1;
      qPainter.drawLine(pos0.x, pos0.y, x1, pos0.y);
    };

Now, I consider the result, at least, as consistent (if not satisfying):

Snapshot of testQThickEllipse (fixed for border width)

The innerRadius appears now as exclusive. If this is not intended, a resp. pre-adjustment of the parameters at the begin of the midpointEllipseThick() could be applied.

Scheff's Cat
  • 19,528
  • 6
  • 28
  • 56
  • 1
    This looks promising but it’s getting very late in my timezone so I’ll have to try this out tomorrow. – Indiana Kernick May 04 '19 at 13:43
  • @Kerndog73 Oh, Australia. I see... ;-) – Scheff's Cat May 04 '19 at 13:44
  • This looks really good but I think there might be an off-by-one error in here somewhere. The thickness at the top and bottom part of the ellipse seems to be one pixel too small. Try drawing something with `outerRadius == innerRadius` and you'll see what I mean. The thickness of the blue and green parts seems to be OK but the magenta and yellow parts are a little thin. I honestly have no idea how this code works so I have no hope of fixing it on my own! – Indiana Kernick May 04 '19 at 23:52
  • Also, those commented error calculations worry me (`//errOuter = errOuterXY();`). Could you double check that they aren't required and remove them, please? BTW, I really like that this only draws horizontal lines. That makes it decently fast for really thick ellipses. Other than those two things, I'm really happy with what you've done! If you fix those problems, I might even open a bounty and award you some extra internet points! – Indiana Kernick May 04 '19 at 23:55
  • @Kerndog73 The commented error calculations denote IMHO useless re-initializations of error when the corresponding curve has not flipped before the resp. loop. Uncommenting them shouldn't change anything (except some needless extra calculations). – Scheff's Cat May 05 '19 at 06:36
  • Oh good so I can safely remove them. How about the off-by-one error? Do you know what might be causing it? – Indiana Kernick May 05 '19 at 06:46
  • @Kerndog73 Concerning the off-by-one error I will have a look after breakfast. ;-) Prohibiting `outerRadius == innerRadius` is not an option? (The `midpointEllipse()` could be used as fall-back in this case.) – Scheff's Cat May 05 '19 at 06:54
  • It doesn’t just affect 1px thick ellipses. It effects all ellipses. Have a close look at `innerRadius = {30, 22}; outerRadius = {32, 24}`. You’ll see there are 3 pixels on the left and right (good) but 2 pixels on the top and bottom (bad). – Indiana Kernick May 05 '19 at 06:59
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/192854/discussion-between-scheff-and-kerndog73). – Scheff's Cat May 05 '19 at 07:00
  • @Kerndog73 Just realized another gap in the final phase and made another small bug fix. To fight through all the edge cases, can be tedious and bothering... – Scheff's Cat May 05 '19 at 11:47
  • @Scheff Have you edited your answer with the latest fixes? If yes, for clarity, you should add a note telling the reader that the published code is now considered as finished. If not, please do it and add a note saying it is finished now. – fpiette Nov 19 '19 at 08:35
  • @fpiette The most bottom source of `midpointEllipseThick()` and the most bottom source of `horiLine()` represent the final state, I found in my local sources. I used varying tests in `paintEvent()` (given in snippets) to explore the various edge cases. TL;DR: You can use the complete sample from the top but, please, replace `midpointEllipseThick()` and `horiLine()` with the latest versions. – Scheff's Cat Nov 19 '19 at 10:44
  • @fpiette Unfortunately, I cannot append the complete sample again. (I tried but exceeded the limit of 30000 lines, sorry.) I could make it a project on github but I've to find free time for this. (It's currently near end of year, the time where everybody remembers what still should be done this year...) ;-) – Scheff's Cat Nov 19 '19 at 10:50
  • @Scheff I would love to have the unit. You can mail it to me. You can easily find my email address by googling "francois piette ics". The top link is likely to be me. (Sorry to not post my email here, this is to prevent spam). Thank you. For your information, I already ported your code (as it appears here) to Delphi and to HLSL. – fpiette Nov 19 '19 at 20:05
1

The problem you are facing is that the outlines of a thick ellipse of constant thickness are not ellipses, they are higher degree curves ! Filling between two ellipses can only give approximations.

On the picture, the red curve corresponds to a constant thickness.

enter image description here

A correct solution is to draw with a thick pen, i.e. sweep a disk of the desired radius by letting its center follow an ellipse, using the standard algorithm.

As such, this is an inefficient procedure, as the successive disks overlap and the pixels will be drawn several times. A solution is to consider the new pixels that are covered by the disk for the eight directions of displacement. These sets of pixels must be computed and tabulated in advance, for the given radius.

To establish the tables, draw a disk and erase it with a disk shifted by one pixel in one of the eight cardinal directions; repeat for all directions.

  • That might actually be a viable solution. Drawing 1px stroked circles in a straight line is only slightly slower than whatever QPainter is doing when it draws thick lines. 1px circles do leave gaps when moving diagonally though. The table of bitmaps sounds like an interesting idea. I’m sceptical on the performance because you’d have to iterate each pixel in a bitmap and most of the pixels aren’t set. I haven’t tried it though. I’ll still award the bounty to Scheff due to the amount of effort he put into it. I might consider changing the accepted answer though. – Indiana Kernick May 06 '19 at 22:27
  • I could do some bit twiddling to make the bitmap loop body branchless. I'm still sceptical because you'd still be looking at a lot of pixels that you don't need to change. – Indiana Kernick May 06 '19 at 22:29
  • 1
    Of course you store a list of the pixels to be set ! –  May 06 '19 at 22:35
  • Ah, that would be smart! – Indiana Kernick May 06 '19 at 22:37
  • @Kerndog73: If I am right, you can build these "delta" lists by means of a modified circle drawing algorithm: every time you draw a pixel, you check which neighbors it has in the cardinal directions. –  May 07 '19 at 05:56
  • @Kerndog73: also note that the method extends to any continuous curve and any brush shape. –  May 07 '19 at 05:57
  • You've lost me. Why are we checking neighboring pixels? – Indiana Kernick May 07 '19 at 06:40
  • @Kerndog73: to find the deltas in the cardinal directions. –  May 07 '19 at 06:41
  • We check the neighboring pixels while drawing the circle? Why? I don't get it. I understood "draw a disk and erase it with a disk shifted by one pixel in one of the eight cardinal directions". At what point would you need to check neighboring pixels? – Indiana Kernick May 07 '19 at 06:46
  • @Kerndog73: this is to avoid drawing the disks, outlines are enough. –  May 07 '19 at 06:49