Possible nouveau module déconvolution

[aurelienpierre]

De bonnes nouvelles pour vous tenir en haleine :



L'image centrale a été défloutée en 10 min à partir de celle de gauche, elle-même générée par un flou gaussien de 5 px à partir de celle de droite. La zore en rouge est le masque à l'intérieur duquel est calculé le flou (pour ne pas le calculer sur toute l'image, ce qui prendrait trop de temps).

Notez en particulier les fines rayures du col de la chemise, qui sont récupérées alors qu'elles ont presque entièrement disparu de l'image floue, et la relative absence d'artefacts (anneaux, double bordures et bruit)

10 minutes c'est ta version Python ou C? Si c'est en C c'est pas très encourageant non? Ou alors tu as plein d'idée pour améliorer les performances.... ?

C'est du Python, avec 3 threads (1/canal RGB), et 450 itérations × 3 canaux.

J'ai modifié l'algorithme source des chercheurs en utilisant 3 kernels de flou (1/canal) au lieu d'un seul, ce qui permet de tenir compte de la variance de la PSF en fonction de la longueur d'onde (et possiblement régler les aberrations chromatiques), mais aussi de parallèliser complètement le traitement des 3 canaux. Le problème est qu'utiliser 3 threads divise le temps d'exécution seulement par 2 en comparaison avec 1 thread, ce qui est probablement dû au Global Interpreter Lock (GIL) de Python, qui ne permet de traiter qu'une opération à la fois (donc c'est du multitasking, pas du vrai multithreading).

J'ai perdu un temps infini à essayer d'optimiser mon code Python, pour un gain infime, ce qui laisse apparaître que Python est une catastrophe pour le calcul à haute performances. Par exemple, j'utilise une librairie mathématique (Numpy) qui se base sur LAPACK pour les calculs matriciels (en Fortran, je pense). Mais certaines fonctions Numpy ressortent des types de données figés peu importe le type de l'entrée (la FFT ressort obligatoirement un type `complex` 64 ou 128 bits, d'autres des flottants 64 bits alors que toutes mes inputs sont 32 bits, ce qui est déjà 2 × trop). Du coup, cette espèce de non-cohérence du typage créée des copies de copies de variables quand tu pourrais juste les remplacer inplace et l'optimisation des accès mémoire est quasi-impossible.

Tout ce que j'ai réussi à faire, c'est utiliser un compilateur JIT qui permet à peu de frais d'assurer un semblant de cohérence des types et de prémâcher les fonctions matricielles de base (avec un support hypothétique des boucles parallélisée mais dans des conditions que je ne remplis pas). Ça permet de gagner genre 3 s/20 itérations. Mais du coup, j'ai abandonné l'idée d'être plus performant. En C, on pourrait gérer des vrais threads avec openmp, et en plus tuiler la FFT (on en calcule 4/itération/canal) bas niveau avec FFTW.

Évidemment, le nombre d'itérations requis dépend de la taille du flou et de la précision attendue. Pour compenser un flou naturel d'objectif, avec une PSF 3×3 ou 5×5 px, on peut s'en sortir avec 50-150 itérations, soit de l'ordre d'une centaine de secondes en Python sur ma photo de 2.8 Mpx. Pour un flou de mouvement plus important, on atteint 2h de calcul, sauf qu'en l'état actuel des choses, c'est ça ou la photo part à la poubelle. Donc d'un point de vue utilisation photo (arrête moi si tu n'es pas d'accord), je ne trouve ça réaliste (quand je fais des panoramas par recomposition et fusion d'exposition avec Hugin, c'est l'ordre de grandeur).

---

Effectivement c'est assez impressionnant !

Par contre, je me demande si le flou gaussien est vraiment un bon test. Si le cas d'utilisation est de sauver une photo pour laquelle la mise au point est mauvaise, c'est plutôt un flou de type "convolution avec un disque" (ou plus précisément avec la forme du diaphragme) (il doit y avoir un mot simple pour désigner ce genre de flou mais je ne trouve pas). L'autre cas d'utilisation possible est d'aller au delà de la résolution de l'appareil (compenser la diffraction et l'éventuel filtre passe-bas devant le capteur), là ça doit être plus proche d'un flou gaussien mais avec un rayon plus petit. Bon, bref, je pense à haute voix, si mes interrogations t'aident tant mieux, sinon tu peux ignorer ;-).

Pas sûr d'avoir compris le principe de la « zone rouge ». Si je comprends bien, c'est la zone sur laquelle tu évalues les paramètres de la déconvolution, et après tu appliques la déconvolution sur toute l'image. C'est quelque chose comme ça ?

En tous cas, encore une fois félicitations pour le boulot accompli ! Continue, et continue à nous tenir (au courant|en haleine) ;-).

En fait, un vrai flou d'objo est pire qu'un disque, c'est un kernel formé par une gaussienne périodique déphasée : http://yehar.com/blog/?p=1495

J'utilise le flou gaussien parce que c'est la première image que j'ai testée et que ça me permet d'avoir un benchmark coût/qualité cohérent depuis le début de mes expérimentations. J'ai d'autres images de test avec des flous de bouger et d'objo réels. En fait, en déconvolution aveugle, la PSF est initialisée avec des 1/somme(PSF), donc peu importe le type de flou de l'image, on va s'en rapprocher par itérations successives. Tu peux accélérer la convergence en entrant manuellement une PSF approchée (c'est de la déconvolution myope et mon algo le permet aussi), mais dans des cas pourris, si la PSF approchée est trop éloignée de la réalité, l'algo risque de converger plus tard ou la solution va juste dégénérer.

Après, c'est pas idiot de faire du sur-échantillonnage, c'est ce qui est fait en astronomie. Je m'étais déjà amusé à essayer et le résultat était bon mais le coût de calcul explose. En utilisation réelle, les photos actuelles de 24 Mpx à 52 Mpx sont déjà largement sur-échantillonnées de toute façon.

La zone rouge, en fait tu as compris. En déconvolution aveugle, on évalue 4 FFT par itération/canal : 2 sur l'image défloutée, puis 2 sur la PSF. On ne peut pas couper à la FFT sur toute l'image pour les 2 premières, mais les 2 dernières n'ont pas forcément besoin d'être évaluées partout, surtout si la PSF varie spatialement (genre zone de bokeh et zones supposées nettes). Du coup, on économise du temps de calcul en masquant l'image. L'algo permet aussi de ne pas masquer, à tes risques et périls.

Idem, pour une déconvolution non aveugle, on a juste besoin de 2 FFT/itération, donc à la fin de l'algo, je rajoute 10 à 50 itérations de déconvolution de Richardson-Lucy classique à partir de la PSF estimée précédemment, ça prend 10-15 s et ça rajoute un petit boost de netteté pour pas cher.

Le vrai challenge en fait, c'est de savoir quels paramètres j'expose à l'utilisateur non-physicien pour lui permettre d'utiliser le module sans avoir à prendre un cours de traitement de signal, mais sans tomber dans une approche Apple avec une boîte noire taille unique.

[EDIT]

Explication des abbréviations :

PSF : https://fr.wikipedia.org/wiki/Fonction_d...t_du_point
FFT : https://fr.wikipedia.org/wiki/Transforma...ier_rapide
JIT : https://fr.wikipedia.org/wiki/Compilatio...vol%C3%A9e
Déconvolution : https://fr.wikipedia.org/wiki/D%C3%A9convolution