Recalage d’images, PIV et corrélation d’images — Les logiciels

Posté par  . Édité par Davy Defaud, Benoît Sibaud, ZeroHeure et bubar🦥. Modéré par ZeroHeure. Licence CC By‑SA.
37
1
oct.
2017
Science

Le recalage d’images est utilisé dans la communauté de l’analyse d’images médicales depuis très longtemps (Barnea & Silverman, 1972, Ledbetter et al. 1979), on peut même remonter plus loin dans le temps avec les travaux de Sidney Bertram, 1963 en fournissant tous les outils via des descriptions de circuits analogiques (vraiment impressionnant pour l’époque).

Cette dépêche fait suite à Recalage d’images, PIV et corrélation d’images — Les bases théoriques. Les bases théoriques ayant été présentées dans la dépêche précédente, le lecteur souhaitant comprendre plus en détail le fonctionnement de ces logiciels pourra s’y référer.

Sommaire

Logiciels libres de DIC, PIV et recalage

Nous ne présenterons ici que des logiciels libres, il existe beaucoup de solutions propriétaires permettant de faire de la corrélation d’images numériques (DIC), de la Particle image velocimetry (PIV) et du recalage, mais elles ne seront pas abordées.

Logiciels de DIC

Yet Another Digital Images Correlation (YaDICs)

Le laboratoire de Mécanique de Lille, comme beaucoup d’autres laboratoires de mécanique, a développé ses propres outils permettant aux chercheurs en mécanique expérimentale d’avoir accès à la transformation que l’on peut avoir lorsque l’on fait un essai mécanique.

YaDICs Figure 1 : Schéma de principe de la plate‐forme YaDICs. En rouge, l’ajout de la méthode de régularisation par filtre médian

Pour résumer notre plate‐forme :

  • métrique des moindres carrés (SSD) ;
  • interpolateurs : bilinéaires et bicubiques ;
  • transformations :
    • globales (mouvement des corps rigides, déformations homogènes, essais brésiliens),
    • locales (possibilité d’utiliser toutes les cinématiques globales localement),
    • élastiques avec des cinématiques éléments finis en 2D et 3D ;
  • optimiseurs :
    • Gauß‐Newton pour le global et le local (très peu de variables par éléments ⇒ Hessienne pas chères),
    • simple descente au gradient pour les éléments finis : \mu_{k+1}=\mu_k-\alpha_k \gamma_k {\partial \mathcal{C}(\mu,\mathcal{I_F},\mathcal{I_M})}/{\partial \mu} avec  \gamma_k=1\!\!1, \alpha_k = 1.2*\alpha_{k-1} si l’erreur diminue et \alpha_k = -0.5*\alpha_{k-1} si l’erreur augmente ;
  • régularisations :
    • ∅ de régularisation par pénalisation de Tikhonov,
    • régularisation par filtrage médian ou gaussien (le gaussien n’étant quasiment jamais utilisé).

Voici un exemple de fichier de paramètres documenté permettant de mettre en place n’importe quelle séquence de calcul à n’importe quelle échelle :

