AFFIMOFF : UNE VISIONNEUSE 3D AVEC PARAMÉTRISATION ET TEXTURE

// parametrie.cpp

#include "parametrie.h"

int *iar, *jar;   // pour l'arête i : iar[i] < jar[i] 
int *fa1, *fa2;   // fa1[i] : 1-ère face de la i-ième arête
double *us, *vs;  // (u,v) pour chaque sommet

int nouvar(int i, int j) {  // l'arête (i,j) est-elle déjà connue ?
    int k, ia, ja;          // si oui, c'est laquelle ?
    if(j<i) {ia=j; ja=i;}
    else    {ia=i; ja=j;}
    for(k=0; k<na; k++)
        if(iar[k]==ia && jar[k]==ja) return k;
    return -1;        // l'arête (i,j) n'est pas trouvée
}

void xar(int i, int j) {  // échanger les arêtes i et j
    xint(&iar[i], &iar[j]);
    xint(&jar[i], &jar[j]);
    xint(&fa1[i], &fa1[j]);
    xint(&fa2[i], &fa2[j]);
}

void topologie() {
    // à partir des ns sommets xs[],ys[],zs[] et des nf faces ifa[],jfa[],kfa[]
    // on va calculer les na arêtes iar[],jar[],fa1[],fa2[], nab : nb d'arêtes
    // de bordure, naf : nb d'arêtes de la frontière et ranger les arêtes de la
    // la frontière au début des tableaux iar[],jar[],fa1[],fa2[]
    int nan, i, j, k, iq, nb, ni, nm, nn, ii, ia, ja;
    if(!etpl) {
        // relation d'Euler : na = ns + nf + c
        // c = -2 : genre 0, c = 0 : genre 1, c = 2 : genre 2, ...
        nan = ns + nf + 30;   // pour prévoir des modèles de genre élevé
        allocint(&iar, nan);
        allocint(&jar, nan);
        allocint(&fa1, nan);
        allocint(&fa2, nan);
        allocdouble(&us, ns);
        for(i=0; i<nf; i++) {
            if(ifa[i]<0 || ifa[i]>=ns) quitter("Une face est erronée");
            if(jfa[i]<0 || jfa[i]>=ns) quitter("Une face est erronée");
            if(kfa[i]<0 || kfa[i]>=ns) quitter("Une face est erronée");
            us[ifa[i]]=1.0;
            us[jfa[i]]=1.0;  // pour marquer les sommets utilisés
            us[kfa[i]]=1.0;
        }
        nn=ns;
        for(i=0; i<nn; i++) 
            if(us[i]<0.5) {  // on élimine les sommets isolés et inutiles
                ns=ns-1;
                for(j=i; j<ns; j++) {xs[j]=xs[j+1]; ys[j]=ys[j+1]; zs[j]=zs[j+1];}
                for(j=0; j<nf; j++) {
                    if(ifa[j]>i) ifa[j]=ifa[j]-1;   // dans : if(ifa[j]>i)
                    if(jfa[j]>i) jfa[j]=jfa[j]-1;   // ifa[j]=i est impossible
                    if(kfa[j]>i) kfa[j]=kfa[j]-1;   // le sommet i n'est pas utilisé
                }
            }
        delete [] us;
        // on renseigne : na, fa1[i] et fa2[i], i=0 à na ( fa2[i]=-1 en bordure )
        iar[0] = ifa[0]; jar[0] = jfa[0];
        if(jar[0]<iar[0]) xint(&iar[0],&jar[0]);
        fa1[0] = 0;
        iar[1] = jfa[0]; jar[1] = kfa[0];
        if(jar[1]<iar[1]) xint(&iar[1],&jar[1]);
        fa1[1] = 0;
        iar[2] = kfa[0]; jar[2] = ifa[0];
        if(jar[2]<iar[2]) xint(&iar[2],&jar[2]);
        fa1[2] = 0;
        na = 3;  
        for(i=0; i<nan; i++) fa2[i]=-1;
        for(i=1; i<nf; i++) {
            if(na+3>nan) quitter("Problème de mémoire pour les arêtes");
            iq=nouvar(ifa[i], jfa[i]);
            if(iq>=0) if(fa2[iq]<0) fa2[iq] = i; else quitter("Une arête est erronée");
            else {
                iar[na] = ifa[i]; jar[na] = jfa[i];
                if(jar[na]<iar[na]) xint(&iar[na],&jar[na]);
                fa1[na] = i;
                na=na+1;
            }
            iq=nouvar(jfa[i], kfa[i]);
            if(iq>=0) if(fa2[iq]<0) fa2[iq] = i; else quitter("Une arête est erronée");
            else {
                iar[na] = jfa[i]; jar[na] = kfa[i];
                if(jar[na]<iar[na]) xint(&iar[na],&jar[na]);
                fa1[na] = i;
                na=na+1;
            }
            iq=nouvar(kfa[i], ifa[i]);
            if(iq>=0) if(fa2[iq]<0) fa2[iq] = i; else quitter("Une arête est erronée");
            else {
                iar[na] = kfa[i]; jar[na] = ifa[i];
                if(jar[na]<iar[na]) xint(&iar[na],&jar[na]);
                fa1[na] = i;
                na=na+1;
            }
        }
        nab=0;      // nombre d'arêtes de bordure pour toutes les ouvertures
        for(i=0; i<na; i++) if(fa2[i]==-1) nab=nab+1;
        j=0;
        for(i=0; i<nab; i++) {        // pour chaque i-ième arête de bordure
            while(fa2[j]!=-1) j=j+1;
            xar(i,j);                 
            j=j+1;    // les nab premières arêtes sont des arêtes de bordure
        }         
        nn = 0;
        if(nab>0) {
            nb = 0;                  // nb : nombre d'ouvertures
            nm = 0;
            do {
                ni = nb;             // y a-t-il une nouvelle ouverture à explorer ?
                for(i=0; i<nab; i++) if(fa2[i]==-1) {nb=nb+1; fa2[i]=nb; ii=i; ia=iar[i]; ja=jar[i]; break;}
                if(nb!=ni) {                   // il y a encore une nouvelle ouverture
                    k=0;
                    do {                         // les arêtes de la nb-ième ouverture
                        for(i=0; i<nab; i++) {   // sont marquées avec fa2[i] = nb
                            if(fa2[i]!=-1) continue;
                            if(iar[i]==ja) {fa2[i]=nb; ja=jar[i]; break;}
                            if(jar[i]==ja) {fa2[i]=nb; ja=iar[i]; break;}
                        }
                        k=k+1;
                        if(k>nab) {
                            quitter("La topologie du modèle n'est pas admissible.");
                        }
                    } while(ja!=ia);     // est-on revenu au début de cette ouverture ?
                    nn=0;
                    for(i=0; i<nab; i++) if(fa2[i]==nb) nn=nn+1;
                    if(nn>nm) {
                        nm=nn;           // la frontière sera l'ouverture ayant le plus d'arêtes
                        for(i=0; i<nab; i++) if(fa2[i]==nb) fa2[i]=-1;  // annulation temporaire
                        fa2[ii]=nb;
                        xar(0,ii);
                        ja=jar[0];
                        nn=1;
                        do {                // pour mettre la frontière au début et dans l'ordre
                            for(i=0; i<nab; i++) {
                                if(fa2[i]!=-1) continue;
                                if(iar[i]==ja) {fa2[i]=nb; ja=jar[i]; xar(nn,i); nn=nn+1; break;}
                                if(jar[i]==ja) {fa2[i]=nb; ja=iar[i]; xar(nn,i); nn=nn+1; break;}
                            }
                        } while(ja!=ia);
                        naf=nn;                       // nombre d'arêtes de la frontière
                    }                  // naf = 0 : aucune ouverture et pas de frontière
                }
            } while(nb!=ni);
            for(i=0; i<nab; i++) fa2[i]=-1;      // pour toutes les arêtes de bordure
        }
        etpl=true;
    }
}

