perf(vulkan): evaluate SIFT and RootSIFT matching

README.md

@@ -27,7 +27,7 @@ Lardon3D ne vise pas simplement "dossier de photos → objet 3D", mais "ensemble
 - **Image View** : vues triées et filtrées pour la TUI
 - **Task** : moteur de tâches avec pause/reprise, annulation et séquences
 - **Task Checkpoint v1** : snapshot durable, fichier atomique et reprise sûre
-- **Project Database v7** : tâches, catalogue, Feature Store multi-type et support intra-image
+- **Project Database v11** : tâches, catalogue, Feature Store, matching et tâches Matcher durables
 - **Task Kind Registry** : identité métier durable et reconstruction runtime explicite
 - **Recovery projet** : reprise automatique sélective et bornée des imports récupérables
 - **Task Queue** : file FIFO avec sélection adaptative et backpressure
@@ -154,7 +154,8 @@ Lardon3D est en développement actif. La persistance des tâches, le catalogue,
 le Feature Store multipasse, le Visual Index ORB, Candidate Pair Generator
 et Matcher v1 sont implémentés. Le runtime Feature + Matcher emploie des tâches
 durables, de petits lots, le Resource Governor interactif et un hot path Vulkan
-ORB exact avec fallback CPU. DAG générique,
+ORB exact avec fallback CPU. La feasibility Vulkan SIFT/RootSIFT a été rejetée ;
+ces deux matchers restent sur OpenCV L2. DAG générique,
 vérification géométrique, SfM et viewer restent des tickets séparés planifiés.
 
 ## Licence

docs/architecture/matcher.md

@@ -77,8 +77,11 @@ Le Match File complet est sérialisé dans un buffer heap borné à 98336 octets
 écrit par un unique `write_exact`, puis synchronisé une fois. Les mesures locales
 restent dans `.opencode/work/current_ticket.md`, pas dans ce contrat canonique.
 À 8192 features, le coût CPU dominant reste l'évaluation exacte des distances
-dans `cv::BFMatcher::knnMatch`. ORB peut remplacer ce seul hot path par Vulkan ;
+dans `cv::BFMatcher::knnMatch`. ORB peut remplacer ce seul hot path par Vulkan.
-SIFT et RootSIFT restent intégralement sur BFMatcher CPU.
+Une feasibility réelle a rejeté Vulkan pour SIFT et RootSIFT : leur accumulation
+flottante ne garantit pas le top-2 OpenCV sur les égalités adversariales et le
+gain n'est présent que sur les grandes paires carrées. Ils restent intégralement
+sur BFMatcher CPU.
 
 `matcher.run` v1 orchestre le Matcher sans connaître son backend interne. La
 tâche persiste uniquement la configuration, l'identité des Feature Sets à

docs/architecture/reconstruction_pipeline.md

@@ -92,7 +92,7 @@ Index, persistance, canonicalisation, idempotence et tâche durable.
 
 | Aspect | Description |
 |--------|-------------|
-| **Matcher v1** | ORB/Hamming exact CPU ou Vulkan, SIFT/RootSIFT L2 CPU, KNN k=2 et Lowe, sans vérification géométrique. |
+| **Matcher v1** | ORB CPU/Vulkan ; SIFT/RootSIFT L2 CPU (Vulkan rejeté) ; k=2 + Lowe ; sans géométrie. |
 | **Persistance** | Match Result NO_MATCH/MATCHED et Match File content-addressed validé. |
 | **Orchestration** | `matcher.run` traite les Candidate Pairs par pages et lots durables de 1/2/4/8. |
 
@@ -228,8 +228,9 @@ Image Catalog (B) ──► Feature Store (C)
 
 ## Statut du pipeline
 
-Import, Image Catalog, Feature Extraction, Feature Store, Visual Index et
+Import, Image Catalog, Feature Extraction, Feature Store, Visual Index,
-Candidate Pair sont **IMPLEMENTED**. Matching, Tracks et SfM sont **PLANNED**.
+Candidate Pair et Matching v1 sont **IMPLEMENTED**. Geometric Verification,
+Tracks et SfM sont **PLANNED**.
 
 Ce document décrit la vision architecturale cible du pipeline de
 reconstruction. Les modules listés ici ne sont pas tous implémentés.

docs/architecture/resource_aware_pipeline.md

@@ -64,6 +64,11 @@ backend ORB top-2 de production possède un contexte lazy réutilisable, 640 Kio
 de buffers bornés et un fallback CPU exact. Le CPU reste le fallback portable
 si Vulkan est absent, incompatible ou désactivé pour la session.
 
+La feasibility SIFT/RootSIFT a borné son prototype Vulkan à 8,125 Mio de
+payload lazy, mais n'a pas franchi la Gate de production. Le Governor ne réserve
+donc aucun slot ni budget GPU pour SIFT/RootSIFT ; leur estimation CPU publiée
+reste inchangée.
+
 ## Profil interactif 8845HS mesuré
 
 - budget CPU Lardon3D : 12 threads logiques, 4 réservés au desktop ;

docs/architecture/vulkan_matcher.md

@@ -1,4 +1,15 @@
-# Backend Vulkan ORB v1
+# Backends Vulkan du Matcher
+
+## Statut
+
+- ORB Vulkan v1 : production, gelé ;
+- SIFT Vulkan : rejeté après feasibility sur Radeon 780M ;
+- RootSIFT Vulkan : rejeté après feasibility sur Radeon 780M.
+
+Le rejet SIFT/RootSIFT ne modifie ni leur contrat CPU, ni leur identité
+persistante. Ils restent sur OpenCV BFMatcher L2.
+
+## Backend Vulkan ORB v1
 
 ## Frontière de correction
 
@@ -66,3 +77,43 @@ lorsque l'outil est présent. Sans cette chaîne,
 le même code compile avec un stub indisponible et tous les Matchers restent CPU.
 L'option Meson `-Dvulkan_orb=disabled` force ce build CPU-only ; `auto` est le
 défaut portable et `enabled` exige explicitement le loader et `glslc`.
+
+## Feasibility SIFT et RootSIFT rejetée
+
+Le prototype, compilé uniquement dans les exécutables de faisabilité, partage
+le contexte ORB et utilise un workgroup par query. Les
+lanes parcourent des tranches fixes des candidats B, accumulent 128 différences
+carrées et la lane 0 fusionne les top-2 locaux dans l'ordre total
+`(distance carrée, train_idx)`. La sortie reste O(A) et aucune matrice A×B
+n'existe. Les buffers F32 maximaux représentent 4 Mio par entrée et 128 Kio de
+sortie, soit 8,125 Mio de payload lazy en plus des 640 Kio ORB.
+
+Sur 24 160 requêtes SIFT et 24 161 requêtes RootSIFT contrôlées, y compris
+des descriptors produits par OpenCV SIFT et les frontières `nextafter` autour
+de Lowe 0,7, aucune divergence Lowe n'a été observée. Ce résultat de corpus
+n'est pas une garantie universelle. FP32 ne reproduit pas les distances bit à
+bit et un corpus d'égalités adversariales démontre
+une divergence d'indices reproductible : OpenCV choisit `(0, 1)` et FP32
+Vulkan `(7, 14)`. FP64 choisit `(0, 128)` et ne restaure donc pas le contrat
+OpenCV. Le recalcul CPU des distances des deux candidats ne peut corriger une
+mauvaise sélection d'indices. Cent répétitions ont produit les mêmes octets GPU
+sur cette Radeon 780M et cette pile RADV ; aucune garantie cross-GPU ou
+cross-driver n'en est déduite.
+
+Le meilleur workgroup mesuré est 64 pour les grandes paires carrées. Après
+fermeture d'une lecture vidéo susceptible d'avoir perturbé la première campagne,
+cinq campagnes contrôlées confirment à 8192² un total médian de 75,243 ms contre
+115,030 ms CPU pour SIFT (1,53×), et 74,744 ms contre 117,484 ms pour RootSIFT
+(1,57×). À 4096², les gains ne sont que 1,14× et 1,11×. Le recalcul CPU par
+`cv::norm` ne reproduit pas BFMatcher bit à bit. Les
+quatre formes asymétriques 256×8192, 1024×4096, 4096×1024 et 8192×256 restent
+toutes plus lentes que le CPU. FP64 mesure environ 271 ms à 8192², soit environ
+0,44× le CPU. Le gain FP32 ne compense pas une identité backend obligatoire,
+une politique de fallback durable distincte et une sélection plus complexe.
+
+Ce rejet concerne Lardon3D v1 sur la cible mesurée, pas Vulkan ou SIFT en
+général. Conclusion : aucune Gate B, aucun sélecteur, aucun fallback SIFT Vulkan, aucune
+modification du fingerprint et aucune migration DB. CPU OpenCV demeure le seul
+backend SIFT/RootSIFT de production. Le shader et ses tests restent des preuves
+reproductibles ; ils ne sont ni exposés par l'API normale, ni embarqués dans les
+binaires de production.

