CS 535 NOTES
CYRUS-BECK CLIPPING ALGORITHM

Cyrus Beck Development

The Cyrus-Beck line clipping algorithm uses vector ideas to compute a line clipped to a convex polygon. The input would be the line (L) and the polygon (a 2D array of vertices). The parametric vector form for a straight line is
P(t) = P0 + t(P1 - P0) which is a "mathematical" shorthand for the pair of equations:
x(t) = x0 + t(x1 - x0)
y(t) = y0 + t(y1 - y0)
Which could be turned into computer code easily.
Consider any edge E from the polygon (basically, an edge is formed from 2 adjacent vertices V[i] and V[i+1]). Let Q be any point on the edge E (we will use one of the vertices in the code below) and consider the vector from Q to any point on the line L. Since the point is on the line L, it has to satisify the vector form of the line equation. Thus, there is a value t with P(t) producing the point, so the vector then becomes P(t) - Q. Using the OUTER NORMAL, N, to the edge E, we can decide inside/outside for the point P(t) by examining the sign of the dot product of N and P(t)-Q.
Remark: The "classic" version of Cyrus-Beck used the INNER normal (Ninner = -Nouter). This results in a change of sign and thus slightly different comparisons in some of the if statements in the code below.
The outer normal is used here because in 3D we always use the outer normal as the surface normal, so to be consistent with the 3D case, the outer normal is used in the development here. However, the code typically found in a reference book was probably developed from the inner normal case - it still will work and can be used as is. Note that INSIDE and OUTSIDE relative to the polygon make sense since the polygon is required to be CONVEX. See below for a discussion of INSIDE and OUTSIDE. The outer normal is the one which lies in the outside. In function terms we would write this as dotProduct (N, P(t)-Q). The sign will be based on the cos(A) where A = angle between N and P(t)-Q and this angle will tell us inside/outside. Remember that N is the outer normal
  
if (dotProduct(N,P(t)-Q) > 0)
  P(t) OUTSIDE since A < 90 degrees
else if (dotProduct(N,P(t)-Q) < 0)
  P(t) INSIDE since A > 90 degrees
else
  P(t) SHOULD BE ON EDGE E 
  this is because N would be perpendicular
  to both E and P(t)-Q
Since we want the intersection of L and E, this last case is the one that we want (dotProduct (N,P(t)-Q) = 0). Unsing the above parametric vector equation for a line --- substituting in P0 + t(P1 - P0) (for the point on L) for P(t). Then we solve for t by using various vector and scalar/vector and dot product rules we get that
t = - dotProduct(N,P0-Q)/dotProduct(N,P1-P0)

A value of 0 for the denominator is bad because we cannot divide by 0, so we consider the possible ways in which the denominator could be 0

  1. N = 0. This should never happen (the outer normal should never be the 0 vector)
  2. P1 - P0 = 0. This means that we have a degenerate line - a single point. Clipping a point does not make sense - we could show the point (inside the region) or not show the point (outside). The numerator would tell us this because we have created a vector from Q to P0.
  3. dotProduct = 0. This means that edge E and line L are parallel. The question becomes is L "inside" or "outside" of the polygon. Since we have a point Q on the edge E and a point P0 on the line L and in the calculation for t we have taken the dot product of the outer normal and the vector between a point on E and a point on L, we can test for inside/outside. If the numerator > 0 then outside (N is the outer normal) and convexity means that L misses the clipping area entirely. If the numerator is < 0 then L is on the inside BUT there will be NO intersection with the edge E to compute so we can go on. If numerator == 0 then E and L are the same line and the problem is finished.
If we look at the parametric form of the line equation, we have that t < 0 or t > 1 are outside of the line segment we are working with. thus, the t value must be between 0 and 1 to be interesting (on the line segment). We can test these via the dot product of the line with the normal - the denominator in the calculation for t above.
if (dotProduct(N,P1-P0) > 0)
  EXITING - since angle A < 90 degrees
else if (dotProduct (N,P1-P0) < 0)
  ENTERING - since angle A > 90 degrees
If we repeat this for each edge in the polygon -
for each (edge E in polygon)
{ generate t value 
  test for range 0 to 1
  if (in range)
    classify as entering or exiting
}
We logically would generate two lists of t's - entering t's and exiting t's. Note that we have required that the polygon be convex, here is where it is used again. Since the polygon is CONVEX only two of these values are "real". One entering and one exiting t will generate the clipped line segment. The values of t that we want are
   t1 = MAX {t | t entering}
   t2 = MIN {t | t exiting}
Finding the MAX and the MIN is easy - just keep "current MAX" and "current MIN" in the above for loop. Thus we only need to have two t values at any time in the above for loop. If we end with tentering > texiting then the line L misses the polygon completely.

Outer Normal Calculation

This leaves us with the problem of determing the outer normal. We know that if a line has slope m then the slope of the perpendicular line is -1/m. We can use this to get the normal - the normal vector would have components (-1/m,1). Here we examine each edge of the polygon and compute the slope m for that edge. Computing m is simple given that we have two points on the line. This gives us a "candidate" normal. Because the polygon is convex, we know that all of the other vertices of the polygon must be on the INSIDE half-plane formed by the edge (as a line). Thus if we take a vector from one of the vertices on the edge with one of the other vertices of the polygon, the sign of the dot product with the candidate normal will tell us inner or outer. If inner then just multiply the normal by -1 to get the outer normal.

