Apparent bug in GeomFill_ConstrainedFilling

I recently ran into a problem where I was constructing a face from its four bounding edges using the GeomFill_ConstrainedFilling package. After getting many mysterious incorrect results, I have tracked it down to a bug in the sortbounds function. The error occurs when the endpoints of the edges have certain relationships, and results in edges being mis-sorted.

Below is the function with my correction. If this is useful, please let me know. If anyone thinks I've made an error, tell me that too!
Regards,
Myles Baker

static void sortbounds(const Standard_Integer nb,
Handle(GeomFill_Boundary)* bound,
Standard_Boolean* rev,
GeomFill_CornerState* stat)
{
// trier les bords (facon bourinos),
// flaguer ceux a renverser,
// flaguer les baillements au coins.
Standard_Integer i,j;
Handle(GeomFill_Boundary) temp;
rev[0] = 0;
gp_Pnt pf,pl;
gp_Pnt qf,ql;
for (i = 0; i if(!rev[i]) bound[i]->Points(pf,pl);
else bound[i]->Points(pl,pf);
for (j = i+1; j bound[j]->Points(qf,ql);
/***********************************
Old code (bug suspected)
************************************
Standard_Real d = qf.Distance(pl);
if(d temp = bound[i+1];
bound[i+1] = bound[j];
bound[j] = temp;
stat[i+1].Gap(d);
rev[i+1] = Standard_False;
}
d = ql.Distance(pl);
if(d temp = bound[i+1];
bound[i+1] = bound[j];
bound[j] = temp;
stat[i+1].Gap(d);
rev[i+1] = Standard_True;
}
*/
/**********************************
New code (Myles Baker, 6/15/04)
**********************************/
Standard_Real df = qf.Distance(pl);
Standard_Real dl = ql.Distance(pl);
if (df

if(df temp = bound[i+1];
bound[i+1] = bound[j];
bound[j] = temp;
stat[i+1].Gap(df);
rev[i+1] = Standard_False;
}
} else {
if(dl temp = bound[i+1];
bound[i+1] = bound[j];
bound[j] = temp;
stat[i+1].Gap(dl);
rev[i+1] = Standard_True;
}
}
/* end of modification */
}
}
if(!rev[nb-1]) bound[nb-1]->Points(pf,pl);
else bound[nb-1]->Points(pl,pf);
bound[0]->Points(qf,ql);
stat[0].Gap(pl.Distance(qf));

// flaguer les angles entre tangentes au coins et entre les normales au
// coins pour les bords contraints.
gp_Pnt pbid;
gp_Vec tgi, nori, tgn, norn;
Standard_Real fi, fn, li, ln;
for (i = 0; i Standard_Integer next = (i+1)%nb;
if(!rev[i]) bound[i]->Bounds(fi,li);
else bound[i]->Bounds(li,fi);
bound[i]->D1(li,pbid,tgi);
if(rev[i]) tgi.Reverse();
if(!rev[next]) bound[next]->Bounds(fn,ln);
else bound[next]->Bounds(ln,fn);
bound[next]->D1(fn,pbid,tgn);
if(rev[next]) tgn.Reverse();
Standard_Real ang = PI - tgi.Angle(tgn);
stat[next].TgtAng(ang);
if(bound[i]->HasNormals() && bound[next]->HasNormals()){
stat[next].Constraint();
nori = bound[i]->Norm(li);
norn = bound[next]->Norm(fn);
ang = nori.Angle(norn);
stat[next].NorAng(ang);
}
}
}