netcdf DIC_parameter_file {
dimensions:
  dim_sample = UNLIMITED;// Allow adding automatically, correlation field along a "time" or "sample" dimension, function of image_list size. This option doesn't work thus dim_tmp (temporarly) has been added (by default Unlimited - don't change)
  dim_tmp = 4;// Number of image couple to correlate. If image_list is larger, only dim_tmp correlation will be done. If image_list is smaller a crash will occur after the last correlation.
  dim_system = 1;// number of image composing the observed field (by default 1 - don't change)
  dim_camera = 1;// 1 for correlation or more for stereocorrelation (by default 1 - don't change)
  dim_image = 2;// couple of image to correlate (by default 2 - don't change)
  dim_string = 180;// image name (by default 180 - don't change except for larger image name)
//  change only if mod_image
    dim_mode = 1;//number of image modes plus 1 (1=default)
    dim_space = 2;//associated dimensions

variables:
//input_outpout format definition. Different format exist: double_double, float_float, char_float
  byte format;
       format:type = "float_float";
//---------------------------Mesh------------------------------
// Only used Meshes (in correlation method) are taken into account
  // mesh template for local transformation (Block Matching)
  byte mesh_IC;// give the name you want
       mesh_IC:path = "none";// for external mesh (not activated yet)
       mesh_IC:type = "generate_centroid";// mesh are automatically generated on images (at least one node at each image corner) and displacement is calculated at element centroid. (by default "generate_centroid" - don't change)
       mesh_IC:pitch = 32, 32;// number of pixels between nodes or centroids. (Warning: pitch can be automatically reduced or increased at boundaries depending on image dimensions)
       mesh_IC:window = 64, 64;// dimension of Elements. Each one is centered around its centroid (warning: according to pitch and window dimensions an overlap can exist: the element overlap is here of 50%)
       mesh_IC:component = "x,y";// spatial dimensions. (Warning: if 3D, "z" is required as well as a 3rd component on pitch and window)
       mesh_IC:store_peak = 1;// number of intercorrelation peak to store (usually 3 in PIV)
  // mesh template for Elastic transformation (local with unity partition - Finite Element Method)
  byte mesh_OFFEM;// give the name you want
       mesh_OFFEM:path = "none";// for external mesh (not activated yet)
       mesh_OFFEM:type = "generate_regularFem";// regular mesh is automatically generated on image (at least one node at each image corner)
       mesh_OFFEM:pitch = 64, 64;// number of pixels between nodes or element size. In the case of finite element, due to unity partition, pitch = window, contrary to mesh_IC. 
       mesh_OFFEM:component = "x,y";// spatial dimensions. (Warning: if 3D, "z" is required as well as a 3rd component on pitch and window)
  // mesh template for global transformation (rigid, affine...)
  byte mesh_OFI;
       mesh_OFI:path = "none";// for external mesh (not activated yet)
       mesh_OFI:type = "global";// there is no mesh for global transformation (by default "global" - don't change)
       mesh_OFI:component = "x,y";// spatial dimensions. (Warning: if 3D, "z" is required)
// mesh template for Local transformation (rigid, affine...)
 byte mesh_OFIB;
    mesh_OFIB:path = "none";// for external mesh (not activated yet)
    mesh_OFIB:type = "generate_centroid_integrated";// default for OFIB (don't change)
    mesh_OFIB:pitch = 16,16;// pitch between two windows
    mesh_OFIB:window = 16,16;// windows size (in this example overlap = 0%)
    mesh_OFIB:component = "x,y";// space component (x,y and z-if 3D)
//---------------------------Mesh------------------------------ 
//---------------------------Correlation Method----------------
 //example of "Inter-Correlation" method definition
 byte IC;// Give the name you want
  IC:type = "IC_fftw";// define the correlation method: "IC_direct", "IC_fftw", "IC_fftw_phase", "IC_fft", "IC_fft_phase"
  IC:mesh = "mesh_IC";// associate the correlation method and the mesh
  IC:peak = "interpP";// define the sub-pixel precision strategy:  "max", "interpP", "interpG", "interpBiLin", "barycenter"
  IC:median_filter = 3;// box for median filter [+-(mf-1)/2] in elements
  IC:store_correlogram = 0;// allow storing the correlogram : 0 = no storage, otherwise correlogram are stored in netcdf format

//example of "Optical Flow Integrated" method definition for "Homogenous strain"
  byte OFI_H; // Give the name you want
    OFI_H:type = "OFI";// Optical-Flow Integrated
    OFI_H:mesh = "mesh_OFI";// associated mesh
    OFI_H:shape = "homogeneous";// assumed kinematics [rigid_body, homogeneous, blade,brasilian,image]
    OFI_H:modes = "Tx,Ty,Rz,Exx,Eyy,Exy";// desired modes
    OFI_H:refX = -1;//shape function X reference, -1 = automatic, i.e. system reference at the center of image, otherwise provide location in pix
    OFI_H:refY = -1;//shape function Y reference, -1 = automatic, i.e. system reference a the center of image, otherwise provide location in pix
    OFI_H:refZ = -1;//shape function Z reference, -1 = automatic, i.e. system reference a the center of image, otherwise provide location in pix
    OFI_H:convergence_speed = 1e-3;//convergence criterion (Epsilon_(i+1) - Epsilon(i)<convergence_speed)
    OFI_H:interpolation_loop = "cubic";//Subpixel interpolation : cubic or linear by default - during iterations
    OFI_H:interpolation_final = "cubic";//Subpixel interpolation : cubic or linear by default - at the final iteration

//example of "Optical Flow Integrated" method definition for "rigid body motion"
  byte OFI_RB;// Give the name you want
    OFI_RB:type = "OFI";// Optical-Flow Integrated
    OFI_RB:mesh = "mesh_OFI";// associated mesh
    OFI_RB:shape = "rigid_body";// assumed kinematics [rigid_body, homogeneous, blade,brasilian,image]
    OFI_RB:modes = "Tx,Ty,Rz";// desired modes
    OFI_RB:refX = -1;//shape function X reference    OFI_RB:refY = -1;//shape function Y reference
    OFI_RB:refZ = -1;//shape function Z reference
    OFI_RB:convergence_speed = 1e-3;//convergence criterion (Epsilon_(i+1) - Epsilon(i)<convergence_speed)
    OFI_RB:interpolation_loop = "cubic";
    OFI_RB:interpolation_final = "cubic";

//example of "Optical Flow Integrated per block" method definition
  byte OFIB;// Give the name you want
    OFIB:type = "OFIB";// Optical-Flow Integrated per Block
    OFIB:mesh = "mesh_OFIB";// associated mesh
    OFIB:shape = "rigid_body";// assumed kinematics [rigid_body, homogeneous, blade,brasilian,image]
    OFIB:modes = "Tx,Ty,Rz";// desired modes
    OFIB:refX = -1;//shape function X reference
    OFIB:refY = -1;//shape function Y reference
    OFIB:refZ = -1;//shape function Z reference
    OFIB:convergence_speed = 1e-3;//convergence criterion (Epsilon_(i+1) - Epsilon(i)<convergence_speed)
    OFIB:median_filter = 3;// box for median filter [+-(mf-1)/2] in elements
    OFIB:interpolation_loop = "cubic";//cubic or linear by default - during iterations
    OFIB:interpolation_final = "cubic";//cubic or linear by default - at the final iteration

//example of "Optical Flow Finite Element Method" method definition
  byte OFFEM;// Give the name you want
    OFFEM:type = "OFFEM";// Optical-Flow Finite Element Method
    OFFEM:mesh = "mesh_OFFEM";// associated mesh
    OFFEM:method = "gradient";// optimization method - newton (warning : not suitable for large 2D images and 3D images) or gradient
    OFFEM:alpha0 = 5e-1;// init coef for gradient method
    OFFEM:shape = "Q4";// assumed kinematics : C8 for 3D
    OFFEM:convergence_speed = 1e-3;//convergence criterion (Epsilon_(i+1) - Epsilon(i)<convergence_speed)
    OFFEM:median_filter = 3;// box for median filter [+-(mf-1)/2] in elements
    OFFEM:iter_max = 100;//max number of authorized iterations : thus calcul stop either if iter_max is reached, if convergence_speed is reached or if the result diverges
    OFFEM:interpolation_loop = "cubic";//cubic or linear by default - during iterations
    OFFEM:interpolation_final = "cubic";//cubic or linear by default - at the final iteration
//---------------------------Correlation Method----------------

//------------Correlation sequence per image couple------------
  byte sequence; 
    sequence:analysis = "OFI_H,OFI_H,OFFEM,OFFEM,OFFEM,OFFEM,OFFEM";//correlation sequence list
    sequence:filter = 6,5,4,3,2,1,0;//associated pyramidal filter: super-pixel = 2^i x 2^i pix, e.g. i=0 -> 1x1, i=6 -> 64x64 
//------------Correlation sequence per image couple------------


//------------------------image list---------------------------
  char image_list(dim_sample, dim_tmp, dim_image, dim_string);

      image_list:type = "File";//image capture from File or on-fly Grab(not activated yet)
      image_list:mask = "mask.png";// mask name (must be in "path" directory)
      image_list:mode = 0;// 0:dependent mode (Image sequence) - 1:separate mode (PIV)
      image_list:path = "./img/";// image location


  //declare desired outputs and physical pixel size
  float output;
      output:display_mode = 0;//Allow a real time field display (not activated yet)
      output:store_img = 0;//store input data in netcdf format
      output:store_def = 1;//store the deformed image after deformation
      output:store_field = 1;//store converged displacement fields on image resolution (Warining: avoid in 3D -> in every cases displacement fields are stored using mesh resolution)
      output:store_res = 1;//store converged residue fields
      output:store_strain = 0;//store converged strain field (not activated yet)
      output:ratio_pxmm = 65;//Allow to output fields in physical units (Warning bug : only put integers, i.e. not 6.5 um but 650 nm for example)
      output:format = "";//tiff or netcdf default
//------------------------image list---------------------------

  char mode_img(dim_sample, dim_mode, dim_space, dim_string);  

//------------------------image list--------------------------- 
// allocate variables
data:

  format = 1;

  mesh_IC = 1;
  mesh_OFIB = 1;
  mesh_OFI = 1;
  mesh_OFFEM = 1;

  IC = 1;
  OFI_H = 1;
  OFI_RB = 1;
  OFIB = 1;
  OFFEM = 1;

  sequence =1;

  image_list =
"00000001.jpg", "00000100.jpg",
"00000001.jpg", "00000150.jpg",
"00000001.jpg", "00000200.jpg",
"00000001.jpg", "00000250.jpg";

//only for mod_image: i.e. you have created your own image modes as follows, for example:
  mode_img =  
"T.png",    "zero.png",
"zero.png", "T.png",
"Exx.png",  "zero.png",
"zero.png", "Eyy.png",
"Eyy.png",  "Exx.png";

}

