{"cells":[{"cell_type":"markdown","metadata":{"id":"I-9uxGnBMZoN"},"source":["# Installation\n","\n","- Install conda or mini-conda: https://docs.conda.io/en/latest/miniconda.html\n","- Install dependencies: conda install -c conda-forge -y opencv notebook ipywidgets matplotlib\n","- Place the file \"utils.py\" in the same location or add it to PYTHONPATH"]},{"cell_type":"code","source":["from google.colab import drive\n","drive.mount('/content/drive/')\n","\n","# Adapter ce chemin à votre hiérarchie\n","%cd /content/drive/MyDrive/uni-ete\n","%ls"],"metadata":{"id":"XLWOKNCuMcJh"},"execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"g7Pxow1HRRin"},"outputs":[],"source":["import math\n","import numpy as np\n","import cv2 as cv\n","import matplotlib.pyplot as plt\n","from utils import *\n","from ipywidgets import interact"]},{"cell_type":"markdown","metadata":{"id":"QgMXnN9wSgrS"},"source":["# Step 1: Fingerprint segmentation"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"E2gtenrUMZoP"},"outputs":[],"source":["fingerprint = cv.imread('samples/sample_1_1.png', cv.IMREAD_GRAYSCALE)\n","show(fingerprint, f'Fingerprint with size (w,h): {fingerprint.shape[::-1]}')"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"wpgIAkrBsnvw"},"outputs":[],"source":["# Calculate the local gradient (using Sobel filters)\n","\"\"\"\n"," On applique le filtre de sobel pour mettre en évidence les contours de l'empreinte digitale.\n","\"\"\"\n","gx, gy = cv.Sobel(fingerprint, cv.CV_64F, 1, 0), cv.Sobel(fingerprint, cv.CV_64F, 0, 1)\n","show((gx, 'Gx'), (gy, 'Gy'))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"LQ1s0yjIs8j4"},"outputs":[],"source":["# Calculate the magnitude of the gradient for each pixel\n","gx2, gy2 = gx**2, gy**2\n","gm = np.sqrt(gx2 + gy2)\n","show((gx2, 'Gx**2'), (gy2, 'Gy**2'), (gm, 'Gradient magnitude'))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"GPhP0Ubks_nA"},"outputs":[],"source":["# Integral over a square window\n","window = (25,25)\n","sum_gm = cv.boxFilter(gm, -1, window)\n","show(sum_gm, 'Sum of the gradient magnitude')"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"5YafVxEOtC_4"},"outputs":[],"source":["# Use a simple threshold for segmenting the fingerprint pattern\n","\n","\"\"\"\n"," On calcule un masque qui permettra de déterminer les contours extérieurs de l'empreinte\n","\"\"\"\n","\n","thr = np.max(sum_gm) * 0.2\n","mask = cv.threshold(sum_gm, thr, 255, cv.THRESH_BINARY)[1].astype(np.uint8)\n","show(fingerprint, mask, cv.merge((mask, fingerprint, fingerprint)))"]},{"cell_type":"markdown","metadata":{"id":"jJ2D5OHU6inW"},"source":["# Step 2: Estimation of local ridge orientation"]},{"cell_type":"markdown","metadata":{"id":"x_ajHsYAYwCD"},"source":["The ridge orientation is estimated as ortoghonal to the gradient orientation, averaged over a window $W$. \n","\n","$G_{xx}=\\sum_W{G_x^2}$\n","\n","$G_{yy}=\\sum_W{G_y^2}$\n","\n","$G_{xy}=\\sum_W{G_xG_y}$\n","\n","$\\theta=\\frac{\\pi}{2} + \\frac{phase(G_{xx}-G_{yy}, 2G_{xy})}{2}$\n","\n","For each orientation, we will also calculate a confidence value (strength), which measures how much all gradients in $W$ share the same orientation. \n","\n","$strength = \\frac{\\sqrt{(G_{xx}-G_{yy})^2+(2G_{xy})^2}}{G_{xx}+G_{yy}}$"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"jnyEfap-tFSw"},"outputs":[],"source":["W = (23, 23)\n","\n","# gxx = cv.GaussianBlur(gx2, W, 0)\n","# gyy = cv.GaussianBlur(gy2, W, 0)\n","\n","# gxy = cv.GaussianBlur(gx * gy, W, 0)\n","\n","gxx = cv.blur(gx2, W)\n","gyy = cv.blur(gy2, W)\n","\n","gxy = cv.blur(gx * gy, W)\n","\n","orientations = (-1) * 0.5 * np.arctan2(2 * gxy, gxx - gyy) + np.pi / 2\n","strengths = np.sqrt((gxx - gyy) ** 2 + (2 * gxy) ** 2) / (gxx + gyy + 1e-5)\n","\n","show(draw_orientations(fingerprint, orientations, strengths, mask, 1, 16), 'Orientation image')"]},{"cell_type":"markdown","metadata":{"id":"SDHYsNmJ7TaV"},"source":["# Step 3: Estimation of local ridge frequency"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"ezc3YDWXCTEm"},"outputs":[],"source":["# region = cv.bitwise_and(fingerprint, fingerprint, mask=mask\n","region = fingerprint[80:200, 100:220]\n","show(region)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Pvv7fS0KQpJ3"},"outputs":[],"source":["# smoothed = cv.medianBlur(region, 5)\n","smoothed = cv.GaussianBlur(region, (5, 5), 0)\n","\n","xs = np.sum(smoothed, axis=1)\n","print(xs)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"8-hQrph0DSZO"},"outputs":[],"source":["x = np.arange(region.shape[0])\n","f, axarr = plt.subplots(1,2, sharey = True)\n","axarr[0].imshow(region,cmap='gray')\n","axarr[1].plot(xs, x)\n","# axarr[1].plot(x, xs)\n","axarr[1].set_ylim(region.shape[0]-1,0)\n","plt.show()"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"zrM-B9JKRXth"},"outputs":[],"source":["# local_maxima = (np.diff(np.sign(np.diff(xs))) < 0).nonzero()[0] + 1\n","local_maxima = np.where((xs[:-2] < xs[1:-1]) & (xs[1:-1] > xs[2:]))[0] + 1\n","print(local_maxima)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"pasX1AILRQdq"},"outputs":[],"source":["x = np.arange(region.shape[0])\n","plt.plot(x, xs)\n","plt.xticks(local_maxima)\n","plt.grid(True, axis='x')\n","plt.show()"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"w6CPgUQgtYCo"},"outputs":[],"source":["# sum = 0\n","# deltas = []\n","\n","# for i in range(len(local_maxima) - 1):\n","# deltas.append(local_maxima[i + 1] - local_maxima[i])\n","# sum += local_maxima[i + 1] - local_maxima[i]\n","\n","# distances = sum / len(deltas)\n","# print(distances)\n","\n","distances = local_maxima[1:] - local_maxima[:-1]\n","print(distances)"]},{"cell_type":"markdown","source":["<!-- -->"],"metadata":{"id":"SaobQYWyE_Oc"}},{"cell_type":"code","execution_count":null,"metadata":{"id":"QlfcrCgPVozS"},"outputs":[],"source":["ridge_period = distances.mean()\n","print(ridge_period)"]},{"cell_type":"markdown","metadata":{"id":"QShcuAYY7wd0"},"source":["# Step 4: Fingerprint enhancement"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"htZDESSYtesf"},"outputs":[],"source":["# Create the filter bank\n","or_count = 8\n","gabor_bank = [gabor_kernel(ridge_period, o) for o in np.arange(0, np.pi, np.pi/or_count)]"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"l7Ldq1VNtgEQ"},"outputs":[],"source":["show(*gabor_bank)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"gmeTXLIvthI_"},"outputs":[],"source":["\n","nf = 255-fingerprint\n","all_filtered = np.array([cv.filter2D(nf, cv.CV_32F, f) for f in gabor_bank])\n","show(nf, *all_filtered)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"usS93oHDt7VQ"},"outputs":[],"source":["y_coords, x_coords = np.indices(fingerprint.shape)\n","\n","orientation_idx = np.round(((orientations % np.pi) / np.pi) * or_count).astype(np.int32) % or_count\n","filtered = all_filtered[orientation_idx, y_coords, x_coords]\n","enhanced = mask & np.clip(filtered, 0, 255).astype(np.uint8)\n","show(fingerprint, filtered, enhanced)"]},{"cell_type":"markdown","metadata":{"id":"LQvHUusjyBfR"},"source":["# Step 5: Detection of minutiae positions"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"jmCPCdMFuCMe"},"outputs":[],"source":["# Binarization\n","_, ridge_lines = cv.threshold(enhanced, 32, 255, cv.THRESH_BINARY)\n","show(enhanced, ridge_lines)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"OtEgLkx7uE_e"},"outputs":[],"source":["# Thinning\n","skeleton = cv.ximgproc.thinning(ridge_lines, thinningType = cv.ximgproc.THINNING_GUOHALL)\n","show(ridge_lines, skeleton)"]},{"cell_type":"code","source":["def compute_crossing_number(values):\n"," return np.count_nonzero(values < np.roll(values, -1))"],"metadata":{"id":"pIIvgsPueUj3"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["cn_filter = np.array([[ 1, 2, 4],\n"," [128, 0, 8],\n"," [ 64, 32, 16]\n"," ])"],"metadata":{"id":"5NYXa5ybeZbP"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["all_8_neighborhoods = [np.array([int(d) for d in f'{x:08b}'])[::-1] for x in range(256)]\n","cn_lut = np.array([compute_crossing_number(x) for x in all_8_neighborhoods]).astype(np.uint8)"],"metadata":{"id":"wwmjIaSUecxO"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Skeleton: from 0/255 to 0/1 values\n","skeleton01 = np.where(skeleton!=0, 1, 0).astype(np.uint8)\n","# Apply the filter to encode the 8-neighborhood of each pixel into a byte [0,255]\n","neighborhood_values = cv.filter2D(skeleton01, -1, cn_filter, borderType = cv.BORDER_CONSTANT)\n","# Apply the lookup table to obtain the crossing number of each pixel from the byte value of its neighborhood\n","cn = cv.LUT(neighborhood_values, cn_lut)\n","# Keep only crossing numbers on the skeleton\n","cn[skeleton==0] = 0"],"metadata":{"id":"nz7u55ZOef_r"},"execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"4NzceQaJuMk_"},"outputs":[],"source":["# Find minutiae.\n","# List format: [(x: int, y: int, type_terminaison: bool), ...]\n","minutiae = [(x,y,cn[y,x]==1) for y, x in zip(*np.where(np.isin(cn, [1,3])))]"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"_NUdi_JAuPVH"},"outputs":[],"source":["show(draw_minutiae(fingerprint, minutiae), skeleton, draw_minutiae(skeleton, minutiae))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"VLWwRQBZuQo_"},"outputs":[],"source":["# Create mask.\n","mask_edges = cv.distanceTransform(cv.copyMakeBorder(mask, 1, 1, 1, 1, cv.BORDER_CONSTANT), cv.DIST_C, 3)[1:-1,1:-1]\n","show(mask, mask_edges)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"4Sb6iyRFuUF_"},"outputs":[],"source":["# Filter minutiae.\n","filtered_minutiae = list(filter(lambda m: mask_edges[m[1], m[0]]>10, minutiae))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"99BiKiH2uV5H"},"outputs":[],"source":["show(draw_minutiae(fingerprint, filtered_minutiae), skeleton, draw_minutiae(skeleton, filtered_minutiae))"]},{"cell_type":"markdown","metadata":{"id":"_ZikeJ62zp-9"},"source":["# Step 6: Estimation of minutiae directions"]},{"cell_type":"code","source":["def compute_next_ridge_following_directions(previous_direction, values):\n"," next_positions = np.argwhere(values!=0).ravel().tolist()\n"," if len(next_positions) > 0 and previous_direction != 8:\n"," # There is a previous direction: return all the next directions, sorted according to the distance from it,\n"," # except the direction, if any, that corresponds to the previous position\n"," next_positions.sort(key = lambda d: 4 - abs(abs(d - previous_direction) - 4))\n"," if next_positions[-1] == (previous_direction + 4) % 8: # the direction of the previous position is the opposite one\n"," next_positions = next_positions[:-1] # removes it\n"," return next_positions"],"metadata":{"id":"0--AIiPl1tjo"},"execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"X24IjohPzp--"},"outputs":[],"source":["r2 = 2**0.5 # sqrt(2)\n","\n","\"\"\"\n","Si nous sommes en p et que nous voulons aller en 7, nous bougeons que d'un pixel,\n","donc la valeur sera (-1, 0, 1).\n","Si nous sommes en p et que nous voulons aller en 4, nous bougeons de 2 pixel via diagonale,\n","donc la valeur sera (1,1, sqrt(2)).\n","\"\"\"\n","# The eight possible (x, y) offsets with each corresponding Euclidean distance\n","xy_steps = [(-1,-1,r2),( 0,-1,1),( 1,-1,r2),( 1, 0,1),( 1, 1,r2),( 0, 1,1),(-1, 1,r2),(-1, 0,1)]\n","\n","# TODO\n","# LUT: for each 8-neighborhood and each previous direction [0,8],\n","# where 8 means \"none\", provides the list of possible directions\n","nd_lut = [[compute_next_ridge_following_directions(pd, x) for pd in range(9)] for x in all_8_neighborhoods]"]},{"cell_type":"code","source":["def follow_ridge_and_compute_angle(x, y, d = 8):\n"," px, py = x, y\n"," length = 0.0\n"," while length < 20: # max length followed\n"," next_directions = nd_lut[neighborhood_values[py,px]][d]\n"," if len(next_directions) == 0:\n"," break\n"," # Need to check ALL possible next directions\n"," if (any(cn[py + xy_steps[nd][1], px + xy_steps[nd][0]] != 2 for nd in next_directions)):\n"," break # another minutia found: we stop here\n"," # Only the first direction has to be followed\n"," d = next_directions[0]\n"," ox, oy, l = xy_steps[d]\n"," px += ox ; py += oy ; length += l\n"," # check if the minimum length for a valid direction has been reached\n"," return math.atan2(-py+y, px-x) if length >= 10 else None"],"metadata":{"id":"sXnx4-WE14hC"},"execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"4T5fvHFMzp--"},"outputs":[],"source":["# List format: [(x: int, y: int, type_terminaison: bool, direction: int), ...]\n","valid_minutiae = []\n","for x, y, term in filtered_minutiae:\n"," d = None\n"," if term: # termination: simply follow and compute the direction\n"," d = follow_ridge_and_compute_angle(x, y)\n"," else: # bifurcation: follow each of the three branches\n"," dirs = nd_lut[neighborhood_values[y,x]][8] # 8 means: no previous direction\n"," if len(dirs)==3: # only if there are exactly three branches\n"," angles = [follow_ridge_and_compute_angle(x+xy_steps[d][0], y+xy_steps[d][1], d) for d in dirs]\n"," if all(a is not None for a in angles):\n"," a1, a2 = min(((angles[i], angles[(i+1)%3]) for i in range(3)), key=lambda t: angle_abs_difference(t[0], t[1]))\n"," d = angle_mean(a1, a2)\n"," if d is not None:\n"," valid_minutiae.append( (x, y, term, d) )"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"8TccPjw_zp--"},"outputs":[],"source":["show(draw_minutiae(fingerprint, valid_minutiae))"]},{"cell_type":"markdown","metadata":{"id":"gAnKIS9pzp-_"},"source":["# Step 7: Creation of local structures"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"0dm6Jlvizp-_"},"outputs":[],"source":["# Compute the cell coordinates of a generic local structure\n","mcc_radius = 70\n","mcc_size = 16\n","\n","g = 2 * mcc_radius / mcc_size\n","x = np.arange(mcc_size)*g - (mcc_size/2)*g + g/2\n","y = x[..., np.newaxis]\n","iy, ix = np.nonzero(x**2 + y**2 <= mcc_radius**2)\n","ref_cell_coords = np.column_stack((x[ix], x[iy]))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"g0-NqFYezp-_"},"outputs":[],"source":["mcc_sigma_s = 7.0\n","mcc_tau_psi = 400.0\n","mcc_mu_psi = 1e-2\n","\n","def Gs(t_sqr):\n"," \"\"\"Gaussian function with zero mean and mcc_sigma_s standard deviation, see eq. (7) in MCC paper\"\"\"\n"," return np.exp(-0.5 * t_sqr / (mcc_sigma_s**2)) / (math.tau**0.5 * mcc_sigma_s)\n","\n","def Psi(v):\n"," \"\"\"Sigmoid function that limits the contribution of dense minutiae clusters, see eq. (4)-(5) in MCC paper\"\"\"\n"," return 1. / (1. + np.exp(-mcc_tau_psi * (v - mcc_mu_psi)))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"zP32k1Drzp-_"},"outputs":[],"source":["# n: number of minutiae\n","# c: number of cells in a local structure\n","\n","xyd = np.array([(x,y,d) for x,y,_,d in valid_minutiae]) # matrix with all minutiae coordinates and directions (n x 3)\n","\n","# rot: n x 2 x 2 (rotation matrix for each minutia)\n","d_cos, d_sin = np.cos(xyd[:,2]).reshape((-1,1,1)), np.sin(xyd[:,2]).reshape((-1,1,1))\n","rot = np.block([[d_cos, d_sin], [-d_sin, d_cos]])\n","\n","# rot@ref_cell_coords.T : n x 2 x c\n","# xy : n x 2\n","xy = xyd[:,:2]\n","# cell_coords: n x c x 2 (cell coordinates for each local structure)\n","cell_coords = np.transpose(rot@ref_cell_coords.T + xy[:,:,np.newaxis],[0,2,1])\n","\n","# cell_coords[:,:,np.newaxis,:] : n x c x 1 x 2\n","# xy : (1 x 1) x n x 2\n","# cell_coords[:,:,np.newaxis,:] - xy : n x c x n x 2\n","# dists: n x c x n (for each cell of each local structure, the distance from all minutiae)\n","dists = np.sum((cell_coords[:,:,np.newaxis,:] - xy)**2, -1)\n","\n","# cs : n x c x n (the spatial contribution of each minutia to each cell of each local structure)\n","cs = Gs(dists)\n","diag_indices = np.arange(cs.shape[0])\n","cs[diag_indices,:,diag_indices] = 0 # remove the contribution of each minutia to its own cells\n","\n","# local_structures : n x c (cell values for each local structure)\n","local_structures = Psi(np.sum(cs, -1))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"GCn07G_0zp-_"},"outputs":[],"source":["@interact(i=(0,len(valid_minutiae)-1))\n","def test(i=0):\n"," show(draw_minutiae_and_cylinder(fingerprint, ref_cell_coords, valid_minutiae, local_structures, i))"]},{"cell_type":"markdown","metadata":{"id":"h9G3SosVzp-_"},"source":["# Step 8: Fingerprint comparison"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"jPULw0vMzp_A"},"outputs":[],"source":["print(f\"\"\"Fingerprint image: {fingerprint.shape[1]}x{fingerprint.shape[0]} pixels\n","Minutiae: {len(valid_minutiae)}\n","Local structures: {local_structures.shape}\"\"\")"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"GEENB_yvzp_A"},"outputs":[],"source":["f1, m1, ls1 = fingerprint, valid_minutiae, local_structures"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"_n5L3L22zp_A"},"outputs":[],"source":["ofn = 'samples/sample_1_2' # Fingerprint of the same finger\n","#ofn = 'samples/sample_2' # Fingerprint of a different finger\n","# f2 = cv.imread(f'{ofn}.png', cv.IMREAD_GRAYSCALE)\n","# m2 = ???\n","# ls2 = ???\n","f2, (m2, ls2) = cv.imread(f'{ofn}.png', cv.IMREAD_GRAYSCALE), np.load(f'{ofn}.npz', allow_pickle=True).values()"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"hoYOwyAJzp_B"},"outputs":[],"source":["# Compute\n","# Compute all pairwise normalized Euclidean distances between local structures in v1 and v2\n","# ls1 : n1 x c\n","# ls1[:,np.newaxis,:] : n1 x 1 x c\n","# ls2 : (1 x) n2 x c\n","# ls1[:,np.newaxis,:] - ls2 : n1 x n2 x c\n","# dists : n1 x n2\n","dists = np.linalg.norm(ls1[:,np.newaxis,:] - ls2, axis = -1)\n","dists /= np.linalg.norm(ls1, axis = 1)[:,np.newaxis] + np.linalg.norm(ls2, axis = 1) # Normalize as in eq. (17) of MCC paper"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"lYSeRY7-zp_B"},"outputs":[],"source":["num_p = 5 # For simplicity: a fixed number of pairs\n","pairs = np.unravel_index(np.argpartition(dists, num_p, None)[:num_p], dists.shape)\n","score = 1 - np.mean(dists[pairs[0], pairs[1]]) # See eq. (23) in MCC paper\n","print(f'Comparison score: {score:.2f}')"]},{"cell_type":"code","source":["@interact(i = (0,len(pairs[0])-1), show_local_structures = False)\n","def show_pairs(i=0, show_local_structures = False):\n"," show(draw_match_pairs(f1, m1, ls1, f2, m2, ls2, ref_cell_coords, pairs, i, show_local_structures))"],"metadata":{"id":"QUO78cF3fwDs"},"execution_count":null,"outputs":[]}],"metadata":{"colab":{"private_outputs":true,"provenance":[]},"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.12.3"},"vscode":{"interpreter":{"hash":"ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"}}},"nbformat":4,"nbformat_minor":0}
\ No newline at end of file
%% Cell type:markdown id: tags:
# Installation
- Install conda or mini-conda: https://docs.conda.io/en/latest/miniconda.html