Quadrics are the only visible surfaces implemented in Boar, but constructive solid geometry requires us to handle two invisible primitive objects as well. The classic An Introduction to Ray Tracing mentions several researchers who solved that particular problem long ago.

There are a few pitfalls, such as the necessity to create pruned CSG subtrees on the fly, and pass them back to the caller. Here is another illustrative fragment of code from Boar (please kindly ignore the ugly memory leak; it's just a small matter of programming ...):

/*
  Intersecting CSG composites with an Axis Aligned Box and simplifying them
  (an illustrative code fragment)

  Feel free to send comments, improvements and questions to
  "hobold@vectorizer.org".

  This code fragment is Copyright 2012 by Holger Bettag, licensed under GPLv2.

  If the GPLv2 is too restrictive for you, contact me. I just want to learn
  about improvements you might find. I am sure we can agree on something that
  won't cost you a dime, and doesn't force all the obligations of the GPLv2 on
  you. The GPLv2 is merely the default.
*/


// potentially simplified objects are returned by reference
VoxClass Intersection::intersectVoxel(const GeObject*& Obj,
                                      const Vec3& Min, const Vec3& Max) const {
  VoxClass status1, status2;

  const GeObject* obj1;
  const GeObject* obj2;
  
  status1 = child1->intersectVoxel(obj1, Min, Max);
  status2 = child2->intersectVoxel(obj2, Min, Max);

  Obj = this;  // default return object is this CSG intersection
  
  // outside of one means outside of intersection
  if ((status1 == voxOutside) || (status2 == voxOutside)) {
    return voxOutside;
  }

  // both surfaces means we need to test against CSG combination
  if ((status1 == voxSurface) && (status2 == voxSurface)) {
    // there might be an opportunity for pruning here, in case a child ...
    // ... could simplify itself
    if ((child1 != obj1) || (child2 != obj2)) {
      // TODO: fix downstream memory leak ...
      Obj = new Intersection((GeObject*)obj1, (GeObject*)obj2);
      if (Obj == NULL) {
        fprintf(stderr,
                "Intersection::intersectVoxel(): no memory for pruned CSG.\n");
        exit(0);
      }
    }
    return voxSurface;
  }

  // surface of only one means we can prune CSG tree
  if (status1 == voxSurface) {
    Obj = obj1;
    return voxSurface;
  }
  if (status2 == voxSurface) {
    Obj = obj2;
    return voxSurface;
  }

  // both inside means inside of intersection
  return voxInside;
}


// potentially simplified objects are returned by reference
VoxClass Union::intersectVoxel(const GeObject*& Obj,
                               const Vec3& Min, const Vec3& Max) const {
  VoxClass status1, status2;

  const GeObject* obj1;
  const GeObject* obj2;
  
  status1 = child1->intersectVoxel(obj1, Min, Max);
  status2 = child2->intersectVoxel(obj2, Min, Max);

  Obj = this;  // default return object is this CSG union
  
  // inside of one means inside union
  if ((status1 == voxInside) || (status2 == voxInside))
    return voxInside;

  // both surfaces means we need to test against CSG combination
  if ((status1 == voxSurface) && (status2 == voxSurface)) {
    // there might be an opportunity for pruning here, in case a child ...
    // ... could simplify itself
    if ((child1 != obj1) || (child2 != obj2)) {
      // TODO: fix downstream memory leak ...
      Obj = new Union((GeObject*)obj1, (GeObject*)obj2);
      if (Obj == NULL) {
        fprintf(stderr,
                "Union::intersectVoxel(): no memory for pruned CSG.\n");
        exit(0);
      }
    }
    return voxSurface;
  }

  // surface of only one means we can prune CSG tree
  if (status1 == voxSurface) {
    Obj = obj1;
    return voxSurface;
  }
  if (status2 == voxSurface) {
    Obj = obj2;
    return voxSurface;
  }

  // outside of both means outside of union
  return voxOutside;
}