From c3177f3a48c499de73fd3843de37d9d9154b4f15 Mon Sep 17 00:00:00 2001 From: Carole Amiot Date: Fri, 1 Mar 2019 08:29:32 +0000 Subject: [PATCH 1/7] Add matching function to library matching call with ctypes --- c/sift/sift4ctypes.cpp | 84 ++++++++++++++++++++++++++++++++++++------ s2plib/sift.py | 78 +++++++++++++++++++++++++++++++++------ 2 files changed, 139 insertions(+), 23 deletions(-) diff --git a/c/sift/sift4ctypes.cpp b/c/sift/sift4ctypes.cpp index bb7359e5..d8aeb216 100644 --- a/c/sift/sift4ctypes.cpp +++ b/c/sift/sift4ctypes.cpp @@ -1,30 +1,32 @@ // Copyright (C) 2019, Julien Michel (CNES) #include +#include #include "Utilities/Parameters.h" #include "LibImages/LibImages.h" #include "LibSift/LibSift.h" +#include "sift_anatomy_20141201/lib_util.h" extern "C"{ /** * This function is meant to be mapped to python using ctypes. - * + * * It computes sifts points of input_buffer which is interpreted as a w x h image. - * Keypoints are returned as a linear buffer of float of size recordSize * nbRecords. - * + * Keypoints are returned as a linear buffer of float of size recordSize * nbRecords. + * * This buffer is the responsibiliy of the caller and should be freed by her. */ - float * sift(const float * input_buffer, const size_t w, const size_t h, - const float thresh_dog, - const unsigned int ss_noct, + float * sift(const float * input_buffer, const size_t w, const size_t h, + const float thresh_dog, + const unsigned int ss_noct, const unsigned int ss_nspo, unsigned int & recordSize, unsigned int & nbRecords) { // Derive Image from buffer Image im(input_buffer, (const size_t) w, (const size_t) h, 1); - + // prepare params object Parameters params; params.setDefaultValues(); @@ -35,7 +37,7 @@ extern "C"{ // run sift Sift sift(params); sift.computeKeyPoints(im); - + // Compute the number of records nbRecords = sift.m_keyPoints->size(); @@ -47,7 +49,7 @@ extern "C"{ descriptorSize = firstPoint->getNbOri() * firstPoint->getNbHist() * firstPoint->getNbHist(); } recordSize = descriptorSize + 4; - + // Allocate output buffer float * out = new float[recordSize*nbRecords]; @@ -64,11 +66,69 @@ extern "C"{ for(unsigned int i = 0; i < descriptorSize;++i) { out[currentIndex+4+i] = (*key)->getPtrDescr()[i]; - } - } + } + } return out; } - /** + + float * matching(const float * k1, + const float * k2, + const unsigned int length_desc, + const unsigned int offset_desc, + const unsigned int nb_sift_k1, + const unsigned int nb_sift_k2, + const float sift_thresh, + const float epi_thresh, + const bool use_relative_method, + unsigned int & nb_match){ + // Structure of k1 and k2 is supposed to be the following one : + // 4 first numbers : pos_y pos_x scale orientation + // length_desc floats representing the descriptors + nb_match = 0; + const float sift_thresh_square = sift_thresh*sift_thresh; + float * matches = new float[nb_sift_k1 * 4]; + // TODO : add the fundamental matrix + // Compute the Euclidian distance for every pairs of points + for (int i = 0; i < nb_sift_k1; i++) { + float distA = INFINITY; + float distB = INFINITY; + int indexA = -1; + + const float * curr_k1_desc = &k1[i*(length_desc+offset_desc) + offset_desc]; + for (int j = 0; j < nb_sift_k2; j++) { + const float * curr_k2_desc = &k2[j*(length_desc+offset_desc) + offset_desc]; + //float dist = euclidean_distance_square(curr_k1_desc, curr_k2_desc, length_desc); + float dist = 0.0; + for (int ll = 0; ll < length_desc; ll++) { + float t = (curr_k1_desc[ll] - curr_k2_desc[ll]); + dist += t*t; + } + // find_the_two_nearest_keys + if (dist < distA) { + distB = distA; + distA = dist; + indexA = j; + } else if (dist < distB) + distB = dist; + } + + float val = distA; + if (use_relative_method){ + val = distA / distB; + } + if (val < sift_thresh_square) { + matches[nb_match*4] = k1[i*(length_desc+offset_desc)]; + matches[nb_match*4+1] = k1[i*(length_desc+offset_desc)+1]; + matches[nb_match*4+2] = k2[indexA*(length_desc+offset_desc)]; + matches[nb_match*4+3] = k2[indexA*(length_desc+offset_desc)+1]; + nb_match += 1; + } + } + return matches; + + } + + /** * This function allows to free the float buffer from python. */ void delete_buffer(float * buffer) diff --git a/s2plib/sift.py b/s2plib/sift.py index 3b321ba9..7f8662a2 100644 --- a/s2plib/sift.py +++ b/s2plib/sift.py @@ -45,7 +45,7 @@ def keypoints_from_nparray(arr, thresh_dog=0.0133, nb_octaves=8, nb_scales=3, of nb_scales (optional): Number of scales offset (optional): offset to apply to sift position in case arr is an extract of a bigger image Returns: - A numpy array of shape (nb_points,132) containing for each row (y,x,scale,orientation, sift_descriptor) + A numpy array of shape (nb_points,132) containing for each row (y,x,scale,orientation, sift_descriptor) """ # retrieve numpy buffer dimensions (h, w) = arr.shape @@ -157,17 +157,27 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, if any, a numpy 2D array containing the list of inliers matches. """ # compute matches + # Open k1 and k2 files containing the two lists of SIFTs + k1_keys_descriptors = list() + with open(k1) as f: + for line in f: + ll = [float(val) for val in line.split(' ')[:-1]] + k1_keys_descriptors.append(ll) + k2_keys_descriptors = list() + with open(k2) as f: + for line in f: + ll = [float(val) for val in line.split(' ')[:-1]] + k2_keys_descriptors.append(ll) + k1_keys_descriptors = np.asarray(k1_keys_descriptors, dtype=np.float32) + k2_keys_descriptors = np.asarray(k2_keys_descriptors, dtype=np.float32) + + matches = keypoints_match_from_nparray(k1_keys_descriptors, k2_keys_descriptors, method, sift_thresh, + epipolar_threshold) + + # Write this to file mfile = common.tmpfile('.txt') - cmd = "matching %s %s -o %s --sift-threshold %f" % ( - k1, k2, mfile, sift_thresh) - if method == 'absolute': - cmd += " --absolute" - if F is not None: - fij = ' '.join(str(x) for x in [F[0, 2], F[1, 2], F[2, 0], - F[2, 1], F[2, 2]]) - cmd = "%s -f \"%s\"" % (cmd, fij) - cmd += " --epipolar-threshold {}".format(epipolar_threshold) - common.run(cmd) + with open(mfile, 'w') as f: + np.savetxt(f, matches, delimiter=' ') matches = np.loadtxt(mfile) if matches.ndim == 2: # filter outliers with ransac @@ -185,6 +195,52 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, return np.loadtxt(mfile) +def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold ): + # load shared library + lib = ctypes.CDLL(sift4ctypes_library) + + # Set expected args and return types + lib.matching.argtypes = (ndpointer(dtype=ctypes.c_float, shape=k1.shape), + ndpointer(dtype=ctypes.c_float, shape=k2.shape), + ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, + ctypes.c_float, ctypes.c_float, + ctypes.c_bool, + ctypes.POINTER(ctypes.c_uint)) + lib.matching.restype = ctypes.POINTER(ctypes.c_float) + + # Get info of descriptor size + nb_sift_k1, descr = k1.shape + sift_offset = 4 + length_descr = descr - sift_offset + + # Transform information of method into boolean + use_relative_method = False + if method == 'relative': + use_relative_method = True + + # Create variables to be updated by function call + nb_matches = ctypes.c_uint() + + # Call sift fonction from sift4ctypes.so + matches_ptr = lib.matching(k1, k2, length_descr, sift_offset, nb_sift_k1, len(k2), + sift_threshold, epi_threshold, + use_relative_method, + ctypes.byref(nb_matches) + ) + + # Transform result into a numpy array + matches = np.asarray([matches_ptr[i] for i in range(0, nb_matches.value * 4)]) + + # Delete results to release memory + lib.delete_buffer.argtypes = (ctypes.POINTER(ctypes.c_float)), + lib.delete_buffer(matches_ptr) + + # Reshape keypoints array + matches = matches.reshape((nb_matches.value, 4)) + + return matches + + def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h): """ Compute a list of SIFT matches between two images on a given roi. From b61a8218b3f5289afbd323f960a95d14e31b4c05 Mon Sep 17 00:00:00 2001 From: Carole Amiot Date: Fri, 1 Mar 2019 10:30:46 +0000 Subject: [PATCH 2/7] Add epipolar verification in keypoints distance computation This module is rewritten in order to get rid of the keypoints structures --- c/sift/Makefile | 6 +-- c/sift/sift4ctypes.cpp | 85 ++++++++++++++++++++++++++++++++++++++---- makefile | 2 - s2plib/sift.py | 13 ++++++- 4 files changed, 92 insertions(+), 14 deletions(-) diff --git a/c/sift/Makefile b/c/sift/Makefile index 89035265..a989eee5 100644 --- a/c/sift/Makefile +++ b/c/sift/Makefile @@ -12,15 +12,15 @@ OBJ = LibImages/LibImages.o \ Utilities/Memory.o \ Utilities/Parameters.o \ Utilities/Time.o \ - Utilities/Utilities.o\ + Utilities/Utilities.o \ sift4ctypes.o \ LIB = libsift4ctypes.so -all: $(LIB) matching +all: $(LIB) $(LIB): $(OBJ) $(CXX) -fPIC -shared -o $@ $^ clean: - rm -f $(LIB) matching $(OBJ) + rm -f $(LIB) $(OBJ) diff --git a/c/sift/sift4ctypes.cpp b/c/sift/sift4ctypes.cpp index d8aeb216..e2263153 100644 --- a/c/sift/sift4ctypes.cpp +++ b/c/sift/sift4ctypes.cpp @@ -6,7 +6,58 @@ #include "Utilities/Parameters.h" #include "LibImages/LibImages.h" #include "LibSift/LibSift.h" -#include "sift_anatomy_20141201/lib_util.h" +#include "linalg.c" + + +// Compute the SQUARE Euclidean distance. +float euclidean_distance_square(const float* x, const float* y, int length) +{ + float d = 0.0; + for (int i = 0; i < length; i++) { + float t = (x[i] - y[i]); + d += t*t; + } + return d; +} + +float distance_epipolar(const float* k1, const float* k2, int length, const float* s1, const float * s2, + const float epi_thresh, const unsigned int offset_desc) +{ + float a = s1[0]; + float b = s1[1]; + float c = s1[2]; + float d = s2[0]; + float e = s2[1]; + float f = s2[2]; + + // Naming follows Ives' convention in which x is the row index + float x1 = k1[0]; + float y1 = k1[1]; + float x2 = k2[0]; + float y2 = k2[1]; + + // rectified x coordinates (in Ives' conventions x is the row index) + float xx1 = b * y1 + a * x1 + c; // scalar product of (b, a, c) and (x1, y1, 1) + float xx2 = e * y2 + d * x2 + f; // scalar product of (e, d, f) and (x2, y2, 1) + + // points satisfy the epipolar constraint when the rectified x are equal + if (fabs(xx1 - xx2) < epi_thresh) + return euclidean_distance_square(&k1[offset_desc], &k2[offset_desc], length); + else + return INFINITY; +} + +float distance(const float* k1, const float* k2, int length, const float* s1, const float * s2, + const float epi_thresh, const unsigned int offset_desc) +{ + // parameters used in distance_epipolar + (void) s1; + (void) s2; + (void) epi_thresh; + + return euclidean_distance_square(&k1[offset_desc], &k2[offset_desc], length); + +} extern "C"{ /** @@ -79,6 +130,8 @@ extern "C"{ const unsigned int nb_sift_k2, const float sift_thresh, const float epi_thresh, + double * fund_mat, + const bool use_fundamental_mat, const bool use_relative_method, unsigned int & nb_match){ // Structure of k1 and k2 is supposed to be the following one : @@ -87,22 +140,40 @@ extern "C"{ nb_match = 0; const float sift_thresh_square = sift_thresh*sift_thresh; float * matches = new float[nb_sift_k1 * 4]; - // TODO : add the fundamental matrix + + float (*keypoints_distance_overloaded)(const float *, + const float *, + int, + const float*, const float*, + const float, + const unsigned int); + + // Use fundamental matrix if given + float s1[3]; float s2[3]; + if (use_fundamental_mat){ + rectifying_similarities_from_affine_fundamental_matrix(s1, s2, fund_mat); + keypoints_distance_overloaded = distance_epipolar; + } + else{ + keypoints_distance_overloaded = distance; + } + // Compute the Euclidian distance for every pairs of points for (int i = 0; i < nb_sift_k1; i++) { float distA = INFINITY; float distB = INFINITY; int indexA = -1; - const float * curr_k1_desc = &k1[i*(length_desc+offset_desc) + offset_desc]; + const float * curr_k1_desc = &k1[i*(length_desc+offset_desc)]; for (int j = 0; j < nb_sift_k2; j++) { - const float * curr_k2_desc = &k2[j*(length_desc+offset_desc) + offset_desc]; - //float dist = euclidean_distance_square(curr_k1_desc, curr_k2_desc, length_desc); - float dist = 0.0; + const float * curr_k2_desc = &k2[j*(length_desc+offset_desc)]; + float dist = keypoints_distance_overloaded(curr_k1_desc, curr_k2_desc, length_desc, + s1, s2, epi_thresh, offset_desc); + /*float dist = 0.0; for (int ll = 0; ll < length_desc; ll++) { float t = (curr_k1_desc[ll] - curr_k2_desc[ll]); dist += t*t; - } + }*/ // find_the_two_nearest_keys if (dist < distA) { distB = distA; diff --git a/makefile b/makefile index ec0ff82f..5ed4886c 100644 --- a/makefile +++ b/makefile @@ -50,7 +50,6 @@ homography: $(BINDIR) sift: $(BINDIR) $(MAKE) -j -C c/sift cp c/sift/libsift4ctypes.so $(LIBDIR) - cp c/sift/matching ${BINDIR} mgm: $(MAKE) -C 3rdparty/mgm #cp 3rdparty/mgm/mgm $(BINDIR) @@ -198,7 +197,6 @@ clean_homography: clean_sift: $(MAKE) -C c/sift clean $(RM) $(LIBDIR)/libsift4ctypes.so - $(RM) $(BINDIR)/matching clean_asift: $(RM) -r $(BINDIR)/build_asift diff --git a/s2plib/sift.py b/s2plib/sift.py index 7f8662a2..35d77775 100644 --- a/s2plib/sift.py +++ b/s2plib/sift.py @@ -172,7 +172,7 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, k2_keys_descriptors = np.asarray(k2_keys_descriptors, dtype=np.float32) matches = keypoints_match_from_nparray(k1_keys_descriptors, k2_keys_descriptors, method, sift_thresh, - epipolar_threshold) + epipolar_threshold, F) # Write this to file mfile = common.tmpfile('.txt') @@ -195,7 +195,7 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, return np.loadtxt(mfile) -def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold ): +def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold, F): # load shared library lib = ctypes.CDLL(sift4ctypes_library) @@ -204,6 +204,7 @@ def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold ) ndpointer(dtype=ctypes.c_float, shape=k2.shape), ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, ctypes.c_float, ctypes.c_float, + ctypes.POINTER(ctypes.c_double), ctypes.c_bool, ctypes.c_bool, ctypes.POINTER(ctypes.c_uint)) lib.matching.restype = ctypes.POINTER(ctypes.c_float) @@ -218,12 +219,20 @@ def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold ) if method == 'relative': use_relative_method = True + # Format fundamental matrix + use_fundamental_matrix = False + coeff_mat = None + if F: + coeff_mat = np.asarray([F[0, 2], F[1, 2], F[2, 0], F[2, 1], F[2, 2]]) + use_fundamental_matrix = True + # Create variables to be updated by function call nb_matches = ctypes.c_uint() # Call sift fonction from sift4ctypes.so matches_ptr = lib.matching(k1, k2, length_descr, sift_offset, nb_sift_k1, len(k2), sift_threshold, epi_threshold, + coeff_mat, use_fundamental_matrix, use_relative_method, ctypes.byref(nb_matches) ) From 1b5ec522524f42eb2c7652b678d6f9bedfa03695 Mon Sep 17 00:00:00 2001 From: Carole Amiot Date: Fri, 1 Mar 2019 14:03:40 +0000 Subject: [PATCH 3/7] [refac] Suppress useless commented code --- c/sift/sift4ctypes.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/c/sift/sift4ctypes.cpp b/c/sift/sift4ctypes.cpp index e2263153..169703f1 100644 --- a/c/sift/sift4ctypes.cpp +++ b/c/sift/sift4ctypes.cpp @@ -169,11 +169,6 @@ extern "C"{ const float * curr_k2_desc = &k2[j*(length_desc+offset_desc)]; float dist = keypoints_distance_overloaded(curr_k1_desc, curr_k2_desc, length_desc, s1, s2, epi_thresh, offset_desc); - /*float dist = 0.0; - for (int ll = 0; ll < length_desc; ll++) { - float t = (curr_k1_desc[ll] - curr_k2_desc[ll]); - dist += t*t; - }*/ // find_the_two_nearest_keys if (dist < distA) { distB = distA; From 0364263adac695fb30d6ff1adc7127c631d5d4f1 Mon Sep 17 00:00:00 2001 From: Carlo de Franchis Date: Tue, 26 Mar 2019 11:07:02 +0100 Subject: [PATCH 4/7] [refactor] remove unused files --- c/sift/matching.c | 108 ------- c/sift/sift_anatomy_20141201/lib_keypoint.c | 313 -------------------- c/sift/sift_anatomy_20141201/lib_keypoint.h | 155 ---------- c/sift/sift_anatomy_20141201/lib_matching.c | 155 ---------- c/sift/sift_anatomy_20141201/lib_matching.h | 31 -- c/sift/sift_anatomy_20141201/lib_util.c | 210 ------------- c/sift/sift_anatomy_20141201/lib_util.h | 71 ----- 7 files changed, 1043 deletions(-) delete mode 100644 c/sift/matching.c delete mode 100644 c/sift/sift_anatomy_20141201/lib_keypoint.c delete mode 100644 c/sift/sift_anatomy_20141201/lib_keypoint.h delete mode 100644 c/sift/sift_anatomy_20141201/lib_matching.c delete mode 100644 c/sift/sift_anatomy_20141201/lib_matching.h delete mode 100644 c/sift/sift_anatomy_20141201/lib_util.c delete mode 100644 c/sift/sift_anatomy_20141201/lib_util.h diff --git a/c/sift/matching.c b/c/sift/matching.c deleted file mode 100644 index 431d8132..00000000 --- a/c/sift/matching.c +++ /dev/null @@ -1,108 +0,0 @@ -/** - * @file matching.c - * @brief [[MAIN]] Matching of keypoints - * - * @li Method 1 : threshold on the ratio of distances to the two nearest keypoints - * @li Method 2 : threshold on the distance to the nearest keypoint - * - * @author Ives Rey-Otero (original) - * @author Carlo de Franchis (modified) - */ - - -#include -#include -#include - -#ifdef USE_TIMING -# include "timing.c" -#endif//USE_TIMING -#include "linalg.c" -#include "pickopt.c" -#include "sift_anatomy_20141201/lib_util.c" -#include "sift_anatomy_20141201/lib_keypoint.c" -#include "sift_anatomy_20141201/lib_matching.c" - - -void print_help(char *v[]) -{ - fprintf(stderr, "usage:\n\t%s file1.txt file2.txt [-o file]" - // 0 1 2 3 4 5 - " [--verbose] [-f \"a b c d e\"]" - " [--absolute] [--sift-threshold t (0.6)]" - " [--epipolar-threshold t (10)]\n", *v); -} - - -int main(int c, char *v[]) -{ - // set default parameters - int n_hist = 4; - int n_ori = 8; - int n_bins = 36; - - // parse arguments - const char *output_file = pick_option(&c, &v, "o", "/dev/stdout"); - bool verbose = (bool) pick_option(&c, &v, "-verbose", NULL); - float sift_thresh = (float) atof(pick_option(&c, &v, "-sift-threshold", "0.6")); - bool meth_flag = ((bool) pick_option(&c, &v, "-absolute", NULL)) == false; - float epi_thresh = (float) atof(pick_option(&c, &v, "-epipolar-threshold", "10")); - - // read the (optional) fundamental matrix - const char *fund_mat_str = pick_option(&c, &v, "f", ""); - double *fund_mat = NULL; - if (strcmp(fund_mat_str, "")) { - // the string is not empty - int n_fund; - fund_mat = alloc_parse_doubles(5, fund_mat_str, &n_fund); - if (n_fund != 5) { - fprintf(stderr, "can't read affine fundamental matrix from \"%s\"\n", fund_mat_str); - return EXIT_FAILURE; - } - } - - // check remaining args - if (c < 3) { - print_help(v); - return EXIT_FAILURE; - } - -#ifdef USE_TIMING - // initialise timer - struct timespec ts; portable_gettime(&ts); -#endif//USE_TIMING - - // Read input keypoint ASCII files - struct sift_keypoints* k1 = sift_malloc_keypoints(); - struct sift_keypoints* k2 = sift_malloc_keypoints(); - sift_read_keypoints(k1, v[1], n_hist, n_ori, n_bins, 1); - sift_read_keypoints(k2, v[2], n_hist, n_ori, n_bins, 1); -#ifdef USE_TIMING - if (verbose) print_elapsed_time(&ts, "read input keypoints:", 35); -#endif//USE_TIMING - - // matching - struct sift_keypoints* out_k1 = sift_malloc_keypoints(); - struct sift_keypoints* out_k2 = sift_malloc_keypoints(); - matching(k1, k2, out_k1, out_k2, sift_thresh, meth_flag, fund_mat, - epi_thresh); -#ifdef USE_TIMING - if (verbose) print_elapsed_time(&ts, "compute matches:", 35); -#endif//USE_TIMING - - // print - fprintf_pairs(output_file, out_k1, out_k2); - -#ifdef USE_TIMING - if (verbose) print_elapsed_time(&ts, "print output:", 35); -#endif//USE_TIMING - fprintf(stderr, "%d matches\n", out_k1->size); - - // cleanup - sift_free_keypoints(k1); - sift_free_keypoints(k2); - free(out_k1); - free(out_k2); - - return EXIT_SUCCESS; -} diff --git a/c/sift/sift_anatomy_20141201/lib_keypoint.c b/c/sift/sift_anatomy_20141201/lib_keypoint.c deleted file mode 100644 index 623eea85..00000000 --- a/c/sift/sift_anatomy_20141201/lib_keypoint.c +++ /dev/null @@ -1,313 +0,0 @@ -/* -IPOL SIFT -Copyright (C) 2014, Ives Rey-Otero, CMLA ENS Cachan - - -Version 20140911 (September 11th, 2014) - -This C ANSI source code is related to the IPOL publication - - [1] "Anatomy of the SIFT Method." - I. Rey Otero and M. Delbracio - Image Processing Online, 2013. - http://www.ipol.im/pub/algo/rd_anatomy_sift/ - -An IPOL demo is available at - http://www.ipol.im/pub/demo/rd_anatomy_sift/ - - - - - -== Patent Warning and License ================================================= - -The SIFT method is patented - - [2] "Method and apparatus for identifying scale invariant features - in an image." - David G. Lowe - Patent number: 6711293 - Filing date: Mar 6, 2000 - Issue date: Mar 23, 2004 - Application number: 09/519,89 - - These source codes are made available for the exclusive aim of serving as - scientific tool to verify the soundness and completeness of the algorithm - description. Compilation, execution and redistribution of this file may - violate patents rights in certain countries. The situation being different - for every country and changing over time, it is your responsibility to - determine which patent rights restrictions apply to you before you compile, - use, modify, or redistribute this file. A patent lawyer is qualified to make - this determination. If and only if they don't conflict with any patent terms, - you can benefit from the following license terms attached to this file. - - -This program is free software: you can use, modify and/or -redistribute it under the terms of the simplified BSD -License. You should have received a copy of this license along -this program. If not, see -. - -*/ -/** - @file sift_keypoint.c - * @brief data structures to store information relative to keypoint - * - * @li struct keypoint keypoint data structure. - * @li struct sift_keypoints list of keypoints with variable length. - * @li print,save, read for lists of keypoints. - * - * @author Ives Rey-Otero - */ - - -#include -#include -#include -#include - -#include "lib_keypoint.h" -#include "lib_util.h" - - -struct keypoint* sift_malloc_keypoint(int n_ori, int n_hist, int n_bins) -{ - struct keypoint* key = xmalloc(sizeof(struct keypoint)); - key->n_ori = n_ori; - key->n_hist = n_hist; - key->n_bins = n_bins; - key->descr = xmalloc(n_ori*n_hist*n_hist*sizeof(float)); - key->orihist = xmalloc(n_bins*sizeof(float)); - // set default value - key->x = 0.; - key->y = 0.; - key->sigma = 0.; - key->theta = 0.; - key->val = 0; - key->o = 0; - key->s = 0; - key->i = 0; - key->j = 0; - for(int n=0;ndescr[n] = 0.; - } - for(int i = 0; i < n_bins; i++){ - key->orihist[i] = 0.; - } - // return pointer - return key; -} - - -static void copy_keypoint(const struct keypoint* kA, struct keypoint* kB) -{ - int l = kB->n_hist * kB->n_hist * kB->n_ori; // length of the feature vector - assert(l == (kA->n_hist * kA->n_hist * kA->n_ori)); - // copy struct - kB->i = kA->i; - kB->j = kA->j; - kB->s = kA->s; - kB->o = kA->o; - kB->x = kA->x; - kB->y = kA->y; - kB->sigma = kA->sigma; - kB->val = kA->val; - kB->theta = kA->theta; - kB->edgeResp = kA->edgeResp; - for(int n = 0; n < l; n++){ - kB->descr[n] = kA->descr[n]; - } - for(int n = 0; n < kB->n_bins; n++){ - kB->orihist[n] = kA->orihist[n]; - } -} - -struct keypoint* sift_malloc_keypoint_from_model_and_copy(const struct keypoint* key) -{ - int n_hist = key->n_hist; /* number of histograms per direction in the feature vector */ - int n_ori = key->n_ori; /* number of bins in each orientation histogram */ - int n_bins = key->n_bins; /* number of bins in the orientation histogram for orientation attribution */ - struct keypoint* keycp = sift_malloc_keypoint(n_ori,n_hist,n_bins); - copy_keypoint(key, keycp); - return keycp; -} - - -struct sift_keypoints* sift_malloc_keypoints() -{ - struct sift_keypoints *keys = xmalloc(sizeof(struct sift_keypoints)); - keys->list = xmalloc(100 * sizeof(struct keypoint*)); - keys->capacity = 100; - keys->size = 0; - return keys; -} - - -static void realloc_sift_keypoints(struct sift_keypoints* keys) -{ - keys->list = (struct keypoint **) xrealloc(keys->list, - 2 * keys->capacity * sizeof(struct keypoint*)); - keys->capacity *= 2; -} - -void sift_add_keypoint_to_list(struct keypoint* key, - struct sift_keypoints* keys) -{ - if (keys->size > keys->capacity - 1) // checks the memory limit - realloc_sift_keypoints(keys); - keys->list[keys->size] = key; - keys->size += 1; -} - - -void sift_free_keypoint(struct keypoint* key) -{ - xfree(key->orihist); - xfree(key->descr); - xfree(key); -} - - -void sift_free_keypoints(struct sift_keypoints* keys) -{ - for (int k = 0; k < keys->size; k++) - sift_free_keypoint(keys->list[k]); - xfree(keys->list); - xfree(keys); -} - -void fprintf_one_keypoint(FILE* f, const struct keypoint* k, int n_descr, int n_bins, int flag) -{ - // coordinates - fprintf(f,"%f %f ", k->y, k->x); - - // scale and orientation - if (flag > -1) - fprintf(f,"%f %f ", k->sigma, k->theta); - - // descriptor - if (flag > 0) - for (int n = 0; n < n_descr; n++) - fprintf(f,"%i ", (int) k->descr[n]); - - // orientation histogram - if (flag > 1) - for (int n = 0; n < n_bins; n++) - fprintf(f,"%f ", k->orihist[n]); -} - - -// dump coordinates, scale, orientation and descriptor to a binary file -void raw_dump_keypoint(FILE* f, const struct keypoint* k, int n_descr) -{ - char buf[4 * sizeof(float) + n_descr * sizeof(uint8_t)]; - float *keypoint = (float *) buf; - uint8_t *descriptor = (uint8_t *) (buf + 4 * sizeof(float)) ; - - // copy the data to a buffer - keypoint[0] = k->y; - keypoint[1] = k->x; - keypoint[2] = k->sigma; - keypoint[3] = k->theta; - for (int j = 0; j < n_descr; j++) - descriptor[j] = (k->descr)[j]; - - // write the buffer to the output file - fwrite(buf, sizeof(char), 4 * sizeof(float) + n_descr, f); -} - - -void fprintf_keypoints(FILE* f, const struct sift_keypoints* keys, int nb, - bool binary, int flag) -{ - // use at most nb keypoints unless nb is zero - int n = nb ? min_2(keys->size, nb) : keys->size; - if (n > 0) { - int n_hist = keys->list[0]->n_hist; - int n_ori = keys->list[0]->n_ori; - int n_descr = n_hist * n_hist * n_ori; - int n_bins = keys->list[0]->n_bins; - for (int i = 0; i < n; i++) { - if (binary) - raw_dump_keypoint(f, keys->list[i], n_descr); - else { - fprintf_one_keypoint(f, keys->list[i], n_descr, n_bins, flag); - fprintf(f,"\n"); - } - } - } -} - - -void sift_save_keypoints(const struct sift_keypoints* keys, const char* name, int flag) -{ - FILE* f = fopen(name,"w"); - if (!f){ - fatal_error("Failed to open %s for writing\n", name); - } - fprintf_keypoints(f, keys, 0, false, flag); - fclose(f); -} - - -void sift_print_keypoints(const struct sift_keypoints* keys, int flag) -{ - fprintf_keypoints(stdout, keys, 0, false, flag); -} - - - - - - -/** @brief read list of oriented keypoints from txt file - * - * @param keys = output list of keypoints - * @param name = input filename - * - * @param n_hist the descriptor has (n_hist*n_hist) weighted coefficients. - * @param n_ori the descriptor has (n_hist*n_hist*n_ori) coefficients - * @param n_bins the size of the orientation histogram - * - */ -void sift_read_keypoints(struct sift_keypoints* keys, - const char* name, - int n_hist, - int n_ori, - int n_bins, - int flag) -{ - size_t buffer_size = 1024 * 1024; // 1MB buffer for long lines. - char* buffer = xmalloc(buffer_size); - FILE* stream = fopen(name,"r"); - if ( !stream) - fatal_error("File \"%s\" not found.", name); - while(fgets(buffer, buffer_size, stream) != NULL){ - int pos = 0; - int read = 0; - struct keypoint* key = sift_malloc_keypoint(n_ori, n_hist, n_bins); - // read coordinates - sscanf(buffer+pos,"%f %f %f %f %n", &key->y - , &key->x - , &key->sigma - , &key->theta - , &read); - pos+=read; - if (flag > 0){ - // read descriptor - for(int i = 0; i < n_hist*n_hist*n_ori; i++){ - sscanf(buffer+pos, "%f %n",&(key->descr[i]),&read); - pos +=read; - } - } - if (flag > 1){ - // read orientation histogram - for(int i=0;iorihist[i]),&read); - pos += read; - } - } - sift_add_keypoint_to_list(key,keys); - } - xfree(buffer); -} diff --git a/c/sift/sift_anatomy_20141201/lib_keypoint.h b/c/sift/sift_anatomy_20141201/lib_keypoint.h deleted file mode 100644 index 9fa627d3..00000000 --- a/c/sift/sift_anatomy_20141201/lib_keypoint.h +++ /dev/null @@ -1,155 +0,0 @@ -/* -IPOL SIFT -Copyright (C) 2014, Ives Rey-Otero, CMLA ENS Cachan - - -Version 20140911 (September 11th, 2014) - -This C ANSI source code is related to the IPOL publication - - [1] "Anatomy of the SIFT Method." - I. Rey Otero and M. Delbracio - Image Processing Online, 2013. - http://www.ipol.im/pub/algo/rd_anatomy_sift/ - -An IPOL demo is available at - http://www.ipol.im/pub/demo/rd_anatomy_sift/ - - - - - -== Patent Warning and License ================================================= - -The SIFT method is patented - - [2] "Method and apparatus for identifying scale invariant features - in an image." - David G. Lowe - Patent number: 6711293 - Filing date: Mar 6, 2000 - Issue date: Mar 23, 2004 - Application number: 09/519,89 - - These source codes are made available for the exclusive aim of serving as - scientific tool to verify the soundness and completeness of the algorithm - description. Compilation, execution and redistribution of this file may - violate patents rights in certain countries. The situation being different - for every country and changing over time, it is your responsibility to - determine which patent rights restrictions apply to you before you compile, - use, modify, or redistribute this file. A patent lawyer is qualified to make - this determination. If and only if they don't conflict with any patent terms, - you can benefit from the following license terms attached to this file. - - -This program is free software: you can use, modify and/or -redistribute it under the terms of the simplified BSD -License. You should have received a copy of this license along -this program. If not, see -. - -*/ -/** - * @file lib_keypoint.h - * @brief data structures to store information relative to keypoint - * - * @li struct keypoint keypoint data structure. - * @li struct sift_keypoints list of keypoints with variable length. - * - * @author Ives Rey-Otero - */ - - - - -#ifndef _LIB_KEYPOINT_H_ -#define _LIB_KEYPOINT_H_ - -#include -#include - - -/** @brief keypoint structure, related to a keypoint - * - * stores SIFT output (keypoint position, scale, orientation, feature vector...) - * stores intermediary results (orientation histogram, summarizes...) - * - * - */ -struct keypoint -{ - float x; // coordinates - float y; - float sigma; // level of blur (it includes the assumed image blur) - float theta; // orientation - - int o; // discrete coordinates - int s; - int i; - int j; - - float val; // normalized operator value (independant of the scalespace sampling) - float edgeResp; // edge response - int n_hist; // number of histograms in each direction - int n_ori; // number of bins per histogram - float* descr; - - int n_bins; // number of bins in the orientation histogram - float* orihist; // gradient orientation histogram -}; - -struct keypoint* sift_malloc_keypoint(int n_ori, int n_hist, int n_bins); - -struct keypoint* sift_malloc_keypoint_from_model_and_copy(const struct keypoint* key); - - -/** @brief list of pointers to keypoint structures (variable size list) - * "vector" of keypoint struct (list of variable size) - */ -struct sift_keypoints{ - int size; // number of elements in the list - int capacity; // current max number of elements in the list (size<=capacity) - struct keypoint** list; // array of pointers to keypoint structure -}; - -/** @brief allocate list of keypoints */ -struct sift_keypoints* sift_malloc_keypoints(); - -/** @brief add element in list of keypoints */ -void sift_add_keypoint_to_list(struct keypoint* key, struct sift_keypoints* keys); - -struct keypoint* sift_add_keypoint_to_listPASC(struct sift_keypoints* keypoints, int n_ori, int n_hist, int n_bins); - -/** @brief free list of keypoints */ -void sift_free_keypoints(struct sift_keypoints *keys); - - -/** @brief Save list of keys on a txt filea (3 formats) - * flag values - * flag >= 0 -> prints coordinates (x,y,sigma,theta) - * flag >= 1 -> prints descriptor - * flag >= 2 -> prints scalespace sample coordinates - * + orientation histogram. - */ -void sift_save_keypoints(const struct sift_keypoints* keys, const char* name, int flag); - -/** @brief Save list of keys on the standard output - */ -void sift_print_keypoints(const struct sift_keypoints* keys, int flag); - - -void fprintf_keypoints(FILE* f, const struct sift_keypoints* keys, int nb, - bool binary, int flag); - -/** @brief Read list of keys from a txt file */ -void sift_read_keypoints(struct sift_keypoints* keys, - const char* name, - int n_hist, - int n_ori, - int n_bins, - int flag); - -// Print a single keypoint -void fprintf_one_keypoint(FILE* f, const struct keypoint *k, int n_descr, int n_bins, int flag); - -#endif // _LIB_KEYPOINT_H_ diff --git a/c/sift/sift_anatomy_20141201/lib_matching.c b/c/sift/sift_anatomy_20141201/lib_matching.c deleted file mode 100644 index b354b31c..00000000 --- a/c/sift/sift_anatomy_20141201/lib_matching.c +++ /dev/null @@ -1,155 +0,0 @@ -/** - * @file sift_matching.c - * @brief data structures to store information relative to a pair of keypoints - * - * @li struct keypointPr : Pair of keypoint data structure. - * @li struct keypointPr_list : List of pairs. - * @li print,save, read for lists of pairs. - * - * @author (original) Ives Rey-Otero - * @author (modified) Carlo de Franchis - * @author (modified) David Youssefi - */ - - -#include -#include -#include - -#include "lib_keypoint.h" -#include "lib_matching.h" -#include "lib_util.h" - -static float keypoints_distance_epipolar(struct keypoint* k1, struct keypoint* k2, - float s1[3], float s2[3], - float epi_thresh, int n) -{ - float a = s1[0]; - float b = s1[1]; - float c = s1[2]; - float d = s2[0]; - float e = s2[1]; - float f = s2[2]; - - // keypoints x, y coordinates - float x1 = k1->x; - float y1 = k1->y; - float x2 = k2->x; - float y2 = k2->y; - - // rectified x coordinates (in Ives' conventions x is the row index) - float xx1 = b * y1 + a * x1 + c; // scalar product of (b, a, c) and (x1, y1, 1) - float xx2 = e * y2 + d * x2 + f; // scalar product of (e, d, f) and (x2, y2, 1) - - // points satisfy the epipolar constraint when the rectified x are equal - if (fabs(xx1 - xx2) < epi_thresh) - return euclidean_distance(k1->descr, k2->descr, n); - else - return INFINITY; -} - - -static float keypoints_distance(struct keypoint* k1, struct keypoint* k2, - float s1[3], float s2[3], float epi_thresh, - int n) -{ - // parameters used in keypoints_distance_epipolar - (void) s1; - (void) s2; - (void) epi_thresh; - - return euclidean_distance(k1->descr, k2->descr, n); -} - - - -void matching(struct sift_keypoints *k1, struct sift_keypoints *k2, - struct sift_keypoints *out_k1, struct sift_keypoints *out_k2, - float sift_thresh, int flag, double fund_mat[5], float epi_thresh) -{ - float s1[3]; float s2[3]; - - float (*keypoints_distance_overloaded)(struct keypoint*, - struct keypoint*, - float*, float*, - float, int); - - if (fund_mat) { - rectifying_similarities_from_affine_fundamental_matrix(s1, s2, fund_mat); - keypoints_distance_overloaded = keypoints_distance_epipolar; - } else - keypoints_distance_overloaded = keypoints_distance; - - if ((k1->size == 0) || (k2->size == 0)) - return; - - int n_hist = k1->list[0]->n_hist; - int n_ori = k1->list[0]->n_ori; - int n = n_hist * n_hist * n_ori; - - for (int i = 0; i < k1->size; i++) { - float distA = INFINITY; - float distB = INFINITY; - int indexA = -1; - - for (int j = 0; j < k2->size; j++) { - float dist = keypoints_distance_overloaded(k1->list[i], k2->list[j], - s1, s2, epi_thresh, n); - // find_the_two_nearest_keys - if (dist < distA) { - distB = distA; - distA = dist; - indexA = j; - } else if (dist < distB) - distB = dist; - } - - float val = (flag == 1 ? distA / distB : distA); - if (val < sift_thresh) { - sift_add_keypoint_to_list(k1->list[i], out_k1); - sift_add_keypoint_to_list(k2->list[indexA], out_k2); - } - } -} - -void fprintf_pairs(const char *filename, const struct sift_keypoints *k1, - const struct sift_keypoints *k2) -{ - FILE *f = fopen(filename, "w"); - for (int i = 0; i < k1->size; i++) { - fprintf_one_keypoint(f, k1->list[i], 0, 0, -1); - fprintf_one_keypoint(f, k2->list[i], 0, 0, -1); - fprintf(f, "\n"); - } - fclose(f); -} - - -void print_pairs(const struct sift_keypoints *k1, - const struct sift_keypoints *k2) -{ - fprintf_pairs("/dev/stdout", k1, k2); -} - - -void save_pairs_extra(const char* name, - const struct sift_keypoints *k1, - const struct sift_keypoints *k2A, - const struct sift_keypoints *k2B) -{ - FILE* f = fopen(name,"w"); - - if (k1->size > 0) { - int n_hist = k1->list[0]->n_hist; - int n_ori = k1->list[0]->n_ori; - int dim = n_hist * n_hist * n_ori; - int n_bins = k1->list[0]->n_bins; - for (int i = 0; i < k1->size; i++) { - fprintf_one_keypoint(f, k1->list[i], dim, n_bins, 2); - fprintf_one_keypoint(f, k2A->list[i], dim, n_bins, 2); - fprintf_one_keypoint(f, k2B->list[i], dim, n_bins, 2); - fprintf(f, "\n"); - } - } - fclose(f); -} diff --git a/c/sift/sift_anatomy_20141201/lib_matching.h b/c/sift/sift_anatomy_20141201/lib_matching.h deleted file mode 100644 index a9046be2..00000000 --- a/c/sift/sift_anatomy_20141201/lib_matching.h +++ /dev/null @@ -1,31 +0,0 @@ -/** - * @file lib_matching.h - * @brief data structures to store information relative to a pair of keypoints - * - * @li struct keypointPr : Pair of keypoint data structure. - * @li struct keypointPr_list : List of pairs. - * @li print,save, read for lists of pairs. - * - * @author Ives Rey-Otero - */ -#ifndef _LIB_MATCHING_H_ -#define _LIB_MATCHING_H_ - -void matching(struct sift_keypoints *k1, struct sift_keypoints *k2, - struct sift_keypoints *out_k1, struct sift_keypoints *out_k2, - float sift_thresh, int flag, - double fund_mat[5], float epi_thresh); - -void fprintf_pairs(const char *filename, - const struct sift_keypoints *k1, - const struct sift_keypoints *k2); - -void print_pairs(const struct sift_keypoints *k1, - const struct sift_keypoints *k2); - -void save_pairs_extra(const char* name, - const struct sift_keypoints *k1, - const struct sift_keypoints *k2A, - const struct sift_keypoints *k2B); - -#endif diff --git a/c/sift/sift_anatomy_20141201/lib_util.c b/c/sift/sift_anatomy_20141201/lib_util.c deleted file mode 100644 index 7b662c0b..00000000 --- a/c/sift/sift_anatomy_20141201/lib_util.c +++ /dev/null @@ -1,210 +0,0 @@ -#include -#include -#include -#include -#include - -// Find the minimum among two int -int min_2(int a, int b) {return b < a ? b : a;} - -// Allocate memory or abord on failure -void* xmalloc(size_t size) -{ - if (size == 0) - fprintf(stderr,"xmalloc: zero size"); - void *p = malloc(size); - if (!p) - { - double sm = size / (0x100000 * 1.0); - fprintf(stderr,"xmalloc: out of memory when requesting " - "%zu bytes (%gMB)",//:\"%s\"", - size, sm);//, strerror(errno)); - } - return p; -} - -// Reallocate memory of abort on failure -void* xrealloc(void* p, size_t size) -{ - void *r = realloc(p, size); - if (!r) fprintf(stderr,"realloc failed"); - return r; -} - -// Free memory allocated by xmalloc or xrealloc. -void xfree(void* p) -{ - if (!p) - fprintf(stderr,"trying to free a null pointer"); - free(p); -} - -// Write a formatted message to standard error and abort the program, for -// example, fatal_error("Failed to open file %s for writing", filename); -void fatal_error(const char* format, ...) -{ - va_list args; - va_start(args, format); - fprintf(stderr, "Fatal error: "); - vfprintf(stderr, format, args); - fprintf(stderr, "\n"); - va_end(args); - exit(1); -} - -// Write formatted message to standard error for debugging. -void debug(const char* format, ...) -{ - va_list args; - va_start(args, format); - fprintf(stderr, "Debug: "); - vfprintf(stderr, format, args); - fprintf(stderr, "\n"); - va_end(args); -} - -// Find the maximum value of an array and store its index in the array. -float find_array_max(const float* array, int length, int *position) -{ - double max = -FLT_MAX; - for(int i = 0; i < length; i++){ - if(array[i] > max){ - *position = i; - max = array[i]; - } - } - return max; -} - -// Find the minimum value of an array and store its index in the array. -float find_array_min(const float* array, int length, int *position) -{ - double min = FLT_MAX; - for(int i = 0; i < length; i++){ - if(array[i] < min){ - min = array[i]; - *position = i; - } - } - return min; -} - -// Find the maximum value of an array. -float array_max(const float* array, int length) -{ - int i; - float max = find_array_max(array, length, &i); - (void)i; - return max; -} - -// Find the minimum value of an array. -float array_min(const float* array, int length) -{ - int i; - float min = find_array_min(array, length, &i); - (void)i; - return min; -} - -// Find the two minimal values of an array. -void find_array_two_min(const float* array, int length, float* minA, float* minB, int* iA, int* iB) -{ - if(array[0] < array[1]){ - *iA = 0; - *iB = 1; - *minA = array[0]; - *minB = array[1]; - } - else{ - *iA = 1; - *iB = 0; - *minA = array[1]; - *minB = array[0]; - } - for(int i = 2; i < length; i++){ - if(array[i] < *minA){ - *minB = *minA; - *minA = array[i]; - *iB = *iA; - *iA = i; - } - else if( array[i] < *minB){ - *minB = array[i]; - *iB = i; - } - } -} - -// Compute the SQUARE Euclidean distance. -float euclidean_distance_square(const float* x, const float* y, int length) -{ - float d = 0.0; - for (int i = 0; i < length; i++) { - float t = (x[i] - y[i]); - d += t*t; - } - return d; -} - -// Compute the Euclidean distance. -float euclidean_distance(const float* x, const float* y, int length) -{ - return sqrt(euclidean_distance_square(x, y, length)); -} - -// L2 norm of a vector -float array_l2norm(const float* array, int length) -{ - float l2norm = 0; - for(int i = 0; i < length; i++){ - l2norm += array[i]*array[i]; - } - l2norm = sqrt(l2norm); - return l2norm; -} - -// Linearly rescale input such that its values are in the range [0, 1]. -void linear_conversion(const float *in, float *out, int l) -{ - float a, b; - float min = array_min(in, l); - float max = array_max(in, l); - - // skip the normalization if max = min - if( max > min + 0.00001){ - a = (1.0 - 0.0) / (max - min); - b = -a * min; - for (int i = 0; i < l; i++) - out[i] = a * in[i] + b; - } - else{ - for (int i = 0; i < l; i++) - out[i] = in[i]; - } -} - -// Compute the x modulus y -float modulus(float x, float y) -{ - float z = x; - int n; - if(z < 0){ - n = (int)((-z)/y)+1; - z += n*y; - } - n = (int)(z/y); - z -= n*y; - return z; -} - -// Multiply the rotation matric R_alpha with [x,y]^T -void apply_rotation(float x, float y, float *rx, float *ry, float alpha) -{ - float c = cos(alpha); - float s = sin(alpha); - float tx = c * x - s * y; - float ty = s * x + c * y; - *rx = tx; - *ry = ty; -} diff --git a/c/sift/sift_anatomy_20141201/lib_util.h b/c/sift/sift_anatomy_20141201/lib_util.h deleted file mode 100644 index 471be0ad..00000000 --- a/c/sift/sift_anatomy_20141201/lib_util.h +++ /dev/null @@ -1,71 +0,0 @@ -#ifndef _LIB_UTIL_H_ -#define _LIB_UTIL_H_ - -#include -#include - -#ifndef M_PI - #define M_PI 3.14159265358979323846 -#endif - -#ifndef M_SQRT2 - #define M_SQRT2 1.41421356237309504880 -#endif - -#ifndef M_LN2 - #define M_LN2 0.693147180559945309417 -#endif - -// Find the minimum among two int -int min_2(int a, int b); - -// Allocate memory or abord on failure -void* xmalloc(size_t size); - -// Reallocate memory of abort on failure -void* xrealloc(void* p, size_t size); - -// Free memory allocated by xmalloc or xrealloc. -void xfree(void* p); - -// Write a formatted message to standard error and abort the program, for -// example, fatal_error("Failed to open file %s for writing", filename); -void fatal_error(const char* format, ...); - -// Write formatted message to standard error for debugging. -void debug(const char* format, ...); - -// Linearly rescale input such that its values are in the range [0, 250]. -void linear_conversion(const float *in, float *out, int l); - -// Find the maximum value of an array. -float array_max(const float* array, int length); - -// Find the minimum value of an array. -float array_min(const float* array, int length); - -// Find the maximum value of an array. -float find_array_max(const float* array, int length, int *position); - -// Find the minimum value of an array. -float find_array_min(const float* array, int length, int *position); - -// Find the two minimal values in an array. -void find_array_two_min(const float* array, int length, float* minA, float* minB, int* iA, int* iB); - -// L2 norm of an array. -float array_l2norm(const float* array, int length); - -// Compute the SQUARE of the euclidean distance -float euclidean_distance_square(const float* a1, const float* a2, int length); - -// Compute the euclidean distance -float euclidean_distance(const float* a1, const float* a2, int length); - -// Compute the x modulus y -float modulus(float x, float y); - -// Multiply the rotation matric R_alpha with [x,y]^T -void apply_rotation(float x, float y, float *rx, float *ry, float alpha); - -#endif // _LIB_UTIL_H_ From e6ecd6939965161176b8665ac64bb07879285d05 Mon Sep 17 00:00:00 2001 From: Carlo de Franchis Date: Tue, 26 Mar 2019 11:26:25 +0100 Subject: [PATCH 5/7] [refactor] pass sift descriptors through memory (possible since merge of PR #193) --- s2p/sift.py | 71 ++++++++++++++++------------------------------- tests/test_s2p.py | 8 +++--- 2 files changed, 28 insertions(+), 51 deletions(-) diff --git a/s2p/sift.py b/s2p/sift.py index b1b166a4..e024c64e 100644 --- a/s2p/sift.py +++ b/s2p/sift.py @@ -23,9 +23,10 @@ # TODO: This is kind of ugly. Cleaner way to do this is to update # LD_LIBRARY_PATH, which we should do once we have a proper config file -sift4ctypes_library = os.path.join(os.path.dirname( - os.path.abspath(__file__)), '../lib/libsift4ctypes.so') -ctypes.CDLL(sift4ctypes_library) +here = os.path.dirname(os.path.abspath(__file__)) +sift4ctypes = os.path.join(os.path.dirname(here), 'lib', 'libsift4ctypes.so') +lib = ctypes.CDLL(sift4ctypes) + # Filter warnings from rasterio reading files wihtout georeferencing warnings.filterwarnings("ignore", category=rio.errors.NotGeoreferencedWarning) @@ -51,9 +52,6 @@ def keypoints_from_nparray(arr, thresh_dog=0.0133, nb_octaves=8, nb_scales=3, of # retrieve numpy buffer dimensions h, w = arr.shape - # load shared library - lib = ctypes.CDLL(sift4ctypes_library) - # Set expected args and return types lib.sift.argtypes = (ndpointer(dtype=ctypes.c_float, shape=(h, w)), ctypes.c_uint, ctypes.c_uint, ctypes.c_float, ctypes.c_uint, ctypes.c_uint, ctypes.POINTER(ctypes.c_uint), ctypes.POINTER(ctypes.c_uint)) @@ -69,7 +67,7 @@ def keypoints_from_nparray(arr, thresh_dog=0.0133, nb_octaves=8, nb_scales=3, of # Transform result into a numpy array keypoints = np.asarray([keypoints_ptr[i] - for i in range(0, nb_points.value*desc_size.value)]) + for i in range(nb_points.value*desc_size.value)]) # Delete results to release memory lib.delete_buffer.argtypes = (ctypes.POINTER(ctypes.c_float)), @@ -94,12 +92,12 @@ def image_keypoints(im, x, y, w, h, max_nb=None, thresh_dog=0.0133, nb_octaves=8 http://www.ipol.im/pub/pre/82/ Args: - im: path to the input image + im (str): path to the input image max_nb (optional): maximal number of keypoints. If more keypoints are detected, those at smallest scales are discarded Returns: - path to the file containing the list of descriptors + numpy array of shape (n, 132) containing, on each row: (y, x, s, o, 128-descriptor) """ # Read file with rasterio with rio.open(im) as ds: @@ -124,10 +122,7 @@ def image_keypoints(im, x, y, w, h, max_nb=None, thresh_dog=0.0133, nb_octaves=8 if max_nb is not None: keypoints = keypoints[:max_nb] - keyfile = common.tmpfile('.txt') - np.savetxt(keyfile, keypoints, delimiter=' ', fmt='%.3f') - - return keyfile + return keypoints def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, @@ -136,7 +131,10 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, Find matches among two lists of sift keypoints. Args: - k1, k2: paths to text files containing the lists of sift descriptors + k1 (array): numpy array of shape (n, 132), where each row represents a + sift keypoint with (y, x, scale, orientation, 128-descriptor) + k2 (array): numpy array of shape (m, 132), where each row represents a + sift keypoint method (optional, default is 'relative'): flag ('relative' or 'absolute') indicating wether to use absolute distance or relative distance @@ -157,21 +155,7 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, if any, a numpy 2D array containing the list of inliers matches. """ # compute matches - # Open k1 and k2 files containing the two lists of SIFTs - k1_keys_descriptors = list() - with open(k1) as f: - for line in f: - ll = [float(val) for val in line.split(' ')[:-1]] - k1_keys_descriptors.append(ll) - k2_keys_descriptors = list() - with open(k2) as f: - for line in f: - ll = [float(val) for val in line.split(' ')[:-1]] - k2_keys_descriptors.append(ll) - k1_keys_descriptors = np.asarray(k1_keys_descriptors, dtype=np.float32) - k2_keys_descriptors = np.asarray(k2_keys_descriptors, dtype=np.float32) - - matches = keypoints_match_from_nparray(k1_keys_descriptors, k2_keys_descriptors, method, sift_thresh, + matches = keypoints_match_from_nparray(k1, k2, method, sift_thresh, epipolar_threshold, F) # Write this to file @@ -179,7 +163,6 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, with open(mfile, 'w') as f: np.savetxt(f, matches, delimiter=' ') - matches = np.loadtxt(mfile) if matches.ndim == 2: # filter outliers with ransac if model == 'fundamental' and len(matches) >= 7: common.run("ransac fmn 1000 .3 7 %s < %s" % (mfile, mfile)) @@ -196,9 +179,8 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold, F): - # load shared library - lib = ctypes.CDLL(sift4ctypes_library) - + """ + """ # Set expected args and return types lib.matching.argtypes = (ndpointer(dtype=ctypes.c_float, shape=k1.shape), ndpointer(dtype=ctypes.c_float, shape=k2.shape), @@ -215,9 +197,7 @@ def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold, length_descr = descr - sift_offset # Transform information of method into boolean - use_relative_method = False - if method == 'relative': - use_relative_method = True + use_relative_method = (method == 'relative') # Format fundamental matrix use_fundamental_matrix = False @@ -230,24 +210,21 @@ def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold, nb_matches = ctypes.c_uint() # Call sift fonction from sift4ctypes.so - matches_ptr = lib.matching(k1, k2, length_descr, sift_offset, nb_sift_k1, len(k2), - sift_threshold, epi_threshold, - coeff_mat, use_fundamental_matrix, - use_relative_method, - ctypes.byref(nb_matches) - ) + matches_ptr = lib.matching(k1.astype('float32'), k2.astype('float32'), + length_descr, sift_offset, len(k1), len(k2), + sift_threshold, epi_threshold, coeff_mat, + use_fundamental_matrix, use_relative_method, + ctypes.byref(nb_matches)) # Transform result into a numpy array - matches = np.asarray([matches_ptr[i] for i in range(0, nb_matches.value * 4)]) + matches = np.asarray([matches_ptr[i] for i in range(nb_matches.value * 4)]) # Delete results to release memory - lib.delete_buffer.argtypes = (ctypes.POINTER(ctypes.c_float)), + lib.delete_buffer.argtypes = ctypes.POINTER(ctypes.c_float), lib.delete_buffer(matches_ptr) # Reshape keypoints array - matches = matches.reshape((nb_matches.value, 4)) - - return matches + return matches.reshape((nb_matches.value, 4)) def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h): diff --git a/tests/test_s2p.py b/tests/test_s2p.py index 5075ea09..713c6138 100644 --- a/tests/test_s2p.py +++ b/tests/test_s2p.py @@ -40,9 +40,7 @@ def unit_gdal_version(): def unit_image_keypoints(): - kpts = s2p.sift.image_keypoints('tests/data/input_triplet/img_02.tif', 100, 100, 200, 200) - - test_kpts = np.loadtxt(kpts) + test_kpts = s2p.sift.image_keypoints('tests/data/input_triplet/img_02.tif', 100, 100, 200, 200) ref_kpts = np.loadtxt('tests/data/expected_output/units/unit_image_keypoints.txt') test_set = set(map(tuple,test_kpts[:,0:2])) @@ -90,7 +88,9 @@ def unit_image_keypoints(): def unit_matching(): - test_matches = s2p.sift.keypoints_match('tests/data/units/sift1.txt','tests/data/units/sift2.txt') + k1 = np.loadtxt('tests/data/units/sift1.txt') + k2 = np.loadtxt('tests/data/units/sift2.txt') + test_matches = s2p.sift.keypoints_match(k1, k2) expected_matches = np.loadtxt('tests/data/expected_output/units/unit_keypoints_match.txt') # Check that numbers of matches are the same From 016d1223aa5ee13ce189022fe6d4314f1ea561bb Mon Sep 17 00:00:00 2001 From: Carlo de Franchis Date: Mon, 8 Apr 2019 15:06:25 +0200 Subject: [PATCH 6/7] [bugfix] replace ctypes.POINTER with ndpointer + handle cases where F is not provided --- s2p/sift.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/s2p/sift.py b/s2p/sift.py index e024c64e..99703988 100644 --- a/s2p/sift.py +++ b/s2p/sift.py @@ -178,16 +178,18 @@ def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, return np.loadtxt(mfile) -def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold, F): +def keypoints_match_from_nparray(k1, k2, method, sift_threshold, + epi_threshold=10, F=None): """ + Wrapper for the sift keypoints matching function of libsift4ctypes.so. """ # Set expected args and return types lib.matching.argtypes = (ndpointer(dtype=ctypes.c_float, shape=k1.shape), ndpointer(dtype=ctypes.c_float, shape=k2.shape), - ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, - ctypes.c_float, ctypes.c_float, - ctypes.POINTER(ctypes.c_double), ctypes.c_bool, - ctypes.c_bool, + ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, + ctypes.c_uint, ctypes.c_float, ctypes.c_float, + ndpointer(dtype=ctypes.c_double, shape=(5,)), + ctypes.c_bool, ctypes.c_bool, ctypes.POINTER(ctypes.c_uint)) lib.matching.restype = ctypes.POINTER(ctypes.c_float) @@ -201,8 +203,8 @@ def keypoints_match_from_nparray(k1, k2, method, sift_threshold, epi_threshold, # Format fundamental matrix use_fundamental_matrix = False - coeff_mat = None - if F: + coeff_mat = np.zeros(5) + if F is not None: coeff_mat = np.asarray([F[0, 2], F[1, 2], F[2, 0], F[2, 1], F[2, 2]]) use_fundamental_matrix = True From 8666c2db619b1bf04c9c1564bfce83c84ce87823 Mon Sep 17 00:00:00 2001 From: Carlo de Franchis Date: Mon, 8 Apr 2019 15:07:27 +0200 Subject: [PATCH 7/7] [bugfix] x y swap --- c/sift/sift4ctypes.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/c/sift/sift4ctypes.cpp b/c/sift/sift4ctypes.cpp index 169703f1..f02d8be7 100644 --- a/c/sift/sift4ctypes.cpp +++ b/c/sift/sift4ctypes.cpp @@ -31,10 +31,10 @@ float distance_epipolar(const float* k1, const float* k2, int length, const floa float f = s2[2]; // Naming follows Ives' convention in which x is the row index - float x1 = k1[0]; - float y1 = k1[1]; - float x2 = k2[0]; - float y2 = k2[1]; + float x1 = k1[1]; + float y1 = k1[0]; + float x2 = k2[1]; + float y2 = k2[0]; // rectified x coordinates (in Ives' conventions x is the row index) float xx1 = b * y1 + a * x1 + c; // scalar product of (b, a, c) and (x1, y1, 1)