DIMANCHE 16 SEPTEMBRE 2012
Régression circulaire

Ce billet est consacré à la régression d'un ensemble de points (x; y) par un cercle défini par son centre (a; b) et son rayon R. Il s'agit plus d'un mémo pour moi-même qu'autre chose. J'en profite pour tester également les iframe sur TinyMCE ainsi que les formules mathématiques en MathML, rédigées sous LibreOffice. Le document html en question est disponible à l'adresse suivante :

http://floriansella.free.fr/images/blogs/RegressionCirculaire/RegressionCirculaire.html

Revenons à notre régression circulaire. N'étant pas mathématicien, mon document manque probablement de rigueur. Cependant, il est amplement suffisant pour des gens comme moi, c'est-à-dire des ingénieurs en traitement d'images.

A noter que vous trouverez quelque chose de plus sérieux à l'adresse suivante :

www.scribd.com/doc/14819165/Regressions-coniques-quadriques-circulaire-spherique

Dans ce document, l'étude se sert notamment d'un changement de variables que j'ai supprimé par simplicité.

Voilà un petit exemple pour mettre en application tout ça. Il s'agit de calculer un ensemble de points autour d'un cercle défini par son centre (20.5; -3.7) et son rayon 8. Les coordonnées de chaque point sont bruitées par le code suivant :

double const xnoise = 2.0 * (double)(rand() % 200 - 100) / 100.0;
double const ynoise = 2.0 * (double)(rand() % 200 - 100) / 100.0;

Voici la sortie des résultats :

  ---------------------------------------------------------------------
||Expected circle :
||         x center = 20.5
||         y center = -3.7
||           radius = 8
  ---------------------------------------------------------------------
||x original|  x noised|     noise   ||y original|  y noised|     noise||
  ---------------------------------------------------------------------
||      28.5|     27.32|     -1.18   ||      -3.7|     -4.36|     -0.66||
||   28.3252|   29.0052|      0.68   ||  -2.03675|  -2.03675|         0||
||   27.8084|   29.1884|      1.38   || -0.446197| 0.0338029|      0.48||
||   26.9722|   26.5322|     -0.44   ||   1.00216|   2.16216|      1.16||
||   25.8532|   27.0932|      1.24   ||   2.24503|   1.52503|     -0.72||
||   24.5002|   24.6002|       0.1   ||   3.22808|   4.12808|       0.9||
||   22.9724|   22.5924|     -0.38   ||   3.90836|   2.44836|     -1.46||
||   21.3366|   22.5566|      1.22   ||   4.25614|   4.07614|     -0.18||
||   19.6642|   21.5642|       1.9   ||   4.25622|   5.09622|      0.84||
||   18.0283|   16.5683|     -1.46   ||   3.90859|   2.62859|     -1.28||
||   16.5004|   18.3204|      1.82   ||   3.22845|   1.30845|     -1.92||
||   15.1474|   15.1874|      0.04   ||   2.24552|   3.30552|      1.06||
||   14.0282|   13.8682|     -0.16   ||   1.00276|   2.64276|      1.64||
||   13.1919|   11.6119|     -1.58   ||  -0.44552|  -0.12552|      0.32||
||    12.675|    13.035|      0.36   ||  -2.03603|  -2.13603|      -0.1||
||      12.5|     11.44|     -1.06   ||  -3.69926|  -3.17926|      0.52||
||   12.6747|   14.0947|      1.42   ||  -5.36252|  -4.60252|      0.76||
||   13.1913|   12.5713|     -0.62   ||  -6.95313|  -6.71313|      0.24||
||   14.0273|   13.3673|     -0.66   ||  -8.40156|  -8.42156|     -0.02||
||   15.1463|   13.8463|      -1.3   ||  -9.64453|  -9.76453|     -0.12||
||   16.4991|   16.5591|      0.06   ||  -10.6277|  -12.4077|     -1.78||
||   18.0269|   18.4669|      0.44   ||  -11.3081|  -10.6481|      0.66||
||   19.6627|   19.1227|     -0.54   ||  -11.6561|  -12.3761|     -0.72||
||   21.3351|   22.1551|      0.82   ||  -11.6563|  -11.4363|      0.22||
||    22.971|    22.031|     -0.94   ||  -11.3088|  -11.9488|     -0.64||
||   24.4989|   25.4389|      0.94   ||  -10.6288|  -11.7488|     -1.12||
||   25.8521|   25.0921|     -0.76   ||  -9.64602|  -8.50602|      1.14||
||   26.9714|   25.7114|     -1.26   ||  -8.40336|  -9.22336|     -0.82||
||   27.8078|   28.2678|      0.46   ||  -6.95516|  -6.13516|      0.82||
||   28.3249|   28.9049|      0.58   ||   -5.3647|   -3.8047|      1.56||
  ---------------------------------------------------------------------
||Computed circle :
||         x center = 20.3688
||         y center = -3.76401
||           radius = 8.25942
||             rmse = 1.01834
  ---------------------------------------------------------------------

