''fast'' Isosurface using polynom-extrapolation

The algorithm fits 2nd order polynoms with maximal 2nd derivative as specified and uses polynom extrapolation to approximate places where a root cannot be. For the exact searching of a root (once it has found a place where a root has to be) it uses a hybrid algorithm (secant method combined with bisection).

There is also a way to report if the specified 2nd derivative is too small and for finding an appriate value for it.

Example scene

//C++ code

class MyFunction: public Isosurface::IsosurfaceFunction {

public:
    MyFunction():evalcount("Function","Evaluations","") {}

    DBL evaluateAt(const Vector3& v) const {

        evalcount.eval();
        DBL x = v.a[0]; DBL y = v.a[1]; DBL z=v.a[2];
        DBL k=5, a=0.95, b=0.8;
        return pow( (x*x+y*y+z*z-a*k*k),2 )-b*( pow(z-k,2)-2*x*x )*( pow(z+k,2)-2*y*y );

    }

private:
    mutable HitStatistics evalcount;
};


Pigment3D* setUpScene_Isosurface() {

    Scene* s = new Scene();

    Pigment3DMapped check = Pigment3DMapped(new Checker());
    check.insertCol(0.5,Colour3(1,0,0));
    check.insertCol(.5,Colour3(0,0,0));

    Pigment3DMapped imap = Pigment3DMapped(new Gradient(Vector3(1,0,0)));
    imap.insertCol(0,Colour3(1,1,1));
    imap.insert(1,check.copy());

    Plane *pl = new Plane(Vector3(0,1,0));
    pl->setPigment(imap);
    pl->setFinish(Finish(0.1,0.6,0.2,1,1));
    pl->translate(Vector3(0,-2,0));
    s->addObject(pl);

    pl = new Plane(Vector3(0,0,1));
    pl->setPigment(Pigment3DConstant(Colour5(1,1,1,0,0)));
    pl->translate(Vector3(0,0,15));
    s->addObject(pl);

    MyFunction* myf = new MyFunction();
    Isosurface* myiso = new Isosurface(myf,new Sphere(6),160,true);
    myiso->scale(Vector3(.5,.5,.5));
    myiso->rotate(Vector3(0,35,0));
    myiso->translate(Vector3(0,0,5));
    myiso->setPigment(Pigment3DConstant(Colour5(1,1,1,0,0.65)));
    myiso->setInterior(Interior(1.33,0));
    myiso->setFinish(Finish(0.05,0.6,0.3,1.5,10));
    s->addObject(myiso);

    PointLight* p = new PointLight(Vector3(0,10,-8),Colour3(1),s);
    s->addLight(p);

    Material3D m = Material3D(new Texture3D(new Pigment3DConstant(Colour5(0)), new Finish(0.1,0.6,0,0,0)), new Interior() );
    s->setDefaultMaterial(m);

    s->object.prepareToRender();

    CameraPerspective* c = new CameraPerspective(s);

    return c;
};

// equvalient pov scene code
#declare check = pigment{checker rgb <1,0,0>,rgb 0};

plane{y,0

  pigment{
    gradient x
    pigment_map{
      [0 rgb 1]
      [1 check]
    }
  }
  finish{ambient 0.1 diffuse .6 reflection .2 phong 1 phong_size 1}

  translate<0,-2,0>

}

plane{z,0
  pigment{rgb 1}
  translate<0,0,15>
}

#declare k=5;
#declare a=0.95;
#declare b=0.8;

isosurface{
  function{
    pow( (x*x+y*y+z*z-a*k*k),2 )-b*( pow(z-k,2)-2*x*x )*( pow(z+k,2)-2*y*y )
  }
  contained_by{
    sphere{0,6}
  }
  max_gradient 900
  scale .5
  rotate<0,35,0>
  translate<0,0,5>
  pigment{rgbt <1,1,1,.65>}
  interior{
    ior 1.33
  }
  finish{ambient 0.05 diffuse .6 reflection .3 phong 1 phong_size 10}
}

light_source{
  <0,10,-8>,1
}

global_settings{
  max_trace_level 100 adc_bailout 1/100
}