L’idée est en fait assez simple :

  • on commence par définir un maillage en fonction d’un type associé à une méthode ;
  • on associe ce maillage à une méthode de corrélation ;
  • on définit quelle méthode on utilise à quelle échelle du processus pyramidale ;
  • on définit les données que l’on veut sortir (déplacements aux pixels, résidus…) ;
  • on définit la séquence d’images.

On peut ensuite compiler le programme avec ncgen3, puis lancer l’exécutable. Les résultats sont visualisables avec le logiciel ncview car tout est fait en netcdf, il est donc simple de charger les données en Matlab et Python. On pourra trouver ce fichier très verbeux, mais c’est le prix à payer lorsque l’on veut pouvoir fabriquer comme on le souhaite ses séquences. En général, on part d’un fichier déjà fait et on le modifie. Il n’y a donc pas forcément besoin de comprendre tout le fichier.

Pydic

Pydic est une interface python très simplifiée à la bibliothèque OpenCV pour des applications plutôt dédiées à des problèmes de mécanique. Pydic permet, de façon très simple, de visualiser les résultats (champ de déformation) grâce à diverses animations graphiques, puis d’exploiter ces résultats grâce au langage Python. La première vocation de pydic est avant tout d’être très simple et pédagogique.

Figure pydic Figure 2 : Schéma de principe de la corrélation d’images numériques dans pydic : prise d’images d’un essai (ici, un essai de flexion), traitement de la séquence d’images pour avoir la transformation (champs de déplacements), calcul et affichage des déformations

