Calcul formel

[FrenchLeaf]

Bonjour,

J'essaye de traiter un problème de vision par ordinateur. (Expliqué ci-apres http://www.cs.unc.edu/~marc/tutorial/node55.html).
Ce problème consiste à trouver une matrice de contrainte entre deux images représentant une meme scène sous deux points de vues (appellé courament Fundamental matrix).

Le problème que j'ai est pour identifier les facteur a3,a2,a1,a0 sur la formule suivante :

Det(alpha*F1 +(1 - alpha)*F2) == a3*lambda + a2*lambda + a1 *lambda +a0 = 0

F1 et F2 sont des matrices 3*3. Alpha est un scalaire.
Le problème est donc comment extraite les valeurs symbolique de a3,a2,a1,a0 pour qu'ensuite je puisse utiliser un solveur cubique sur le polynome obtenu (et ainsi trouvez lambda).

Je sais que c'est possible de le faire sous Maple (mais je n'ai pas maple), Donc j'ai découvert Xcas et j'essaye de trouvez la solution ...

Je pense que la méthode à utiliser est SOLVE. Mais je ne sais pas trop comment l'utiliser pour résoudre ce problème :

Je suis parti comme suit (mais apres je ne sais pas comment continuer. Pourriez vous m'aider s'il vous plait ? )
//- Definition de la fonction F :
f = det( (alpha* [ [a1,b1,c1] , [d1,e1,f1] , [g1,h1,i1] ]) * ( (1-alpha) * [ [a2,b2,c2] , [d2,e2,f2] , [g2,h2,i2] ] ) );
//- Definition de la fonction G (les facteurs que je recherche)
g = a3*lambda + a2*lambda + a1 *lambda +a0;


Merci d'avance,
Pierre.

[parisse]

Je pense qu'il y a 2 erreurs de formulation dans les données de votre problème: alpha et lambda sont identiques, et on doit trouver a_3*lambda^3+..., autrement dit a_3, etc sont les coefficients d'un polynome de degré 3 en lambda. Donc, si les matrices F1 et F2 sont quelconques
f := det( (lambda* [ [a1,b1,c1] , [d1,e1,f1] , [g1,h1,i1] ]) +( (1-lambda) * [ [a2,b2,c2] , [d2,e2,f2] , [g2,h2,i2] ] ) );
ensuite coeff(f,lambda,3) donne le coefficient de degré 3 donc a_3, etc.
Ca donne bien sur une formule contenant beaucoup de termes dans le cas où les coeffs a1, b1, etc. sont des inconnues et même si on sait calculer les solutions d'une équation de degré 3 par une formule, celle-ci serait extrêmement compliquée et sans intérêt. Il vaudra beaucoup mieux calculer les solutions numériques de l'équation en fonction des valeurs numériques des coefficients (par exemple via fsolve dans Xcas). Et à mon avis, le mieux est d'écrire directement une fonction Xcas (ou C++ utilisant giac) qui trouve 1 ou 3 racines réelles à partir des 2 matrices F1 et F2 où les coefficients sont numériques
f(F1,F2):=fsolve(det(x*F1+(1-x)*F2)=0,x);
devrait faire l'affaire (et marcher pour des matrices F1 et F2 de taille différente de 3 par la même occasion).
Sauf erreur, on peut aussi remarquer qu'en multipliant par F1^(-1) lorsque F1 est inversible, on obtient
det(x*I+(1-x)*F2*F1^-1)=0
et comme x=1 n'est pas solution, alors det(F2*F1^(-1)-x/(x-1)I)=0 donc x/(x-1) est une valeur propre de F2*F1^(-1), on peut donc trouver les valeurs de x à coup de solve(x/(x-1)=y,x) pour y dans egvl(F2*F1^(-1)).

[FrenchLeaf]

Bonjour,

Tout d'abord merci pour votre réponse. (rapide qui plus est et bien explicite).

J'ai en effet fait une erreur dans les notations. Je m'en excuse.

Mon problème est que les matrices F1 et F2 sont différente à chaque itération de mon programme. Je ne peux donc pas allez manuellement sous Xcas pour le résoudre à chaque fois.

Je voudrai une solution générique avec les termes de F1 et F2 exprimé sous forme de scalaire.
Les coefficients a_3, a_2, a_1, a_0 obtenu avec coeff(f,lambda,3) devrait me suffire. Meme si comme vous l'avez dit cela n'est pas très digeste

Concernant l'inversibilité de F1 et F2 je viens de regarder et leur determinant sont très proche de 0. Je ne peux donc pas utiliser la solution 2 que vous avez proposé.

Autre question : Pensez vous que je peux parametriser le problème ainsi ? :

Det(F1 +lambda*F2) == a_3*lambda + a_2*lambda + a_1 *lambda +a_0 = 0

Les coefficients a_3,a_2,a_1,a_0 sont alors bcp plus digeste :

Code : Tout sélectionner

