Fixing the code for Line Clipping Algorithm

Question

Midpoint subdivision algorithm [page-93(104)]works on the basis of dividing a line into smaller segments and tests each segment to find whether they are within the visible boundary of the clipping region or not.

In the Binary search algorithm, we find the middle element and then either choose right hand side or left hand side.

But, here, as the following image shows, after first segmentation, we find that both of the subsections are actually disputed. So, both of them are candidates for further subdivisions. So, we can not proceed like Binary Search.

I am using iterative method. But, the following code doesn't work:

    Line2d GetClippedLine()
    {
        Line2d clippingCandidate = this->line;

        std::vector<Line2d> lines = clippingCandidate.GetMidpointSubLines();

        while(lines[0] != lines[1])
        {
            lines = clippingCandidate.GetMidpointSubLines();

            Line2d one = lines[0];
            Line2d two = lines[1]; 

            if(one.IsClippingCandidate(rectangle))
            {
                clippingCandidate = one;
            }
            if(two.IsClippingCandidate(rectangle))
            {           
                clippingCandidate = two;
            }

            if(one.IsVisible(rectangle))
            {
                Coordinates2d::Draw(one, Yellow);
            }
            if(two.IsVisible(rectangle))
            {
                Coordinates2d::Draw(two, Yellow);
            }

            clippingCandidate.Show();
            //std::cout<<"++";
            //two.Show();
            std::cout<<"\n";
        }

        return clippingCandidate;
    }

Answer 1

You're exactly right to ask. Explanations of midpoint subdivision have arisen that are very sloppy or just wrong. It looks like your code is based on one of these bad sources.

M-S is only useful for finding intersections when you already know that a segment straddles a clipping boundary (one end point on each side), and it's generally implemented with integers. It was originally used as a subroutine in a variation of the complete clipping algorithm of Cohen and Sutherland.

See the Wikipedia article if you're not familiar with C-S. The "out codes" guide successive clipping against infinite lines containing the viewport boundaries. In the pseudocode there, you'd replace the floating point math with M-S.

Say you are clipping against the left boundary at x=C, and the line segment that straddles it is P0(x0,y0)---P1(x1,y1). Say also x0<C<=x1, so P0 is known to be outside the boundary. Then the M-S algorithm is:

tx1 = x1; // don't modify P1; it's inside the boundary
ty1 = y1;
while (x0 < C) {
  xm = (x0 + tx1 + 1) >> 1;
  ym = (y0 + ty1 + 1) >> 1;
  if (xm <= C) { // the midpoint is on or outside the boundary
    x0 = xm;     //   move P0
    y0 = ym;
  } else {       // the midpoint is strictly inside
    tx1 = xm;    //   move P1 
    ty1 = ym;
  }
}
// The clipped segment is (x0,x1)--(y0,y1).

You need 3 other minor variations of this for the other 3 boundaries.

Termination conditions are tricky. The + 1s are necessary to avoid looping forever in the case x0 = C-1 and tx1 = C: (C + C - 1 + 1) >> 1 == C, so the next iteration will terminate.

