Subsections

Proprietà varie

Oltre all'algoritmo di approssimazione di una curva assegnata con una curva di Bézier (che spiegherò più avanti), ho pensato di realizzare delle funzioni che risolvono alcuni interessanti problemi.

Suddivisione di una curva

Una curva di Bézier di solito è definita con $0\leq t \leq 1$. È possibile tracciare solo una parte di curva se cambio l'intervallo di definizione del parametro: in generale posso considerare $l\leq t\leq m$, con $0\leq l\leq m
\leq1$.

Inoltre, utilizzando la tecnica del Blossom (fioritura di punti) in un caso particolare, è possibile risalire ai punti di controllo del frammento di curva (punti che ho chiamato sotto-punti di controllo).

La tecnica del Blossom consiste nell'eseguire l'algoritmo di de Casteljau variando ad ogni iterazione il valore del parametro t: la curva viene così a dipendere da un numero maggiore di variabili (in pratica diventa una superficie di $\mathbb{R}^k$ con k pari al numero di parametri diversi usati.

Se voglio calcolare un punto di controllo cj della curva ridotta con con $0\leq l\leq t\leq m\leq 1$, se i punti di controllo sono n, basta che eseguo un Blossom considerando l come valore del parametro nelle prime n-j iterazioni e m nelle ultime i. In simboli:

cj=b[l<n-i>, m<i>].

A questo servono gli ulteriori parametri della funzione Castel(...): l corrisponde a sta_t, m a end_t e l'ultimo parametro indica j (cioè quale punto di controllo si vuole calcolare).
void paper::findSubCP(punto* cp,int ncp,double sta_t, double end_t,
                      punto* out){
  for (int i=0;i<ncp;++i)
    out[i]=Castel(cp,ncp,sta_t,end_t,i);
}

Intersezione di una curva con una retta

Una proprietà importante delle curve di Bézier è che il Minimax, ossia quel rettangolo che circoscrive il poligono di controllo, contiene sempre la curva di Bézier relativa a quel poligono. Dunque calcolarlo è molto semplice (basta un ciclo sui punti di controllo alla ricerca di quelli ``estremi'').

Questa proprietà si può sfruttare per costruire un algoritmo che cerca l'intersezione della curva assegnata con un segmento che può essere introdotto dall'utente:

$\diamond$
si esegue un semplice test se il segmento è contenuto o interseca il Minimax;

$\diamond$
se il test fallisce sicuramente il segmento non interseca la curva;

$\diamond$
altrimenti si divide la curva in due parti rispetto al parametro t e si richiama ricorsivamente l'algoritmo sulle due metà;

$\diamond$
l'algoritmo termina quando il test è vero e il Minimax è molto piccolo (tipicamente un rettangolo di pochi pixel di lato); questo rettangolo individua con buona precisione il punto di intersezione del segmento con la retta.
Questo algoritmo trova tutti i punti di intersezione del segmento con la curva, in quanto nella divisione riesegue il test su entrambe le metà della curva, scartandone una solo se non ci sono intersezioni.

La funzione che effettua il test tra il segmento e il minimax è riportata qui sotto: i parametri di input sono le coordinate dell'angolo inferiore destro e superiore sinistro del minimax. I vertici del segmento sono contenuti nella variabile p_line[].

Il funzionamento è molto semplice: si controlla dapprima se uno dei vertici del segmento è interno al Minimax, altrimenti si cerca una eventuale intersezione del segmento con uno dei lati del rettangolo.

bool paper::lineAndMinimax(double min_x, double min_y, double max_x,
                           double max_y){
  double x1=p_line[0].x;
  double x2=p_line[1].x;
  double y1=p_line[0].y;
  double y2=p_line[1].y;
  // Primo vertice dentro
  if ( (x1>min_x) && (x1<max_x) && (y1>min_y) && (y1<max_y) )
     return true;
  // Secondo vertice dentro
  if ( (x2>min_x) && (x2<max_x) && (y2>min_y) && (y2<max_y) )
     return true;
  // I vertici sono fuori... 
  // La retta la vedo con una rappresentazione parametrica :
  //         x=(1-t)*x1+t*x2
  //         y=(1-t)*y1+t*y2
  // Valore del parametro che da' l'intersezione
  // con x=min_x (per x)
  double t=(min_x - x1) / (x2-x1);
  for (int i=0;i<2;++i){
    // Calcola la y dell'intersezione
    double inter_y=t*(y2-y1)+y1;
    // Il SEGMENTO interseca la retta verticale
    // x=min_x (cioe' 0<t<1) e rimane
    // all'interno del lato del rettangolo? 
    if ((t>=0) && (t<=1) && (inter_y <=max_y) && (inter_y >=min_y)) 
      return true;
    // Ripeti tutto per x=max_x
    t=(max_x - x1) / (x2-x1);
  }

  // Valore del parametro che da' l'intersezione con y=min_y (per y)
  t=(min_y - y1) / (y2-y1);
  for (int i=0;i<2;++i){
    // Calcola la x dell'intersezione
    double inter_x=t*(x2-x1)+x1;
    // Il SEGMENTO interseca la retta verticale
    // x=min_x (cioe' 0<t<1) e rimane
    // all'interno del lato del rettangolo? 
    if ((t>=0) && (t<=1) && (inter_x <=max_x) && (inter_x >=min_x)) 
      return true;
    // Ripeti tutto per y=max_y
    t=(max_y - y1) / (y2-y1);
  }
  return false;
}

La funzione che realizza l'algoritmo completo è riportata qui. Si nota che per la suddivisione della curva si usa la funzione findSubCP(...) definita prima per trovare i sotto-punti di controllo. La funzione richiama se stessa ricorsivamente ed usa il vettore int_p dichiarato all'esterno per memorizzare i punti di intersezione trovati, e l'indice int_p_i che memorizza sempre la prima posizione disponibile per l'inserimento:

void paper::xintersect(double st_t, double en_t, punto* int_p,
                       int& int_p_i){
  double min_x,max_x;
  double max_y,min_y;
  punto *subcp=new punto[ncp];
  findSubCP(cp,ncp,st_t,en_t,subcp);
  findMinimax(subcp,ncp,min_x,min_y,max_x,max_y);
  delete subcp;
  if(lineAndMinimax(min_x,min_y,max_x,max_y)){
    // If Minimax is little enough ...( 4x4 )
    if ( ((int(  (max_x-min_x)*xscale)) < 4) && 
         ((int( (-min_y+max_y)*yscale)) < 4) ) { 
        // ... then add this point ...
      displayMinimax(min_x,min_y,max_x,max_y);    
      //... but if there is another point near this? 
      for (int i=0;i<int_p_i;++i)

        if((abs(((max_x+min_x)/2 - int_p[i].x)*xscale)+ 
             abs(((max_y+min_y)/2 - int_p[i].y)*yscale) )<6)
          return; // Don't add it !

      int_p[int_p_i].x=(max_x+min_x)/2;
      int_p[int_p_i++].y=(min_y+max_y)/2;
    }
    else{ // reduce minimax
      double eps=1e-6;
      //first half    
      xintersect(st_t,(st_t+en_t+eps)/2, int_p,int_p_i); 
      //second half
      xintersect((st_t+en_t-eps)/2,en_t, int_p,int_p_i); 
    }
  }
  return;// does not intersects
}

Riduzione del grado

Nel programma è presente una funzione per ridurre il numero di punti di controllo di una curva. In generale la curva ottenuta dopo la riduzione non sarà identica all'originale, dato che la riduzione procede per estrapolazioni: la riduzione anche di un solo punto fa perdere alcune ``informazioni'' sulla curva, ma se è effettuata con un certo criterio, la curva ``ridotta'' non differirà molto dall'originale.

La tecnica di base della riduzione consiste nel considerare una curva come proveniente da una processo contrario di elevazione del grado (parto da n-1 punti di controllo e costruisco la stessa curva che ha n punti di controllo), processo che al contrario non muta la forma della curva. Nel processo di elevazione, detti $\tilde{b}_i$ i vecchi punti di controllo e bi i nuovi si ha:

\begin{displaymath}
b_i =\frac{i}{n}\tilde{b}_{i-1} + \frac{n-i}{n}\tilde{b}_i\qquad i=0,1,
\ldots, n
\end{displaymath}

e questo dà origine a due formule per la riduzione:

\begin{displaymath}
\tilde{b}_i= \frac{ nb_i-i\tilde{b}_{i-1} }{n-i} \qquad i=0, 1,\ldots, n-1
\end{displaymath}

e

\begin{displaymath}
\tilde{b}_{i-1}= \frac{ nb_i - (n-i)\tilde{b}_i} {i} \qquad i=n, n-1,\ldots, 1
\end{displaymath}

La prima delle due formule è più precisa nella prima parte del poligono di controllo, mentre la seconda e più precisa nell'ultima. Per questo nel mio programma io le uso entrambe, applicando la prima nella prima metà dei punti di controllo e la seconda nella seconda metà; se poi i punti sono dispari, il punto centrale viene calcolato come la media dei due punti che si ottengono applicando prima la prima formula e poi la seconda. Riporto qui il codice di questa funzione che in input accetta il vettore di punti di controllo old il numero di punti n ed un (eventuale) vettore (ret) dove memorizzare i punti che risultano dalla riduzione.

inline punto* degreeReduction(punto *old, int n, punto *ret=NULL){
  if (ret==NULL)
    ret=new punto[n-1];
  ret[0]=old[0];             //primo
  ret[n-2]=old[n-1];         //ultimo
  for (int i=1;i<n/2;++i){
    // parte iniziale
    ret[i].x=(n*old[i].x-i*ret[i-1].x)/(n-i);
    ret[i].y=(n*old[i].y-i*ret[i-1].y)/(n-i);
    // parte finale
    ret[n-i-2].x=(n*old[n-i-1].x-(i)*ret[n-i-1].x)/(n-i);
    ret[n-i-2].y=(n*old[n-i-1].y-(i)*ret[n-i-1].y)/(n-i);
  }
  // n dispari; fai la media dei punti centrali
  if (dispari(n-1)){
    int i=n/2-1;
  ret[i].x=( ((n*old[i].x-i*ret[i-1].x)/(n-i)) + 
             ((n*old[n-i-1].x-(i)*ret[n-i-1].x)/(n-i)) )/2;
  ret[i].y=( ((n*old[i].y-i*ret[i-1].y)/(n-i)) + 
             ((n*old[n-i-1].y-(i)*ret[n-i-1].y)/(n-i)) )/2;
  }
  return ret;
}

Se però si prova ad effettuare più di una riduzione applicando più volte questo metodo, si nota che ben presto la curva comincia ad allontanarsi molto da quella originale.

Per questo ho effettuato una piccola aggiunta che funziona bene quando si hanno molti punti di controllo: divido l'insieme dei punti in sotto-gruppi da 4 punti ed applico la riduzione a questi quattro, come se costituissero gli unici punti di una piccola curva; alla fine riunisco tutti i sottoinsiemi da 3 punti ottenuti. Così facendo uso l`algoritmo di riduzione in modo da perdere meno precisione possibile: passando da 4 punti a 3, il primo e l'ultimo punto rimangono fissi, mentre il secondo è ottenuto come media delle due formule di riduzione. La Fig. 1 mostra come l'algoritmo riesca a ridurre in un colpo circa un quarto dei punti di controllo, mantenendo la curva molto vicina all'originale.

Riporto qui il codice:

void paper::degRe(){ 
  int punti_gruppo=4;
  if (ncp<punti_gruppo+1) return; // pochi punti!
  punto *tmp=new punto[(ncp/punti_gruppo)*(punti_gruppo-1)];
  for (int i=0;i<ncp/punti_gruppo;++i)
    degreeReduction(&cp[i*(punti_gruppo)], punti_gruppo,
                    &tmp[(i*(punti_gruppo-1))]); 
  
  int i=0;
  for (;i<(ncp/punti_gruppo)*(punti_gruppo-1);++i)
    cp[i]=tmp[i];
  // qui i e' gia' stato incrementato...
  for (int j=0;j<ncp%(punti_gruppo);++j) 
     // gli ultmi [1 2 3..n-1] punti vengono ricopiati 
    cp[i++]=cp[ncp-(ncp%punti_gruppo)+j];
  delete tmp;
  ncp=i; 
  repaint(TRUE);
  ++riduzioni;
}

Tocci Giovanni 2001-09-17