AUnion.cpp

Go to the documentation of this file.

/**
 * @file AUnion.cpp
 * Surface Modeling Library
 * CS497
 *
 * @author William Nagel
 */

#include "AUnion.h"

REGISTER_IMPLICIT(AUnion,"BinaryOp:Algebraic:Union");

/**
 * Blending region is limited to f < r1, g < r2, f > 0, g > 0.
 */
bool AUnion::indomain(double f, double g)
{
  return f <= m_r1 && g <= m_r2 && f > 0 && g > 0;
}

bool AUnion::indomain(Intervald f, Intervald g)
{
  Intervald r1(0.0,m_r1);
  Intervald r2(0.0,m_r2);
  return (f.overlaps(r1) && g.overlaps(r2));
}

/** 
 * Computes an elliptic union blend.
 * h(f,g) = -( (f - r1)^d/r1^d + (g - r2)^d/r2^d - 1).
 * Negated such that surface is negative inside.
 * Returns smaller of f or g if either is outside its
 * blending radius.
 *
 * When f = 0 and as g approaches r2, returns zero.
 * When g = 0 and as f approaches r1, returns zero.
 */
double AUnion::h(double f, double g)
{
  if (indomain(f,g))
    return 1.0 - pow(1.0 - f/m_r1, m_d) - 
                 pow(1.0 - g/m_r2, m_d);
  else
    return (f < g) ? f : g;
}

Intervald AUnion::h(Intervald f, Intervald g)
{
  if (indomain(f,g))
    return Intervald(1.0) - 
          (Intervald(1.0) - f/Intervald(m_r1)).pow(Intervald(m_d)) - 
          (Intervald(1.0) - g/Intervald(m_r2)).pow(Intervald(m_d));
  else
    return (f < g) ? f : g;
}

/**
 * dh(f,g)/df = d f^(d-1)/r1
 */
double AUnion::hf(double f, double g)
{
  if (indomain(f,g))
    return m_d*pow(1.0 - f/m_r1,m_d-1.0)/m_r1;
  else
    return (f < g) ? 1.0 : 0.0;
}

Intervald AUnion::hf(Intervald f, Intervald g)
{
  if (indomain(f,g))
    return Intervald(m_d) * 
          (Intervald(1.0) - f/Intervald(m_r1)).pow(Intervald(m_d-1.0)) / 
           Intervald(m_r1);
  else
    return (f < g) ? Intervald(1.0) : Intervald(0.0);
}

/**
 * dh(f,g)/dg = d g^(d-1)/r2
 */
double AUnion::hg(double f, double g)
{
  if (indomain(f,g))
    return m_d*pow(1.0 - g/m_r2,m_d-1.0)/m_r2;
  else
    return (f < g) ? 0.0 : 1.0;
}

Intervald AUnion::hg(Intervald f, Intervald g)
{
  if (indomain(f,g))
    return Intervald(m_d) *
          (Intervald(1.0) - g/Intervald(m_r2)).pow(Intervald(m_d-1.0)) / 
           Intervald(m_r2);
  else
    return (f < g) ? Intervald(0.0) : Intervald(1.0);
}

/**
 * d2h(f,g)/df2 = d(d-1) f^(d-2)/r1.
 * m_d should be at least two.
 */
double AUnion::hff(double f, double g)
{
  if (indomain(f,g))
    return m_d*(m_d-1.0)*pow(1.0 - f/m_r1,m_d-2.0)/(m_r1*m_r1);
  else
    return 0.0;
}

Intervald AUnion::hff(Intervald f, Intervald g)
{
  if (indomain(f,g))
    return Intervald(m_d) * Intervald(m_d-1.0) * 
          (Intervald(1.0) - f/Intervald(m_r1)).pow(Intervald(m_d-2.0)) / 
           Intervald(m_r1*m_r1);
  else
    return Intervald(0.0);
}

/**
 * d2h(f,g)/dg2 = d(d-1) g^(d-2)/r2
 */
double AUnion::hgg(double f, double g)
{
  if (indomain(f,g))
    return m_d*(m_d-1.0)*pow(1.0 - g/m_r2,m_d-2.0)/(m_r2*m_r2);
  else
    return 0.0;
}

Intervald AUnion::hgg(Intervald f, Intervald g)
{
  if (indomain(f,g))
    return Intervald(m_d) * Intervald(m_d-1.0) *
          (Intervald(1.0) - g/Intervald(m_r2)).pow(Intervald(m_d-2.0)) / 
           Intervald(m_r2*m_r2);
  else
    return Intervald(0.0);
}

/** No mixed second derivatives.
 */
double AUnion::hfg(double f, double g)
{
  return 0.0;
}

Intervald AUnion::hfg(Intervald f, Intervald g)
{
  return Intervald(0.0);
}

---