Par défaut, pydic utilise :

  • le calcul du flot optique par la méthode de Lucas‐Kanade (méthode locale utilisant une métrique de type SSD et une descente de Gauß‐Newton) avec des imagettes de référence alignées sur une grille ;
  • un lissage des déplacements obtenus par la méthode des splines (hypothèse cinématique de continuité).

Il est également possible de changer ce comportement en utilisant des imagettes de référence non alignées sur une grille. Dans ce cas, les positions des imagettes sont choisies grâce à l’algorithme de Shi‐Tomasi, qui permet d’optimiser ce choix. Il est alors nécessaire de recalculer ces déplacements sur une grille. Pour cela, les différentes méthodes d’interpolation proposées sont celles disponibles par la méthode griddata (linéaire, plus proche voisins ou cubique) de la bibliothèque scipy.

Une limitation de cette approche est que l’hypothèse cinématique utilisée par pydic (lissage par spline) ne permet pas de prendre en compte les discontinuités fortes pouvant apparaître lors de fissurations. Aussi, des travaux en cours consistent à intégrer l’algorithme du flot optique dense qui devrait permettre de s’affranchir de cette hypothèse.

Digital Corelation Engine (Dice)

Les logiciels de DIC local et global ont comme originalité d’avoir mis en place un algorithme du simplexe pour la minimisation de la fonctionnelle qui évite d’avoir à calculer le gradient de l’image. Le logiciel est assez récent, je n’ai pas eu l’occasion de le tester (mes tentatives de compilation n’ont pas réussi…).

dolphin_dic

C’est certainement d’un point de vue du mécanicien des solides, un peu versé dans le numérique, la solution libre la plus prometteuse. Elle est basée sur le logiciel FEniCS pour représenter les éléments finis (fonction de paramétrisation de la transformation) et déterminer les opérateurs de régularisation. FEniCS est certainement la bibliothèque éléments finis la plus ambitieuse du moment, Martin Genet a fait un travail très intéressant en le liant avec VTK pour une interpolation des images rapides et en faisant en sorte que l’on puisse facilement intervenir dans le code via du Python.

Logiciels de recalage d’images

Elastix et SimpleElastix

