Stefan Röttger

// get overlap ratio of range a with range b
double compute_overlap(double a1,double a2,double b1,double b2)
   {
   if (a1>=a2) return(0.0); // malformed range a
   if (b1>=b2) return(0.0); // malformed range b
   if (a1>=b2 || a2<=b1) return(0.0); // no overlap of a and b
   if (a1>=b1 && a2<=b2) return(1.0); // full overlap of a with b
   if (b1>=a1 && b2<=a2) return((b2-b1)/(a2-a1)); // full overlap of b with a
   if (a1<b1) return((a2-b1)/(a2-a1)); // partial left-side overlap of a with b
   if (a2>b2) return((b2-a1)/(a2-a1)); // partial right-side overlap of a with b
   return(0.0);
   }

// compute contour overlap
double compute_contour(double ds,double dt,     // pixel size
                       double value,            // scalar function value
                       double dvds,double dvdt, // partial derivative of scalar function
                       double isovalue)         // iso value to be extracted
{
   double slope;
   double range2;
   double overlap;

   // compute slope
   slope=sqrt(dvds*dvds+dvdt*dvdt);

   // map foot prints
   range2=0.25*(ds+dt)*slope;

   // compute overlap of cell with contour
   overlap=compute_overlap(value-range2,value+range2,        // foot print range of pixel,
                           isovalue-range2,isovalue+range2); // foot print range of iso line

   return(overlap);
}