import math import numpy as np import scipy.ndimage as ndi from .._shared.utils import check_nD, _supported_float_type from ..feature.util import DescriptorExtractor, FeatureDetector from .._shared.filters import gaussian from ..transform import rescale from ..util import img_as_float from ._sift import _local_max, _ori_distances, _update_histogram def _edgeness(hxx, hyy, hxy): """Compute edgeness (eq. 18 of Otero et. al. IPOL paper)""" trace = hxx + hyy determinant = hxx * hyy - hxy * hxy return (trace * trace) / determinant def _sparse_gradient(vol, positions): """Gradient of a 3D volume at the provided `positions`. For SIFT we only need the gradient at specific positions and do not need the gradient at the edge positions, so can just use this simple implementation instead of numpy.gradient. """ p0 = positions[..., 0] p1 = positions[..., 1] p2 = positions[..., 2] g0 = vol[p0 + 1, p1, p2] - vol[p0 - 1, p1, p2] g0 *= 0.5 g1 = vol[p0, p1 + 1, p2] - vol[p0, p1 - 1, p2] g1 *= 0.5 g2 = vol[p0, p1, p2 + 1] - vol[p0, p1, p2 - 1] g2 *= 0.5 return g0, g1, g2 def _hessian(d, positions): """Compute the non-redundant 3D Hessian terms at the requested positions. Source: "Anatomy of the SIFT Method" p.380 (13) """ p0 = positions[..., 0] p1 = positions[..., 1] p2 = positions[..., 2] two_d0 = 2 * d[p0, p1, p2] # 0 = row, 1 = col, 2 = octave h00 = d[p0 - 1, p1, p2] + d[p0 + 1, p1, p2] - two_d0 h11 = d[p0, p1 - 1, p2] + d[p0, p1 + 1, p2] - two_d0 h22 = d[p0, p1, p2 - 1] + d[p0, p1, p2 + 1] - two_d0 h01 = 0.25 * ( d[p0 + 1, p1 + 1, p2] - d[p0 - 1, p1 + 1, p2] - d[p0 + 1, p1 - 1, p2] + d[p0 - 1, p1 - 1, p2] ) h02 = 0.25 * ( d[p0 + 1, p1, p2 + 1] - d[p0 + 1, p1, p2 - 1] + d[p0 - 1, p1, p2 - 1] - d[p0 - 1, p1, p2 + 1] ) h12 = 0.25 * ( d[p0, p1 + 1, p2 + 1] - d[p0, p1 + 1, p2 - 1] + d[p0, p1 - 1, p2 - 1] - d[p0, p1 - 1, p2 + 1] ) return (h00, h11, h22, h01, h02, h12) def _offsets(grad, hess): """Compute position refinement offsets from gradient and Hessian. This is equivalent to np.linalg.solve(-H, J) where H is the Hessian matrix and J is the gradient (Jacobian). This analytical solution is adapted from (BSD-licensed) C code by Otero et. al (see SIFT docstring References). """ h00, h11, h22, h01, h02, h12 = hess g0, g1, g2 = grad det = h00 * h11 * h22 det -= h00 * h12 * h12 det -= h01 * h01 * h22 det += 2 * h01 * h02 * h12 det -= h02 * h02 * h11 aa = (h11 * h22 - h12 * h12) / det ab = (h02 * h12 - h01 * h22) / det ac = (h01 * h12 - h02 * h11) / det bb = (h00 * h22 - h02 * h02) / det bc = (h01 * h02 - h00 * h12) / det cc = (h00 * h11 - h01 * h01) / det offset0 = -aa * g0 - ab * g1 - ac * g2 offset1 = -ab * g0 - bb * g1 - bc * g2 offset2 = -ac * g0 - bc * g1 - cc * g2 return np.stack((offset0, offset1, offset2), axis=-1) class SIFT(FeatureDetector, DescriptorExtractor): """SIFT feature detection and descriptor extraction. Parameters ---------- upsampling : int, optional Prior to the feature detection the image is upscaled by a factor of 1 (no upscaling), 2 or 4. Method: Bi-cubic interpolation. n_octaves : int, optional Maximum number of octaves. With every octave the image size is halved and the sigma doubled. The number of octaves will be reduced as needed to keep at least 12 pixels along each dimension at the smallest scale. n_scales : int, optional Maximum number of scales in every octave. sigma_min : float, optional The blur level of the seed image. If upsampling is enabled sigma_min is scaled by factor 1/upsampling sigma_in : float, optional The assumed blur level of the input image. c_dog : float, optional Threshold to discard low contrast extrema in the DoG. It's final value is dependent on n_scales by the relation: final_c_dog = (2^(1/n_scales)-1) / (2^(1/3)-1) * c_dog c_edge : float, optional Threshold to discard extrema that lie in edges. If H is the Hessian of an extremum, its "edgeness" is described by tr(H)²/det(H). If the edgeness is higher than (c_edge + 1)²/c_edge, the extremum is discarded. n_bins : int, optional Number of bins in the histogram that describes the gradient orientations around keypoint. lambda_ori : float, optional The window used to find the reference orientation of a keypoint has a width of 6 * lambda_ori * sigma and is weighted by a standard deviation of 2 * lambda_ori * sigma. c_max : float, optional The threshold at which a secondary peak in the orientation histogram is accepted as orientation lambda_descr : float, optional The window used to define the descriptor of a keypoint has a width of 2 * lambda_descr * sigma * (n_hist+1)/n_hist and is weighted by a standard deviation of lambda_descr * sigma. n_hist : int, optional The window used to define the descriptor of a keypoint consists of n_hist * n_hist histograms. n_ori : int, optional The number of bins in the histograms of the descriptor patch. Attributes ---------- delta_min : float The sampling distance of the first octave. It's final value is 1/upsampling. float_dtype : type The datatype of the image. scalespace_sigmas : (n_octaves, n_scales + 3) array The sigma value of all scales in all octaves. keypoints : (N, 2) array Keypoint coordinates as ``(row, col)``. positions : (N, 2) array Subpixel-precision keypoint coordinates as ``(row, col)``. sigmas : (N,) array The corresponding sigma (blur) value of a keypoint. scales : (N,) array The corresponding scale of a keypoint. orientations : (N,) array The orientations of the gradient around every keypoint. octaves : (N,) array The corresponding octave of a keypoint. descriptors : (N, n_hist*n_hist*n_ori) array The descriptors of a keypoint. Notes ----- The SIFT algorithm was developed by David Lowe [1]_, [2]_ and later patented by the University of British Columbia. Since the patent expired in 2020 it's free to use. The implementation here closely follows the detailed description in [3]_, including use of the same default parameters. References ---------- .. [1] D.G. Lowe. "Object recognition from local scale-invariant features", Proceedings of the Seventh IEEE International Conference on Computer Vision, 1999, vol.2, pp. 1150-1157. :DOI:`10.1109/ICCV.1999.790410` .. [2] D.G. Lowe. "Distinctive Image Features from Scale-Invariant Keypoints", International Journal of Computer Vision, 2004, vol. 60, pp. 91–110. :DOI:`10.1023/B:VISI.0000029664.99615.94` .. [3] I. R. Otero and M. Delbracio. "Anatomy of the SIFT Method", Image Processing On Line, 4 (2014), pp. 370–396. :DOI:`10.5201/ipol.2014.82` Examples -------- >>> from skimage.feature import SIFT, match_descriptors >>> from skimage.data import camera >>> from skimage.transform import rotate >>> img1 = camera() >>> img2 = rotate(camera(), 90) >>> detector_extractor1 = SIFT() >>> detector_extractor2 = SIFT() >>> detector_extractor1.detect_and_extract(img1) >>> detector_extractor2.detect_and_extract(img2) >>> matches = match_descriptors(detector_extractor1.descriptors, ... detector_extractor2.descriptors, ... max_ratio=0.6) >>> matches[10:15] array([[ 10, 412], [ 11, 417], [ 12, 407], [ 13, 411], [ 14, 406]]) >>> detector_extractor1.keypoints[matches[10:15, 0]] array([[ 95, 214], [ 97, 211], [ 97, 218], [102, 215], [104, 218]]) >>> detector_extractor2.keypoints[matches[10:15, 1]] array([[297, 95], [301, 97], [294, 97], [297, 102], [293, 104]]) """ def __init__( self, upsampling=2, n_octaves=8, n_scales=3, sigma_min=1.6, sigma_in=0.5, c_dog=0.04 / 3, c_edge=10, n_bins=36, lambda_ori=1.5, c_max=0.8, lambda_descr=6, n_hist=4, n_ori=8, ): if upsampling in [1, 2, 4]: self.upsampling = upsampling else: raise ValueError("upsampling must be 1, 2 or 4") self.n_octaves = n_octaves self.n_scales = n_scales self.sigma_min = sigma_min / upsampling self.sigma_in = sigma_in self.c_dog = (2 ** (1 / n_scales) - 1) / (2 ** (1 / 3) - 1) * c_dog self.c_edge = c_edge self.n_bins = n_bins self.lambda_ori = lambda_ori self.c_max = c_max self.lambda_descr = lambda_descr self.n_hist = n_hist self.n_ori = n_ori self.delta_min = 1 / upsampling self.float_dtype = None self.scalespace_sigmas = None self.keypoints = None self.positions = None self.sigmas = None self.scales = None self.orientations = None self.octaves = None self.descriptors = None @property def deltas(self): """The sampling distances of all octaves""" deltas = self.delta_min * np.power( 2, np.arange(self.n_octaves), dtype=self.float_dtype ) return deltas def _set_number_of_octaves(self, image_shape): size_min = 12 # minimum size of last octave s0 = min(image_shape) * self.upsampling max_octaves = int(math.log2(s0 / size_min) + 1) if max_octaves < self.n_octaves: self.n_octaves = max_octaves def _create_scalespace(self, image): """Source: "Anatomy of the SIFT Method" Alg. 1 Construction of the scalespace by gradually blurring (scales) and downscaling (octaves) the image. """ scalespace = [] if self.upsampling > 1: image = rescale(image, self.upsampling, order=1) # smooth to sigma_min, assuming sigma_in image = gaussian( image, sigma=self.upsampling * math.sqrt(self.sigma_min**2 - self.sigma_in**2), mode='reflect', ) # Eq. 10: sigmas.shape = (n_octaves, n_scales + 3). # The three extra scales are: # One for the differences needed for DoG and two auxiliary # images (one at either end) for peak_local_max with exclude # border = True (see Fig. 5) # The smoothing doubles after n_scales steps. tmp = np.power(2, np.arange(self.n_scales + 3) / self.n_scales) tmp *= self.sigma_min # all sigmas for the gaussian scalespace sigmas = self.deltas[:, np.newaxis] / self.deltas[0] * tmp[np.newaxis, :] self.scalespace_sigmas = sigmas # Eq. 7: Gaussian smoothing depends on difference with previous sigma # gaussian_sigmas.shape = (n_octaves, n_scales + 2) var_diff = np.diff(sigmas * sigmas, axis=1) gaussian_sigmas = np.sqrt(var_diff) / self.deltas[:, np.newaxis] # one octave is represented by a 3D image with depth (n_scales+x) for o in range(self.n_octaves): # Temporarily put scales axis first so octave[i] is C-contiguous # (this makes Gaussian filtering faster). octave = np.empty( (self.n_scales + 3,) + image.shape, dtype=self.float_dtype, order='C' ) octave[0] = image for s in range(1, self.n_scales + 3): # blur new scale assuming sigma of the last one gaussian( octave[s - 1], sigma=gaussian_sigmas[o, s - 1], mode='reflect', out=octave[s], ) # move scales to last axis as expected by other methods scalespace.append(np.moveaxis(octave, 0, -1)) if o < self.n_octaves - 1: # downscale the image by taking every second pixel image = octave[self.n_scales][::2, ::2] return scalespace def _inrange(self, a, dim): return ( (a[:, 0] > 0) & (a[:, 0] < dim[0] - 1) & (a[:, 1] > 0) & (a[:, 1] < dim[1] - 1) ) def _find_localize_evaluate(self, dogspace, img_shape): """Source: "Anatomy of the SIFT Method" Alg. 4-9 1) first find all extrema of a (3, 3, 3) neighborhood 2) use second order Taylor development to refine the positions to sub-pixel precision 3) filter out extrema that have low contrast and lie on edges or close to the image borders """ extrema_pos = [] extrema_scales = [] extrema_sigmas = [] threshold = self.c_dog * 0.8 for o, (octave, delta) in enumerate(zip(dogspace, self.deltas)): # find extrema keys = _local_max(np.ascontiguousarray(octave), threshold) if keys.size == 0: extrema_pos.append(np.empty((0, 2))) continue # localize extrema oshape = octave.shape refinement_iterations = 5 offset_max = 0.6 for i in range(refinement_iterations): if i > 0: # exclude any keys that have moved out of bounds keys = keys[self._inrange(keys, oshape), :] # Jacobian and Hessian of all extrema grad = _sparse_gradient(octave, keys) hess = _hessian(octave, keys) # solve for offset of the extremum off = _offsets(grad, hess) if i == refinement_iterations - 1: break # offset is too big and an increase would not bring us out of # bounds wrong_position_pos = np.logical_and( off > offset_max, keys + 1 < tuple([a - 1 for a in oshape]) ) wrong_position_neg = np.logical_and(off < -offset_max, keys - 1 > 0) if not np.any(np.logical_or(wrong_position_neg, wrong_position_pos)): break keys[wrong_position_pos] += 1 keys[wrong_position_neg] -= 1 # mask for all extrema that have been localized successfully finished = np.all(np.abs(off) < offset_max, axis=1) keys = keys[finished] off = off[finished] grad = [g[finished] for g in grad] # value of extremum in octave vals = octave[keys[:, 0], keys[:, 1], keys[:, 2]] # values at interpolated point w = vals for i in range(3): w += 0.5 * grad[i] * off[:, i] h00, h11, h01 = hess[0][finished], hess[1][finished], hess[3][finished] sigmaratio = self.scalespace_sigmas[0, 1] / self.scalespace_sigmas[0, 0] # filter for contrast, edgeness and borders contrast_threshold = self.c_dog contrast_filter = np.abs(w) > contrast_threshold edge_threshold = np.square(self.c_edge + 1) / self.c_edge edge_response = _edgeness( h00[contrast_filter], h11[contrast_filter], h01[contrast_filter] ) edge_filter = np.abs(edge_response) <= edge_threshold keys = keys[contrast_filter][edge_filter] off = off[contrast_filter][edge_filter] yx = ((keys[:, :2] + off[:, :2]) * delta).astype(self.float_dtype) sigmas = self.scalespace_sigmas[o, keys[:, 2]] * np.power( sigmaratio, off[:, 2] ) border_filter = np.all( np.logical_and( (yx - sigmas[:, np.newaxis]) > 0.0, (yx + sigmas[:, np.newaxis]) < img_shape, ), axis=1, ) extrema_pos.append(yx[border_filter]) extrema_scales.append(keys[border_filter, 2]) extrema_sigmas.append(sigmas[border_filter]) octave_indices = np.concatenate( [np.full(len(p), i) for i, p in enumerate(extrema_pos)] ) if len(octave_indices) == 0: raise RuntimeError( "SIFT found no features. Try passing in an image containing " "greater intensity contrasts between adjacent pixels." ) extrema_pos = np.concatenate(extrema_pos) extrema_scales = np.concatenate(extrema_scales) extrema_sigmas = np.concatenate(extrema_sigmas) return extrema_pos, extrema_scales, extrema_sigmas, octave_indices def _fit(self, h): """Refine the position of the peak by fitting it to a parabola""" return (h[0] - h[2]) / (2 * (h[0] + h[2] - 2 * h[1])) def _compute_orientation( self, positions_oct, scales_oct, sigmas_oct, octaves, gaussian_scalespace ): """Source: "Anatomy of the SIFT Method" Alg. 11 Calculates the orientation of the gradient around every keypoint """ gradient_space = [] # list for keypoints that have more than one reference orientation keypoint_indices = [] keypoint_angles = [] keypoint_octave = [] orientations = np.zeros_like(sigmas_oct, dtype=self.float_dtype) key_count = 0 for o, (octave, delta) in enumerate(zip(gaussian_scalespace, self.deltas)): gradient_space.append(np.gradient(octave)) in_oct = octaves == o if not np.any(in_oct): continue positions = positions_oct[in_oct] scales = scales_oct[in_oct] sigmas = sigmas_oct[in_oct] oshape = octave.shape[:2] # convert to octave's dimensions yx = positions / delta sigma = sigmas / delta # dimensions of the patch radius = 3 * self.lambda_ori * sigma p_min = np.maximum(0, yx - radius[:, np.newaxis] + 0.5).astype(int) p_max = np.minimum( yx + radius[:, np.newaxis] + 0.5, (oshape[0] - 1, oshape[1] - 1) ).astype(int) # orientation histogram hist = np.empty(self.n_bins, dtype=self.float_dtype) avg_kernel = np.full((3,), 1 / 3, dtype=self.float_dtype) for k in range(len(yx)): hist[:] = 0 # use the patch coordinates to get the gradient and then # normalize them r, c = np.meshgrid( np.arange(p_min[k, 0], p_max[k, 0] + 1), np.arange(p_min[k, 1], p_max[k, 1] + 1), indexing='ij', sparse=True, ) gradient_row = gradient_space[o][0][r, c, scales[k]] gradient_col = gradient_space[o][1][r, c, scales[k]] r = r.astype(self.float_dtype, copy=False) c = c.astype(self.float_dtype, copy=False) r -= yx[k, 0] c -= yx[k, 1] # gradient magnitude and angles magnitude = np.sqrt(np.square(gradient_row) + np.square(gradient_col)) theta = np.mod(np.arctan2(gradient_col, gradient_row), 2 * np.pi) # more weight to center values kernel = np.exp( np.divide(r * r + c * c, -2 * (self.lambda_ori * sigma[k]) ** 2) ) # fill the histogram bins = np.floor( (theta / (2 * np.pi) * self.n_bins + 0.5) % self.n_bins ).astype(int) np.add.at(hist, bins, kernel * magnitude) # smooth the histogram and find the maximum hist = np.concatenate((hist[-6:], hist, hist[:6])) for _ in range(6): # number of smoothings hist = np.convolve(hist, avg_kernel, mode='same') hist = hist[6:-6] max_filter = ndi.maximum_filter(hist, [3], mode='wrap') # if an angle is in 80% percent range of the maximum, a # new keypoint is created for it maxima = np.nonzero( np.logical_and( hist >= (self.c_max * np.max(hist)), max_filter == hist ) ) # save the angles for c, m in enumerate(maxima[0]): neigh = np.arange(m - 1, m + 2) % len(hist) # use neighbors to fit a parabola, to get more accurate # result ori = (m + self._fit(hist[neigh]) + 0.5) * 2 * np.pi / self.n_bins if ori > np.pi: ori -= 2 * np.pi if c == 0: orientations[key_count] = ori else: keypoint_indices.append(key_count) keypoint_angles.append(ori) keypoint_octave.append(o) key_count += 1 self.positions = np.concatenate( (positions_oct, positions_oct[keypoint_indices]) ) self.scales = np.concatenate((scales_oct, scales_oct[keypoint_indices])) self.sigmas = np.concatenate((sigmas_oct, sigmas_oct[keypoint_indices])) self.orientations = np.concatenate((orientations, keypoint_angles)) self.octaves = np.concatenate((octaves, keypoint_octave)) # return the gradient_space to reuse it to find the descriptor return gradient_space def _rotate(self, row, col, angle): c = math.cos(angle) s = math.sin(angle) rot_row = c * row + s * col rot_col = -s * row + c * col return rot_row, rot_col def _compute_descriptor(self, gradient_space): """Source: "Anatomy of the SIFT Method" Alg. 12 Calculates the descriptor for every keypoint """ n_key = len(self.scales) self.descriptors = np.empty( (n_key, self.n_hist**2 * self.n_ori), dtype=np.uint8 ) # indices of the histograms hists = np.arange(1, self.n_hist + 1, dtype=self.float_dtype) # indices of the bins bins = np.arange(1, self.n_ori + 1, dtype=self.float_dtype) key_numbers = np.arange(n_key) for o, (gradient, delta) in enumerate(zip(gradient_space, self.deltas)): in_oct = self.octaves == o if not np.any(in_oct): continue positions = self.positions[in_oct] scales = self.scales[in_oct] sigmas = self.sigmas[in_oct] orientations = self.orientations[in_oct] numbers = key_numbers[in_oct] dim = gradient[0].shape[:2] center_pos = positions / delta sigma = sigmas / delta # dimensions of the patch radius = self.lambda_descr * (1 + 1 / self.n_hist) * sigma radius_patch = math.sqrt(2) * radius p_min = np.asarray( np.maximum(0, center_pos - radius_patch[:, np.newaxis] + 0.5), dtype=int ) p_max = np.asarray( np.minimum( center_pos + radius_patch[:, np.newaxis] + 0.5, (dim[0] - 1, dim[1] - 1), ), dtype=int, ) for k in range(len(p_max)): rad_k = float(radius[k]) ori = float(orientations[k]) histograms = np.zeros( (self.n_hist, self.n_hist, self.n_ori), dtype=self.float_dtype ) # the patch r, c = np.meshgrid( np.arange(p_min[k, 0], p_max[k, 0]), np.arange(p_min[k, 1], p_max[k, 1]), indexing='ij', sparse=True, ) # normalized coordinates r_norm = np.subtract(r, center_pos[k, 0], dtype=self.float_dtype) c_norm = np.subtract(c, center_pos[k, 1], dtype=self.float_dtype) r_norm, c_norm = self._rotate(r_norm, c_norm, ori) # select coordinates and gradient values within the patch inside = np.maximum(np.abs(r_norm), np.abs(c_norm)) < rad_k r_norm, c_norm = r_norm[inside], c_norm[inside] r_idx, c_idx = np.nonzero(inside) r = r[r_idx, 0] c = c[0, c_idx] gradient_row = gradient[0][r, c, scales[k]] gradient_col = gradient[1][r, c, scales[k]] # compute the (relative) gradient orientation theta = np.arctan2(gradient_col, gradient_row) - ori lam_sig = self.lambda_descr * float(sigma[k]) # Gaussian weighted kernel magnitude kernel = np.exp((r_norm * r_norm + c_norm * c_norm) / (-2 * lam_sig**2)) magnitude = ( np.sqrt(gradient_row * gradient_row + gradient_col * gradient_col) * kernel ) lam_sig_ratio = 2 * lam_sig / self.n_hist rc_bins = (hists - (1 + self.n_hist) / 2) * lam_sig_ratio rc_bin_spacing = lam_sig_ratio ori_bins = (2 * np.pi * bins) / self.n_ori # distances to the histograms and bins dist_r = np.abs(np.subtract.outer(rc_bins, r_norm)) dist_c = np.abs(np.subtract.outer(rc_bins, c_norm)) # the orientation histograms/bins that get the contribution near_t, near_t_val = _ori_distances(ori_bins, theta) # create the histogram _update_histogram( histograms, near_t, near_t_val, magnitude, dist_r, dist_c, rc_bin_spacing, ) # convert the histograms to a 1d descriptor histograms = histograms.reshape(-1) # saturate the descriptor histograms = np.minimum(histograms, 0.2 * np.linalg.norm(histograms)) # normalize the descriptor descriptor = (512 * histograms) / np.linalg.norm(histograms) # quantize the descriptor descriptor = np.minimum(np.floor(descriptor), 255) self.descriptors[numbers[k], :] = descriptor def _preprocess(self, image): check_nD(image, 2) image = img_as_float(image) self.float_dtype = _supported_float_type(image.dtype) image = image.astype(self.float_dtype, copy=False) self._set_number_of_octaves(image.shape) return image def detect(self, image): """Detect the keypoints. Parameters ---------- image : 2D array Input image. """ image = self._preprocess(image) gaussian_scalespace = self._create_scalespace(image) dog_scalespace = [np.diff(layer, axis=2) for layer in gaussian_scalespace] positions, scales, sigmas, octaves = self._find_localize_evaluate( dog_scalespace, image.shape ) self._compute_orientation( positions, scales, sigmas, octaves, gaussian_scalespace ) self.keypoints = self.positions.round().astype(int) def extract(self, image): """Extract the descriptors for all keypoints in the image. Parameters ---------- image : 2D array Input image. """ image = self._preprocess(image) gaussian_scalespace = self._create_scalespace(image) gradient_space = [np.gradient(octave) for octave in gaussian_scalespace] self._compute_descriptor(gradient_space) def detect_and_extract(self, image): """Detect the keypoints and extract their descriptors. Parameters ---------- image : 2D array Input image. """ image = self._preprocess(image) gaussian_scalespace = self._create_scalespace(image) dog_scalespace = [np.diff(layer, axis=2) for layer in gaussian_scalespace] positions, scales, sigmas, octaves = self._find_localize_evaluate( dog_scalespace, image.shape ) gradient_space = self._compute_orientation( positions, scales, sigmas, octaves, gaussian_scalespace ) self._compute_descriptor(gradient_space) self.keypoints = self.positions.round().astype(int)