docs/performance/target_hardware.md

@@ -70,6 +70,25 @@ de parité couvrent exactement le top-2 jusqu'à 8192, le Match File complet et
 fallback CPU. Ces nombres décrivent la machine mesurée et ne sont pas un contrat
 portable de latence.
 
+## Feasibility Vulkan SIFT / RootSIFT
+
+Sur le même RADV PHOENIX, `shaderFloat64`, les timestamps compute, un subgroup
+de 64 lanes et les workgroups jusqu'à 1024 invocations sont disponibles. La
+meilleure variante SIFT FP32 utilise 64 lanes. Après fermeture d'une lecture
+vidéo et retour à un PSI avg10 nul, cinq campagnes donnent à 8192² : SIFT
+115,030 ms CPU contre 75,243 ms total potentiel (1,53×), et RootSIFT
+117,484 ms contre 74,744 ms (1,57×). À 4096², les gains totaux tombent à 1,14×
+et 1,11×. Les quatre paires asymétriques testées perdent face au CPU. FP64
+atteint environ 271 ms à 8192².
+
+La campagne numérique trouve une divergence top-2 sur des sommes égales
+adversariales pour FP32 comme FP64, des distances et Match Files différents,
+mais aucune divergence Lowe sur le corpus testé. Le backend devrait donc porter
+une identité scientifique propre et ne pourrait pas employer le fallback CPU
+transparent d'ORB. Cette complexité n'est pas justifiée par le profil de
+performance : SIFT et RootSIFT Vulkan sont rejetés pour Lardon3D v1 sur cette
+cible, OpenCV reste la référence production.
+
 Un run soutenu de 5000 dispatchs mixtes 1024/4096/8192/4096 a traité environ
 1232 paires/s en 4,06 s. Le processus de benchmark complet a culminé à environ
 202 Mio RSS ; les compteurs `pswpin` et `pswpout` sont restés à zéro. Après le

include/lardon3d/orb_vulkan_backend.h

@@ -23,6 +23,16 @@ typedef struct {
   uint32_t second_distance;
 } Lardon3DOrbTop2;
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+typedef struct {
+  uint32_t neighbor_count;
+  uint32_t best_index;
+  float best_squared_distance;
+  uint32_t second_index;
+  float second_squared_distance;
+} Lardon3DSiftTop2;
+#endif
+
 typedef enum {
   LARDON3D_ORB_VULKAN_OK = 0,
   LARDON3D_ORB_VULKAN_UNAVAILABLE,
@@ -53,6 +63,13 @@ Lardon3DOrbVulkanResult lardon3d_orb_vulkan_top2(
     uint32_t feature_count_a, const unsigned char *descriptors_b,
     uint32_t feature_count_b, Lardon3DOrbTop2 *output, size_t output_capacity);
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+Lardon3DOrbVulkanResult lardon3d_sift_vulkan_top2(
+    Lardon3DOrbVulkanBackend *backend, const float *descriptors_a,
+    uint32_t feature_count_a, const float *descriptors_b,
+    uint32_t feature_count_b, Lardon3DSiftTop2 *output, size_t output_capacity);
+#endif
+
 bool lardon3d_orb_vulkan_backend_info(Lardon3DOrbVulkanBackend *backend,
                                       Lardon3DOrbVulkanInfo *info);
 

meson.build

@@ -66,7 +66,42 @@ if matcher_vulkan_enabled
         output: 'orb_top2_spv.h',
         command: [python, files('tools/embed_spirv.py'), '@INPUT@', '@OUTPUT@'],
     )
+    sift_top2_spv = custom_target(
+        'sift-top2-spv',
+        input: 'shaders/sift_top2.comp',
+        output: 'sift_top2.spv',
+        command: [glslc, '-fshader-stage=compute', '@INPUT@', '-o', '@OUTPUT@'],
+    )
+    sift_top2_header = custom_target(
+        'sift-top2-header',
+        input: sift_top2_spv,
+        output: 'sift_top2_spv.h',
+        command: [python, files('tools/embed_spirv.py'), '@INPUT@', '@OUTPUT@'],
+    )
+    sift_top2_fp64_spv = custom_target(
+        'sift-top2-fp64-spv',
+        input: 'shaders/sift_top2.comp',
+        output: 'sift_top2_fp64.spv',
+        command: [
+            glslc,
+            '-DLARDON3D_SIFT_FP64=1',
+            '-fshader-stage=compute',
+            '@INPUT@',
+            '-o',
+            '@OUTPUT@',
+        ],
+    )
+    sift_top2_fp64_header = custom_target(
+        'sift-top2-fp64-header',
+        input: sift_top2_fp64_spv,
+        output: 'sift_top2_fp64_spv.h',
+        command: [python, files('tools/embed_spirv.py'), '@INPUT@', '@OUTPUT@'],
+    )
     matcher_backend_sources += [orb_top2_header]
+    sift_feasibility_sources = matcher_backend_sources + [
+        sift_top2_header,
+        sift_top2_fp64_header,
+    ]
     matcher_backend_dependencies += [vulkan]
 endif
 
@@ -638,6 +673,19 @@ if matcher_vulkan_enabled
     test('orb-vulkan-backend', orb_vulkan_backend_test, timeout: 120,
          env: opencv_test_environment)
 
