Friendlyskies AoI Forum

SlugFiller · October 8, 2009, 3:36 am

Sorry for the cross-post from the SF forum, but this wasn't getting much attention there, and was probably OT for that forum anyway.

I used to think I was pretty good with math. Oh well...

For the past few days I've been working on implementing an object type which supports Doo-Sabin subdivision surfaces. I did initially have some difficulty since it converts edges to faces, which sort of interferes with face-vertex parameter consistency when crossing the edge. I was able to resolve that with a correct split-weld paradigm which keeps separate sides of an edge as separate faces without changing the Doo-Sabin geometry.

I now seem to be having problem with the non-uniform subdivision, though. For one thing, it produces alot more faces than a trimesh set on approximation with comparable geometry, while leaving a surface that is far less smooth, and contains far larger faces (I'm not even sure how that's possible). It also takes alot of time doing that.

What more, as it approaches the tolerance distance, it starts creating many visible creases in what should be a perfectly smooth surface. This I don't get. It's probably some bug in the welding, but I can't even begin to guess what.

Even after adding geometry-based welding, it just ends up all wrong.

So, my questions are:
1. Do you know a good free hoster where I could post the code for review? (~148k uncompressed)
2. Once I take care of the prior, can anyone here help me figure this thing out?

Jep · October 8, 2009, 5:38 am

Can't help with math, don't even understand what you're on about smile but code can be posted right here.
The BB code tag {code}{/code} can be used. Replace {} with [].

Code:

This is some code.

And of course you can link to your code that you've put on some other site.

SlugFiller · October 8, 2009, 7:53 am

And of course you can link to your code that you've put on some other site.

Hence Q1.

The code is ~148k, I believe it's above the PunBB max. I can narrow it down to the subdivision calculations, though:

