openvgfordummies/bacwar/metaballs.cpp

/**
  *
  */


#include "metaballs.h"



MetaballSystem::MetaballSystem(int ballCount0, Ball *balls0, float goo0, float threshold0) : balls(balls0), ballCount(ballCount0), goo(goo0), threshold(threshold0)
{
    minSize=MAXFLOAT;
    for(int i=0;i<ballCount;i++)
    {
        if(balls[i].size<minSize)
        {
            minSize=balls[i].size;
        }
    }
}

/**Return the metaball field's force at point 'pos'.*/

float MetaballSystem::calcForce(QVector2D pos)
{
      float force = 0;
      for(int i=0;i<ballCount;i++)
      {
          //### Formula (1)
          // TODO div = abs(ball.pos - pos)**self.goo

          QVector2D v=(balls[i].pos-pos);
          float div = (goo==2.0 ? v.lengthSquared() : powf(v.length(), goo));
          if(div != 0) //# to prevent division by zero
              force += balls[i].size / div;
          else
              force += 10000; //#"big number"
      }
      return force;
}

/**Return a normalized (magnitude = 1) normal at point 'pos'.*/
  QVector2D MetaballSystem::calcNormal(QVector2D pos)
  {
      QVector2D np(0,0);

      for(int i=0;i<ballCount;i++)
      {
          //### Formula (3)
          //div = abs(ball.pos - pos)**(2 + self.goo)
          //np += -self.goo * ball.size * (ball.pos - pos) / div

          QVector2D v=(balls[i].pos-pos);
          //float div = pow(v.length(), 2+goo);
          float div = v.lengthSquared() * (goo==2.0 ? v.lengthSquared() : powf(v.length(), goo));

          np += -goo * balls[i].size * v / div;

      }
      return np.normalized();

  }

/**Return a normalized (magnitude = 1) tangent at point 'pos'.*/
  QVector2D MetaballSystem::calcTangent(QVector2D pos)
  {
      QVector2D np = calcNormal(pos);
      //### Formula (7)
      return QVector2D(-np.y(), np.x());
  }

  /**Step once towards the border of the metaball field, return
     new coordinates and force at old coordinates.
  */
  float MetaballSystem::stepOnceTowardsBorder(QVector2D &pos)
  {
      float force = calcForce(pos);
      QVector2D np = calcNormal(pos);
      //### Formula (5)
      float stepsize = powf(minSize / threshold, 1 / goo) -
                 powf(minSize / force, 1 / goo) + 0.01;
      pos += np * stepsize;
      return force;
  }

  /**Track the border of the metaball field and return new
     coordinates.
  */
  QVector2D MetaballSystem::trackTheBorder(QVector2D pos)
  {
      float force = 9999999;
      //# loop until force is weaker than the desired threshold

      int foo=0;
      while (force > threshold)
      {
          foo++;
          if(foo>100) break;
          force = stepOnceTowardsBorder(pos);
      }
      return pos;
  }



/**Euler's method.
   The most simple way to solve differential systems numerically.
*/
QVector2D MetaballSystem::eulerTangent(QVector2D pos, float h)
{
    return pos + h * calcTangent(pos);
}


/**Runge-Kutta 2 (=mid-point).
   This is only a little more complex than the Euler's method,
   but significantly better.
*/
QVector2D MetaballSystem::rungeKutta2Tangent(QVector2D pos, float h)
{
    return pos + h * calcTangent(pos + calcTangent(pos) * h / 2);
}