Elastix est probablement le logiciel de recalage d’images le plus polyvalent, il est développé par deux équipes néerlandaises (Erasmus MC, Leiden University Medical Center), elles ont en permanence des doctorants et ingénieurs qui travaillent sur cette plate‐forme. C’est donc une plate‐forme pérenne et de par leur culture en informatique, elle dispose des dernières avancées comme l’utilisation d’accélération des calculs par OpenCL.
C’est vraiment un outil intéressant, car ils ont implémenté beaucoup de métriques telle que la SSD, la NCC, plusieurs informations mutuelles. Il est possible de pénaliser la métrique par différentes énergies « mécaniques » et il y a un grand nombre d’optimiseurs. C’est donc une boîte à outils remarquable qui permet bien souvent d’arriver à des résultats très satisfaisants. Ils n’ont en revanche pas du tout orienté leur cadriciel vers le traitement temps réel ou la métrologie, le mécanicien devra donc bien analyser et comprendre ces outils avant de parvenir à les utiliser sereinement. Le LML a utilisé Elastix jusqu’à ce qu’il dispose de ses propres outils, il nous arrive encore d’utiliser Elastix lorsqu’il arrive que cela ne fonctionne pas avec nos outils. Il n’est pas possible de mettre en place la régularisation par filtrage telle que celle présentée dans les méthodes de type démon ou dans YaDICs.

Voici un exemple de code Python permettant de faire du recalage d’image non rigide avec SimpleElastix :

import SimpleITK as sitk

elastixImageFilter = sitk.ElastixImageFilter()
elastixImageFilter.SetFixedImage(sitk.ReadImage("fixedImage.png"))
elastixImageFilter.SetMovingImage(sitk.ReadImage("movingImage.png"))

parameterMapVector = sitk.VectorOfParameterMap()
parameterMapVector.append(sitk.GetDefaultParameterMap("affine"))
parameterMapVector.append(sitk.GetDefaultParameterMap("bspline"))
elastixImageFilter.SetParameterMap(parameterMapVector)

elastixImageFilter.Execute()
sitk.WriteImage(elastixImageFilter.GetResultImage())

On constate tout de suite que ce code est extrêmement compact, il est ainsi extrêmement simple de mettre en place un recalage d’image via cette interface développée par Kasper Marstal, un doctorant de Stefan Klein.

Il est aussi très simple de mettre en place ces algorithmes sur une séquence de par l’aspect orienté objet de leur code, comme montré dans l’exemple ci‐dessous :

transformixImageFilter = sitk.TransformixImageFilter()
transformixImageFilter.SetTransformParameterMap(transformParameterMap)

population = ['image1.hdr', 'image2.hdr', ... , 'imageN.hdr']

for filename in population:
    transformixImageFilter.SetMovingImage(sitk.ReadImage(filename))
    transformixImageFilter.Execute()
    sitk.WriteImage(transformixImageFilter.GetResultImage(), "result_"+filename)

OpenCV

Bien que cette bibliothèque ne soit pas dédiée au recalage d’images, elle implémente de nombreux algorithmes, qui sont au niveau de l’état de l’art de ce qu’il se fait. La plupart sont accessibles via Python, si ce n’est le calcul via le processeur graphique. On pourra par exemple noter le développement de l’algorithme DeepFlow qui se sert de l’estimation par DeepMatching du champ de déplacement dense pour pénaliser un algorithme de flot optique régularisé par Tikhonov. cela permet de rajouter une notion de distance qui place cette méthode dans les approches hybrides iconiques et intensité. Bien que plus lent que les algorithmes habituels, il a une plus grande robustesse que les algorithmes de flots optiques traditionnels.
OpenCV dispose également de méthodes permettant de trouver des champs épars (telles que Lucas‐Kanade) avec les méthodes éprouvées, puis de projeter sur un champ dense.
OpenCV se distingue aussi par la mise à disposition d’outils permettant l’identification des paramètres intrinsèques d’un système de stéréo‐vision, ce qui permet de mettre en place simplement une chaîne de stéréo‐corrélation.
De plus, OpenCV permet de mettre en place des méthodes de type scale‐invariant feature transform (disponible dans contrib, de par les brevets associés à cette méthode), qui permettent de trouver des correspondances pour des changements d’échelles bien plus grands que ce que l’on peut faire avec le flot optique où même le recalage d’images. Cela peut faire également une bonne pénalisation iconique d’une fonctionnelle de flot optique, comme cela a été fait dans SIFTFlow (non implémenté dans OpenCV à l’heure actuelle).
J’ai testé cette bibliothèque pour des opérations de base, comme l’interpolation d’images par une transformation non rigide, elle est extrêmement rapide. En Python, elle va généralement \approx 6 \times plus vite que Scikit Image, sur des opérations profitant des instructions de vectorisations et accélérations diverses du processeur.

ImageJ