Le cercle prédit est localisé en (20.36; -3.76) avec un rayon de 8.25 et un rmse (Root Mean Square Error) de 1.01834.

Le code est disponible à l'adresse suivante :

http://floriansella.free.fr/images/blogs/RegressionCirculaire/CircularRegression.zip

Aucune dépendance (boost, opencv, ou je sais pas quoi d'autre ... ) n'est requise juste un "g++ main.cpp matrix.cpp -o test" pour compiler. Vous pouvez utiliser ce code sans restriction aucune. Son utilisation est très simple :

// variables
std::vector<double> xlist; // ... The list of the x coordinates of the points that match a circle.
std::vector<double> ylist; // ... The list of the y coordinates of the points that match a circle.
double xcenter = 0.0; // To store the x coordinate of the center of the circle computed by the function.
double ycenter = 0.0; // To store the y coordinate of the center of the circle computed by the function.
double radius = 0.0; // To store the radius of the circle computed by the function.
double rmse = 0.0; // To store the root mean square error resulted from the regression.
 
try
{
   // regression of the points by a circle
   Matrix::circularRegression(xlist, ylist, xcenter, ycenter, radius, rmse);
}
catch(std::exception const & error)
{
   std::cout << error.what() << std::endl;
}
catch(...)
{
   std::cout << "An undefined error has occured" << std::endl;
}

Un petit merci si ça vous a servi à quelque chose. Sinon pas grave.

A noter que ceci a fait l'objet d'une contribution de ma part sous developpez.com. Je dis pas ça pour me la péter mais le fil des messages associé à cette contribution aura tôt fait d'être plus étoffé que ce pauvre blog vide de visiteurs. La contribution en question est disponible à l'adresse suivante :

http://www.developpez.net/forums/d1262000/autres-langages/algorithmes/contribuez/cpp-regression-circulaire/

PS pour les mathématiciens : pas taper !

(Dernière édition 2012-09-17 à 21:18:51)

 

double const xnoise = 2.0 * (double)(rand() % 200 - 100) / 100.0;

double const ynoise = 2.0 * (double)(rand() % 200 - 100) / 100.0;


PUBLIE PAR FLORIAN A 17:39:05
Enregistré dans la ou les catégories maths 
0 commentaires / Poster un commentaire

VENDREDI 07 MAI 2010
PSHUFB pour une conversion RGB / Niveaux de gris

Récemment, j'ai découvert l'instruction PSHUFB (Packed Shuffle Bytes) qui se trouve dans le 4ième jeu d'instruction SSE qui s'appelle SSSE3 (Supplemental Streaming SIMD Extension 3). On trouve un petit blabla à ce sujet sur Wikipedia.

Cela faisait déjà quelques années, que de temps en temps, j'essayais de me faire une conversion RGB/Niveaux de gris en xmmx sans y arriver efficacement pour cause d'instructions peu intéressantes pour le faire. Il est possible que je m'y prenais mal aussi. L'instruction PSHUFB permet désormais une conversion d'une image RGB en niveaux de gris 5 pixels par 5 pixels soit au moins 5 fois plus vite qu'en C/C++ traditionnel.

Par rapport aux instructions PSHUFW (SSE) et PSHUFLW / PSHUFHW / PSHUFD (SSE2), l'adressage des octets à remanier se fait différemment.

Pour chacune de ces 4 instructions (PSHUFW / PSHUFLW / PSHUFHW / PSHUFD), il y a 3 opérandes SRC, DST et un octet ORDER. L'octet ORDER permet de coder l'adresse des octets à copier de SRC vers DST. Les adresses, dans l'octet ORDER, sont codées par paire de bits. Par exemple, si la valeur de l'octet ORDER est de 141, 10001101 en binaire, alors les opérations suivantes seront effectuéees :


DST[0] = SRC [1] car les bits n°0 (1) et n°1 (0) de ORDER codent 1
DST[1] = SRC [3] car les bits n°2 (1) et n°3 (1)
de ORDER codent 3
DST[2] = SRC [0] car les bits n°4 (0) et n°5 (0)
de ORDER codent 0
DST[3] = SRC [2] car les bits n°6 (0) et n°7 (1)
de ORDER codent 2

A noter également que l'instruction PSHUFW ne fonctionne que sur les registres MMX 64 bits. Les instructions PSHUFLW et PSHUFHW fonctionnent uniquement sur 8 octets, respectivement le QUAD WORD de poids faible et le QUAD WORD de poids fort. Ces 2 instructions là permettent de faire un simili PSHUFW sur des registres XMMX 128 bits. En fait, les 4 instructions PSHUFW (SSE) et PSHUFLW / PSHUFHW / PSHUFD (SSE2) ne permettent de manipuler que 4 données 16 bits (WORD) ou 32 bits (DOUBLE WORD) à la fois et uniquement à l'intérieur d'un 64 bits (QUAD WORD).

Rien qu'en ça l'instruction PSHUFB est une révolution. On peut directement manipuler les 16 octets d'un registre 128 bits et les redisposer n'importe où dans le registre 128 bits de destination.

Différence appréciable également : ce n'est pas un octet ORDER qui définit la manipulation mais directement un registre 128 bits (le registre source). PSHUFB n'a que 2 opérandes. Ce registre source charge 16 bytes signés, un byte négatif permettant d'affecter à l'octet de destination associé un 0. On remarquera que PSHUFB donne la possibilité de faire l'équivalent d'un PSHUFW pour des registres 128 bits.

Un exemple d'utilisation, en syntaxe Intel :


vector addresses(16)
addresses[ 0] = 4;    // code l'opération DST[ 0] = SRC[ 4]
addresses[ 1] = 12;   // code l'opération DST[ 1] = SRC[12]
addresses[ 2] = -1;   // code l'opération DST[ 2] = 0
...
addresses[15] = 15;   // code l'opération DST[15] = SRC[15]

char * ptr = &addresses[0];
__asm  mov      eax, ptr   ;    
__asm  movdqu  xmm0, [eax] ; // chargement des adresses dans xmm0
__asm  pshufb  xmm1, xmm0  ; // manipulation de xmm1 par xmm0

 

Concernant, la conversion RGB en niveaux de gris, le code commenté est disponible ici. Il s'agit d'un code C/C++/XMMX qui prend un buffer RGB 24 bits, démultiplexe les 3 plans R, G et B (démultiplexation en 16 bits par paresse) et calcule en 16 bits le produit R*Cr + G*Cg + B*Cb ensuite ramené à 8 bits. Les constantes Cr, Cg et Cb sont les coefficients de pondération des 3 composantes. Le code, à valeur de test, fonctionne uniquement sur 5 pixels RGB uniquement et n'est pas optimisé. Il s'agit là de tester l'instruction PSHUFB dans le cadre d'une conversion RGB - niveaux de gris.

Le résultat est le suivant ;


  RGB =   20  30  40  50  60  70  80  90 100 110 120 130 140 150 160 170 (données en entrée)
    R =   20   0  50   0  80   0 110   0 140   0 170   0   0   0   0   0
(données en sortie)
    G =   30   0  60   0  90   0 120   0 150   0   0   0   0   0   0   0 (données en sortie)
    B =   40   0  70   0 100   0 130   0 160   0   0   0   0   0   0   0 (données en sortie)
 grey =   29  59  89 119 149  56   0   0   0   0   0   0   0   0   0   0 (données en sortie)
 

Voilà quelques considérations qui sortent un peu de l'ordinaire ...


PUBLIE PAR FLORIAN A 19:01:11
Enregistré dans la ou les catégories xmmx 
0 commentaires / Poster un commentaire

MERCREDI 30 SEPTEMBRE 2009
Google sketchup

Dans le cadre de mon activité professionnelle, pour une présentation à une conférence du Vision Show, j'ai voulu utiliser Google Sketchup pour illustrer les différents systèmes de vision. Je suis parti d'un projet vide. Pour le premier système (triangulation laser), j'y ai passé un peu moins d'une journée : il y a quand même pas mal de choses à assimiler. Pour chacun des autres, ça m'a pris entre 1 heure et 2 heures d'autant que je réutilisais les composants précédemment créés (caméra, champ de vision, objets, etc).

Voilà le résultat :

Google Sketchup : triangulation laser

Dans l'ordre,

  • capteur ponctuel en z
  • stéréovision
  • stéréovision par projection de franges
  • triangulation laser
  • projection de franges
  • temps de vol

Au début, le logiciel semble d'apparence simple (on a l'impression d'avoir peu de fonctionnalités lors du premier lancement avec la configuration par défaut de l'IHM) puis on progresse et on voit l'étendue des possibilités. Je considère ce logiciel pour la 3D comme l'est Paint.NET pour la 2D. En gros, Google Sketchup est à AutoCAD ce que Paint.NET est à Photoshop. Gratuit, fonctionnel et relativement puissant sans le côté professionnel qui peut rendre le logiciel impénétrable pour des non initiés. De plus, une énorme base de modèles 3D est disponible en libre téléchargement (http://sketchup.google.com/3dwarehouse/?hl=fr&ct=lc).

Aller une deuxième image en affichage volumique :

Google Sketchup : triangulation laser

Le fichier source à télécharger ici.


PUBLIE PAR FLORIAN A 20:21:20
Enregistré dans la ou les catégories vision 
0 commentaires / Poster un commentaire



FLORIAN SELLA

Bonjour je suis ingénieur en traitement d'images et j'aime bien dessiner à mes heures perdues. Vous trouverez donc sur ce blog des essais et réflexions sur ces 2 sujets. Bonne visite.


ARCHIVES DU BLOG

maths(1)
xmmx(1)
vision(1)



BLOGS

Blog tout confondu
Blog dessins
Blog traitement d'images
Blog divers


S'ABONNER

Fill RSS