// Fit a label to a path3.
// Author: Jens Schwaiger

import three;
private real eps=100*realEpsilon;

triple nextnormal(triple p, triple q)
{
 triple nw=p-(dot(p,q)*q);
 return abs(nw) < 0.0001 ? p : unit(nw);
}

triple[] firstframe(path3 p, triple optional=O)
{
 triple[] start=new triple[3];
 start[0]=dir(p,reltime(p,0));
 start[1]=(abs(cross(start[0],optional)) < eps) ? perp(start[0]) :
   unit(cross(start[0],optional));
 start[2]=cross(start[0],start[1]);
 return start;
}

// Modification of the bishop frame construction contained in
// space_tube.asy (from Philippe Ivaldi's modules).
// For noncyclic path3s only
triple[] nextframe(path3 p, real reltimestart, triple[] start, real
                  reltimeend, int subdiv=20)
{
 triple[][] bf=new triple[subdiv+1][3];
 real lg=reltimeend-reltimestart;
 if(lg <= 0) return start;
 bf[0]=start;
 int n=subdiv+1;
 for(int i=1; i < n; ++i)
   bf[i][0]=dir(p,reltime(p,reltimestart+(i/subdiv)*lg));

 for(int i=1; i < n; ++i) {
   bf[i][1]=nextnormal(bf[i-1][1],bf[i][0]);
   bf[i][2]=cross(bf[i][0],bf[i][1]);
 }
 return bf[subdiv];
}

surface labelpath(string s, path3 p, real angle=90, triple optional=O)
{
 real Cos=Cos(angle);
 real Sin=Sin(angle);
 path[] text=texpath(Label(s,(0,0),Align,basealign));
 text=scale(1/(max(text).x-min(text).x))*text;
 path[][] decompose=containmentTree(text);

 real[][] xpos=new real[decompose.length][2];
 surface sf;
 for(int i=0; i < decompose.length; ++i) {// Identify positions along x-axis
   xpos[i][1]=i;
   real pos0=0.5(max(decompose[i]).x+min(decompose[i]).x);
   xpos[i][0]=pos0;
 }
 xpos=sort(xpos); // sort by distance from 0;
 triple[] pos=new triple[decompose.length];
 real lg=arclength(p);
 //create frames;
 triple[] first=firstframe(p,optional);
 triple[] t0=first;
 real tm0=0;
 triple[][] bfr=new triple[decompose.length][3];
 for(int j=0; j < decompose.length; ++j) {
   bfr[j]=nextframe(p,tm0,t0,xpos[j][0]);
   tm0=xpos[j][0]; t0=bfr[j];
 }
 transform3[] mt=new transform3[bfr.length];
 for(int j=0; j < bfr.length; ++j) {
   triple f2=Cos*bfr[j][1]+Sin*bfr[j][2];
   triple f3=Sin*bfr[j][1]+Cos*bfr[j][2];
   mt[j]=shift(relpoint(p,xpos[j][0]))*transform3(bfr[j][0],f2,f3);
 }
 for(int j=0; j < bfr.length; ++j) {
   path[] dc=decompose[(int) xpos[j][1]];
   pair pos0=(0.5(max(dc).x+min(dc).x),0);
   sf.append(mt[j]*surface(scale(lg)*shift(-pos0)*dc));
 }
 return sf;
}