ImageJ est une solution d’analyse d’images assez généraliste, c’est un peu le couteau suisse de l’analyse d’images multi‐dimensionnelle et multi‐composante. Il existe une multitude de greffons, des gens ont développé des outils pour le recalage d’images, il en existe bien d’autres (liste des algorithmes) ! Il y a aussi des outils plus spécifiques tels que le greffon PIV.

CImg

Je pense qu’il est normal de citer David Tschumperlé qui a souvent contribué à ce site par des dépêches sur sa bibliothèque Cimg. C’est d’ailleurs cette bibliothèque que nous avons utilisée pour faire YaDICs ! Un grand merci à lui pour cette bibliothèque.

Un exemple assez simple de flot optique multi‐échelle peut être trouvé dans son dépôt de code.

Logiciels de PIV

UVMat, Mpiv, PIVmat, OpenPIV

Il s’agit de logiciels de PIV développés en MATLab, permettant d’obtenir des champs de vitesses par PIV en blockmatching par intercorrélation. Elles proposent à peu près toutes la même approche.

Si l’on regarde la description des utilitaires disponibles dans OpenPIV, on constate que ce sont en fait plusieurs projets qui réimplémentent l’état de l’art des algorithmes de PIV. Il y a aussi des outils de post‐traitement permettant d’avoir un ordre d’idée des pressions et des outils de décomposition de type POD facilitant l’analyse énergétique des écoulements. Ce sont des outils métier permettant de faire l’analyse de l’image jusqu’à la compréhension des phénomènes physiques. Les algorithmes disponibles dans OpenPIV sont moins sophistiqués que ceux du recalage d’images. Dans cette suite logicielle, ils ne résolvent pas le flot optique mais calcul l’inter‐corrélation à l’aide d’une transformation de Fourier rapide, car le passage dans Fourier permet de transformer la convolution en multiplication, c’est ainsi une simple recherche de maximum, les coordonnées du maximum donnant alors la valeur du champ de déplacement sur l’imagette (ZOI). Cette méthode est assez robuste au bruit mais moins précise que les méthodes de flot optique, c’est pour cela qu’ils règlent le temps entre deux images de manière à avoir une amplitude de 10 pixels de déplacement, ils contrebalancent ainsi la moins bonne sensibilité de la méthode par sa plus grande robustesse. Cette approche est géométrique, elle cherche une distance dans un espace particulier, elle ne fait donc pas partie des approches basées sur l’intensité du niveau de gris. Il existe plein de méthodes qui arrive à augmenter sa sensibilité, mais la quantification de l’image en imagette limite beaucoup les possibilités d’amélioration.
D’autres chercheurs tels que Thomas Coppertti pensent qu’il est tout à fait possible d’utiliser les flots optiques dans ces applications, ce qu’il a fait avec succès dans différents articles. Il propose des approches classiques des régularisations par un opérateur \Gamma adapté à la mécanique des fluides (divergent ou laplacien).
Il n’existe personne à ma connaissance, utilisant les approches de type DeepFlow, couplant flot optique régularisé et approche géométrique, alors que celles‐ci semblent particulièrement adaptées à leur problème.

JPIV

JPIV est un logiciel de PIV. Le programme est libre, indépendant de la plate‐forme et disponible sous les termes de la licence GNU General Public License.

JPIV est utile pour :

  • évaluation PIV multi‐passe, multi-grille, basée sur la transformation de Fourier rapide ;
  • évaluation d’ensemble d’un pixel ;
  • filtrage des champs vectoriels, détection des valeurs aberrantes ;
  • calcul de vorticité et autres dérivés ;
  • reconstruction du champ de vitesse 3D ;
  • affichage de l’image (PNG, TIFF, PGM, IMX, IM7) ;
  • affichage du champ vectoriel ;
  • exportation de graphiques de champs vectoriels (EPS, PDF, EMF, SVG, PNG, JPEG…) ;
  • extraction de profil de vitesse, estimation de débit ;
  • statistiques sur les champs vectoriels ;
  • contrôle de ligne de commande et scripts via BeanShell ;
  • traitement par lots d’une liste de fichiers ainsi que de répertoires.

Autres codes

Il est aujourd’hui bien plus simple de voir ce que les gens développent avec des sites comme GitHub.

Une recherche de Digital Image Correlation donne une vingtaine de résultats, c’est à peu près la même chose pour une recherche sur Particule Image velocimetry. En revanche, lorsque l’on cherche Image Registration, on trouve alors plus de 300 projets, ce qui montre bien qu’il y a bien plus de développement libre en informatique dans la section analyse d’images (en tout cas sur GitHub) que dans le domaine de la mécanique.