void weights(int k, double *wij, double *wji) {
/*
// autre choix possible : wij=wji=distance(sommet i, sommet j)
    double p1[3], p2[3];
    int i, j;
    i = iar[k];
    j = jar[k];
    p1[0]=xg(i); p1[1]=yg(i); p1[2]=zg(i);
    p2[0]=xg(j); p2[1]=yg(j); p2[2]=zg(j);
    *wij=dist2p(p1,p2);
    *wji=dist2p(p1,p2);
// autre choix possible : wij=wji=1.0
    *wij=1.0;
    *wji=1.0;
*/
// calcul des poids des arêtes par la méthode "Mean Value"
// "Mean Value" est dite : 2d reproductible
    int i, j, i1, i2, i3;
    double gamma, delta, p1[3], p2[3], p3[3], xij, xji;
    i = iar[k];
    j = jar[k];
    i1=ifa[fa1[k]]; i2=jfa[fa1[k]]; i3=kfa[fa1[k]]; 
    p1[0]=xg(i1); p1[1]=yg(i1); p1[2]=zg(i1);
    p2[0]=xg(i2); p2[1]=yg(i2); p2[2]=zg(i2);
    p3[0]=xg(i3); p3[1]=yg(i3); p3[2]=zg(i3);
    if(i==i1) gamma=angletri(p1,p2,p3);
    if(i==i2) gamma=angletri(p2,p3,p1);
    if(i==i3) gamma=angletri(p3,p1,p2);
    xij=tan(gamma/2.0);
    if(j==i1) gamma=angletri(p1,p2,p3);
    if(j==i2) gamma=angletri(p2,p3,p1);
    if(j==i3) gamma=angletri(p3,p1,p2);
    xji=tan(gamma/2.0);
    if(fa2[k]>=0) {
        i1=ifa[fa2[k]]; i2=jfa[fa2[k]]; i3=kfa[fa2[k]]; 
        p1[0]=xg(i1); p1[1]=yg(i1); p1[2]=zg(i1);
        p2[0]=xg(i2); p2[1]=yg(i2); p2[2]=zg(i2);
        p3[0]=xg(i3); p3[1]=yg(i3); p3[2]=zg(i3);
        if(i==i1) delta=angletri(p1,p2,p3);
        if(i==i2) delta=angletri(p2,p3,p1);
        if(i==i3) delta=angletri(p3,p1,p2);
        xij=xij+tan(delta/2.0);
        if(j==i1) delta=angletri(p1,p2,p3);
        if(j==i2) delta=angletri(p2,p3,p1);
        if(j==i3) delta=angletri(p3,p1,p2);
        xji=xji+tan(delta/2.0);
    }
    p1[0]=xg(i); p1[1]=yg(i); p1[2]=zg(i);
    p2[0]=xg(j); p2[1]=yg(j); p2[2]=zg(j);
    *wij=xij/(2.0*dist2p(p1,p2));
    *wji=xji/(2.0*dist2p(p1,p2));
}

