====== Searching in large image databases -- Image retrieval ======
In the previous labs, we were searching for correspondences between image pairs. A disadvantage of this method is its quadratic time complexity in the number of images. In this lab, we'll look at image (or object) retrieval – a common problem in many computer vision applications. The goal is to find images corresponding to our query image in an extensive image database. A naive method – trying to match all pairs – is too slow for large databases. Moreover, the set of similar images can be only a tiny fraction of the whole database. Therefore, faster and more efficient methods have been developed...
====== 1. Image Representation with a Set of Visual Words. TF-IDF Weighting ======

One of the fastest methods for image retrieval is based on the so-called bag-of-words model. The basic idea is to represent each image in a database by a set of visual words from a visual vocabulary. 
The visual word can be imagined as a representative of an often occurring image patch. With a visual vocabulary, images can be represented by visual words similarly to documents in a natural language.

In the previous labs, the images were represented by a set of descriptions of normalized patches. Visual words can be chosen from a space of such descriptions. Each description can be understood as a point in this space. We'll choose the visual words by vector quantization - one visual word for each cluster of descriptors. We will implement one method commonly used for vector quantization (clustering): K-means.











===== Vector Quantization - K-means =====

One of the simplest method for vector quantization is k-means. You should already know this algorithm from the [[https://cw.fel.cvut.cz/old/courses/ae4b33rpz/labs/11_k-means/start|Pattern Recognition]] course:
  - Initialization: Choose a defined number $k$ of clusters centers (the red crosses in the image). Random points can be chosen from the set of descriptors, but a more sophisticated method exists. 
  - Assign all points to the nearest center -- points are assigned to<latex>k</latex> clusters.
  - Shift each center into the center of gravity (mean) of the assigned points (compute the center of gravity from points in the same cluster). In the case when the mean has no assigned points a new position has to be chosen for this mean.

Repeat steps 2 and 3 until the assignments in step 2 do not change, or until the global change of point-to-mean distances decreases below some threshold.

{{https://cw.fel.cvut.cz/old/_media/courses/ae4m33mpv/labs/3_indexing/kmeans_en.png?800}}

Implement k-means.
  * Write a function ''nearest(means, data)'', which finds the nearest vector from the set of means (matrix with size KxD, where K is the number of means and D is the dimension of the descriptor space) for each row vector in matrix data (with size NxD, where N is the number of points). The output of this function is indices idxs (vector of length N) of the nearest means and distances dists (vector of length N) to the nearest mean.
  * Write a function [means, err]=''kmeans(K, data)'', which finds K //means// in the vector space //data// and return the sum //err// of distances between //data// points and means to which they are assigned.

We prepared the //data// for you (in the test archive at the bottom of this assignment) to make this assignment a little bit easier. (So the next paragraph is just to put things into a context. You don't have to do these steps in first assignment.) 

/*==== Image indexing ====*/

If we would like to use the functions ''kmeans'' and ''nearest'' to represent image by visual words, we would do so as follows:
  - use ''detect_and_describe'' function to compute the image descriptors using the same setting for all images and choosing parameters such that there is ca. 2000-3000 descriptors per image in average.
  - then we apply ''kmeans'' function to all descriptors to find e.g. 5000 centers of clusters (means). The cluster centers -- vectors in descriptor scape -- and their indices are used as visual words for image representation. The set of clusters (descriptor vectors) is called visual vocabulary.
  - for image indexing, we assign a visual word index to each image descriptor using ''nearest'' function and we store the vector of visual words representing the image. Similarly, we store information about geometry with parameters of frame [x;y;a11;a12;a21;a22].

/*With function ''detect_and_describe'' and the one set of parameters, compute the descriptions of images. Set the parameters in such way that the average number of description will be from 2000 to 3000 for one image. Apply the function ''kmeans'' to all descriptions and find 5000 means - the visual words. Use function ''nearest'' to assign descriptors in each image to visual words.*/
===== Inverted file =====

After the previous step, every image is represented by a vector of visual words. By summing the occurrences of the same words we get a vector representation of the image with counts of visual words (see rows <latex>\alpha</latex>,<latex>\beta</latex>,<latex>\gamma</latex> with occurences of words A,B,C,D on the picture below, left). To be able to search effectively in the image database, we need to estimate the similarity. The standard way is to sum the distances of corresponding descriptions. In the bag of words method the vector quantization is used to approximate the description distance. The distance between descriptions is 0 if they are assigned to the same visual word and infinity otherwise. For images represented as vectors of visual words we define similarity as:
 
<latex>\displaystyle\mathrm{score}(\mathbf{x},\mathbf{y}) = \cos{\varphi} = \frac{\mathbf{x}\cdot \mathbf{y}}{\|\mathbf{x}\|\|\mathbf{y}\|} = \frac{1}{\|\mathbf{x}\|\|\mathbf{y}\|}\sum_{i=1}^D x_i y_i</latex>

where //x// and //y// are vectors of visual words (bags of words). We can take advantage of the fact that a part of the similarity can be computed ahead and normalize the size of vectors. After that, the similarity can be computed as a simple dot product of two vectors.

The visual vocabulary can be very big, often with a few thousands or even millions of words. To simplify the distance estimation between such long vectors, the so called inverted file is used. Inverted file is a structure, which for each visual word (A, B, C, D in the picture) contains a list of images <latex>\alpha</latex>,<latex>\beta</latex>,<latex>\gamma</latex> in which the word appears together with the multiplicity.

{{  https://cw.fel.cvut.cz/old/_media/courses/a4m33mpv/cviceni/3_indexace/inverted.png  }}

We will represent the database as a sparse matrix. A sparse matrix is a structure for the memory-efficient representation of huge matrices where most of the entries are zeros. We will use scipy implementation ([[https://docs.scipy.org/doc/scipy/reference/sparse.html | see scipy documentation for details]]). Our inverted file will be a 2D sparse matrix, where columns are images, and rows are visual words. For instance, the third row of the matrix will contain weights of the third word in all images. We compute the weight of the word from its frequency divided by the length of the vector of visual words (i.e., square root of the sum of squared frequencies).
  * Implement a function ''create_db(vw, num_words)'', which builds the inverted file //DB// in the form described above (sparse matrix NxM, where N=''num_words'' is number of visual words in the vocabulary and M is the number of images. Parameter //vw// is a np.ndarray of objects (object here is another np.ndarray), where each array contains the list of words in the i<sup>th</sup> image. Normalize the resulting vectors, so each column of our database sums to one. (For the sparse matrix, use the scipy function [[https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html#scipy.sparse.csr_matrix|csr_matrix]]).

===== TF-IDF weighting =====

In real images, visual words tend to appear with different frequencies. Similar to words in the text, some words are more frequent than others. The number of visual words in one image is changing depending on the scene's complexity and the detector. To deal with these differences, we have to introduce a suitable weighting of visual words. One weighting scheme, widely used also in text analysis and document retrieval, is called [[http://en.wikipedia.org/wiki/Tf–idf|TF-IDF]] (term frequency–inverse document frequency). Weight of each word $i$ in the image $j$ consists of two values:

$$ \mathrm{tf_{i,j}} = \frac{n_{i,j}}{\sum_k n_{k,j}},$$

where $n_{i,j}$ is the number of occurences of word $i$ in image $j$; and

$$ \mathrm{idf_{i}} = \log \frac{|D|}{|\{d: t_{i} \in d\}|},$$

where $|D|$ is the number of documents and  $|\{d: t_{i} \in d\}|$ is the number of documents containing visual word $i$. For words that are not present in any image, <latex>\mathrm{idf_{i}}=0</latex>. The resulting weight is:

$$ \mathrm{(tf\mbox{-}idf)_{i,j}} = \mathrm{tf_{i,j}} \cdotp \mathrm{idf_{i}}$$

We need two functions
  * a function ''get_idf(vw, num_words)'', which computes the IDF of all words based on the lists of words in all documents. The result will be a ''num_words''$\times$1 matrix.
  * adjust a function ''create_db'' to function ''create_db_tfidf(vw, num_words, idf)'', which instead of word frequencies computes their weights according to TF-IDF (frequencies of each word multiplied by the IDF of the word). Normalize the resulting vectors, so each **column** of our database sums to one. This normalization is the $$ \frac{1}{||x||},$$ part of cosine distance weighting. This way, we save some computation during query time.

===== Image ranking according to similarity with query image =====

Thanks to the inverted file, we can rank images according to their similarity to the query. The query is defined by a bounding box in the image around the object of interest. Take the visual words which are contained inside the rectangle. Compute the length and weights of visual words of the query and multiply it with the corresponding columns of the inverted file, matrix //DB//. The result is a similarity matrix between the query and database images. (You don't have to create visual words selection on your own for now. It is mentioned to show context).
  * write a function ''query(DB, q, idf)'', which computes similarity of the query //q// - list of visual words and images in the inverted file DB. Parameter //idf// are IDF weights of all visual words. The result of this function is an array //score// with similarities ordered in descending order and //img_ids//, indices of images in descending order according to similarities.
===== What you should upload? =====
Upload file image_search.py with the before mentioned implemented methods.

===== Testing =====
To test your code, you can use a jupyter notebook test.ipynb in {{https://gitlab.fel.cvut.cz/mishkdmy/mpv-python-assignment-templates/-/tree/master/assignment_4_5_indexing_template}}. Data can be found {{https://cw.fel.cvut.cz/old/_media/courses/a4m33mpv/cviceni/3_indexace/tfidf_test.zip|here tfidf_test.zip}}. The notebook is also provided with expected results in commentaries. Feel free to use function //get_A_matrix_from_geom// in spatial_verification.py which returns geometry arranged into affine matrix.
====== 2. Fast Spatial Verification. Query Expansion. ======
After an image ranking according to vector similarities, we find the set of images, which have common words with the query. In an ideal case, we will find several views of the same scene. We know that the visual words are representing parts of the scene in the image. Therefore, the coordinates of the visual words in images of the same scene are spatially related. This relation is defined by the mutual position of the scene and the camera at the moment when the images were taken. We can use this relation to verify if the images contain the same scene or object.


===== Spatial Verification =====

The tentative correspondences of visual words will be needed for spatial verification. They will be used to discover a geometric transformation between the coordinates of visual words from the query and words from each relevant image.

After ordering according to tf-idf, we have a set of ''K'' relevant images. These are images that have at least several visual words common with the query. For each pair, query and a similar image, we find tentative correspondences. The same visual words represent features with similar descriptions. Still,  there can be several features with the same visual words in one image and we have to try each pair because we do not have the original description anymore. If there are ''M'' features with the same particular visual word in the query image and ''N'' features with the same word in a similar image, the tentative corresponding pairs (pairs of word indices in the image) will be a Cartesian product of the feature sets. We treat each visual word in the image this way. There may be a lot of features with the same visual word in both images, which produces a large number of tentative pairs (NxM), from which only ''min(N,M)'' can be true correspondences. For this reason we will set two limits on the total number of tentative correspondence in each image pair (''max_tc'') and the maximal number of pairs (''max_MxN''). Moreover, to control the number of tentative correspondences, we will prefer the visual words with smaller product ''MxN''. 
In summary, when computing the tentative correspondences, the following scenario is applied: 1) sort common visual words w.r.t. ''MxN'', 2) keep adding tentative correspondences in ascending order of ''MxN'' until a termination condition is met (either ''MxN'' goes beyond ''max_MxN'' or the number of correspondences goes beyond ''max_tc'').

  * Write a function ''get_tentative_correspondencies(query_vw, vw, relevant_idxs, max_tc, max_MxN)'', which computes tentative correspondences. Correspondences are stored in an np.ndarray of objects (arrays). The //k//-th array will hold an array of //T<sub>k</sub>// tentative correspondences between query and //k//-th relevant image, as an array //Tx2<sub>k</sub>// storing pairs of indices into ''query_vw'' and ''vw'', respectively. The input is a vector of visual words //query_vw//in query image; array of arrays //vw// with arrays of visual words in images in the database; array //relevant// (K, )  of indices of relevant images. Finally, //max_tc// and //max_MxN// defines termination criteria.

We need a geometric model for the verification tentative correspondences. In the same manner, as in the previous task, we can use a homography or an epipolar geometry, but the verification of a huge amount of images can be too time-consuming. It happens quite often that only a small fraction of tentative correspondences are true correspondences. Therefore, it is more feasible to take a simpler model. We can use the fact that our feature points are not only points but affine (resp. similarity) frames. Remember that a frame defines a geometric relation between the canonical frame and an image. If we compose an affine (resp. similarity) transformation from query image into canonical frame and from canonical frame into relevant image from the database, we will get an affine (resp. similarity) transformation of visual word neighborhood from the query image to the image from the database. 

It is clear that geometric transformation between the query and the image from the database is not necessarily affine (resp. similarity). Still from the assumption that we deal with a planar object, we can take it as a local approximation of perspective transformation. With a big enough threshold, we can consider all "roughly" correct points to be consistent. The advantage of this method is its speed. From each tentative correspondence, we get a hypothesis of transformation -- it is not necessary to use sampling. This way, we can verify all hypotheses for each pair query--database image and save the best one. The number of inliers will be the new score, which we will use for the reordering of relevant images. 

  * Write a function ''ransac_affine(q_geom, geometries, correspondencies, relevant, inlier_threshold)'', which estimates the number of inliers and the best hypothesis A of transformation from query to database image. Input is an array //q_geom// (Qx6, where Q is the number of words in the query) of geometries (affine frames [x;y;a11;a12;a21;a22]) of query visual words; array of cells //geometries// with geometries of images from the database; array of cells //corrs// and a list of relevant images //relevant// as in the previous function. The function returns an array //scores// (of size Kx1) with the numbers of inliers and a matrix of transformations //A// (matrix Kx3x3) between query and database image. //Inlier_threshold// is a maximal Euclidean distance of the reprojected points.

  * Join the spatial verification with tf-idf voting into one function ''query_spatial_verification(query_visual_words, query_geometry, bbox, visual_words, geometries, DB, idf, params)'', which orders images from database according to number of inliers. The result will be output of function ''ransac_affine'' joined with output of function ''query'' where query score is higher than parameter //minimum_score//. Only top K (parameter //max_spatial//) images will be processed in spatial verification. Score of relevant images will be added to score from function ''query''.  //img_ids// will be the ordered list of image indices according to the new score. Return also ordered inliers_counts and corresponding transformation matrices. Input parameters will be the same as described above, with one expception -- parameters //query_visual_words// and //q_geom// will contain all visual words from the query and parameter //bbox// (a bounding box in form [xmin, ymin, xmax, ymax]) which will be used for visual word selection. Parameter //params// will contain the union of parameters for functions ''ransac_affine'' and ''get_tentative_correspondencies'', maximal number of images for spatial verification (maximal length of parameter //relevant//) will be set in field //max_spatial//. 



===== Query Expansion =====
**(optional task for 2 bonus points)**

The basic idea of query expansion is to use knowledge of spatial transformation between query and database image, for enriching the set of visual words of the query. In the first step, correct images (results) for the first query have to be chosen. To avoid the expansion of unwanted images, it is necessary to choose only images that are really related to query image. This can be done by choosing only images with high enough score (for instance, at least 10 inliers). For each chosen image, the centers of visual words are projected by transformation A<sup>-1</sup> and words which are projected inside the query bounding box are added into the query. It is a good idea to constrain the total number of visual words in the expanded query. It is desirable to choose visual words, which add information (it is useless to add a visual word from each image with frequency 100, it will only cause a problem during correspondence finding). Therefore, the words are ordered according to ''tf'', and then the first //max_qe_keypoints// words are chosen.
  * add query expansion to the function ''query_spatial_verification''. Query expansion will be used if the parameter //use_query_expansion// is //True//. The other parameters are the same.

===== Testing =====

To test your code, you can use [[https://gitlab.fel.cvut.cz/mishkdmy/mpv-python-assignment-templates/-/blob/9bd445dbf2b9b62bccec260b973d761d86f7cce7/assignment_4_5_indexing_template/spatial_verification.ipynb|a jupyter notebook]] (where you can also see a reference solution). In template you are also provided with utils.py where you can find tools for results visualization. Copy the file with {{ :courses:mpv:labs:3_indexing:mpvdb_haff2.mat |mpvdb_haff2.mat}} database, unpack archive with {{ :courses:mpv:labs:3_indexing:mpvimgs.zip |images}} and put it into data directory in [[https://gitlab.fel.cvut.cz/mishkdmy/mpv-python-assignment-templates/-/tree/master/assignment_4_5_indexing_template|template]].

Upload your image_search.py and spatial_verification.py in a single archive into 06_spatial task. **Do not upload utils.py.**

**For this task, there is no auto evaluation; results will be checked by hand!**

====== 3. Image Retrieval in Big Databases  ======
On the set of 1499 images (225 images by your colleagues from a previous year, 73 images of Pokemons, 506 images from dataset [[https://archive.org/details/ukbench|ukbench]]) and 695 images from the dataset [[http://www.robots.ox.ac.uk/~vgg/data/oxbuildings/|Oxford Buildings]]) we have computed features with detector sshessian with Baumberg iteration (haff2, returns affine covariant points). We have estimated dominant orientation and computed SIFT descriptors. With k-means algorithm (with approximate NN) we have estimated 50000 cluster centers from SIFT descriptions and we have assigned all descriptions. You can find the data in file [[http://cmp.felk.cvut.cz/~sulcmila/teaching/MPV/mpvdb50k_haff2.mat|mpvdb50k_haff2.mat]], where array of cells with the lists of visual words is stored in variable VW, array of cells with geometries in variable GEOM and names of the images in the variable NAMES. Cluster centers are stored in file [[http://cmp.felk.cvut.cz/~sulcmila/teaching/MPV/mpvcx50k_haff2.mat|mpvcx50k_haff2.mat]] in variable CX. Compressed images are in file [[http://cmp.felk.cvut.cz/~sulcmila/teaching/MPV/mpvimgs1499.zip|mpvimgs1499.zip]]. Descriptions to images are in file [[http://cmp.felk.cvut.cz/~sulcmila/teaching/MPV/mpvdesc_haff2.mat|mpvdesc_haff2.mat]], stored as array of cells DESC. Elements of these cells are matrices of type UINT8 of size 128xN<sub>i</sub>, where N<sub>i</sub> is number of points in i<sup>th</sup> image. Compute results on the database of **haff2** detector (it should work better for more difficult perspective transformation).

Choose 7 images:
  * 3 from Oxford Buildings dataset
  * 2 Pokemons: X=(your index in class MPV)-1, your images are those with index from (2*X) to (2*X+1).
  * 2 from dataset ukbench.


On each image chose rectangle with object of query. Choose visual words, which have centers inside the bounding box (including boundary). With function ''query_spatial_verification'' query the database and find results.
Prepare a jupyter notebook where you show your results using
''utils.query_and_visu(q_id, visual_words, geometries, bbox_xyxy, DB, idf, options.copy(), img_names, figsize=(7, 7))'' for visualization (as done in jupyter notebook in previous task).

Then you transform your notebook into file ''results.html'' using
''jupyter nbconvert --to html my_results_on_big_db.ipynb --ExecutePreprocessor.timeout=60  --EmbedImagesPreprocessor.resize=small --output results.html''

Then put results.html into archive and upload file into **07_big_db**. 
<note important>It is important to have file size < 3MB (brute can't show results file higher than 3MB). Using previously shown command, mainly ''figsize=(7, 7)'' in visualisation, shouldn't be a problem.</note>

/*Save their indices into variable //query_id//. Assign decription with function ''nearest''. On each image chose rectangle with object of query and save coordinates (not indices! The top-left corner has coordinates (0,0)) into column of matrix //query_bbx// of size 4x20 of type DOUBLE in order [xmin;ymin;xmax;ymax]. Choose visual words, which have centers inside the bounding box (including boundary). With function ''querysp'' query the database and store the results (score, list of images and transformation A) into structure //query_results// in correct position (1..20) into fields //score//, //img_ids//, //A//. All variables //query_// and variable //opt// save into file ''results.mat'' and submit into upload system as task ''07_results''. */

/*Archive the functions kmeans.m, nearest.m, createdb.m, createdb_tfidf.m, query.m, corrm2m.m, ransacm2m.m and querysp.m together with all used non-standard functions you have created (add the functions from first lab only if they are needed in this task) into one file and submit to as task ''07_spatial''.*/