On constate donc qu’il y a un grand nombre de solutions existantes pour résoudre des problèmes très proches ! Et, je ne parle pas ici de tous les gens que je connais qui développent dans leur coin un petit bout de code pour répondre à leur besoin !

Autres domaines

Outre le médical et la mécanique, il existe d’autres domaines ayant recours à du recalage d’images (ou de vidéos).

Citons par exemple l’exploration spatiale et notamment les images en provenance de sondes, satellites et rovers martiens, qui sont utilisées pour recréer des paysages ou des panoramas. Damia Bouic, qui traite les images brutes des rovers martiens et fournit des tutoriels pour GIMP et Hugin, a été récemment interviewée par Florence Porcel, vidéaste de la chaîne La folle histoire de l’univers ou encore le morphing (transition entre deux images) ou la souris optique (mouvement estimé par le recalage des images acquises).

Il est toutefois important de noter que toutes les méthodes de recalage utilisées pour la création de panoramas ne sont pas basées sur la méthode présentée ici, par exemple Hugin utilise la technologie brevetée SIFT pour trouver les correspondances entre les différentes images, c’est une approche géométrique pure et non basée sur les mesures de similarité ou sur l’approche iconique. Cette approche est généralement plus efficace lors de grandes déformations mais ne permet pas d’avoir des champs denses, il faut donc reconstruire par interpolation de données non structurées une transformation, les deux principales méthodes sont la triangulation ou l’utilisation de splines à base radiale, comme dans l’article sur l’interpolation multivariée.