Having said that, midpoint subdivision is pretty much obsolete. It was useful with processors that had only integer arithmetic (often the case until at least the mid 90's; I implemented it in 8088 assembly language in 1984). Finding the midpoint requires only division by 2, which is integer right shift, so it was possible to clip with no more than ceiling(log_2 n) fast iterations for coordinates of max size n. These days with floating point units operating at gigaflop rates, it's probably faster and certainly easier to clip with floating point.

Addition Just for fun, implemented in C:

#include <stdio.h>
#include <stdlib.h>

typedef unsigned OUTCODE;
typedef int COORD;
typedef int BOOL;
#define TRUE 1
#define FALSE 0

#define XMIN 0
#define YMIN 0
#define XMAX 5000
#define YMAX 3000

// Not strictly portable, but usually fine.
#define SIGN_BIT (~(~0u >> 1))

#define LEFT SIGN_BIT
#define TOP (LEFT >> 1)
#define RIGHT (TOP >> 1)
#define BOTTOM (RIGHT >> 1)
#define ALL (LEFT | BOTTOM | RIGHT | TOP)

// Mask the sign bit.
#define M(X) ((X) & SIGN_BIT)

// Shift previous value and mask in the new sign bit.
#define SM(Prev, New) (((OUTCODE)(Prev) >> 1) | M(New))

__inline OUTCODE outcode(COORD x, COORD y) {
  return SM(SM(SM(M(YMAX - y), XMAX - x), y - YMIN), x - XMIN);
}

// In the S-T coordinate system, pO is outside boundary C and will be moved
// to the boundary while pI doesn't move. I is the termination correction.
#define MOVE_TO_BOUNDARY(SO, TO, SI, TI, C, I, IS_OUTSIDE) do { \
  COORD tsi = SI, tti = TI; \
  while (SO IS_OUTSIDE C) { \
    COORD sm = (tsi + SO + I) >> 1; \
    COORD tm = (tti + TO + I) >> 1; \
    if (sm IS_OUTSIDE ## = C) { \
      SO = sm; \
      TO = tm; \
    } else { \
      tsi = sm; \
      tti = tm; \
    } \
  } \
} while (0)

BOOL clip(COORD *x0p, COORD *y0p, COORD *x1p, COORD *y1p) {
  COORD x0 = *x0p, y0 = *y0p, x1 = *x1p, y1 = *y1p;
  OUTCODE code0 = outcode(x0, y0);
  OUTCODE code1 = outcode(x1, y1);
  for (;;) {
    if ((code0 | code1) == 0) {
      *x0p = x0; *y0p = y0; *x1p = x1; *y1p = y1;
      return TRUE;
    } else if (code0 & code1) {
      return FALSE;
    } else if (code0) {
      if (code0 & BOTTOM)     MOVE_TO_BOUNDARY(y0, x0, y1, x1, YMAX, 0, >);
      else if (code0 & RIGHT) MOVE_TO_BOUNDARY(x0, y0, x1, y1, XMAX, 0, >);
      else if (code0 & TOP)   MOVE_TO_BOUNDARY(y0, x0, y1, x1, YMIN, 1, <);
      else /* LEFT */         MOVE_TO_BOUNDARY(x0, y0, x1, y1, XMIN, 1, <);
      code0 = outcode(x0, y0);
    } else {
      if (code1 & BOTTOM)     MOVE_TO_BOUNDARY(y1, x1, y0, x0, YMAX, 0, >);
      else if (code1 & RIGHT) MOVE_TO_BOUNDARY(x1, y1, x0, y0, XMAX, 0, >);
      else if (code1 & TOP)   MOVE_TO_BOUNDARY(y1, x1, y0, x0, YMIN, 1, <);
      else /* LEFT */         MOVE_TO_BOUNDARY(x1, y1, x0, y0, XMIN, 1, <);
      code1 = outcode(x1, y1);
    }
  }
}

int main(void) {
  int n = 0, margin = 2000;
  for (;;) {
    // Generate some random points around the viewport.
    int x0 = rand() % (2 * margin + XMAX - XMIN) - margin;
    int y0 = rand() % (2 * margin + YMAX - YMIN) - margin;
    int x1 = rand() % (2 * margin + XMAX - XMIN) - margin;
    int y1 = rand() % (2 * margin + YMAX - YMIN) - margin;
    printf("a(%d, %d)--(%d, %d) %x--%x\n", x0, y0, x1, y1,
      outcode(x0,y0) >> 28, outcode(x1,y1) >> 28);
    BOOL r = clip(&x0, &y0, &x1, &y1);
    printf("a(%d, %d)--(%d, %d): %d\n", x0, y0, x1, y1, r);
  }
  return 0;
}

On my MacBook it clips a billion segments in 90 seconds. It would be interesting to see how this compares to floating point C-S.