title: Dual Contouring and Quadratic Error Function
body: <p>I've implemented marching cubes, dual marching cubes and adaptive marching cubes in C#, only to find out that I need dual contouring for my purposes.
  I've read all works about dual contouring and I get all but the core of the dual contouring itself: minimizing the quadratic error function (QEF).</p>
  
  <p>Right now, I'm calculating internal voxel's vertices position simply by finding the mean between all edgePoints sharing that single vertex (3 to 6 edges) and it works well, but it obviously doesn't create the internal vertices in the right places.</p>
  
  <p>Here is the piece of code I'm trying to create. Any help would be very appreciated</p>
  
  <pre><code>/// &lt;summary&gt;
      /// ORIGINAL WORK: Dual Contouring of Hermite Data by Tao Ju (remember me of a MechCommander 2 character)
      /// 2.3 Representing and minimizing QEFs
      /// The function E[x] can be expressed as the inner
      /// product (Ax-b)T (Ax-b) where A is a matrix whose rows are the
      /// normals ni and b is a vector whose entries are ni*pi.  &lt;------------ (dot product?)&gt;
      /// Typically, the quadratic function E[x] is expanded into the form
      /// E[x] = xT AT Ax - 2xT AT b + bT b (2)
      /// where the matrix AT A is a symmetric 3x3 matrix, AT b is a column
      /// vector of length three and bT b is a scalar. The advantage of this expansion
      /// is that only the matrices AT A, AT b and bT b need be stored
      /// (10 floats), as opposed to storing the matrices A and b. Furthermore,
      /// a minimizing value ˆ x for E[x] can be computed by solving
      /// the normal equations AT Aˆ x = AT b.
      /// &lt;/summary&gt;
      public Vector3 GetMinimumError(Vector3 p0, Vector3 p1, Vector3 p2, Vector3 n0, Vector3 n1, Vector3 n2)
      {
          //so, here we are. I'm creating a vector to store the final value.
          Vector3 position = Vector3.Zero;
  
          //Values of b are supposed to b (:P) three floats. The only way i know to find a float value
          //by multiplying 2 vectors is to use dot product.
          Vector3 b = new Vector3(
                 Vector3.Dot(p0, n0),
                 Vector3.Dot(p1, n1),
                 Vector3.Dot(p2, n2));
  
          //What the transpose of a vector is supposed to be?
          //I don't know, but i think should be the vector itself :)
          float bTb = Vector3.Dot(b, b); 
  
          //i create a square matrix 3x3, so i can use c# matrix transformation libraries.
          //i know i will probably have to build bigger matrix later on, but it should fit for now
          Matrix A = new Matrix(
              n0.X, n0.Y, n0.Z, 0,
              n1.X, n1.Y, n1.Z, 0,
              n2.X, n2.Y, n2.Z, 0,
              0, 0, 0, 0);
  
          //easy
          Matrix AT = Matrix.Transpose(A);
  
          //EASY
          Matrix ATA = Matrix.Multiply(AT, A);
  
          //Another intuition. Hope makes sense...
          Vector3 ATb = Vector3.Transform(b, AT);
  
          //...
          // some cool stuff about solving
          // the normal equations AT Aˆ x = AT b
          //...
  
          return position; //profit!
      }
  </code></pre>