Code:

  private void calcMixParams(int[][] dest, double[][] destweights, int ofs, double srcweight1, double srcweight2, int[] sources1, double[] weights1, int[] sources2, double[] weights2)
  {
    TreeMap<Integer, Double> map = new TreeMap<Integer, Double>();
    for (int i = 0; i < sources1.length; i++)
      {
        double toset = weights1[i]*srcweight1;
        Double val = map.get(sources1[i]);
        if (val != null)
          toset += val;
        map.put(sources1[i], toset);
      }
    for (int i = 0; i < sources2.length; i++)
      {
        double toset = weights2[i]*srcweight2;
        Double val = map.get(sources2[i]);
        if (val != null)
          toset += val;
        map.put(sources2[i], toset);
      }
    dest[ofs] = new int[map.size()];
    destweights[ofs] = new double[map.size()];
    int count = 0;
    for (Map.Entry<Integer, Double> entry : map.entrySet())
      {
        dest[ofs][count] = entry.getKey();
        destweights[ofs][count] = entry.getValue();
        count++;
      }
  }

  private CenterGonMesh calcSubdivisions(double tol)
  {
    Vec3[] verts = new Vec3[this.vertices.size()];
    double[][] vertweights = new double[this.vertices.size()][1];
    int[][] vertparamsources = new int[this.vertices.size()][1];
    double[][] vertparamweights = new double[this.vertices.size()][1];
    double[] vertsmooth = new double[this.vertices.size()];
    boolean[] vertpure = new boolean[this.vertices.size()];
    int[] edges1 = new int[this.edges.size()];
    int[] edges2 = new int[this.edges.size()];
    double[] edgesmooth = new double[this.edges.size()];
    boolean[] edgepure = new boolean[this.edges.size()];
    int[][] faces = new int[this.faces.size()][];
    int[] facesources = new int[this.faces.size()];
    Vec3[] centers = new Vec3[this.faces.size()];
    int[][] centerparamsources = new int[this.faces.size()][1];
    double[][] centerparamweights = new double[this.faces.size()][1];
    double tol_edgesplit = 1e-5;//tol*tol;
    double tol_vertsplit = 1e-5;//tol_edgesplit*4.0;
    double tol_edgeweld = 1e-5;//tol_edgesplit*0.25;
    double tol_geom = 1e-5;//tol*tol;
    for (int i = 0; i < verts.length; i++)
      {
        verts[i] = this.vertices.get(i).getPos();
        vertparamsources[i][0] = i;
        vertparamweights[i][0] = 1.0;
        vertsmooth[i] = this.vertices.get(i).getSmoothness();
        vertpure[i] = true;
      }
    for (int i = 0; i < edges1.length; i++)
      {
        Edge e = this.edges.get(i);
        edges1[i] = e.getFrom().getIndex();
        edges2[i] = e.getTo().getIndex();
        edgesmooth[i] = e.getSmoothness();
        edgepure[i] = true;
      }
    for (int i = 0; i < faces.length; i++)
      {
        Face f = this.faces.get(i);
        faces[i] = new int[f.getEdgeList().size()];
        int count = 0;
        for (Edge e : f.getEdgeList())
          faces[i][count++] = e.getIndex();
        facesources[i] = i;
        centers[i] = f.getCenter();
        centerparamsources[i][0] = i+verts.length;
        centerparamweights[i][0] = 1.0;
      }
    if (this.smoothing == QUADRIC_SMOOTHING)
      {
        for (int s = 0; s < /*MAX_SUBDIVISIONS*/3 && ((edges1.length >> MAX_SUBDIVISIONS) == 0); s++)
          {
            int facevertcount = 0;
            for (int i = 0; i < faces.length; i++)
              facevertcount += faces[i].length;
            Vec3[] facevert_newpos = new Vec3[facevertcount];
            int[] facevert_nexts = new int[facevertcount];
            int[] facevert_verts = new int[facevertcount];
            int[] facevert_faces = new int[facevertcount];
            int[][] facevert_edges = new int[facevertcount][4];
            Vec3[][] edge_newpos = new Vec3[edges1.length][2];
            int[][] edge_faceverts = new int[edges1.length][4];
            int[] vert_linkcounts = new int[verts.length];
            Vec3[] vert_newpos = new Vec3[verts.length];
            for (int i = 0; i < verts.length; i++)
              vert_newpos[i] = new Vec3();
            for (int i = 0; i < edges1.length; i++)
              edge_faceverts[i][0] = edge_faceverts[i][1] = edge_faceverts[i][2] = edge_faceverts[i][3] = -1;
            int count = 0;
            for (int i = 0; i < faces.length; i++)
              {
                int first = count;
                int lastedge = faces[i][faces[i].length-1];
                for (int j = 0; j < faces[i].length; j++)
                  {
                    int thisedge = faces[i][j];
                    int common = edges1[lastedge];
                    if (edges1[lastedge] != edges1[thisedge] && edges1[lastedge] != edges2[thisedge])
                      common = edges2[lastedge];
                    facevert_newpos[count] = centers[i].times(QUADRIC_WEIGHT_CENTER).plus(verts[common].times(QUADRIC_WEIGHT_VERTEX));
                    vert_newpos[common].add(facevert_newpos[count]);
                    vert_linkcounts[common]++;
                    facevert_nexts[count] = count+1;
                    facevert_verts[count] = common;
                    facevert_faces[count] = i;
                    facevert_edges[count][0] = lastedge;
                    facevert_edges[count][1] = thisedge;
                    if (edges1[lastedge] == common)
                      {
                        facevert_edges[count][2] = 0;
                        edge_faceverts[lastedge][2] = count;
                      }
                    else
                      {
                        facevert_edges[count][2] = 1;
                        edge_faceverts[lastedge][3] = count;
                      }
                    if (edges1[thisedge] == common)
                      {
                        facevert_edges[count][3] = 0;
                        edge_faceverts[thisedge][0] = count;
                      }
                    else
                      {
                        facevert_edges[count][3] = 1;
                        edge_faceverts[thisedge][1] = count;
                      }
                    count++;
                    lastedge = thisedge;
                  }
                facevert_nexts[count-1] = first;
              }
            for (int i = 0; i < edges1.length; i++)
              {
                edge_newpos[i][0] = verts[edges1[i]].times(QUADRIC_WEIGHT_EDGE_1).plus(verts[edges2[i]].times(QUADRIC_WEIGHT_EDGE_2));
                edge_newpos[i][1] = verts[edges2[i]].times(QUADRIC_WEIGHT_EDGE_1).plus(verts[edges1[i]].times(QUADRIC_WEIGHT_EDGE_2));
                if (edgesmooth[i] > 0.0)
                  {
                    if (edge_faceverts[i][0] >= 0 && edge_faceverts[i][2] >= 0)
                      {
                        edge_newpos[i][0].scale(1.0-edgesmooth[i]);
                        edge_newpos[i][0].add(facevert_newpos[edge_faceverts[i][0]].plus(facevert_newpos[edge_faceverts[i][2]]).times(0.5*edgesmooth[i]));
                      }
                    if (edge_faceverts[i][1] >= 0 && edge_faceverts[i][3] >= 0)
                      {
                        edge_newpos[i][1].scale(1.0-edgesmooth[i]);
                        edge_newpos[i][1].add(facevert_newpos[edge_faceverts[i][1]].plus(facevert_newpos[edge_faceverts[i][3]]).times(0.5*edgesmooth[i]));
                      }
                  }
                vert_newpos[edges1[i]].add(edge_newpos[i][0]);
                vert_linkcounts[edges1[i]]++;
                vert_newpos[edges2[i]].add(edge_newpos[i][1]);
                vert_linkcounts[edges2[i]]++;
              }
            for (int i = 0; i < verts.length; i++)
              {
                if (vert_linkcounts[i] <= 0 || vertsmooth[i] <= 0.0)
                  {
                    vert_newpos[i] = verts[i];
                    continue;
                  }
                vert_newpos[i].scale(vertsmooth[i]/vert_linkcounts[i]);
                vert_newpos[i].add(verts[i].times(1.0-vertsmooth[i]));
              }
            int[] facevert_welds = new int[facevertcount];
            boolean[][] edge_splits = new boolean[edges1.length][2];
            for (int i = 0; i < facevertcount; i++)
              {
                if (facevert_newpos[i].distance2(vert_newpos[facevert_verts[i]]) < tol_vertsplit)
                  {
                    facevert_welds[i] = 4;
                    vert_linkcounts[facevert_verts[i]]--;
                  }
                else
                  facevert_welds[i] = 0;
              }
            for (int i = 0; i < edges1.length; i++)
              {
                if ((edge_faceverts[i][0] >= 0 && facevert_welds[edge_faceverts[i][0]] == 4) || (edge_faceverts[i][2] >= 0 && facevert_welds[edge_faceverts[i][2]] == 4))
                  {
                    edge_splits[i][0] = false;
                    vert_linkcounts[edges1[i]]--;
                  }
                else if (edge_newpos[i][0].distance2(vert_newpos[edges1[i]]) < tol_edgesplit)
                  {
                    edge_splits[i][0] = false;
                    vert_linkcounts[edges1[i]]--;
                  }
                else
                  edge_splits[i][0] = true;
                if ((edge_faceverts[i][1] >= 0 && facevert_welds[edge_faceverts[i][1]] == 4) || (edge_faceverts[i][3] >= 0 && facevert_welds[edge_faceverts[i][3]] == 4))
                  {
                    edge_splits[i][1] = false;
                    vert_linkcounts[edges2[i]]--;
                  }
                else if (edge_newpos[i][1].distance2(vert_newpos[edges2[i]]) < tol_edgesplit)
                  {
                    edge_splits[i][1] = false;
                    vert_linkcounts[edges2[i]]--;
                  }
                else
                  edge_splits[i][1] = true;
              }
            boolean[] vert_welds = new boolean[verts.length];
            for (int i = 0; i < verts.length; i++)
              vert_welds[i] = (verts[i].distance2(vert_newpos[i]) < tol_geom);
            for (int i = 0; i < edges1.length; i++)
              {
                if (edge_splits[i][0] && vert_welds[edges1[i]])
                  {
                    Vec3 proj = verts[edges2[i]].minus(verts[edges1[i]]);
                    Vec3 test = edge_newpos[i][0].minus(verts[edges1[i]]);
                    proj.scale(test.dot(proj)/proj.length2());
                    if (test.distance2(proj) < tol_geom)
                      {
                        edge_splits[i][0] = false;
                        vert_linkcounts[edges1[i]]--;
                      }
                  }
                if (edge_splits[i][1] && vert_welds[edges2[i]])
                  {
                    Vec3 proj = verts[edges1[i]].minus(verts[edges2[i]]);
                    Vec3 test = edge_newpos[i][1].minus(verts[edges2[i]]);
                    proj.scale(test.dot(proj)/proj.length2());
                    if (test.distance2(proj) < tol_geom)
                      {
                        edge_splits[i][0] = false;
                        vert_linkcounts[edges1[i]]--;
                      }
                  }
              }
            for (int i = 0; i < facevertcount; i++)
              {
                if (facevert_welds[i] == 4)
                  continue;
                if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]] || edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                  continue;
                if (vert_welds[facevert_verts[i]])
                  {
                    facevert_welds[i] = 4;
                    vert_linkcounts[facevert_verts[i]]--;
                  }
              }
            boolean canstop = true;
            for (int i = 0; i < verts.length; i++)
              if (vert_linkcounts[i] > 0)
                canstop = false;
            if (canstop)
              break;
            for (int i = 0; i < facevertcount; i++)
              if (facevert_welds[i] == 4)
                vert_newpos[facevert_verts[i]] = verts[facevert_verts[i]];
            for (int i = 0; i < verts.length; i++)
              if (vert_linkcounts[i] <= 0)
                vert_newpos[i] = verts[i];
            for (int i = 0; i < facevertcount; i++)
              {
                if (facevert_welds[i] == 4)
                  continue;
                if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]] && facevert_newpos[i].distance2(edge_newpos[facevert_edges[i][0]][facevert_edges[i][2]]) < tol_edgeweld)
                  facevert_welds[i] |= 1;
                if (edge_splits[facevert_edges[i][1]][facevert_edges[i][3]] && facevert_newpos[i].distance2(edge_newpos[facevert_edges[i][1]][facevert_edges[i][3]]) < tol_edgeweld)
                  facevert_welds[i] |= 2;
              }
            int vertcount = verts.length;
            int edgecount = edges1.length;
            int facecount = faces.length;
            for (int i = 0; i < edges1.length; i++)
              {
                if (edge_splits[i][0])
                  {
                    vertcount++;
                    edgecount++;
                  }
                if (edge_splits[i][1])
                  {
                    vertcount++;
                    edgecount++;
                  }
              }
            for (int i = 0; i < facevertcount; i++)
              {
                if (facevert_welds[i] == 0 || facevert_welds[i] == 1 || facevert_welds[facevert_nexts[i]] == 0 || facevert_welds[facevert_nexts[i]] == 2)
                  {
                    edgecount++;
                    facecount++;
                  }
                if (facevert_welds[i] == 4)
                  continue;
                if (facevert_welds[i] != 0)
                  {
                    if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]] && edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                      {
                        edgecount++;
                        facecount++;
                      }
                    continue;
                  }
                vertcount++;
                edgecount++;
                if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]] || edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                  {
                    edgecount++;
                    facecount++;
                  }
              }
            Vec3[] newverts = new Vec3[vertcount];
            int[][] newvertparamsources = new int[vertcount][];
            double[][] newvertparamweights = new double[vertcount][];
            double[] newvertsmooth = new double[vertcount];
            boolean[] newvertpure = new boolean[vertcount];
            count = 0;
            int[] vert_newverts = new int[verts.length];
            for (int i = 0; i < verts.length; i++)
              {
                newverts[count] = vert_newpos[i];
                newvertparamsources[count] = vertparamsources[i];
                newvertparamweights[count] = vertparamweights[i];
                newvertsmooth[count] = 1.0;
                newvertpure[count] = true;
                vert_newverts[i] = count;
                count++;
              }
            int[][] edge_newverts = new int[edges1.length][2];
            for (int i = 0; i < edges1.length; i++)
              {
                if (edge_splits[i][0])
                  {
                    newverts[count] = edge_newpos[i][0];
                    calcMixParams(newvertparamsources, newvertparamweights, count, QUADRIC_WEIGHT_EDGE_1, QUADRIC_WEIGHT_EDGE_2, vertparamsources[edges1[i]], vertparamweights[edges1[i]], vertparamsources[edges2[i]], vertparamweights[edges2[i]]);
                    newvertsmooth[count] = 1.0;
                    newvertpure[count] = false;
                    edge_newverts[i][0] = count;
                    count++;
                  }
                if (edge_splits[i][1])
                  {
                    newverts[count] = edge_newpos[i][1];
                    calcMixParams(newvertparamsources, newvertparamweights, count, QUADRIC_WEIGHT_EDGE_2, QUADRIC_WEIGHT_EDGE_1, vertparamsources[edges1[i]], vertparamweights[edges1[i]], vertparamsources[edges2[i]], vertparamweights[edges2[i]]);
                    newvertsmooth[count] = 1.0;
                    newvertpure[count] = false;
                    edge_newverts[i][1] = count;
                    count++;
                  }
              }
            int[] facevert_newverts = new int[facevertcount];
            for (int i = 0; i < facevertcount; i++)
              {
                if (facevert_welds[i] == 0)
                  {
                    newverts[count] = facevert_newpos[i];
                    calcMixParams(newvertparamsources, newvertparamweights, count, QUADRIC_WEIGHT_VERTEX, QUADRIC_WEIGHT_CENTER, vertparamsources[facevert_verts[i]], vertparamweights[facevert_verts[i]], centerparamsources[facevert_faces[i]], centerparamweights[facevert_faces[i]]);
                    newvertsmooth[count] = 1.0;
                    newvertpure[count] = false;
                    facevert_newverts[i] = count;
                    count++;
                  }
              }
            int[] newedges1 = new int[edgecount];
            int[] newedges2 = new int[edgecount];
            double[] newedgesmooth = new double[edgecount];
            boolean[] newedgepure = new boolean[edgecount];
            count = 0;
            int[][] edge_newedges = new int[edges1.length][3];
            for (int i = 0; i < edges1.length; i++)
              {
                int v0 = vert_newverts[edges1[i]];
                int v1 = edge_newverts[i][0];
                int v2 = edge_newverts[i][1];
                int v3 = vert_newverts[edges2[i]];
                if (edge_splits[i][0])
                  newedges1[count] = v1;
                else
                  newedges1[count] = v0;
                if (edge_splits[i][1])
                  newedges2[count] = v2;
                else
                  newedges2[count] = v3;
                newedgesmooth[count] = edgesmooth[i];
                newedgepure[count] = edgepure[i];
                edge_newedges[i][0] = count;
                count++;
                if (edge_splits[i][0])
                  {
                    newedges1[count] = v0;
                    newedges2[count] = v1;
                    newedgesmooth[count] = 1.0;
                    newedgepure[count] = edgepure[i];
                    edge_newedges[i][1] = count;
                    count++;
                  }
                if (edge_splits[i][1])
                  {
                    newedges1[count] = v2;
                    newedges2[count] = v3;
                    newedgesmooth[count] = 1.0;
                    newedgepure[count] = edgepure[i];
                    edge_newedges[i][2] = count;
                    count++;
                  }
              }
            int[][] facevert_newedges = new int[facevert_verts.length][4];
            for (int i = 0; i < facevertcount; i++)
              {
                facevert_newedges[i][0] = facevert_newedges[i][1] = facevert_newedges[i][2] = facevert_newedges[i][3] = -1;
                int j = facevert_nexts[i];
                if (facevert_welds[i] == 0 || facevert_welds[i] == 1 || facevert_welds[j] == 0 || facevert_welds[j] == 2)
                  {
                    if (facevert_welds[i] == 0)
                      newedges1[count] = facevert_newverts[i];
                    else if (facevert_welds[i] == 1)
                      newedges1[count] = edge_newverts[facevert_edges[i][0]][facevert_edges[i][2]];
                    else if (facevert_welds[i] == 4)
                      newedges1[count] = vert_newverts[facevert_verts[i]];
                    else
                      newedges1[count] = edge_newverts[facevert_edges[i][1]][facevert_edges[i][3]];
                    if (facevert_welds[j] == 0)
                      newedges2[count] = facevert_newverts[j];
                    else if (facevert_welds[j] == 2)
                      newedges2[count] = edge_newverts[facevert_edges[j][1]][facevert_edges[j][3]];
                    else if (facevert_welds[j] == 4)
                      newedges2[count] = vert_newverts[facevert_verts[j]];
                    else
                      newedges2[count] = edge_newverts[facevert_edges[j][0]][facevert_edges[j][2]];
                    newedgesmooth[count] = 1.0;
                    newedgepure[count] = false;
                    facevert_newedges[i][0] = count;
                    count++;
                    newedgesmooth[edge_newedges[facevert_edges[j][1]][0]] = 1.0;
                  }
                if (facevert_welds[i] == 4) {
                  if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]])
                    facevert_newedges[i][2] = edge_newedges[facevert_edges[i][0]][facevert_edges[i][2]+1];
                  if (edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                    facevert_newedges[i][1] = edge_newedges[facevert_edges[i][1]][facevert_edges[i][3]+1];
                  continue;
                }
                if (facevert_welds[i] != 0)
                  {
                    if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]] && edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                      {
                        newedges1[count] = edge_newverts[facevert_edges[i][0]][facevert_edges[i][2]];
                        newedges2[count] = edge_newverts[facevert_edges[i][1]][facevert_edges[i][3]];
                        newedgesmooth[count] = 1.0;
                        newedgepure[count] = false;
                        facevert_newedges[i][facevert_welds[i]] = count;
                        count++;
                      }
                    else if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]])
                      facevert_newedges[i][1] = edge_newedges[facevert_edges[i][0]][facevert_edges[i][2]+1];
                    else
                      facevert_newedges[i][2] = edge_newedges[facevert_edges[i][1]][facevert_edges[i][3]+1];
                    continue;
                  }
                if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]])
                  newedges1[count] = edge_newverts[facevert_edges[i][0]][facevert_edges[i][2]];
                else
                  newedges1[count] = vert_newverts[facevert_verts[i]];
                newedges2[count] = facevert_newverts[i];
                newedgesmooth[count] = 1.0;
                newedgepure[count] = false;
                facevert_newedges[i][2] = count;
                count++;
                if (edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                  newedges1[count] = edge_newverts[facevert_edges[i][1]][facevert_edges[i][3]];
                else if (!edge_splits[facevert_edges[i][0]][facevert_edges[i][2]])
                  {
                    facevert_newedges[i][1] = facevert_newedges[i][2];
                    continue;
                  }
                else
                  newedges1[count] = vert_newverts[facevert_verts[i]];
                newedges2[count] = facevert_newverts[i];
                newedgesmooth[count] = 1.0;
                newedgepure[count] = false;
                facevert_newedges[i][1] = count;
                count++;
              }
            int[][] newfaces = new int[facecount][];
            int[] newfacesources = new int[facecount];
            Vec3[] newcenters = new Vec3[facecount];
            int[][] newcenterparamsources = new int[facecount][];
            double[][] newcenterparamweights = new double[facecount][];
            count = 0;
            int firstfv = 0;
            for (int i = 0; i < faces.length; i++)
              {
                int counter = faces[i].length;
                for (int j = 0; j < faces[i].length; j++)
                  {
                    if (facevert_newedges[firstfv+j][3] != -1)
                      counter++;
                    if (facevert_newedges[firstfv+j][0] == -1)
                      {
                        if (facevert_welds[firstfv+j] == 4 && facevert_newedges[firstfv+j][1] != -1)
                          counter++;
                        if (facevert_welds[facevert_nexts[firstfv+j]] == 4 && facevert_newedges[facevert_nexts[firstfv+j]][2] != -1)
                          counter++;
                      }
                  }
                newfaces[count] = new int[counter];
                counter = 0;
                for (int j = 0; j < faces[i].length; j++)
                  {
                    if (facevert_newedges[firstfv+j][3] != -1)
                      newfaces[count][counter++] = facevert_newedges[firstfv+j][3];
                    if (facevert_newedges[firstfv+j][0] == -1)
                      {
                        if (facevert_welds[firstfv+j] == 4 && facevert_newedges[firstfv+j][1] != -1)
                          newfaces[count][counter++] = facevert_newedges[firstfv+j][1];
                        newfaces[count][counter++] = edge_newedges[facevert_edges[firstfv+j][1]][0];
                        if (facevert_welds[facevert_nexts[firstfv+j]] == 4 && facevert_newedges[facevert_nexts[firstfv+j]][2] != -1)
                          newfaces[count][counter++] = facevert_newedges[facevert_nexts[firstfv+j]][2];
                      }
                    else
                      newfaces[count][counter++] = facevert_newedges[firstfv+j][0];
                  }
                newfacesources[count] = i;
                newcenters[count] = centers[i];
                newcenterparamsources[count] = centerparamsources[i];
                newcenterparamweights[count] = centerparamweights[i];
                count++;
                firstfv += faces[i].length;
              }
            for (int i = 0; i < facevertcount; i++)
              {
                if (facevert_welds[i] == 3)
                  {
                    if (facevert_newedges[i][3] != -1)
                      {
                        newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][facevert_edges[i][3]+1], facevert_newedges[i][3], edge_newedges[facevert_edges[i][0]][facevert_edges[i][2]+1]};
                        newfacesources[count] = facevert_faces[i];
                        newcenters[count] = vert_newpos[facevert_verts[i]];
                        newcenterparamsources[count] = vertparamsources[facevert_verts[i]];
                        newcenterparamweights[count] = vertparamweights[facevert_verts[i]];
                        count++;
                      }
                  }
                else if (facevert_welds[i] == 1)
                  {
                    if (edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                      {
                        newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][facevert_edges[i][3]+1], facevert_newedges[i][1], edge_newedges[facevert_edges[i][0]][facevert_edges[i][2]+1]};
                        newfacesources[count] = facevert_faces[i];
                        newcenters[count] = vert_newpos[facevert_verts[i]];
                        newcenterparamsources[count] = vertparamsources[facevert_verts[i]];
                        newcenterparamweights[count] = vertparamweights[facevert_verts[i]];
                        count++;
                      }
                  }
                else if (facevert_welds[i] == 2)
                  {
                    if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]])
                      {
                        newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][facevert_edges[i][3]+1], facevert_newedges[i][2], edge_newedges[facevert_edges[i][0]][facevert_edges[i][2]+1]};
                        newfacesources[count] = facevert_faces[i];
                        newcenters[count] = vert_newpos[facevert_verts[i]];
                        newcenterparamsources[count] = vertparamsources[facevert_verts[i]];
                        newcenterparamweights[count] = vertparamweights[facevert_verts[i]];
                        count++;
                      }
                  }
                else if (facevert_welds[i] == 0)
                  {
                    if (edge_splits[facevert_edges[i][0]][facevert_edges[i][2]] || edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                      {
                        if (!edge_splits[facevert_edges[i][0]][facevert_edges[i][2]])
                          newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][facevert_edges[i][3]+1], facevert_newedges[i][1], facevert_newedges[i][2]};
                        else if (!edge_splits[facevert_edges[i][1]][facevert_edges[i][3]])
                          newfaces[count] = new int[] {facevert_newedges[i][1], facevert_newedges[i][2], edge_newedges[facevert_edges[i][0]][facevert_edges[i][2]+1]};
                        else
                          newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][facevert_edges[i][3]+1], facevert_newedges[i][1], facevert_newedges[i][2], edge_newedges[facevert_edges[i][0]][facevert_edges[i][2]+1]};
                        newfacesources[count] = facevert_faces[i];
                        newcenters[count] = vert_newpos[facevert_verts[i]];
                        newcenterparamsources[count] = vertparamsources[facevert_verts[i]];
                        newcenterparamweights[count] = vertparamweights[facevert_verts[i]];
                        count++;
                      }
                  }
                if (facevert_newedges[i][0] == -1)
                  continue;
                int j = facevert_nexts[i];
                if (facevert_newedges[i][1] == -1)
                  newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][0], facevert_newedges[j][2], facevert_newedges[i][0]};
                else if (facevert_newedges[j][2] == -1)
                  newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][0], facevert_newedges[i][0], facevert_newedges[i][1]};
                else
                  newfaces[count] = new int[] {edge_newedges[facevert_edges[i][1]][0], facevert_newedges[j][2], facevert_newedges[i][0], facevert_newedges[i][1]};
                newfacesources[count] = facevert_faces[i];
                if (edge_splits[facevert_edges[i][1]][0])
                  {
                    if (edge_splits[facevert_edges[i][1]][1])
                      {
                        newcenters[count] = edge_newpos[facevert_edges[i][1]][0].plus(edge_newpos[facevert_edges[i][1]][1]);
                        newcenters[count].scale(0.5);
                      }
                    else
                      newcenters[count] = vert_newpos[edges2[facevert_edges[i][1]]];
                  }
                else
                  {
                    if (edge_splits[facevert_edges[i][1]][1])
                      newcenters[count] = vert_newpos[edges1[facevert_edges[i][1]]];
                    else
                      {
                        newcenters[count] = vert_newpos[facevert_verts[i]].plus(vert_newpos[facevert_verts[j]]);
                        newcenters[count].scale(0.5);
                      }
                  }
                calcMixParams(newcenterparamsources, newcenterparamweights, count, 0.5, 0.5, vertparamsources[facevert_verts[i]], vertparamweights[facevert_verts[i]], vertparamsources[facevert_verts[j]], vertparamweights[facevert_verts[j]]);
                count++;
              }
            verts = newverts;
            vertparamsources = newvertparamsources;
            vertparamweights = newvertparamweights;
            vertsmooth = newvertsmooth;
            vertpure = newvertpure;
            edges1 = newedges1;
            edges2 = newedges2;
            edgesmooth = newedgesmooth;
            edgepure = newedgepure;
            faces = newfaces;
            facesources = newfacesources;
            centers = newcenters;
            centerparamsources = newcenterparamsources;
            centerparamweights = newcenterparamweights;
          }
      }
    // TODO: Actually do other subdivisions
    return generateSubdivided(verts, vertparamsources, vertparamweights, vertsmooth, vertpure, edges1, edges2, edgesmooth, edgepure, faces, facesources, centers, centerparamsources, centerparamweights);
  }

  private CenterGonMesh generateSubdivided(Vec3[] verts, int[][] vertparamsources, double[][] vertparamweights, double[] vertsmooth, boolean[] vertpure, int[] edges1, int[] edges2, double[] edgesmooth, boolean[] edgepure, int[][] faces, int[] facesources, Vec3[] centers, int[][] centerparamsources, double[][] centerparamweights)
  {
    CenterGonMesh result = new CenterGonMesh();
    result.smoothing = this.smoothing;
    result.triangulation = this.triangulation;
    Vertex[] vertlist = new Vertex[verts.length];
    Edge[] edgelist = new Edge[edges1.length];
    Face[] facelist = new Face[faces.length];
    for (int i = 0; i < vertlist.length; i++)
      {
        vertlist[i] = result.createVertex();
        vertlist[i].setPos(verts[i]);
        vertlist[i].setSmoothness(vertsmooth[i]);
        if (vertpure[i])
          vertlist[i].setPure();
      }
    for (int i = 0; i < edgelist.length; i++)
      {
        edgelist[i] = result.createEdge(vertlist[edges1[i]], vertlist[edges2[i]]);
        edgelist[i].setSmoothness(edgesmooth[i]);
        if (edgepure[i])
          edgelist[i].setPure();
      }
    for (int i = 0; i < facelist.length; i++)
      {
        Edge[] faceedges = new Edge[faces[i].length];
        for (int j = 0; j < faces[i].length; j++)
          faceedges[j] = edgelist[faces[i][j]];
        facelist[i] = result.createFace(faceedges);
        facelist[i].setCenter(centers[i]);
      }
    result.copyTextureAndMaterial(this);
    ParameterValue[] values = result.getParameterValues();
    if (values != null)
      {
        for (int i = 0; i < values.length; i++)
          {
            if (values[i] instanceof VertexParameterValue)
              {
                VertexParameterValue value = (VertexParameterValue) values[i];
                double oldval[] = value.getValue();
                double newval[] = new double[vertlist.length+facelist.length];
                for (int j = 0; j < vertlist.length; j++)
                  {
                    newval[j] = 0.0;
                    for (int k = 0; k < vertparamsources[j].length; k++)
                      newval[j] += oldval[vertparamsources[j][k]]*vertparamweights[j][k];
                  }
                for (int j = 0; j < facelist.length; j++)
                  {
                    newval[j+vertlist.length] = 0.0;
                    for (int k = 0; k < centerparamsources[j].length; k++)
                      newval[j+vertlist.length] += oldval[centerparamsources[j][k]]*centerparamweights[j][k];
                  }
                values[i] = new VertexParameterValue(newval);
              }
            else if (values[i] instanceof FaceParameterValue)
              {
                FaceParameterValue value = (FaceParameterValue) values[i];
                double oldval[] = value.getValue();
                double newval[] = new double[facelist.length];
                for (int j = 0; j < newval.length; j++)
                  newval[j] = oldval[facesources[j]];
                values[i] = new FaceParameterValue(newval);
              }
            else if (values[i] instanceof FaceVertexParameterValue)
              {
                FaceVertexParameterValue value = (FaceVertexParameterValue) values[i];
                double newval[][] = new double[facelist.length][];
                for (int j = 0; j < newval.length; j++)
                  {
                    newval[j] = new double[faces[j].length];
                    newval[j][0] = 0.0;
                    for (int k = 0; k < centerparamsources[j].length; k++)
                      if (centerparamsources[j][k] < this.vertices.size())
                        {
                          Vertex v = this.vertices.get(centerparamsources[j][k]);
                          for (FaceVertex fv : this.faces.get(facesources[j]).getFaceVertexList())
                            if (fv.getVertex() == v)
                              {
                                newval[j][0] += value.getValue(facesources[j], fv.getIndex()+1)*centerparamweights[j][k];
                                break;
                              }
                        }
                      else
                        newval[j][0] += value.getValue(facesources[j], 0)*centerparamweights[j][k];
                    for (FaceVertex fvout : facelist[j].getFaceVertexList())
                      {
                        int index = fvout.getIndex()+1;
                        int[] sources = vertparamsources[fvout.getVertex().getIndex()];
                        double[] weights = vertparamweights[fvout.getVertex().getIndex()];
                        newval[j][index] = 0.0;
                        for (int k = 0; k < sources.length; k++)
                          if (sources[k] < this.vertices.size())
                            {
                              Vertex v = this.vertices.get(sources[k]);
                              for (FaceVertex fv : this.faces.get(facesources[j]).getFaceVertexList())
                                if (fv.getVertex() == v)
                                  {
                                    newval[j][index] += value.getValue(facesources[j], fv.getIndex()+1)*weights[k];
                                    break;
                                  }
                            }
                          else
                            newval[j][index] += value.getValue(facesources[j], 0)*weights[k];
                      }
                  }
                values[i] = new FaceVertexParameterValue(newval);
              }
          }
        result.setParameterValues(values);
      }
    return result;
  }

