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/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/sift4ctypes.cpp b/c/sift/sift4ctypes.cpp index bb7359e5..f02d8be7 100644 --- a/c/sift/sift4ctypes.cpp +++ b/c/sift/sift4ctypes.cpp @@ -1,30 +1,83 @@ // Copyright (C) 2019, Julien Michel (CNES) #include +#include #include "Utilities/Parameters.h" #include "LibImages/LibImages.h" #include "LibSift/LibSift.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[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) + 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"{ /** * 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 +88,7 @@ extern "C"{ // run sift Sift sift(params); sift.computeKeyPoints(im); - + // Compute the number of records nbRecords = sift.m_keyPoints->size(); @@ -47,7 +100,7 @@ extern "C"{ descriptorSize = firstPoint->getNbOri() * firstPoint->getNbHist() * firstPoint->getNbHist(); } recordSize = descriptorSize + 4; - + // Allocate output buffer float * out = new float[recordSize*nbRecords]; @@ -64,11 +117,84 @@ 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, + 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 : + // 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]; + + 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)]; + for (int j = 0; j < nb_sift_k2; j++) { + 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); + // 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/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_ diff --git a/makefile b/makefile index 09de734c..94428f90 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) @@ -191,7 +190,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/s2p/sift.py b/s2p/sift.py index f1558fa2..99703988 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) @@ -46,14 +47,11 @@ def keypoints_from_nparray(arr, thresh_dog=0.0133, nb_octaves=8, nb_scales=3, of 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 - # 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,19 +155,14 @@ 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 + matches = keypoints_match_from_nparray(k1, k2, method, sift_thresh, + epipolar_threshold, F) + + # 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 if model == 'fundamental' and len(matches) >= 7: common.run("ransac fmn 1000 .3 7 %s < %s" % (mfile, mfile)) @@ -185,6 +178,57 @@ 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=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, + 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) + + # 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 = (method == 'relative') + + # Format fundamental matrix + use_fundamental_matrix = False + 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 + + # Create variables to be updated by function call + nb_matches = ctypes.c_uint() + + # Call sift fonction from sift4ctypes.so + 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(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 + return matches.reshape((nb_matches.value, 4)) + + 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. diff --git a/tests/test_s2p.py b/tests/test_s2p.py index c9a065d3..caf4b96f 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])) @@ -103,8 +101,8 @@ def unit_matching(): """ Test matching of two SIFT keypoints lists. """ - computed = s2p.sift.keypoints_match('tests/data/units/sift1.txt', - 'tests/data/units/sift2.txt') + computed = s2p.sift.keypoints_match(np.loadtxt('tests/data/units/sift1.txt'), + np.loadtxt('tests/data/units/sift2.txt')) expected = np.loadtxt('tests/data/expected_output/units/unit_keypoints_match.txt') assert_arrays_are_equal(computed, expected)