f := det( [ [a1,b1,c1] , [d1,e1,f1] , [g1,h1,i1] ] +( lambda * [ [a2,b2,c2] , [d2,e2,f2] , [g2,h2,i2] ] ) );
coeff(f,lambda,3) = e2*i2*a2-e2*g2*c2-i2*b2*d2-f2*h2*a2+f2*b2*g2+h2*d2*c2
coeff(f,lambda,2) = a1*e2*i2-a1*f2*h2+e2*i1*a2-e2*g2*c1-e2*c2*g1+i2*e1*a2-i2*b2*d1-i2*d2*b1-f2*h1*a2+f2*b2*g1+f2*g2*b1-h2*f1*a2+h2*d2*c1+h2*c2*d1-e1*g2*c2-i1*b2*d2+f1*b2*g2+h1*d2*c2
coeff(f,lambda,1) = a1*e2*i1+a1*i2*e1-a1*f2*h1-a1*h2*f1-e2*g1*c1-i2*b1*d1+f2*b1*g1+h2*d1*c1+e1*i1*a2-e1*g2*c1-e1*c2*g1-i1*b2*d1-i1*d2*b1-f1*h1*a2+f1*b2*g1+f1*g2*b1+h1*d2*c1+h1*c2*d1
coeff(f,lambda,0) = a1*e1*i1-a1*f1*h1-e1*g1*c1-i1*b1*d1+f1*b1*g1+h1*d1*c1
//-- J'applique ensuite un "cubic solveur" et le tour est joué.

Ma question est désormais, puise remplacer automatique des membres dans les expressions précédente a_3... a_0 afin d'avoir une notation plus compacte.
Exemple mettre a1 en facteur pour a_1 réduirait la longeur de l'expression. Je pense Xcas dispose de méthode permettant de remplacer un terme par un nouveau ou bien de dire je souhaite factoriser tel terme.
Factor (a_1) ne donne rien, et simplify de a_1 a fait planté Xcas (la première fois).

Merci encore pour votre réponse très instructive.

Note : Je viens de trouvez une personne qui à obtenu les dit facteur pour le problème posé sous Maple :

Det( lambda*F1 + (1-lambda)*F2) == a_3*lambda + a_2*lambda + a_1 *lambda +a_0 = 0

On obtient cet indigeste code :

Code : Tout sélectionner

double  a=F1(0,0), j=F2(0,0), aa=a-j,
          b=F1(0,1), k=F2(0,1), bb=b-k,
          c=F1(0,2), l=F2(0,2), cc=c-l,
          d=F1(1,0), m=F2(1,0), dd=d-m,
          e=F1(1,1), n=F2(1,1), ee=e-n,
          f=F1(1,2), o=F2(1,2), ff=f-o,
          g=F1(2,0), p=F2(2,0), gg=g-p,
          h=F1(2,1), q=F2(2,1), hh=h-q,
          i=F1(2,2), r=F2(2,2), ii=i-r;

double a1=ee*ii-ff*hh, b1=ee*r+ii*n-ff*q-hh*o, c1=r*n-o*q;
double d1=bb*ii-cc*hh, e1=bb*r+ii*k-cc*q-hh*l, f1=r*k-l*q;
double g1=bb*ff-cc*ee, h1=bb*o+ff*k-cc*n-ee*l, i1=o*k-l*n;

float P[4]={0.f,0.f,0.f,0.f};

  P[0]=(aa*a1-dd*d1+gg*g1);
  P[1]=(aa*b1+a1*j-dd*e1-d1*m+gg*h1+g1*p);
  P[2]=(aa*c1+b1*j-dd*f1-e1*m+gg*i1+h1*p);
  P[3]=(c1*j-f1*m+i1*p);

Pierre.

[parisse]

Pour faire des factorisations partielles dans l'expression, le mieux est de sélectionner une sous-expression de la réponse à la souris, et d'appliquer la commande factor (par le menu ou le raccourci clavier Alt-F). De manière automatique, on peut aussi factoriser le coefficient de a1 en utilisant factor(coeff(...,a1,1)) et en refaisant de même avec le reste coeff(...,a1,0).

Sur la 2ème paramétrisation, la réponse est oui, puisqu'il suffit de remplacer F1,F2 par F2,F1-F2 puisque
x*F1+(1-x)*F2=F2+x*(F1-F2)

Concernant l'inclusion dans un programme C, vous avez 2 options, l'option que vous utilisez (convertir à la main le résultat de Xcas dans votre programme et appeler une routine C pour le solveur) et qui sera la plus rapide a l'execution (par ex. si vous avez 1 milliards de cas à traiter), ou l'option d'utiliser la librairie giac (il faut alors écrire le programme en C++) qui permet d'accéder à toutes les fonctions de xcas directement depuis un source C++, le code compilé sera un peu plus lent, par contre pas de risque d'erreur de recopie de coefficients et cela peut vous intéresser pour utiliser d'autres fonctionnalités. Si cette voie vous intéresse, je peux détailler un peu plus ce que ça donnerait pour votre cas.

[FrenchLeaf]

Bonjour,

Merci beaucoup pour vous explication aussi concise que d'habitude .

Je vais tester demain toutes vos suggestions.

Concernant le côté implémentation je préfère directement insérer à la main les resultats de Xcas. (ce qui toutefois demande de l'attention pour éviter les erreurs de notations entre noms de variables... mais on y arrive, à force de persévérance).

Notes pour information je trouve Xcas excellent.
Il m'a permit ces derniers jours de résoudre le problème dont je viens vous parlez ainsi que de trouvez des erreurs de notations dans un papier scientifique que je relisais (http://www.cs.ubc.ca/~mbrown/minimal/minimal.html).
Je ne manquerai pas autour de moi de montrer Xcas !
C'est dommage qu'il ne soit pas utilisé dans l'école d'ou je viens...

Je ré-écrirai un message à la fin de ce post montrant les différentes étapes pour déboucher à la solution. On ne sait jamais cela pourra peut être aider quelqu'un d'autre.

Encore merci,
Pierre M