La curva quindi deve essere abbastanza ``buona''.
L'algoritmo di campionamento è molto semplice. Partendo da una certa posizione iniziale cerca di ``camminare'' sulla curva, guardando sempre i pixel adiacenti e scegliendo il ``più nero''.
Ovviamente ci sono dei controlli per fare in modo che l'algoritmo vada ``sempre avanti'' e non percorra più volte lo stesso tratto di curva: quando si controllano i pixel adiacenti, vengono esclusi quelli che non sono lungo la direzione di movimento della precedente scelta; quando ciò non è possibile, si cambia direzione, escludendo sempre quella contraria alla direzione attuale (se questa fosse l'unica scelta possibile, vuol dire che si è arrivati alla fine della curva).
Riporto qui il codice:
const int max_points=300; // max numero di punti di campionamento static int passo=25; // in pixel const int soglia_MAX=100; punto *camp; // vettore di punti di campionamento int campiona(QImage *q){ int soglia=soglia_MAX; punto destro; camp=new punto[max_points]; // Cerca il pixel nero piu' a sx for(int c=0;c<q->width();++c){ bool row=false; for (int r=0;r<q->height();++r) if ( qGray(q->pixel(c,r))<soglia){ soglia=qGray(q->pixel(c,r)); camp[0].x=c; camp[0].y=r; row=true; } if (row) break; // rompe il for sulle colonne } // Cerca il pixel nero piu' a dx for(int c=q->width()-1;c>=0;--c){ bool row=false; for (int r=q->height()-1;r>=0;--r) if ( qGray(q->pixel(c,r))<soglia){ soglia=qGray(q->pixel(c,r)); destro.x=c; destro.y=r; row=true; } if (row) break; // rompe il for sulle colonne } /* Adesso cammino sulla curva: guardo gli 8 pixel intorno al mio e scelgo il piu' nero. Questo mi dira' la direzione da seguire (ascisse crescenti o decrescenti.) Quando ad un certo punto le ascisse cercheranno di decrescere, avro' finito di percorrere una catena. */ int x_cur=int(camp[0].x); int y_cur=int(camp[0].y); int x_old1=0, y_old1=0; int tmp_gray,i; int direzione=0; int cambi=0; // cambi di direzione // 3 0 // X <-- direzioni // 2 1 // bool a_destra, in_alto; for (i=passo;i<max_points*passo;++i){ a_destra = (direzione==0 || direzione==1); in_alto = (direzione==0 || direzione==3); soglia=soglia_MAX; if ( a_destra) for (int j=0;j<3;++j){ // colonna dx if ((tmp_gray=qGray(q->pixel(x_cur+1,y_cur-1+j))) <soglia){ soglia=tmp_gray; camp[i/passo].x=x_cur+1; camp[i/passo].y=y_cur-1+j; } } //for j if (!a_destra) for (int j=0;j<3;++j){ // colonna sx if ((tmp_gray=qGray(q->pixel(x_cur-1,y_cur-1+j))) <soglia){ soglia=tmp_gray; camp[i/passo].x=x_cur-1; camp[i/passo].y=y_cur-1+j; } } //for j if ((tmp_gray=qGray(q->pixel(x_cur,y_cur-1))) <soglia) // in alto if ( (y_cur-1!=y_old1) && (in_alto) ){ soglia=tmp_gray; camp[i/passo].x=x_cur; camp[i/passo].y=y_cur-1; } if ((tmp_gray=qGray(q->pixel(x_cur,y_cur+1))) <soglia) // in basso if ( (y_cur+1!=y_old1) && (!in_alto) ){ soglia=tmp_gray; camp[i/passo].x=x_cur; camp[i/passo].y=y_cur+1; } if (soglia==soglia_MAX){ if(++cambi==1) if (direzione==1) direzione=0; else direzione=(direzione+1)%4; if (cambi==2) if (direzione==0) direzione=2; else direzione=(direzione+2)%4; //no -1 altrimenti...! if (cambi>2){ return i/passo+1; } --i; // ripeti il ciclo continue; } cambi=0; x_old1=x_cur; y_old1=y_cur; x_cur=int(camp[i/passo].x); y_cur=int(camp[i/passo].y); if ((x_cur==destro.x)&&(y_cur==destro.y)){ return i/passo+1; // raggiunta la fine della curva } } return i/passo+1; }
A questo punto ho elaborato un algoritmo che sposta i punti in modo da modificare la relativa curva di Bézier in modo che si avvicini sempre di più ai punti che provengono dal campionamento. L'algoritmo è Greedy, cioè cerca sempre il punto di controllo che effettua lo spostamento maggiore e lo sposta.
Dopo aver scelto il punto da spostare, ho pensato che la direzione migliore di spostamento, quella che modifica di più la curva, è quella perpendicolare alla curva stessa, perpendicolare calcolata in un punto della curva prossimo al punto di controllo che deve essere spostato.
Per calcolare la perpendicolare in un punto di una curva di Bézier, è
sufficiente trovare la derivata in un punto della curva(che mi dà la direzione
della tangente): sfruttando la forma polinomiale di Bernstein e posto
si ha:
inline punto Bern_der(punto contr[],int punti,double t){ punto tmp; tmp.x=0; tmp.y=0; for(int j=0;j<=punti-1;++j){ double coef=BernPoly(j,punti-2,t); punto del=delta(contr,1,j); tmp.x+=(del.x)*coef; tmp.y+=(del.y)*coef; } tmp.x*=punti-1; tmp.y*=punti-1; return tmp; }mentre la funzione
delta(...)
usa una formula più generica
per restituire il
.
inline punto delta (punto contr[],int r, int i){ punto ret; ret.x=0; ret.y=0; for (int j=0;j<=r;j++){ ret.x+=binom[r][j]*pow_i(-1.0,r-j)*contr[i+j].x; ret.y+=binom[r][j]*pow_i(-1.0,r-j)*contr[i+j].y; } return ret; }
A questo punto non resta che descrivere come l'algoritmo cerca il punto da spostare:
void paper::guess(){ for (int volte=0;volte<20;volte++){ double l1,l2,maxl1=0,maxl,maxm; int maxj,j; disegna_curva(curva,*t_val,curva->size()); for (unsigned int i=10;i<points;i+=10){ // Trovo la perpendicolare alla curva punto x0y0=Bern(cp,ncp,t_val->at(i)); punto tmp2=Bern_der(cp,ncp,t_val->at(i)); double m=-tmp2.x/tmp2.y; // m_perp=-1/m_tg // Trova i due punti di controllo che // intersecano la perpendicolare // do' una certa lunghezza alla perp (~20 px) double l=20.0/(xscale*m); if (abs(m)<1){ // cout <<"m piccolo!"<<endl; l=30.0/xscale; } l1=0;j=1; // il primo e l'ultimo punto NON SI MUOVONO while (l1==0&&j<num_camp-2){ for (int h=0;h<2;++h){ double p[8]; //test intersezione p[0]=x0y0.x; p[1]=x0y0.y; //A0 p[2]=x0y0.x+l; p[3]= m*(l)+x0y0.y; //A1 p[4]=campioni[j].x/xscale; //B0 p[5]=campioni[j].y/yscale; p[6]=campioni[j+1].x/xscale; //B1 p[7]=campioni[j+1].y/yscale; if (interseca(p,l1,l2)) break; else {l=-l;l1=0;} //prova dalla parte opposta... } ++j; } if (maxl1<l1){ maxl1=l1; maxj=j; maxm=m; maxl=l; } } int y=maxj,x=maxj; // trovo il corrispondente punto "ridotto" for (int h=0;h<riduzioni;++h){ y=(x/4)*3; //blocco di appartenenza if ((x%4==1)||(x%4==2)) ++y; if (x%4==3) y+=2; x=y; } QString t("Moving"+QString::number(y)+"("+QString::number(maxj)+ ") with l1= "+ QString::number(maxl)); emit(message(t)); moveControlPoint(y,int((cp[y].x+maxl)*xscale+.5), int((cp[y].y+maxm*maxl)*yscale+.5) ); } }Non riporto qui il codice relativo alla gestione dell'interfaccia grafica, che per altro non interessa per lo studio vero e proprio delle curve di Bézier.
Nelle prossime pagine riporto le figure citate nel testo e altre che mostrano alcune caratteristiche del programma: intersezione e divisione di una curva).
Infine riporto alcuni esempi di prova dell'algoritmo di approssimazione.
Alcuni commenti finali: le curve di prova le ho disegnate con il mouse, quindi sono un po' irregolari. Le caratteristiche di continuità delle curve di Bézier alcune volte impediscono una corretta approssimazione, soprattutto nei casi in cui c'è una brusca variazione di direzione nella curva (discontinuità nelle derivate).
Tocci Giovanni 2001-09-17