I want to build an octree over a quarter of a sphere (for debugging and testing). The octree generator relies on the AIsosurface interface to compute the density and normal at any given point in space. For example, for a full sphere the corresponding code is:

// returns <0 if the point is inside the solid
virtual float GetDensity( float _x, float _y, float _z ) const override
{
    Float3  P = Float3_Set( _x, _y, _z );
    Float3  v = Float3_Subtract( P, m_origin );
    float   l = Float3_LengthSquared( v );
    float   d = Float_Sqrt(l) - m_radius;
    return d;
}
// estimates the gradient at the given point
virtual Float3 GetNormal( float _x, float _y, float _z ) const override
{       
    Float3  P   = Float3_Set( _x, _y, _z );
    float   d   = this->AIsosurface::GetDensity( P );
    float   Nx  = this->GetDensity( _x + 0.001f, _y, _z ) - d;
    float   Ny  = this->GetDensity( _x, _y + 0.001f, _z ) - d;
    float   Nz  = this->GetDensity( _x, _y, _z + 0.001f ) - d;
    Float3  N = Float3_Normalized( Float3_Set( Nx, Ny, Nz ) );
    return N;
}

What is a nice and fast way to compute those values when the shape is bounded by a low number of half-spaces?