+    executable(
+        'test-sift-vulkan-feasibility',
+        sources: [
+            'tests/test_sift_vulkan_feasibility.cpp',
+            'src/match_file.c',
+        ] + sift_feasibility_sources,
+        cpp_args: ['-DLARDON3D_SIFT_VULKAN_FEASIBILITY'],
+        build_by_default: false,
+        include_directories: include_directories('include'),
+        dependencies: [threads, openssl, opencv_benchmark]
+            + matcher_backend_dependencies,
+    )
+
     executable(
         'benchmark-orb-vulkan',
         sources: [
@@ -646,6 +694,17 @@ if matcher_vulkan_enabled
         include_directories: include_directories('include'),
         dependencies: [opencv] + matcher_backend_dependencies,
     )
+
+    executable(
+        'benchmark-sift-vulkan',
+        sources: [
+            'tests/benchmark_sift_vulkan.cpp',
+        ] + sift_feasibility_sources,
+        cpp_args: ['-DLARDON3D_SIFT_VULKAN_FEASIBILITY'],
+        build_by_default: false,
+        include_directories: include_directories('include'),
+        dependencies: [opencv] + matcher_backend_dependencies,
+    )
 endif
 
 matcher_e2e_test = executable(

shaders/sift_top2.comp

@@ -0,0 +1,117 @@
+#version 450
+
+#ifdef LARDON3D_SIFT_FP64
+#extension GL_ARB_gpu_shader_fp64 : require
+#define Metric double
+#else
+#define Metric float
+#endif
+
+layout(local_size_x_id = 0) in;
+
+layout(set = 0, binding = 0, std430) readonly buffer DescriptorsA {
+  float descriptors_a[];
+};
+
+layout(set = 0, binding = 1, std430) readonly buffer DescriptorsB {
+  float descriptors_b[];
+};
+
+struct Top2 {
+  uint best_index;
+  float best_squared_distance;
+  uint second_index;
+  float second_squared_distance;
+};
+
+layout(set = 0, binding = 2, std430) writeonly buffer Top2Output {
+  Top2 top2[];
+};
+
+layout(push_constant) uniform Counts {
+  uint count_a;
+  uint count_b;
+} counts;
+
+shared uint local_best_index[256];
+shared Metric local_best_distance[256];
+shared uint local_second_index[256];
+shared Metric local_second_distance[256];
+
+bool before(Metric left_distance, uint left_index,
+            Metric right_distance, uint right_index) {
+  return left_distance < right_distance ||
+         (left_distance == right_distance && left_index < right_index);
+}
+
+void insert_candidate(inout uint best_index, inout Metric best_distance,
+                      inout uint second_index, inout Metric second_distance,
+                      uint candidate_index, Metric candidate_distance) {
+  if (candidate_index == 0xffffffffu) {
+    return;
+  }
+  if (before(candidate_distance, candidate_index, best_distance, best_index)) {
+    second_index = best_index;
+    second_distance = best_distance;
+    best_index = candidate_index;
+    best_distance = candidate_distance;
+    return;
+  }
+  if (candidate_index != best_index &&
+      before(candidate_distance, candidate_index,
+             second_distance, second_index)) {
+    second_index = candidate_index;
+    second_distance = candidate_distance;
+  }
+}
+
+void main() {
+  uint query_index = gl_WorkGroupID.x;
+  uint lane = gl_LocalInvocationID.x;
+  if (query_index >= counts.count_a) {
+    return;
+  }
+
+  uint best_index = 0xffffffffu;
+  Metric best_distance = Metric(1.0) / Metric(0.0);
+  uint second_index = 0xffffffffu;
+  Metric second_distance = Metric(1.0) / Metric(0.0);
+
+  for (uint train_index = lane; train_index < counts.count_b;
+       train_index += gl_WorkGroupSize.x) {
+    Metric distance = Metric(0.0);
+    for (uint component = 0; component < 128; ++component) {
+      Metric difference =
+          Metric(descriptors_a[query_index * 128 + component]) -
+          Metric(descriptors_b[train_index * 128 + component]);
+      distance += difference * difference;
+    }
+    insert_candidate(best_index, best_distance, second_index, second_distance,
+                     train_index, distance);
+  }
+
+  local_best_index[lane] = best_index;
+  local_best_distance[lane] = best_distance;
+  local_second_index[lane] = second_index;
+  local_second_distance[lane] = second_distance;
+  barrier();
+
+  // A single lane merges the bounded local results in lane order. This makes
+  // the reduction and its tie-break independent of subgroup scheduling.
+  if (lane == 0) {
+    best_index = 0xffffffffu;
+    best_distance = Metric(1.0) / Metric(0.0);
+    second_index = 0xffffffffu;
+    second_distance = Metric(1.0) / Metric(0.0);
+    for (uint source = 0; source < gl_WorkGroupSize.x; ++source) {
+      insert_candidate(best_index, best_distance, second_index, second_distance,
+                       local_best_index[source], local_best_distance[source]);
+      insert_candidate(best_index, best_distance, second_index, second_distance,
+                       local_second_index[source], local_second_distance[source]);
+    }
+    top2[query_index].best_index = best_index;
+    top2[query_index].best_squared_distance = float(best_distance);
+    top2[query_index].second_index = second_index;
+    top2[query_index].second_squared_distance = float(second_distance);
+  }
+}

src/orb_vulkan_backend.cpp

@@ -19,6 +19,10 @@
 #include <vulkan/vulkan.h>
 
 #include "orb_top2_spv.h"
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+#include "sift_top2_spv.h"
+#include "sift_top2_fp64_spv.h"
+#endif
 
 namespace {
 
@@ -30,6 +34,10 @@ static_assert(kDescriptorBufferBytes * 2 + kOutputBufferBytes ==
               LARDON3D_ORB_VULKAN_PERMANENT_BUFFER_BYTES);
 constexpr uint64_t kDefaultVulkanWorkThreshold = 768ULL * 768ULL;
 constexpr uint32_t kDefaultWorkgroupSize = 32;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+constexpr VkDeviceSize kSiftDescriptorBufferBytes =
+    static_cast<VkDeviceSize>(LARDON3D_FEATURE_MAX_FEATURES) * 128 * sizeof(float);
+#endif
 
 enum class BackendState {
   kUninitialized,
@@ -53,6 +61,15 @@ struct RawTop2 {
   uint32_t second_distance;
 };
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+struct RawSiftTop2 {
+  uint32_t best_index;
+  float best_squared_distance;
+  uint32_t second_index;
+  float second_squared_distance;
+};
+#endif
+
 static uint64_t elapsed_ns(std::chrono::steady_clock::time_point start) {
   auto elapsed = std::chrono::steady_clock::now() - start;
   return static_cast<uint64_t>(
@@ -106,18 +123,32 @@ struct Lardon3DOrbVulkanBackend {
   bool dedicated_compute_queue = false;
   VkPhysicalDeviceProperties properties{};
   VkPhysicalDeviceMemoryProperties memory_properties{};
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  VkPhysicalDeviceFeatures features{};
+#endif
   VkCommandPool command_pool = VK_NULL_HANDLE;
   VkCommandBuffer command_buffer = VK_NULL_HANDLE;
   VkDescriptorSetLayout descriptor_set_layout = VK_NULL_HANDLE;
   VkPipelineLayout pipeline_layout = VK_NULL_HANDLE;
   VkPipeline pipeline = VK_NULL_HANDLE;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  VkPipeline sift_pipeline = VK_NULL_HANDLE;
+#endif
   VkDescriptorPool descriptor_pool = VK_NULL_HANDLE;
   VkDescriptorSet descriptor_set = VK_NULL_HANDLE;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  VkDescriptorSet sift_descriptor_set = VK_NULL_HANDLE;
+#endif
   VkQueryPool query_pool = VK_NULL_HANDLE;
   bool timestamps_available = false;
   Buffer descriptors_a;
   Buffer descriptors_b;
   Buffer output;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  Buffer sift_descriptors_a;
+  Buffer sift_descriptors_b;
+  Buffer sift_output;
+#endif
   uint32_t workgroup_size = kDefaultWorkgroupSize;
   uint64_t initialization_ns = 0;
   uint64_t last_dispatch_ns = 0;
@@ -152,12 +183,22 @@ static void destroy_vulkan(Lardon3DOrbVulkanBackend *backend) {
   destroy_buffer(backend, &backend->descriptors_a);
   destroy_buffer(backend, &backend->descriptors_b);
   destroy_buffer(backend, &backend->output);
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  destroy_buffer(backend, &backend->sift_descriptors_a);
+  destroy_buffer(backend, &backend->sift_descriptors_b);
+  destroy_buffer(backend, &backend->sift_output);
+#endif
   if (backend->query_pool != VK_NULL_HANDLE) {
     vkDestroyQueryPool(backend->device, backend->query_pool, nullptr);
   }
   if (backend->pipeline != VK_NULL_HANDLE) {
     vkDestroyPipeline(backend->device, backend->pipeline, nullptr);
   }
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  if (backend->sift_pipeline != VK_NULL_HANDLE) {
+    vkDestroyPipeline(backend->device, backend->sift_pipeline, nullptr);
+  }
+#endif
   if (backend->pipeline_layout != VK_NULL_HANDLE) {
     vkDestroyPipelineLayout(backend->device, backend->pipeline_layout, nullptr);
   }
@@ -185,8 +226,14 @@ static void destroy_vulkan(Lardon3DOrbVulkanBackend *backend) {
   backend->descriptor_set_layout = VK_NULL_HANDLE;
   backend->pipeline_layout = VK_NULL_HANDLE;
   backend->pipeline = VK_NULL_HANDLE;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  backend->sift_pipeline = VK_NULL_HANDLE;
+#endif
   backend->descriptor_pool = VK_NULL_HANDLE;
   backend->descriptor_set = VK_NULL_HANDLE;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  backend->sift_descriptor_set = VK_NULL_HANDLE;
+#endif
   backend->query_pool = VK_NULL_HANDLE;
   backend->timestamps_available = false;
 }
@@ -296,6 +343,9 @@ static bool select_device(Lardon3DOrbVulkanBackend *backend) {
   vkGetPhysicalDeviceProperties(backend->physical_device, &backend->properties);
   vkGetPhysicalDeviceMemoryProperties(backend->physical_device,
                                       &backend->memory_properties);
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  vkGetPhysicalDeviceFeatures(backend->physical_device, &backend->features);
+#endif
   backend->workgroup_size = configured_workgroup_size();
   return backend->workgroup_size <=
              backend->properties.limits.maxComputeWorkGroupInvocations &&
@@ -314,6 +364,11 @@ static bool create_device_and_commands(Lardon3DOrbVulkanBackend *backend) {
   device_info.sType = VK_STRUCTURE_TYPE_DEVICE_CREATE_INFO;
   device_info.queueCreateInfoCount = 1;
   device_info.pQueueCreateInfos = &queue_info;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  VkPhysicalDeviceFeatures enabled_features{};
+  enabled_features.shaderFloat64 = backend->features.shaderFloat64;
+  device_info.pEnabledFeatures = &enabled_features;
+#endif
   if (vkCreateDevice(backend->physical_device, &device_info, nullptr,
                      &backend->device) != VK_SUCCESS) {
     return false;
@@ -340,6 +395,41 @@ static bool create_device_and_commands(Lardon3DOrbVulkanBackend *backend) {
   return true;
 }
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+static bool create_shader_pipeline(Lardon3DOrbVulkanBackend *backend,
+                                   const uint32_t *code, size_t code_size,
+                                   VkPipeline *pipeline) {
+  VkShaderModuleCreateInfo shader_info{};
+  shader_info.sType = VK_STRUCTURE_TYPE_SHADER_MODULE_CREATE_INFO;
+  shader_info.codeSize = code_size;
+  shader_info.pCode = code;
+  VkShaderModule shader = VK_NULL_HANDLE;
+  if (vkCreateShaderModule(backend->device, &shader_info, nullptr, &shader) != VK_SUCCESS) {
+    return false;
+  }
+  VkSpecializationMapEntry workgroup_entry{0, 0, sizeof(uint32_t)};
+  VkSpecializationInfo specialization{};
+  specialization.mapEntryCount = 1;
+  specialization.pMapEntries = &workgroup_entry;
+  specialization.dataSize = sizeof(backend->workgroup_size);
+  specialization.pData = &backend->workgroup_size;
+  VkPipelineShaderStageCreateInfo stage{};
+  stage.sType = VK_STRUCTURE_TYPE_PIPELINE_SHADER_STAGE_CREATE_INFO;
+  stage.stage = VK_SHADER_STAGE_COMPUTE_BIT;
+  stage.module = shader;
+  stage.pName = "main";
+  stage.pSpecializationInfo = &specialization;
+  VkComputePipelineCreateInfo pipeline_info{};
+  pipeline_info.sType = VK_STRUCTURE_TYPE_COMPUTE_PIPELINE_CREATE_INFO;
+  pipeline_info.stage = stage;
+  pipeline_info.layout = backend->pipeline_layout;
+  VkResult result = vkCreateComputePipelines(backend->device, VK_NULL_HANDLE, 1,
+                                             &pipeline_info, nullptr, pipeline);
+  vkDestroyShaderModule(backend->device, shader, nullptr);
+  return result == VK_SUCCESS;
+}
+#endif
+
 static bool create_pipeline(Lardon3DOrbVulkanBackend *backend) {
   VkDescriptorSetLayoutBinding bindings[3]{};
   for (uint32_t index = 0; index < 3; ++index) {
@@ -468,10 +558,18 @@ static bool create_buffer(Lardon3DOrbVulkanBackend *backend, VkDeviceSize size,
 static bool create_buffers_and_descriptors(Lardon3DOrbVulkanBackend *backend) {
   VkDescriptorPoolSize pool_size{};
   pool_size.type = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  pool_size.descriptorCount = 6;
+#else
   pool_size.descriptorCount = 3;
+#endif
   VkDescriptorPoolCreateInfo pool_info{};
   pool_info.sType = VK_STRUCTURE_TYPE_DESCRIPTOR_POOL_CREATE_INFO;
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+  pool_info.maxSets = 2;
+#else
   pool_info.maxSets = 1;
+#endif
   pool_info.poolSizeCount = 1;
   pool_info.pPoolSizes = &pool_size;
   if (vkCreateDescriptorPool(backend->device, &pool_info, nullptr,
@@ -528,6 +626,54 @@ static bool create_buffers_and_descriptors(Lardon3DOrbVulkanBackend *backend) {
   return true;
 }
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+static bool create_sift_resources_locked(Lardon3DOrbVulkanBackend *backend) {
+  if (backend->sift_pipeline != VK_NULL_HANDLE) {
+    return true;
+  }
+  const char *fp64 = std::getenv("LARDON3D_VULKAN_SIFT_FP64");
+  bool use_fp64 = fp64 && std::strcmp(fp64, "1") == 0;
+  const uint32_t *code = use_fp64 ? lardon3d_sift_top2_fp64_spv
+                                  : lardon3d_sift_top2_spv;
+  size_t code_size = use_fp64 ? lardon3d_sift_top2_fp64_spv_size
+                              : lardon3d_sift_top2_spv_size;
+  if (!create_shader_pipeline(backend, code, code_size,
+                              &backend->sift_pipeline)) {
+    return false;
+  }
+  VkDescriptorSetAllocateInfo set_info{};
+  set_info.sType = VK_STRUCTURE_TYPE_DESCRIPTOR_SET_ALLOCATE_INFO;
+  set_info.descriptorPool = backend->descriptor_pool;
+  set_info.descriptorSetCount = 1;
+  set_info.pSetLayouts = &backend->descriptor_set_layout;
+  if (vkAllocateDescriptorSets(backend->device, &set_info,
+                               &backend->sift_descriptor_set) != VK_SUCCESS ||
+      !create_buffer(backend, kSiftDescriptorBufferBytes,
+                     &backend->sift_descriptors_a) ||
+      !create_buffer(backend, kSiftDescriptorBufferBytes,
+                     &backend->sift_descriptors_b) ||
+      !create_buffer(backend, kOutputBufferBytes, &backend->sift_output)) {
+    return false;
+  }
+  VkDescriptorBufferInfo buffer_info[3] = {
+      {backend->sift_descriptors_a.buffer, 0, backend->sift_descriptors_a.size},
+      {backend->sift_descriptors_b.buffer, 0, backend->sift_descriptors_b.size},
+      {backend->sift_output.buffer, 0, backend->sift_output.size},
+  };
+  VkWriteDescriptorSet writes[3]{};
+  for (uint32_t index = 0; index < 3; ++index) {
+    writes[index].sType = VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET;
+    writes[index].dstSet = backend->sift_descriptor_set;
+    writes[index].dstBinding = index;
+    writes[index].descriptorCount = 1;
+    writes[index].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
+    writes[index].pBufferInfo = &buffer_info[index];
+  }
+  vkUpdateDescriptorSets(backend->device, 3, writes, 0, nullptr);
+  return true;
+}
+#endif
+
 static bool initialize_locked(Lardon3DOrbVulkanBackend *backend) {
   if (backend->state == BackendState::kAvailable) {
     return true;
@@ -587,6 +733,63 @@ static bool synchronize_host_read(Lardon3DOrbVulkanBackend *backend,
   return vkInvalidateMappedMemoryRanges(backend->device, 1, &range) == VK_SUCCESS;
 }
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+static VkResult record_and_submit_sift(Lardon3DOrbVulkanBackend *backend,
+                                  VkPipeline pipeline,
+                                  VkDescriptorSet descriptor_set,
+                                  uint32_t count_a, uint32_t count_b) {
+  VkResult result = vkResetCommandBuffer(backend->command_buffer, 0);
+  if (result != VK_SUCCESS) {
+    return result;
+  }
+  VkCommandBufferBeginInfo begin_info{};
+  begin_info.sType = VK_STRUCTURE_TYPE_COMMAND_BUFFER_BEGIN_INFO;
+  begin_info.flags = VK_COMMAND_BUFFER_USAGE_ONE_TIME_SUBMIT_BIT;
+  result = vkBeginCommandBuffer(backend->command_buffer, &begin_info);
+  if (result != VK_SUCCESS) {
+    return result;
+  }
+  if (backend->timestamps_available) {
+    vkCmdResetQueryPool(backend->command_buffer, backend->query_pool, 0, 2);
+    vkCmdWriteTimestamp(backend->command_buffer, VK_PIPELINE_STAGE_TOP_OF_PIPE_BIT,
+                        backend->query_pool, 0);
+  }
+  vkCmdBindPipeline(backend->command_buffer, VK_PIPELINE_BIND_POINT_COMPUTE,
+                    pipeline);
+  vkCmdBindDescriptorSets(backend->command_buffer, VK_PIPELINE_BIND_POINT_COMPUTE,
+                          backend->pipeline_layout, 0, 1,
+                          &descriptor_set, 0, nullptr);
+  uint32_t counts[2] = {count_a, count_b};
+  vkCmdPushConstants(backend->command_buffer, backend->pipeline_layout,
+                     VK_SHADER_STAGE_COMPUTE_BIT, 0, sizeof(counts), counts);
+  vkCmdDispatch(backend->command_buffer, count_a, 1, 1);
+  VkMemoryBarrier barrier{};
+  barrier.sType = VK_STRUCTURE_TYPE_MEMORY_BARRIER;
+  barrier.srcAccessMask = VK_ACCESS_SHADER_WRITE_BIT;
+  barrier.dstAccessMask = VK_ACCESS_HOST_READ_BIT;
+  vkCmdPipelineBarrier(backend->command_buffer, VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT,
+                       VK_PIPELINE_STAGE_HOST_BIT, 0, 1, &barrier, 0, nullptr, 0,
+                       nullptr);
+  if (backend->timestamps_available) {
+    vkCmdWriteTimestamp(backend->command_buffer, VK_PIPELINE_STAGE_BOTTOM_OF_PIPE_BIT,
+                        backend->query_pool, 1);
+  }
+  result = vkEndCommandBuffer(backend->command_buffer);
+  if (result != VK_SUCCESS) {
+    return result;
+  }
+  VkSubmitInfo submit_info{};
+  submit_info.sType = VK_STRUCTURE_TYPE_SUBMIT_INFO;
+  submit_info.commandBufferCount = 1;
+  submit_info.pCommandBuffers = &backend->command_buffer;
+  result = vkQueueSubmit(backend->queue, 1, &submit_info, VK_NULL_HANDLE);
+  if (result != VK_SUCCESS) {
+    return result;
+  }
+  return vkQueueWaitIdle(backend->queue);
+}
+#endif
+
 static VkResult record_and_submit(Lardon3DOrbVulkanBackend *backend,
                                   uint32_t count_a, uint32_t count_b) {
   VkResult result = vkResetCommandBuffer(backend->command_buffer, 0);
@@ -751,6 +954,74 @@ extern "C" Lardon3DOrbVulkanResult lardon3d_orb_vulkan_top2(
   return LARDON3D_ORB_VULKAN_OK;
 }
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+extern "C" Lardon3DOrbVulkanResult lardon3d_sift_vulkan_top2(
+    Lardon3DOrbVulkanBackend *backend, const float *descriptors_a,
+    uint32_t feature_count_a, const float *descriptors_b,
+    uint32_t feature_count_b, Lardon3DSiftTop2 *output, size_t output_capacity) {
+  if (!backend || feature_count_a > LARDON3D_FEATURE_MAX_FEATURES ||
+      feature_count_b > LARDON3D_FEATURE_MAX_FEATURES ||
+      (feature_count_a > 0 && (!descriptors_a || !output ||
+                              output_capacity < feature_count_a)) ||
+      (feature_count_b > 0 && !descriptors_b)) {
+    return LARDON3D_ORB_VULKAN_INVALID_ARGUMENT;
+  }
+  if (feature_count_a == 0) {
+    return LARDON3D_ORB_VULKAN_OK;
+  }
+  if (feature_count_b == 0) {
+    for (uint32_t index = 0; index < feature_count_a; ++index) {
+      output[index] = Lardon3DSiftTop2{};
+    }
+    return LARDON3D_ORB_VULKAN_OK;
+  }
+
+  std::lock_guard<std::mutex> lock(backend->mutex);
+  if (!initialize_locked(backend)) {
+    return LARDON3D_ORB_VULKAN_UNAVAILABLE;
+  }
+  if (!create_sift_resources_locked(backend)) {
+    return fail_session_locked(backend);
+  }
+  VkDeviceSize bytes_a = static_cast<VkDeviceSize>(feature_count_a) * 128 *
+                         sizeof(float);
+  VkDeviceSize bytes_b = static_cast<VkDeviceSize>(feature_count_b) * 128 *
+                         sizeof(float);
+  std::memcpy(backend->sift_descriptors_a.mapping, descriptors_a,
+              static_cast<size_t>(bytes_a));
+  std::memcpy(backend->sift_descriptors_b.mapping, descriptors_b,
+              static_cast<size_t>(bytes_b));
+  if (!synchronize_host_write(backend, backend->sift_descriptors_a, bytes_a) ||
+      !synchronize_host_write(backend, backend->sift_descriptors_b, bytes_b)) {
+    return fail_session_locked(backend);
+  }
+
+  auto start = std::chrono::steady_clock::now();
+  VkResult dispatch_result = record_and_submit_sift(
+      backend, backend->sift_pipeline, backend->sift_descriptor_set,
+      feature_count_a, feature_count_b);
+  backend->last_dispatch_ns = elapsed_ns(start);
+  if (dispatch_result != VK_SUCCESS ||
+      !synchronize_host_read(backend, backend->sift_output)) {
+    return fail_session_locked(backend);
+  }
+  read_gpu_time(backend);
+
+  const RawSiftTop2 *raw =
+      static_cast<const RawSiftTop2 *>(backend->sift_output.mapping);
+  uint32_t neighbors = std::min(feature_count_b, 2U);
+  for (uint32_t index = 0; index < feature_count_a; ++index) {
+    output[index].neighbor_count = neighbors;
+    output[index].best_index = raw[index].best_index;
+    output[index].best_squared_distance = raw[index].best_squared_distance;
+    output[index].second_index = neighbors == 2 ? raw[index].second_index : 0;
+    output[index].second_squared_distance =
+        neighbors == 2 ? raw[index].second_squared_distance : 0.0F;
+  }
+  return LARDON3D_ORB_VULKAN_OK;
+}
+#endif
+
 extern "C" bool lardon3d_orb_vulkan_backend_info(
     Lardon3DOrbVulkanBackend *backend, Lardon3DOrbVulkanInfo *info) {
   if (!backend || !info) {
@@ -807,6 +1078,24 @@ extern "C" Lardon3DOrbVulkanResult lardon3d_orb_vulkan_top2(
   return LARDON3D_ORB_VULKAN_UNAVAILABLE;
 }
 
+#ifdef LARDON3D_SIFT_VULKAN_FEASIBILITY
+extern "C" Lardon3DOrbVulkanResult lardon3d_sift_vulkan_top2(
+    Lardon3DOrbVulkanBackend *backend, const float *, uint32_t feature_count_a,
+    const float *, uint32_t feature_count_b, Lardon3DSiftTop2 *output,
+    size_t output_capacity) {
+  if (!backend || (feature_count_a > 0 && (!output || output_capacity < feature_count_a))) {
+    return LARDON3D_ORB_VULKAN_INVALID_ARGUMENT;
+  }
+  if (feature_count_a == 0 || feature_count_b == 0) {
+    for (uint32_t index = 0; index < feature_count_a; ++index) {
+      output[index] = Lardon3DSiftTop2{};
+    }
+    return LARDON3D_ORB_VULKAN_OK;
+  }
+  return LARDON3D_ORB_VULKAN_UNAVAILABLE;
+}
+#endif
+
 extern "C" bool lardon3d_orb_vulkan_backend_info(
     Lardon3DOrbVulkanBackend *backend, Lardon3DOrbVulkanInfo *info) {
   if (!backend || !info) {

tests/benchmark_sift_vulkan.cpp

@@ -0,0 +1,188 @@
+#include <lardon3d/orb_vulkan_backend.h>
+
+#include <algorithm>
+#include <chrono>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+
+#include <opencv2/core.hpp>
+#include <opencv2/features2d.hpp>
+
+static uint32_t random_u32(uint32_t *state) {
+  uint32_t value = *state;
+  value ^= value << 13;
+  value ^= value >> 17;
+  value ^= value << 5;
+  *state = value;
+  return value;
+}
+
+static std::vector<float> make_descriptors(uint32_t count, bool rootsift,
+                                           uint32_t seed) {
+  std::vector<float> descriptors(static_cast<size_t>(count) * 128);
+  for (uint32_t row = 0; row < count; ++row) {
+    float sum = 0.0F;
+    for (uint32_t component = 0; component < 128; ++component) {
+      float value = static_cast<float>(random_u32(&seed) & 0xffffU) / 65535.0F;
+      descriptors[static_cast<size_t>(row) * 128 + component] = value;
+      sum += value;
+    }
+    if (rootsift) {
+      for (uint32_t component = 0; component < 128; ++component) {
+        float &value = descriptors[static_cast<size_t>(row) * 128 + component];
+        value = std::sqrt(value / sum);
+      }
+    }
+  }
+  return descriptors;
+}
+
+static double median(std::vector<double> values) {
+  std::sort(values.begin(), values.end());
+  return values[values.size() / 2];
+}
+
+static double elapsed_ms(std::chrono::steady_clock::time_point start) {
+  return std::chrono::duration<double, std::milli>(
+             std::chrono::steady_clock::now() - start)
+      .count();
+}
+
+static bool benchmark(uint32_t count_a, uint32_t count_b, bool rootsift) {
+  std::vector<float> a = make_descriptors(count_a, rootsift, 0x12345678U);
+  std::vector<float> b = make_descriptors(count_b, rootsift, 0x87654321U);
+  cv::Mat matrix_a(static_cast<int>(count_a), 128, CV_32FC1, a.data());
+  cv::Mat matrix_b(static_cast<int>(count_b), 128, CV_32FC1, b.data());
+  cv::BFMatcher matcher(cv::NORM_L2, false);
+  std::vector<std::vector<cv::DMatch>> cpu_output;
+  matcher.knnMatch(matrix_a, matrix_b, cpu_output, 2);
+
+  std::vector<double> cpu_samples;
+  for (int repetition = 0; repetition < 5; ++repetition) {
+    auto start = std::chrono::steady_clock::now();
+    matcher.knnMatch(matrix_a, matrix_b, cpu_output, 2);
+    cpu_samples.push_back(elapsed_ms(start));
+  }
+
+  std::vector<Lardon3DSiftTop2> gpu_output(count_a);
+  Lardon3DOrbVulkanBackend *backend = lardon3d_orb_vulkan_backend_create();
+  if (!backend) {
+    return false;
+  }
+  auto cold_start = std::chrono::steady_clock::now();
+  if (lardon3d_sift_vulkan_top2(backend, a.data(), count_a, b.data(), count_b,
+                                gpu_output.data(), gpu_output.size()) !=
+      LARDON3D_ORB_VULKAN_OK) {
+    lardon3d_orb_vulkan_backend_destroy(backend);
+    return false;
+  }
+  double cold_ms = elapsed_ms(cold_start);
+
+  Lardon3DOrbVulkanBackend *shared_backend =
+      lardon3d_orb_vulkan_backend_create();
+  if (!shared_backend) {
+    lardon3d_orb_vulkan_backend_destroy(backend);
+    return false;
+  }
+  unsigned char orb_descriptor[32]{};
+  Lardon3DOrbTop2 orb_output{};
+  if (lardon3d_orb_vulkan_top2(shared_backend, orb_descriptor, 1,
+                               orb_descriptor, 1, &orb_output, 1) !=
+      LARDON3D_ORB_VULKAN_OK) {
+    lardon3d_orb_vulkan_backend_destroy(shared_backend);
+    lardon3d_orb_vulkan_backend_destroy(backend);
+    return false;
+  }
+  auto after_orb_start = std::chrono::steady_clock::now();
+  if (lardon3d_sift_vulkan_top2(shared_backend, a.data(), count_a, b.data(),
+                                count_b, gpu_output.data(), gpu_output.size()) !=
+      LARDON3D_ORB_VULKAN_OK) {
+    lardon3d_orb_vulkan_backend_destroy(shared_backend);
+    lardon3d_orb_vulkan_backend_destroy(backend);
+    return false;
+  }
+  double after_orb_ms = elapsed_ms(after_orb_start);
+  lardon3d_orb_vulkan_backend_destroy(shared_backend);
+
+  std::vector<double> warm_samples;
+  std::vector<double> gpu_samples;
+  for (int repetition = 0; repetition < 5; ++repetition) {
+    auto start = std::chrono::steady_clock::now();
+    if (lardon3d_sift_vulkan_top2(backend, a.data(), count_a, b.data(), count_b,
+                                  gpu_output.data(), gpu_output.size()) !=
+        LARDON3D_ORB_VULKAN_OK) {
+      lardon3d_orb_vulkan_backend_destroy(backend);
+      return false;
+    }
+    warm_samples.push_back(elapsed_ms(start));
+    Lardon3DOrbVulkanInfo info;
+    if (!lardon3d_orb_vulkan_backend_info(backend, &info)) {
+      lardon3d_orb_vulkan_backend_destroy(backend);
+      return false;
+    }
+    gpu_samples.push_back(static_cast<double>(info.gpu_ns) / 1000000.0);
+  }
+  double cpu_ms = median(cpu_samples);
+  double warm_ms = median(warm_samples);
+  auto finalization_start = std::chrono::steady_clock::now();
+  double finalization_checksum = 0.0;
+  uint32_t exact_final_distances = 0;
+  uint32_t compared_final_distances = 0;
+  for (uint32_t query = 0; query < count_a; ++query) {
+    const Lardon3DSiftTop2 &top2 = gpu_output[query];
+    const uint32_t indices[2] = {top2.best_index, top2.second_index};
+    for (uint32_t neighbor = 0; neighbor < top2.neighbor_count; ++neighbor) {
+      double distance = cv::norm(matrix_a.row(static_cast<int>(query)),
+                                 matrix_b.row(static_cast<int>(indices[neighbor])),
+                                 cv::NORM_L2);
+      finalization_checksum += distance;
+      if (query < cpu_output.size() && neighbor < cpu_output[query].size() &&
+          indices[neighbor] ==
+              static_cast<uint32_t>(cpu_output[query][neighbor].trainIdx)) {
+        ++compared_final_distances;
+        if (static_cast<float>(distance) == cpu_output[query][neighbor].distance) {
+          ++exact_final_distances;
+        }
+      }
+    }
+  }
+  double finalization_ms = elapsed_ms(finalization_start);
+  bool final_exact = compared_final_distances == count_a * 2U &&
+                     exact_final_distances == compared_final_distances;
+  if (!std::isfinite(finalization_checksum)) {
+    lardon3d_orb_vulkan_backend_destroy(backend);
+    return false;
+  }
+  std::printf("%s,%u,%u,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.2f,%s\n",
+              rootsift ? "rootsift" : "sift", count_a, count_b, cpu_ms,
+              cold_ms, after_orb_ms, warm_ms, median(gpu_samples),
+              finalization_ms, cpu_ms / (warm_ms + finalization_ms),
+              final_exact ? "yes" : "no");
+  lardon3d_orb_vulkan_backend_destroy(backend);
+  return true;
+}
+
+int main() {
+  cv::setNumThreads(12);
+  std::printf("kind,count_a,count_b,cpu_ms,cold_ms,after_orb_ms,warm_ms,gpu_ms,"
+              "finalization_ms,gain_with_finalization,final_exact\n");
+  const uint32_t sizes[] = {256, 1024, 4096, 8192};
+  for (bool rootsift : {false, true}) {
+    for (uint32_t size : sizes) {
+      if (!benchmark(size, size, rootsift)) {
+        return EXIT_FAILURE;
+      }
+    }
+  }
+  const uint32_t asymmetric[][2] = {
+      {256, 8192}, {1024, 4096}, {4096, 1024}, {8192, 256},
+  };
+  for (const auto &counts : asymmetric) {
+    if (!benchmark(counts[0], counts[1], false)) {
+      return EXIT_FAILURE;
+    }
+  }
+  return EXIT_SUCCESS;
+}

tests/test_sift_vulkan_feasibility.cpp

@@ -0,0 +1,346 @@
+#include <lardon3d/orb_vulkan_backend.h>
+
+extern "C" {
+#include <lardon3d/match_file.h>
+}
+
+#include <array>
+#include <cmath>
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <limits>
+#include <vector>
+
+#include <fcntl.h>
+#include <openssl/evp.h>
+#include <unistd.h>
+
+#include <opencv2/core.hpp>
+#include <opencv2/features2d.hpp>
+#include <opencv2/imgproc.hpp>
+
+struct Comparison {
+  uint64_t queries = 0;
+  uint64_t index_divergences = 0;
+  uint64_t distance_bit_divergences = 0;
+  uint64_t lowe_divergences = 0;
+};
+
+static uint32_t random_u32(uint32_t *state) {
+  uint32_t value = *state;
+  value ^= value << 13;
+  value ^= value >> 17;
+  value ^= value << 5;
+  *state = value;
+  return value;
+}
+
+static bool fp64_requested() {
+  const char *value = std::getenv("LARDON3D_VULKAN_SIFT_FP64");
+  return value && std::strcmp(value, "1") == 0;
+}
+
+static std::vector<float> make_descriptors(uint32_t count, bool rootsift,
+                                           uint32_t seed) {
+  std::vector<float> descriptors(static_cast<size_t>(count) * 128);
+  for (uint32_t row = 0; row < count; ++row) {
+    float sum = 0.0F;
+    for (uint32_t component = 0; component < 128; ++component) {
+      float value = static_cast<float>(random_u32(&seed) & 0xffffU) / 65535.0F;
+      descriptors[static_cast<size_t>(row) * 128 + component] = value;
+      sum += value;
+    }
+    if (rootsift) {
+      for (uint32_t component = 0; component < 128; ++component) {
+        float &value = descriptors[static_cast<size_t>(row) * 128 + component];
+        value = std::sqrt(value / sum);
+      }
+    }
+  }
+  return descriptors;
+}
+
+static bool compare_case(Lardon3DOrbVulkanBackend *backend,
+                         const std::vector<float> &a,
+                         const std::vector<float> &b,
+                         Comparison *comparison) {
+  uint32_t count_a = static_cast<uint32_t>(a.size() / 128);
+  uint32_t count_b = static_cast<uint32_t>(b.size() / 128);
+  cv::Mat matrix_a(static_cast<int>(count_a), 128, CV_32FC1,
+                   const_cast<float *>(a.data()));
+  cv::Mat matrix_b(static_cast<int>(count_b), 128, CV_32FC1,
+                   const_cast<float *>(b.data()));
+  cv::BFMatcher matcher(cv::NORM_L2, false);
+  std::vector<std::vector<cv::DMatch>> expected;
+  matcher.knnMatch(matrix_a, matrix_b, expected, 2);
+
+  std::vector<Lardon3DSiftTop2> actual(count_a);
+  if (lardon3d_sift_vulkan_top2(backend, a.data(), count_a, b.data(), count_b,
+                                actual.data(), actual.size()) !=
+      LARDON3D_ORB_VULKAN_OK) {
+    return false;
+  }
+  for (uint32_t query = 0; query < count_a; ++query) {
+    ++comparison->queries;
+    const auto &cpu = expected[query];
+    const auto &gpu = actual[query];
+    uint32_t neighbors = static_cast<uint32_t>(cpu.size());
+    bool indices_equal = gpu.neighbor_count == neighbors;
+    if (neighbors > 0) {
+      indices_equal = indices_equal &&
+                      gpu.best_index == static_cast<uint32_t>(cpu[0].trainIdx);
+    }
+    if (neighbors > 1) {
+      indices_equal = indices_equal &&
+                      gpu.second_index == static_cast<uint32_t>(cpu[1].trainIdx);
+    }
+    if (!indices_equal) {
+      ++comparison->index_divergences;
+      if (comparison->index_divergences <= 4) {
+        std::fprintf(stderr,
+                     "index divergence q=%u cpu=(%d,%d) gpu=(%u,%u)\n",
+                     query, cpu.empty() ? -1 : cpu[0].trainIdx,
+                     cpu.size() < 2 ? -1 : cpu[1].trainIdx,
+                     gpu.best_index, gpu.second_index);
+      }
+    }
+    if (neighbors > 1) {
+      float best = std::sqrt(gpu.best_squared_distance);
+      float second = std::sqrt(gpu.second_squared_distance);
+      if (best != cpu[0].distance || second != cpu[1].distance) {
+        ++comparison->distance_bit_divergences;
+      }
+      bool cpu_lowe = cpu[1].distance != 0.0F &&
+                      cpu[0].distance < 0.7F * cpu[1].distance;
+      bool gpu_lowe = second != 0.0F && best < 0.7F * second;
+      if (cpu_lowe != gpu_lowe) {
+        ++comparison->lowe_divergences;
+      }
+    }
+  }
+  return true;
+}
+
+static bool files_and_hashes_differ(const char *path_a, const char *path_b) {
+  std::array<unsigned char, 128> bytes_a{};
+  std::array<unsigned char, 128> bytes_b{};
+  int fd_a = open(path_a, O_RDONLY | O_CLOEXEC);
+  int fd_b = open(path_b, O_RDONLY | O_CLOEXEC);
+  if (fd_a < 0 || fd_b < 0) {
+    if (fd_a >= 0) {
+      close(fd_a);
+    }
+    if (fd_b >= 0) {
+      close(fd_b);
+    }
+    return false;
+  }
+  ssize_t size_a = read(fd_a, bytes_a.data(), bytes_a.size());
+  ssize_t size_b = read(fd_b, bytes_b.data(), bytes_b.size());
+  close(fd_a);
+  close(fd_b);
+  if (size_a <= 0 || size_a != size_b) {
+    return false;
+  }
+
+  unsigned char hash_a[EVP_MAX_MD_SIZE];
+  unsigned char hash_b[EVP_MAX_MD_SIZE];
+  unsigned int hash_size_a = 0;
+  unsigned int hash_size_b = 0;
+  bool hashed_a = EVP_Digest(bytes_a.data(), static_cast<size_t>(size_a), hash_a,
+                             &hash_size_a, EVP_sha256(), nullptr) == 1;
+  bool hashed_b = EVP_Digest(bytes_b.data(), static_cast<size_t>(size_b), hash_b,
+                             &hash_size_b, EVP_sha256(), nullptr) == 1;
+  return hashed_a && hashed_b &&
+         std::memcmp(bytes_a.data(), bytes_b.data(), static_cast<size_t>(size_a)) != 0 &&
+         (hash_size_a != hash_size_b ||
+          std::memcmp(hash_a, hash_b, hash_size_a) != 0);
+}
+
+static bool match_file_difference(Lardon3DOrbVulkanBackend *backend) {
+  std::vector<float> query(128, 0.0F);
+  std::vector<float> candidates(2 * 128);
+  for (size_t component = 0; component < 128; ++component) {
+    float difference = component < 16 ? 0.25F :
+                       component < 64 ? 0.00025F : 0.00000025F;
+    candidates[component] = difference;
+    candidates[128 + component] = difference * 2.0F;
+  }
+  cv::Mat matrix_a(1, 128, CV_32FC1, query.data());
+  cv::Mat matrix_b(2, 128, CV_32FC1, candidates.data());
+  cv::BFMatcher matcher(cv::NORM_L2, false);
+  std::vector<std::vector<cv::DMatch>> cpu;
+  matcher.knnMatch(matrix_a, matrix_b, cpu, 2);
+  Lardon3DSiftTop2 gpu{};
+  if (cpu.size() != 1 || cpu[0].size() != 2 ||
+      lardon3d_sift_vulkan_top2(backend, query.data(), 1, candidates.data(), 2,
+                                &gpu, 1) != LARDON3D_ORB_VULKAN_OK ||
+      cpu[0][0].trainIdx != static_cast<int>(gpu.best_index) ||
+      !(cpu[0][0].distance < 0.7F * cpu[0][1].distance)) {
+    return false;
+  }
+  Lardon3DMatchFileEntry cpu_entry = {
+      0, static_cast<uint32_t>(cpu[0][0].trainIdx), cpu[0][0].distance};
+  Lardon3DMatchFileEntry gpu_entry = {
+      0, gpu.best_index, std::sqrt(gpu.best_squared_distance)};
+  if (cpu_entry.distance == gpu_entry.distance) {
+    return false;
+  }
+
+  char cpu_path[] = "/tmp/lardon3d-sift-cpu-match-XXXXXX";
+  char gpu_path[] = "/tmp/lardon3d-sift-gpu-match-XXXXXX";
+  int cpu_fd = mkstemp(cpu_path);
+  int gpu_fd = mkstemp(gpu_path);
+  bool ok = cpu_fd >= 0 && gpu_fd >= 0;
+  if (ok) {
+    ok = lardon3d_match_file_write(cpu_fd, 2, 128, 1, 2, &cpu_entry, 1) ==
+             LARDON3D_MATCH_FILE_OK &&
+         lardon3d_match_file_write(gpu_fd, 2, 128, 1, 2, &gpu_entry, 1) ==
+             LARDON3D_MATCH_FILE_OK;
+  }
+  if (cpu_fd >= 0) {
+    close(cpu_fd);
+  }
+  if (gpu_fd >= 0) {
+    close(gpu_fd);
+  }
+  ok = ok && files_and_hashes_differ(cpu_path, gpu_path);
+  unlink(cpu_path);
+  unlink(gpu_path);
+  return ok;
+}
+
+static bool run_distribution(bool rootsift) {
+  Lardon3DOrbVulkanBackend *backend = lardon3d_orb_vulkan_backend_create();
+  if (!backend) {
+    return false;
+  }
+  Comparison comparison;
+  bool ok = true;
+  const uint32_t sizes[] = {1, 2, 16, 64, 256, 1024, 8192};
+  for (uint32_t size : sizes) {
+    std::vector<float> a = make_descriptors(size, rootsift, 0x12345678U + size);
+    std::vector<float> b = make_descriptors(size, rootsift, 0x87654321U + size);
+    ok = compare_case(backend, a, b, &comparison) && ok;
+  }
+
+  std::vector<float> ties(8 * 128, 0.25F);
+  ok = compare_case(backend, ties, ties, &comparison) && ok;
+
+  std::vector<float> query(128, 0.0F);
+  std::vector<float> near = make_descriptors(4096, false, 0x31415926U);
+  for (float &value : near) {
+    value = (value - 0.5F) * 0.0001F;
+  }
+  ok = compare_case(backend, query, near, &comparison) && ok;
+
+  std::vector<float> equal_sums(1024 * 128);
+  for (uint32_t row = 0; row < 1024; ++row) {
+    uint32_t offset = (row * 37U) % 128U;
+    for (uint32_t component = 0; component < 128; ++component) {
+      uint32_t source = (component + offset) % 128U;
+      float magnitude = source < 16U ? 0.25F :
+                        source < 64U ? 0.00025F : 0.00000025F;
+      equal_sums[static_cast<size_t>(row) * 128 + component] = magnitude;
+    }
+  }
+  ok = compare_case(backend, query, equal_sums, &comparison) && ok;
+
+  const float ratios[] = {
+      std::nextafter(0.7F, 0.0F),
+      0.7F,
+      std::nextafter(0.7F, std::numeric_limits<float>::infinity()),
+  };
+  for (float ratio : ratios) {
+    std::vector<float> boundary(2 * 128, 0.0F);
+    boundary[0] = ratio;
+    boundary[128] = 1.0F;
+    ok = compare_case(backend, query, boundary, &comparison) && ok;
+  }
+
+  const uint32_t asymmetric[][2] = {
+      {256, 8192}, {1024, 4096}, {4096, 1024}, {8192, 256},
+  };
+  for (const auto &counts : asymmetric) {
+    std::vector<float> asymmetric_a =
+        make_descriptors(counts[0], rootsift, 0x10203040U + counts[0]);
+    std::vector<float> asymmetric_b =
+        make_descriptors(counts[1], rootsift, 0x50607080U + counts[1]);
+    ok = compare_case(backend, asymmetric_a, asymmetric_b, &comparison) && ok;
+  }
+
+  cv::Mat image(768, 768, CV_8UC1);
+  uint32_t image_seed = rootsift ? 0x13572468U : 0x24681357U;
+  for (int row = 0; row < image.rows; ++row) {
+    for (int column = 0; column < image.cols; ++column) {
+      uint32_t noise = random_u32(&image_seed);
+      uint32_t row_term = static_cast<uint32_t>(row) * 17U;
+      uint32_t column_term = static_cast<uint32_t>(column) * 31U;
+      image.at<unsigned char>(row, column) =
+          static_cast<unsigned char>((noise + row_term + column_term) & 0xffU);
+    }
+  }
+  cv::GaussianBlur(image, image, cv::Size(3, 3), 0.8);
+  cv::Ptr<cv::SIFT> sift = cv::SIFT::create(1024);
+  std::vector<cv::KeyPoint> keypoints;
+  cv::Mat real_descriptors;
+  sift->detectAndCompute(image, cv::noArray(), keypoints, real_descriptors);
+  if (real_descriptors.rows < 2 || real_descriptors.cols != 128) {
+    ok = false;
+  } else {
+    std::vector<float> real(
+        real_descriptors.ptr<float>(),
+        real_descriptors.ptr<float>() +
+            static_cast<size_t>(real_descriptors.rows) * 128);
+    if (rootsift) {
+      for (int row = 0; row < real_descriptors.rows; ++row) {
+        float sum = 0.0F;
+        for (size_t component = 0; component < 128; ++component) {
+          sum += real[static_cast<size_t>(row) * 128 + component];
+        }
+        for (size_t component = 0; component < 128; ++component) {
+          float &value = real[static_cast<size_t>(row) * 128 + component];
+          value = std::sqrt(value / sum);
+        }
+      }
+    }
+    ok = compare_case(backend, real, real, &comparison) && ok;
+  }
+
+  std::vector<Lardon3DSiftTop2> reference(near.size() / 128);
+  std::vector<Lardon3DSiftTop2> repeated(reference.size());
+  ok = lardon3d_sift_vulkan_top2(
+           backend, near.data(), static_cast<uint32_t>(reference.size()),
+           near.data(), static_cast<uint32_t>(reference.size()),
+           reference.data(), reference.size()) == LARDON3D_ORB_VULKAN_OK && ok;
+  for (int repetition = 0; repetition < 100 && ok; ++repetition) {
+    ok = lardon3d_sift_vulkan_top2(
+             backend, near.data(), static_cast<uint32_t>(repeated.size()),
+             near.data(), static_cast<uint32_t>(repeated.size()),
+             repeated.data(), repeated.size()) == LARDON3D_ORB_VULKAN_OK &&
+         std::memcmp(reference.data(), repeated.data(),
+                     reference.size() * sizeof(reference[0])) == 0;
+  }
+
+  if (!fp64_requested()) {
+    ok = match_file_difference(backend) && ok;
+  }
+
+  std::printf("%s queries=%llu index_divergences=%llu distance_divergences=%llu "
+              "lowe_divergences=%llu\n",
+              rootsift ? "rootsift" : "sift",
+              static_cast<unsigned long long>(comparison.queries),
+              static_cast<unsigned long long>(comparison.index_divergences),
+              static_cast<unsigned long long>(comparison.distance_bit_divergences),
+              static_cast<unsigned long long>(comparison.lowe_divergences));
+  lardon3d_orb_vulkan_backend_destroy(backend);
+  return ok;
+}
+
+int main() {
+  cv::setNumThreads(12);
+  return run_distribution(false) && run_distribution(true)
+             ? EXIT_SUCCESS
+             : EXIT_FAILURE;
+}

tools/embed_spirv.py

@@ -15,14 +15,17 @@ def main() -> int:
 
     words = struct.unpack(f"<{len(source) // 4}I", source)
     output = pathlib.Path(sys.argv[2])
+    stem = output.stem
+    symbol = f"lardon3d_{stem}"
+    guard = f"LARDON3D_{stem.upper()}_H"
     lines = [
-        "#ifndef LARDON3D_ORB_TOP2_SPV_H",
+        f"#ifndef {guard}",
-        "#define LARDON3D_ORB_TOP2_SPV_H",
+        f"#define {guard}",
         "",
         "#include <stddef.h>",
         "#include <stdint.h>",
         "",
-        "static const uint32_t lardon3d_orb_top2_spv[] = {",
+        f"static const uint32_t {symbol}[] = {{",
     ]
     for start in range(0, len(words), 6):
         chunk = ", ".join(f"0x{word:08x}U" for word in words[start : start + 6])
@@ -30,8 +33,8 @@ def main() -> int:
     lines.extend(
         [
             "};",
-            "static const size_t lardon3d_orb_top2_spv_size =",
+            f"static const size_t {symbol}_size =",
-            "    sizeof(lardon3d_orb_top2_spv);",
+            f"    sizeof({symbol});",
             "",
             "#endif",
             "",