void reallocijv(int *ma, int **ia, int **ja, double **va) {
        reallocint(ia,*ma,*ma+5000);
        reallocint(ja,*ma,*ma+5000);     // pour augmenter la taille de ia, ja, va
        reallocdouble(va,*ma,*ma+5000);
        *ma=*ma+5000;
}

void vadd(int *ma, int *na, int **ia, int **ja, double **va, int i, int j, double v) {
    // pour faire dans (ia, ja, va) l'équivalent de a[i][j]=a[i][j]+v
    for(int k=0; k<*na; k++) if((*ia)[k]==i && (*ja)[k]==j) {(*va)[k]=(*va)[k]+v; return;}
    // si on désactive la boucle for() ci-dessus, cela marche aussi,
    // peut-être un peu plus vite et peut-être avec na final plus grand
    if(*na >= *ma) reallocijv(ma, ia, ja, va);
    (*ia)[*na]=i; (*ja)[*na]=j; (*va)[*na]=v; *na=*na+1;
}

void parametrise() {
    // à partir de la topologie on va calculer la paramétrisation us[],vs[]
    // c'est à dire les ns points homologues des sommets dans le plan (u,v)
    // on choisit ici une méthode linéaire à frontière fixe circulaire et à coordonnées barycentriques
    char info[512];
    int i, j, k, ni;
    int *ks;  // ks[i] : -1 sur la frontière et j pour le j-ième sommet intérieur 
    double alfa, wij, wji, pi=3.14159265358979;
    if(!epuv) {
        if(!etpl) topologie();
        if(naf==0) {
            if(!cdit) { 
                strcpy(info,"Le modèle n'a pas d'ouverture.\n\nLa paramétrisation est impossible.\n\n");
                strcat(info,"Les affichages de la frontière, du plan (u,v)\n et du modèle texturé ");
                strcat(info,"avec une image\n ne sont pas possibles.");
                inform(info);
            }
            cdit=true;
            return;
        }
        /* paramétrisation initiale */
        allocdouble(&us, ns);
        allocdouble(&vs, ns);
        /* placement dans (u,v) des naf sommets de la frontière */
        alfa = (2.0*pi)/(double)naf;
        k=jar[0];
        us[k]=cos(alfa);
        vs[k]=sin(alfa);
        for(i=1; i<naf; i++) {
            if(iar[i]==k) k=jar[i];  else k=iar[i];
            if(us[k]==0.0) us[k]=cos(alfa*(double)(i+1));
            if(vs[k]==0.0) vs[k]=sin(alfa*(double)(i+1));
        }
        /*  on a maintenant naf sommets placés et fixes dans (u,v)  *
         *  et ns-naf sommets intérieurs et mis au centre à placer  *
         *  ce qui fait en tout 2(ns-naf) inconnues us[i] et vs[i]  */
        if(naf==ns) {epuv=true; return;}
        allocint(&ks, ns);
        ks[0]=1;
        do {    // ici on marque les sommets reliés au sommet n° 0
            k=0;
            for(i=0; i<na; i++) {
                if(ks[iar[i]]!=0 && ks[jar[i]]==0){ks[jar[i]]=1; k=1;}
                if(ks[jar[i]]!=0 && ks[iar[i]]==0){ks[iar[i]]=1; k=1;}
            }
        } while(k!=0);
        // y a-t-il des parties isolées ?
        k=0;
        for(i=0; i<ns; i++) if(ks[i]==0) k=1;
        if(k==1) {
            inform("Paramétrisation impossible : \n\nle modèle a plusieurs parties disjointes.");
            delete [] ks;
            return;
        }
        // ks[i]=-1 : sommet fixe, ks[i]=j : j-ième sommet inconnu dans (u,v)
        for(i=0; i<ns; i++) ks[i]=0;
        for(i=0; i<naf; i++) {ks[iar[i]]=-1; ks[jar[i]]=-1;}  // sommets déjà fixés dans (u,v)
        j=0;
        for(i=0; i<ns; i++) if(ks[i]==0) {ks[i]=j; j=j+1;}    // sommets inconnus dans (u,v)
        ni=ns-naf;
        if(j>ni) ni=j;          // Cela ne devrait pas arriver. Mais parfois ...
        // Si cela arrive, doit-on le signaler ? ou le refuser ? ou l'accepter ?
        if(ni==0) {delete [] ks; epuv=true; return;}  // aucun sommet inconnu
        // Il y a ni inconnues parmi les ns u[i] et ni autres parmi les v[i]
        // les ni équations en u et celles en v sont indépendantes        
/*
        // option 1 : système linéaire avec matrice pleine et méthode directe
        if(ni > 600) {
            inform("Paramétrisation impossible : \n\nle modèle est trop grand.");
            delete [] ks;
            return;
        }
        double **aa, *bb;
        aa = new double*[ni];  // allocation dynamique de : a[ni][ni]
        for(j=0; j<ni; j++) aa[j] = new double[ni]; 
        for(i=0; i<ni; i++) {
            for(j=0; j<ni; j++) aa[i][j]=0.0;
        }
        allocdouble(&bb, ni);
        for(k=0; k<na; k++) {        // mise en place du
            weights(k, &wij, &wji);  // premier système linéaire
            i = iar[k];            
            j = jar[k];
            if(ks[i]>=0) {
                aa[ks[i]][ks[i]] = aa[ks[i]][ks[i]] + wij;
                if(ks[j]<0) bb[ks[i]] = bb[ks[i]] + wij*us[j];
                else {aa[ks[i]][ks[j]] = aa[ks[i]][ks[j]] - wij;}
            }
            if(ks[j]>=0) {
                aa[ks[j]][ks[j]] = aa[ks[j]][ks[j]] + wji;
                if(ks[i]<0) bb[ks[j]] = bb[ks[j]] + wji*us[i];
                else {aa[ks[j]][ks[i]] = aa[ks[j]][ks[i]] - wji;}
            }
        } 
        jordan(aa,bb,ni);  // solution pour les us[i] 
        j=0;
        for(i=0; i<ns; i++) {
            if(ks[i]>=0) {
                us[i]=bb[j];
                j=j+1;
            }
        }
        for(i=0; i<ni; i++) {
            for(j=0; j<ni; j++) aa[i][j]=0.0;
            bb[i]=0.0;
        }
        for(k=0; k<na; k++) {        // mise en place du
            weights(k, &wij, &wji);  // second système linéaire
            i = iar[k];
            j = jar[k];
            if(ks[i]>=0) {
                aa[ks[i]][ks[i]] = aa[ks[i]][ks[i]] + wij;
                if(ks[j]<0) bb[ks[i]] = bb[ks[i]] + wij*vs[j];
                else {aa[ks[i]][ks[j]] = aa[ks[i]][ks[j]] - wij;}
            }
            if(ks[j]>=0) {
                aa[ks[j]][ks[j]] = aa[ks[j]][ks[j]] + wji;
                if(ks[i]<0) bb[ks[j]] = bb[ks[j]] + wji*vs[i];
                else {aa[ks[j]][ks[i]] = aa[ks[j]][ks[i]] - wji;}
            }
        }
        jordan(aa,bb,ni);  // solution pour les vs[i] 
        j=0;
        for(i=0; i<ns; i++) {
            if(ks[i]>=0) {
                vs[i]=bb[j];
                j=j+1;
            }
        }
        for(j=0; j<ni; j++) delete [](aa[j]);
        delete [] aa;
        delete [] bb;
        // fin de l'option 1
*/
        // option 2 : système linéaire avec matrice creuse et méthodes itératives
        // solu=1 : GMRES (mgmres_st) et solu=2 : gradient conjugué (ijvsolve)
        int mm, nn, mr, solu, *ia, *ja;
        double rsq, *va, *xx, *yy;
        solu = 1;
        mm=8*ni;                // mm : dimension de ia, ja, va
        allocint(&ia, mm);      // quand a[i][j] est non nul
        allocint(&ja, mm);      // pour k=0,nn-1 on fait :
        allocdouble(&va, mm);   // ia[k]=i, ja[k]=j
        allocdouble(&xx, ni);   // va[k]=a[i][j]
        allocdouble(&yy, ni);
        nn=0;
        for(k=0; k<na; k++) {  // mise en place du premier système linéaire
            weights(k, &wij, &wji);               // pour les us[] inconnus
            i = iar[k];
            j = jar[k];
            if(ks[i]>=0) {
                vadd(&mm,&nn,&ia,&ja,&va,ks[i],ks[i],wij);
                if(ks[j]<0) yy[ks[i]] = yy[ks[i]] + wij*us[j];
                else {
                    vadd(&mm,&nn,&ia,&ja,&va,ks[i],ks[j],(-wij));
                }
            }
            if(ks[j]>=0) {
                vadd(&mm,&nn,&ia,&ja,&va,ks[j],ks[j],wji);
                if(ks[i]<0) yy[ks[j]] = yy[ks[j]] + wji*us[i];
                else {
                    vadd(&mm,&nn,&ia,&ja,&va,ks[j],ks[i],(-wji));
                }
            }
        } 
        for(i=0; i<ni; i++) xx[i]=0.0;    // pour la première itération
        if(ni==1) xx[0]=yy[0]/va[0];
        else {
            mr = ni-1; if(mr > 25) mr=25;
            if(solu==1) mgmres_st(ni, nn, ia, ja, va, xx, yy, 5, mr, 1.0E-9, 1.0E-9);
            if(solu==2) ijvsolve(ni, nn, ia, ja, va, yy, xx, &rsq);
        }
        j=0;
        for(i=0; i<ns; i++) {
            if(ks[i]>=0) {
                us[i]=xx[j];  // rangement de la solution pour les us[]
                j=j+1;
            }
        }
        for(i=0; i<ni; i++) yy[i]=0.0;
        for(k=0; k<na; k++) {  // mise en place du second système linéaire
            weights(k, &wij, &wji);              // pour les vs[] inconnus
            i = iar[k];
            j = jar[k];
            if(ks[i]>=0) {
                if(ks[j]<0) yy[ks[i]] = yy[ks[i]] + wij*vs[j];
            }
            if(ks[j]>=0) {
                if(ks[i]<0) yy[ks[j]] = yy[ks[j]] + wji*vs[i];
            }
        } 
        for(i=0; i<ni; i++) xx[i]=0.0;
        if(ni==1) xx[0]=yy[0]/va[0];
        else {
            mr = ni-1; if(mr > 25) mr=25;
            if(solu==1) mgmres_st(ni, nn, ia, ja, va, xx, yy, 5, mr, 1.0E-9, 1.0E-9);
            if(solu==2) ijvsolve(ni, nn, ia, ja, va, yy, xx, &rsq);
        }
        j=0;
        for(i=0; i<ns; i++) {
            if(ks[i]>=0) {
                vs[i]=xx[j];  // rangement de la solution pour les vs[]
                j=j+1;
            }
        }
        delete [] ia;
        delete [] ja;
        delete [] va;
        delete [] xx;
        delete [] yy;
        // fin de l'option 2

        delete [] ks;
        epuv=true;
    }
}