Aller plus loin

  • # ZNCC

    Posté par  . Évalué à 1.

    Merci pour cette dépêche très intéressante.

    Au boulot je me sers pas mal de ZNCC pour la corrélation dense d'images satellites. Ce n'est pas du tout le même besoin que ce que tu décris, puisque dans notre cas il y a très peu de bruit et les images sont très texturées.

    Au passage je cherche une méthode de dérivation de ZNCC, pour optimiser la convergence vers un minimum local, si ça te dit quelque chose :)

    • [^] # Re: ZNCC

      Posté par  . Évalué à 2.

      Salut,

      Si tu as des images très définies et texturé, tu peux essayer de jouer avec Elastix. Il n'y a pas la ZNCC mais la NCC donne souvent de bons résultats. Si tu as des exemples où la NCC serait en défaut par rapport à la ZNCC, cela m'intéresse.

      En ce qui concerne la méthode de dérivation, je ne comprends pas trop ta question. La dérivation d'une métrique est unique, il n'y a pas plusieurs façon de le faire. Ce qui change c'est les techniques de descente au gradient que tu vas faire à partir de ta dérivée. Je te renvois vers la précédente dépêche ou je détail tout cela.

      En ce qui concerne la l’optimisation de la convergence, on cherche en général à converger vers le minimum global, si tu converges vers un minimum local différent du minimum global ta solution est fausse. Ce qu'il est systématiquement fait pour converger vers un minimum global est de faire un processus pyramidal car la métrique est moins chahutée lorsqu'on l'applique à une image moyennée (ou lissée). Si ce qui t'intéresse est de converger rapidement vers la solution alors il n'y a pas vraiment de réponse universelle. Tu peux choisir de n'évaluer qu'un nombre de points tirés aléatoirement ce qui réduit la taille de ton problème mais augmente le nombre d'itérations ou choisir de prendre la descentes au gradient intelligentes de type quasi-Newton ou tu vas chercher la meilleure estimation de la hessienne afin de minimiser ton nombre d'itérations. En mécanique, en général on va plutôt évaluer la métrique sur tous les échantillons et utiliser des méthodes au gradient de type Gauss-Newton. En analyse d'image médicales l'équipe d'Elastix utilise plutôt le gradient stochastique avec une simple descente au gradient (steppest). Si tu décortiques la précedente dépêche tu trouvera en lien tous les articles, ce qui te permettra de comprendre ce dont je te parle.

      Je pense que pour ton cas une descente au gradient avec un Gauss-Newton devrait bien faire l'affaire.

      • [^] # Re: ZNCC

        Posté par  . Évalué à 1.

        En théorie ZNCC est plus robuste aux variations d'intensité, mais je t'avoue que je n'ai pas comparé.

        La dérivation en x et y de ZNCC n'est pas évidente il me semble, si tu as la formule sous la main je veux bien voir. Et encore plus si tu as la hessienne :)

        Non je parle bien de minimum local, dans mon cas d'usage je peux considérer l'estimation initiale suffisamment bonne pour être au voisinage du minimum global. Mais c'est vraiment très spécifique !

        Merci pour toutes tes remarques en tout cas.

        • [^] # Re: ZNCC

          Posté par  . Évalué à 2.

          En théorie ZNCC est plus robuste aux variations d'intensité, mais je t'avoue que je n'ai pas comparé.

          Du coup, ça se regarde ! En général complexifier la métrique si ce n'est pas nécessaire n'est pas une bonne idée …

          La dérivation en x et y de ZNCC n'est pas évidente il me semble, si tu as la formule sous la main je veux bien voir. Et encore plus si tu as la hessienne :)

          Marius Staring explique comment faire la dérivation de la ZNCC dans cet article. Ils appellent NCC la ZNCC. Tu peux trouver l'implémentation de la ZNCC dans Elastix, regardes page 6 du manuel. Au passage, Gauss-Newton ne nécessite pas le calcul de la hessienne mais l'estime à partir du jacobien. La page wikipedia de Gauss-Newton l'explique très bien à la section "Dérivation à partir de la méthode de Newton" en prenant l'exemple d'une forme quadratique.

          Non je parle bien de minimum local, dans mon cas d'usage je peux considérer l'estimation initiale suffisamment bonne pour être au voisinage du minimum global. Mais c'est vraiment très spécifique !

          Non, tu parles bien du minimum global ! Ce que tu dis, c'est justement qu'il n'y a plus de minimum local autour de la solution à laquelle tu te trouves et que tu veux trouver une solution qui s'approche encore plus du minimum global.

  • # Brevet pas valide en Europe

    Posté par  . Évalué à 3.

    Merci pour la contribution.

    J'ai mois même bossé sur une implémentation open-source de SIFT et j'ai pris contact avec les ayant droits: pas de problème à ce que le code soit en licence MIT, les contraintes ne sont que pour les utilisateurs, pas pour le développeur. C'est a dire qu'il suffit pour les utilisateurs de ne pas être sur le sol américain. Et même dans ce cas là, il leur suffit de ne pas "faire d'argent avec" pour ne pas avoir besoin de licence.

    Si cela intéresse du monde, le code est ici:
    https://github.com/silx-kit/silx/blob/master/doc/source/Tutorials/Sift/sift.ipynb
    et je l'utilise pour recaler des time-lap (entre autre).

    • [^] # Re: Brevet pas valide en Europe

      Posté par  . Évalué à 2.

      Ta contribution bien sympa !

      Je dois t'avouer qu'il y a encore des choses qui m'échappent dans les SIFT, si j'ai bien compris tu crées des descripteurs de dimension 128 pilotés la Différence de Gaussiennes (DoG) et des invariants du tenseur de structure de l'image aux différetnes échelle des DoG. Cette partie là ne me pose pas trop de problème. Ce que j'ai plus de mal à comprendre (et que j'aimerais bien mieux comprendre), c'est les techniques permettant d'appairer les descripteurs entre les deux images. Visiblement c'est loin d'être simple. Si tu as bien compris le truc, j’apprécierais vraiment que tu l'expliques, si tu arrive à en trouver le temps.

      • [^] # Re: Brevet pas valide en Europe

        Posté par  . Évalué à 1.

        Salut,

        Concernant la création des descripteurs, tu as l'air au point. Ce dont je me souviens (c'était il y a quelques années) c'est qu'après avoir trouvé l'orientation principale du gradient, on fait la même chose sur un voisinage (4x4 si je me souviens bien) avec des histogrammes sur les (8) orientations dans ces voisinages. C'est la dessus que se trouve le brevet je crois (en plus de la DoG).

        Le matching est assez simple en fait: il est basé sur norme L1 entre deux descripteurs, cette distance doit être la plus petite possible. L'idée un peu astucieuse est de dire qu'un descripteur A marche avec B d'un l'ensemble de keypoints si C, le deuxième meilleur de cet ensemble match beaucoup moins bien.
        soit: dAB < k.dAC avec k=0.53 (0.73*0.73, je sais pas d'où cela viennent ces chiffres).

        Je ne penses pas qu'il y ait de brevet non plus sur le matching. Il faut souvent valider le matching par un RANSAC.

        Voici une implémentation python (pour référence) du matching.
        https://github.com/silx-kit/silx/blob/master/silx/opencl/sift/match.py#L337
        Le reste du fichier c'est l'implémentation GPU.

Suivre le flux des commentaires

Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.