Cyrus Beck Code

// Cyrus-Beck 2-D Line Clipping algorithm
// for ease of coding we will treat a 2D point and a 2D vector
// as the same 
struct Point2D
{ float x,y;
}
// for simplicity we set an upper bound on the number of 
// points allowed to define a polygon - by moving to a class
// with a constructor we could make the array any size we wanted
const int MAXP = 100;
struct Polygon
{ int nPoints;
  Point2D v[MAXP];
}
const int MAXN = 100;
typedef Point2D Normal[MAXN];

// compute the outer normals.  
// note that this requires that the polygon be convex
// to always work
void CalcNormals (Polygon p, Normal & n)
{ int i,j,k;
  point2D v;
  for (i = 0; i < p.nPoints; i++)
  { j = (i+1)%p.nPoints;
    k = (i+2)%p.nPoints;
    // make vector be -1/mI + 1J
    n[i].x = -(p.v[j].y - p.v[i].y)/(p.v[j].x - p.v[i].x);
    n[i].y = 1.0;
    v.x = p.v[k].x - p.v[i].x;
    v.y = p.v[k].y - p.v[i].y;
    if (DotProduct (n[i],v) > 0)    // inner normal
    { n[i].x *= -1;
      n[i].y  = -1;
    }
  }
}
    
float DotProduct (Point2D v1, Point2D v2)
{ 
   return v1.x*v2.x + v1.y*v2*y;
}

void CBClip (Point2D p1, Point2D p2, Normal n, Polygon p, Boolean & visible,
                     Point2D & rp, Point2D & q)
{ float t1,t2,t,num,den;
  Point2D dirV,F;          // vectors
  int I;
  
    // start largest at smallest legal value and smallest 
    // at largest legal value
  t1 = 0.0;
  t2 = 1.0;
   // compute the direction vector
  dirV.x = p2.x - p1.x;
  dirV.y = p2.y - p1.y;

  visible = TRUE;
  i = 0;
  while ( (i < p.nPoints) && visible)
  { F.x = p1.x - p.v[i].x;
    F.y = p1.y - p.v[i].y;
    
    num  = DotProduct (n[i],F);
    den   =  DotProduct (n[i],dirV);

    if (den == 0.0)          // Parallel or Point
    { // parallel - if outside then forget the line; if inside then there are no 
      // intersections with this side 
      // but there may be with other edges, so in this case just keep going
       if (num > 0.0)
        visible = FALSE;   //   Parallel and outside or point (p1 == p2) and outside
    }
    else 
    { t = -(num/den);
      if (den < 0.0)        // entering
      { if (t <= 1.0)
          if (t > t1)
            t1 = t;
      }
      else if ( t >= 0.0)    //exiting
        if (t < t2)
          t2 = t;
    }
    i++;
  }
  if ( t1 <= t2)
  { rp.x = p1.x + t1*dirV.x;
    rp.y = p1.y + t1*dirV.y;
    q.x = p1.x + t2.dirV.x
    q.y = p1.y + t2*dirV.y
  }
  else
    visible = FALSE;
}

Liang Barsky Specialization

The Liang Barsky clipping algorithm is Cyrus Beck specialized to the "standard" rectangle. In that case the outer normals are simple - they are the standard unit vectors i,-i,j, and -j (i = (1,0), j = (0,1)). This also makes the polygon edges have simple equations and the dot product values needed for the numerator and denominator values for the four sides in the above equation for t would be very simple to compute by hand and "hard code" and so the above program can be rewritten. That rewriting results in code that looks quite different from Cyrus Beck.
void LBLineClipping (float &x0, float &y0, float &x1, float &y1, bool &visible)
{ float dx = x1 - x0;
   float dy = y1 - y0;
   float tE = 0.0, tL = 1.0;
   visible = FALSE;
   if ( dx == 0 && dy == 0 && ClipPoint(x0,y0))
     visible = TRUE;
   else 
   { if (CLIP (dx,xmin-x0,tE,tL))          // Inside relative to left side
        if (CLIP (-dx,x0-xmax,tE,tL))      // Inside relative to right side
          if (CLIP (dy,y0-ymax,tE,tL))     // Inside relative to bottom
            if (CLIP(-dy,y0-ymin,tE,tL))   // Inside relative to top
            { visible = TRUE;
               if (tL < 1.0)
               { x1 = x0 + tL*dx;
                  y1 = y0 + tL*dy;
               }
               if (tE > 0.0)
               { x0 += tE*dx;
                  y0 += tE*dy;
               }
            }
   }
}

bool ClipPoint (float x, float y)
{ // True if the point is inside the clip area .
  // you write this one - xmin <= x <= xmax and 
  // ymin <= y <= ymax

}

bool CLIP (float denom, float numer, float &tE, float &tL)
{ float t;
   if (denom > 0)
   { t = numer/denom;
      if (t > tL)
        return FALSE;
      else if (t > tE)
        tE = t;
   }
   else if (denom < 0)
   { t = numer / denom;
      if ( t < tE)
        return FALSE;
      else 
        rL = t;
   }
   else if ( numer > 0)
     return FALSE;
   return TRUE;
}