Vidiot · October 8, 2009, 11:51 am

I´m not sure this will be helpful, but anyway:

-> Take a look at PME source code
and in case that doesn´t help:
-> Ask Francois G. (But I guess he´ll see this here sooner or later)

Good luck
Harald

SlugFiller · October 8, 2009, 1:57 pm

-> Take a look at PME source code

PME uses Catmull-Clark, while I am trying to achieve Doo-Sabin.
I did look at PME for some concepts, but the final geometry is too different between the two. For one thing, Catmull-Clark transforms the mesh into quads even at a single step, while Doo-Sabin maintains most of the same topology.

Vidiot · October 9, 2009, 3:36 pm

I know!
I have them all in Hexagon (Catmull, Doo-S, Bezier Interpolation, Loop Subdivison, Butterfly...), and read about the differences (not that *I* do understand everything about it).

Francois has had something in mind in the active phase of developenent of PME,
something like implementing another subdivision scheme so the user has the choice.
But that was long ago...
Maybe he´s working on the same stuff right now for "clay app".

Cheers
Harald

fpgui · October 9, 2009, 3:59 pm

Several considerations (please note I'm not familiar with Doo-Sabin):

-it seems to take place on quad or ngon meshes. Moreover, since even on quadmeshes it produces triangular faces, you need an ngon face structure. Of course I suggest you use PolyMesh or a half-edge structure of your own which is very easy to work with.

-Subdivision is tricky to implement. Set up test cases, like a single face, 4 faces making a square, a tetrahedron or a cube. Dump the result from first subdivision, print it out, draw the structure on a sheet of paper along with edge/vertices indices. Check the whole structure manually. Check subdiv vertices coordinates. Compare with coordinates you computed manually. Set vertices coordinates to reasonable (eg unsmoothed) positions to first check subdivision topology. Check everything. Write a method which checks if a mesh is sound, run it after each subdivision. It may sound tedious but guess what : it truly is. I have a notebook with print outs of mesh structures, coordinates, checks and everything.

-I don't know how you display the subdivided mesh, but use shade smoothing. All subdivided trimeshes are shade smoothed in the end and that helps a lot. Othewise you need an unreasonable number of faces for the mesh to appear smooth.

-Quad mesh subdivision generally ends up with more faces than trimeshes subdivs for a given surface precision (at least Catmull-Clark). Once a face reaches surface precision you shouldn't subdivide it again. This is adaptive subdivision and it is much more complex with Catmull-Clark quad mesh subdiv than loop or butterfly trimesh subdiv, which one of the reasons why Peter chose trimesh modeling rather than quadmesh. With Doo-Sabin, though, I don't have a clue.

-You're (almost) on your own. Your code is far too long to be read and understood unless we spend a lot of time studying it. Working on a mesh structure is something one get used to through practice, like adding vertices, faces or edges. Subdivision combines these difficulties with other complicated stuff. You should first test the algorithms you use to modify mesh topology, like add a vertex, draw mesh. Add a face, draw mesh, etc. Once you're confident in these algorithms, you can use them safely.

I certainly wouldn't like to pull your interest away from this topic. It's certainly very rewarding when it works. But it's always helpful to split a complicated algorithm into manageable pieces that can be checked once at a time. I can tell you I spent a lot of time only to [i]recode[\i] mesh subdivision for Clay.

Even if subdivision code is too long/complex for us to read and criticize doesn't mean we're not willing to help. Feel free to post advances/problems on this thread.

fpgui · October 9, 2009, 4:44 pm

Vidiot wrote:
I know!
Francois has had something in mind in the active phase of developenent of PME,
something like implementing another subdivision scheme so the user has the choice.

I don't remember what it was about... Or was it the stuff about setting some vertices normals? After some second thought, it turned out to be a lot more complex than standard subdivision for IMO marginal user benefit.

SlugFiller · October 9, 2009, 6:33 pm

-it seems to take place on quad or ngon meshes.

N-gons, or rather, "center" gons. This would make more sense if I could post the entire plugin somewhere, but I just don't happen to have any free webspace lying around.
Still, if you'll note either the structure enumeration at the beginning of calcSubdivisions, or the generation process during generateSubdivided, you'll note the use of multi-edge faces, and get/set of centers.

Set up test cases, like a single face, 4 faces making a square, a tetrahedron or a cube.

The full plugin creates a test-cube. The first subdivision works "as advertised", but the second produces non-planar quads.
In re-reading the Doo-Sabin paper I noticed they have an adjusted version for correct deplanarization (is that a word?), but they are very unspecific regarding their formulas, and I don't think I can implement it in my shifted centers paradigm.
The problem is, their system is intended for untextured models, and I had to make some adjustments to make it work with face-vertex mapping on the original mesh. It is equivalent to the linear formulation, but I need to apply some extra smoothing, and I just can't figure out how to make that work.

-I don't know how you display the subdivided mesh, but use shade smoothing. All subdivided trimeshes are shade smoothed in the end and that helps a lot. Othewise you need an unreasonable number of faces for the mesh to appear smooth.

I do, although the code for that is in a function called getRawRenderingMesh which is not in the above code (due to size constraints). Again, it would be alot easier to explain if you could play with what I have, but I just don't have a host.

Once a face reaches surface precision you shouldn't subdivide it again. This is adaptive subdivision

I have non-uniform subdivision, and this is precisely where things get tricky. This is governed by tol_geom in the code above, and once I set it, the whole subdivision geometry goes to hell. In that much, I can't tell if it's the underlying math, or if I have some bug in the welding process. I do know that at least with non-uniform turned off, the mesh has the right topology, but still has visible creases in areas where it should be smooth.

You're (almost) on your own. Your code is far too long to be read and understood unless we spend a lot of time studying it.

And that's just the subdivision part... I suppose I could have spent more time adding comments, but it would make more sense if you could see how it looks compiled. Hmm... Maybe I should check if my photobucket account is still active.

You should first test the algorithms you use to modify mesh topology, like add a vertex, draw mesh. Add a face, draw mesh, etc. Once you're confident in these algorithms, you can use them safely.

Already tested and retested these. I even use some graph-theory stuff to ensure faces have consistent orientation, and I've tested it by attempting to create the cube with some faces backwards (it correctly auto-flipped the faces that needed flipping).

My problem is with the actual math of the subdivision. It should produce smooth quadric beziers between any 4 centroids (though it doesn't create quads, it does make sure all vertices have a connectivity of 4 after one iteration), but it simply doesn't do that. I could "cheat" by dividing a single step, then outputting beziers for the renderer (although there is no RTBezier, go figure), but that would prevent more interesting uses like subdividing the mesh for editing purposes.

Even if subdivision code is too long/complex for us to read and criticize doesn't mean we're not willing to help. Feel free to post advances/problems on this thread.

Quite frankly I don't know what else to post. I could really use someone that has some understanding of the underlying math. I've been calculating pretty much until my head caught on fire, and I can't seem to get it to work again.
The best I could come up with is that it isn't working now because the resulting surface is not consistent with quadric beziers as predicted, but I'm just not sure how to correct this.

Maybe I'll at least try posting some images of my problem.

SlugFiller · October 9, 2009, 6:39 pm

Okay, my PhotoBucket account is down. So... Recommendations for an image hoster?

fpgui · October 10, 2009, 3:42 pm

When I encountered problems with adaptive subdivision, I used a 3x3 grid mesh and marked the center face for subdivision, leaving the other faces undivided. Then I could see how the adaptive subdivision would work and if it as free of bugs. This is probably an obvious advice and you don't seem to be in need of these, but that's all I can suggest until you can post pictures (picasa ?)

SlugFiller · October 11, 2009, 9:15 am

Yup, Picasa Web is good. Good thing I've got an active gmail account.

Well, here's a pic of various uniform subdivision steps with errors pointed out:
http://picasaweb.google.com/lh/photo/1u … directlink

The errors are less visible in shaded from the angle taken, but they become quite apparent when you spin the object around. The errors cause visible bumps which become more notable with each additional subdivision.

For the plannarization I've considered using least mean plane or average normals, but these are CPU-expensive, and would not result in a quadric bezier.

With non-uniform subdivision activated, I get this for 2nd and 3rd steps:
http://picasaweb.google.com/lh/photo/xl … directlink

Now here there are visible wrinkles which become worse with every step. I can only hope this is due to a programatic error, and not a mathematical one.

Last edited by SlugFiller (October 11, 2009, 9:16 am)

SlugFiller · October 11, 2009, 9:38 am

Well, I've figured one thing out. Since almost each face is center-cut for texturing, my mesh requires 4 times the geometry for the same amount of detail. This means an 80000 triangles centergon mesh would be required to match the details of a 10000 triangles trimesh.

This is a bit disappointing, but I don't see any workaround for that. I suppose I could work with that...

fpgui · October 11, 2009, 12:39 pm

Concerning subdivision error, keep the minimum amount f faces to reproduce the error then look in detail at the whole process.

As for why your mesh would require more faces then a trimesh, I don't understand what you mean (beginning with that stuff on center-cut and center-gon...)

SlugFiller · October 11, 2009, 2:15 pm

Concerning subdivision error, keep the minimum amount f faces to reproduce the error then look in detail at the whole process.

I'm not sure what you mean by that. Are you talking about what is seen in the first image, or the second?

As for why your mesh would require more faces then a trimesh, I don't understand what you mean (beginning with that stuff on center-cut and center-gon...)

If you compare the wireframes to the flat shade, you'll notice that the wireframe appears to have double the "resolution" in quads.
This is partially because of the Doo-Sabin geometry. Unlike Catmull-Clark, which simply splits existing faces, Doo-Sabin actually adds faces where edges and vertices are in the original mesh. For a textured mesh, this is a problem, because a new face created from an edge would have to inherit half of its texture from each of the faces to which the edge was adjacent.
I resolved the problem by creating two faces for an edge instead of one. Each is a half of the face that should be created, set to be textured as one of the original two faces. This additional split is linear - the two new faces are planar (assuming the intended single face was planar, which it should be).
This additional split, causes most quads in the subdivided surface to be further divided linearly into 4 quads each. In other words, I use nearly quadruple the geometry for the same shape.
And of course, each quad is two triangles when triangulated, so that's roughly double the triangles for the same edge-length.
So, to get comparable apparent face size in the flat shade, I would need 8 times the number of triangles than a trimesh set on "approximating". Mind you that at 20000 triangles, you can't tell the difference between shaded and flat view for a trimesh (due to angle granularity). I would need 160000 triangles to get that sort of performance.

Friendlyskies AoI Forum

Announcement

#1 October 8, 2009, 3:36 am

Need mathematical help [Repost]

#2 October 8, 2009, 5:38 am

Re: Need mathematical help [Repost]

Code:

#3 October 8, 2009, 7:53 am

Re: Need mathematical help [Repost]

Code:

#4 October 8, 2009, 11:51 am

Re: Need mathematical help [Repost]

#5 October 8, 2009, 1:57 pm

Re: Need mathematical help [Repost]

#6 October 9, 2009, 3:36 pm

Re: Need mathematical help [Repost]

#7 October 9, 2009, 3:59 pm

Re: Need mathematical help [Repost]

#8 October 9, 2009, 4:44 pm

Re: Need mathematical help [Repost]

Vidiot wrote:

#9 October 9, 2009, 6:33 pm

Re: Need mathematical help [Repost]

#10 October 9, 2009, 6:39 pm

Re: Need mathematical help [Repost]

#11 October 10, 2009, 3:42 pm

Re: Need mathematical help [Repost]

#12 October 11, 2009, 9:15 am

Re: Need mathematical help [Repost]

#13 October 11, 2009, 9:38 am

Re: Need mathematical help [Repost]

#14 October 11, 2009, 12:39 pm

Re: Need mathematical help [Repost]

#15 October 11, 2009, 2:15 pm

Re: Need mathematical help [Repost]

Board footer