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.

Answer 2

I think the reason for this is hw implementations. Don't think about it as clipping. Think about it as an exercise on how to sort lines (or later projected triangles) into a quadtree. In the end you can keep clipping to half node subdivisions until every node is just a pixel wide. Do not treat your outside parts as special. They are just really big quads.

Answer 3

Working in geometric modelling for several years, I have always been implementing all subdivision algorithms in recursive way. It is much simpler to code and understand, and I'm pretty sure it is not slower than an iterative solution.

General case

If the clipping is concave, then the answer may consist of arbitrarily many segments. In this case you need a stack (replacing the code stack from recursion). Here is the example of an implementation:

vector<Line> result;  //ordered array of line segments = clipping result
stack<Line> pieces;  //stack of pieces that are not processed yet
pieces.push(inputLine);

while (!pieces.empty()) {
    auto curr = pieces.top();
    pieces.pop();
    if (GetLength(curr) < eps)  //hard stop criterion
      continue;
    auto relative = GetSegmentPositionRelativeToRegion(curr, clippingRegion);
    if (relative == FULLY_INSIDE)
        result.push_back(curr);
    else if (relative == INTERSECTS) { //not fully inside, not fully outside
        Line left, right;   //halves of the curr line segment
        SplitAtMiddle(curr, left, right);
        pieces.push(left);
        pieces.push(right);
    }
}

Convex case

If your clipping region is guaranteed to be convex, then there is exactly one line segment in the result. In this case the subdivision process has very simple structure. Suppose that the input segment is not fully inside or outside clipping region (i.e. it intersects).

At first you subdivide the segment so that one half is outside and the other one intersects with clipping. It means that you can always remember only a single working line segment.

At some moment this rule breaks, and afterwards there are two branches of subdivision. At the left branch there may be two cases of segment subdivision: (intersects + inside) or (outside + intersects). On the right side it is reflected: (inside + intersects) or (intersects + outside). In all the cases only one subsegment intersects, so only one working segment remains. You can write separate code for each of the branches in a loop.

P.S. Of course this description ignores singular cases, when the clipping region passes exactly through the middle point. As for me, implementing this algorithm is not worth the trouble.

Answer 4

You just have to subdivide any piece of the line that is not completely inside or outside the rectangle.

What I don't get is what this algorithm is good for. At least looking at your code, it would be much faster to just calculate the intersections with the rectangle.