Shape from focus¶
Avez-vous déjà observé une fourmi au microscope ? Le problème est que la profondeur de champ d'un microscope est en général assez faible, limitant la possibilité d'observer un objet en 3d. Seule une faible portion de l'image est alors nette. En revanche, le réglage du focus est souvent très fin. Il est donc possible de scanner verticalement l'objet et d'observer successivement les zones nettes de l'objet, correspondant à différentes altitudes.
L'idée est de prendre une série de photos à des hauteur différentes pour ensuite les combiner numériquement et ainsi créer une image finale ayant une grande profondeur de champ.
Remarque: C'est une technique existante et déjà étudiée. Il existe plusieurs articles concernant ce sujet (google scholar).
J'ai voulu tester cette idée sur une fleur, étant donné qu'il est plus difficile d'attraper une fourmi. La taille de la fleur est inférieure au millimètre.
Avec un peu de patience on obtient un cube d'une vingtaine d'images correspondant à différente coupe de l'objet.
Nous allons voir comment détecter numériquement les zones nettes à l'intérieur de chaque image. Cette mesure est obtenue en évaluant la dérivée seconde (Laplacien en 2d) des valeurs de gris de l'image. Une zone nette présente des variations rapides, donc des courbures importantes, alors qu'une zone floue est lisse, avec des courbures faibles.
On utilise la librairie ndimage de Scipy (doc):
from scipy import ndimage
Chargement des photos¶
n = 20 # nombre de photos
images = []
for i in range(1, n+1):
filename = './data_shapefromfocus/'+'focus%02d.jpg'%i
images.append( ndimage.imread(filename) )
Exemple de 3 images succesives: (voir aussi l'image n°5)
On travail sur les valeurs de gris seulement. La fonction suivante permet de convertir les valeurs RGB en niveau de gris:
def rgb2gray(rgb):
return np.dot(rgb[...,:3], [0.299, 0.587, 0.114])
plt.figure( figsize=(10, 12) );
plt.imshow( rgb2gray( images[6] ), cmap=plt.cm.plasma);
plt.axis('off');
On souhaite évaluer le degrée de netteté de chaque pixels. Pour cela l'idée est de calculer la dérivée seconde (le Laplacien) en chaque point des images. Le filtre gaussian_laplace(I, sigma) permet de réaliser cela (doc). Comme la courbure peut être négative ou positive on calcul le carré du laplacien.
Le paramètre important est la largeur du filtre $\sigma$. On le trouve par essais-erreur.
sigma = 4
laplaciens = []
for I in images:
Igray = rgb2gray( I )
curv = ndimage.filters.gaussian_laplace( Igray , sigma)
curv2 = curv**2
laplaciens.append( curv2 )
Effet du filtre¶
En appliquant le filtre à une image ayant un unique pixel "noir" (=1), on obtient la fonction du filtre. L'application du filtre consiste au produit de convolution de cette fonction avec l'image. La courbe suivante est le profil passant par le centre.
plt.figure( figsize=(10, 12) );
plt.imshow( laplaciens[4], cmap=plt.cm.plasma, vmin=0, vmax=1 );
plt.axis('off');
Profils transverses¶
On regarde maintenant des profils dans la profondeur du cube, pour quelques pixels au hazard. Ceci permet de vérifier que les courbes obtenues onts bien un unique maximum, indiquant l'image la plus nette pour ce pixel.
Lx, Ly, nColors = np.shape( images[0] )
points = [ (np.random.randint(Lx), np.random.randint(Ly)) for i in range(6) ]
plt.figure(figsize=(12, 4))
for pixel in points:
profil = [ L[pixel] for L in laplaciens]
plt.plot( profil ); plt.xlabel('n° image');
Remarque: Le pic du maximum correspond souvent à une image unique, l'image suivante et précedente étant complètement flou... il n'y a donc pas de recouvrement entre les images. Pour obtenir un meilleur résultat, il faudrait un recouvrement de nettété entre les images, c'est-à-dire plus d'images...
On cherche ensuite l'indice du maximum pour chaque pixel:
cArgMax = np.argmax( np.array(laplaciens), 0 )
plt.figure( figsize=(10, 12) );
cmap = plt.get_cmap('jet', len(images));
cax = plt.imshow( cArgMax, cmap=cmap);
plt.colorbar(cax, orientation='horizontal')
plt.axis('off');
On filtre ensuite par la mediane. Là encore la taille du filtre est ajustée par essais.
n = 8
cMax = ndimage.filters.median_filter(cArgMax, footprint=np.ones((n, n)) )
cMax = np.round( cMax )
plt.figure( figsize=(10, 12) );
plt.imshow( cMax, cmap=cmap);
plt.axis('off');
Puis on forme une image unique en choisiant le pixel le plus nette dans chaque image:
I_merged = images[0]
for k, I in enumerate( images ):
I_merged[ cMax==k ] = images[k][ cMax==k ]
plt.figure( figsize=(10, 12) );
plt.imshow( I_merged );
plt.axis('off');
Le microscope électronique à balayage est quand même beaucoup mieux pour faire des photos de fourmis: https://commons.wikimedia.org/wiki/File:Ant_with_gear.jpg
