Image registration establishes a geometric mapping between two or more images capturing the same physical scene. The objective is to align corresponding structural elements by minimizing spatial discrepancies. The process typically relies on spatial transformation models, which are categorized based on their mathematical complexity.
Spatial Transformation Models
Linear transformations handle rigid body movements, including translation, rotation, and scaling. Affine transformations preserve parallelism but allow shearing. For scenes with perspective distortion, projective transformations are employed. A projective model involves eight degrees of freedom, requiring a minimum of four control point pairs for solution.
Nonlinear transformations address complex distortions. Polynomial models are widely used, where the order of the polynomial dictates the fitting accuracy. A second-order polynomial requires solving for twelve coefficients, necessitating at least six control points. Higher orders improve local adaptation but increase computational load and risk overfitting.
Multi-Scale Feature Detection
Robust registration depends on identifying stable keypoints across varying scales. The Harris corner detector evaluates intensity variations within a local neighborhood. When sliding a window in any direction produces significant intensity changes, the center pixel is classified as a corner.
Custom Filtering Kernels
Before gradient computation, image smoothing is applied using directional filters. Below are implementations for a ROEWA exponential decay filter and a circular Gaussian kernel.
#include <opencv2/opencv.hpp>
#include <cmath>
cv::Mat generate_roewa_filter(int radius, float decay_factor) {
int diameter = 2 * radius + 1;
cv::Mat filter = cv::Mat::zeros(diameter, diameter, CV_32F);
for (int y = -radius; y <= radius; ++y) {
float* row_ptr = filter.ptr<float>(y + radius);
for (int x = -radius; x <= radius; ++x) {
row_ptr[x + radius] = std::exp(-std::abs(static_cast<float>(x + y)) / decay_factor);
}
}
return filter;
}
cv::Mat create_circular_gaussian(int radius, float sigma) {
int diameter = 2 * radius + 1;
cv::Mat kernel = cv::Mat::zeros(diameter, diameter, CV_32F);
float inv_2sigma2 = -1.0f / (2.0f * sigma * sigma);
float total_weight = 0.0f;
for (int y = -radius; y <= radius; ++y) {
float* row_ptr = kernel.ptr<float>(y + radius);
for (int x = -radius; x <= radius; ++x) {
if ((x * x + y * y) <= radius * radius) {
float val = std::exp((x * x + y * y) * inv_2sigma2);
row_ptr[x + radius] = val;
total_weight += val;
}
}
}
kernel /= total_weight;
return kernel;
}
Scale-Space Gradient Computation
To detect features invariant to scale, the algorithm constructs a pyramid of responses. Directional filtering decomposes the smoothing kernel to compute gradients efficiently.
void compute_multi_scale_features(const cv::Mat& source_img, std::vector<cv::Mat>& response_maps,
std::vector<cv::Mat>& magnitudes, std::vector<cv::Mat>& orientations) {
cv::Mat gray;
if (source_img.channels() > 1) cv::cvtColor(source_img, gray, cv::COLOR_RGB2GRAY);
else gray = source_img;
gray.convertTo(gray, CV_32F, 1.0 / 255.0);
const int scales = 5;
response_maps.resize(scales);
magnitudes.resize(scales);
orientations.resize(scales);
for (int s = 0; s < scales; ++s) {
float current_sigma = static_cast<float>(0.8f * std::pow(1.2f, s));
int k_radius = cvRound(2.0f * current_sigma);
cv::Mat base_kernel = generate_roewa_filter(k_radius, current_sigma);
cv::Mat q_bottom_right, q_top_left, q_right, q_left;
base_kernel(cv::Range(k_radius + 1, 2 * k_radius + 1), cv::Range::all())
.copyTo(q_bottom_right(cv::Range(k_radius + 1, 2 * k_radius + 1), cv::Range::all()));
base_kernel(cv::Range(0, k_radius), cv::Range::all())
.copyTo(q_top_left(cv::Range(0, k_radius), cv::Range::all()));
base_kernel(cv::Range::all(), cv::Range(k_radius + 1, 2 * k_radius + 1))
.copyTo(q_right(cv::Range::all(), cv::Range(k_radius + 1, 2 * k_radius + 1)));
base_kernel(cv::Range::all(), cv::Range(0, k_radius))
.copyTo(q_left(cv::Range::all(), cv::Range(0, k_radius)));
cv::Mat m_br, m_tl, m_r, m_l;
cv::filter2D(gray, m_br, -1, q_bottom_right);
cv::filter2D(gray, m_tl, -1, q_top_left);
cv::filter2D(gray, m_r, -1, q_right);
cv::filter2D(gray, m_l, -1, q_left);
cv::Mat gx, gy;
cv::log(m_r / m_l, gx);
cv::log(m_br / m_tl, gy);
cv::magnitude(gx, gy, magnitudes[s]);
cv::phase(gx, gy, orientations[s], true);
cv::Mat tensor_11 = (current_sigma * current_sigma) * gx.mul(gx);
cv::Mat tensor_12 = (current_sigma * current_sigma) * gx.mul(gy);
cv::Mat tensor_22 = (current_sigma * current_sigma) * gy.mul(gy);
float blur_sigma = std::sqrt(2.0f) * current_sigma;
int blur_size = cvRound(3.0f * blur_sigma) * 2 + 1;
cv::GaussianBlur(tensor_11, tensor_11, cv::Size(blur_size, blur_size), blur_sigma);
cv::GaussianBlur(tensor_12, tensor_12, cv::Size(blur_size, blur_size), blur_sigma);
cv::GaussianBlur(tensor_22, tensor_22, cv::Size(blur_size, blur_size), blur_sigma);
float det = tensor_11.mul(tensor_22) - tensor_12.mul(tensor_12);
float trace_sq = (tensor_11 + tensor_22).mul(tensor_11 + tensor_22);
response_maps[s] = det - 0.04f * trace_sq;
}
}
Keypoint Localization & Orientation Assignment
Features are identified by locating local maxima across scale and space. A 3x3x3 neighborhood comparison ensures stability.
std::vector<cv::KeyPoint> locate_keypoints(const std::vector<cv::Mat>& responses,
const std::vector<cv::Mat>& mags,
const std::vector<cv::Mat>& orients,
float threshold) {
std::vector<cv::KeyPoint> keypoints;
int num_scales = responses.size();
int rows = responses[0].rows;
int cols = responses[0].cols;
for (int s = 1; s < num_scales - 1; ++s) {
const cv::Mat& curr_resp = responses[s];
for (int r = 1; r < rows - 1; ++r) {
const float* up = responses[s-1].ptr<float>(r);
const float* mid = curr_resp.ptr<float>(r);
const float* down = responses[s+1].ptr<float>(r);
const float* up_prev = responses[s-1].ptr<float>(r-1);
const float* up_next = responses[s-1].ptr<float>(r+1);
const float* mid_prev = curr_resp.ptr<float>(r-1);
const float* mid_next = curr_resp.ptr<float>(r+1);
const float* down_prev = responses[s+1].ptr<float>(r-1);
const float* down_next = responses[s+1].ptr<float>(r+1);
for (int c = 1; c < cols - 1; ++c) {
float center_val = mid[c];
if (center_val > threshold && center_val > mid[c-1] && center_val > mid[c+1] &&
center_val > mid_prev[c] && center_val > mid_next[c] &&
center_val > up_prev[c-1] && center_val > up_prev[c] && center_val > up_prev[c+1] &&
center_val > up[c-1] && center_val > up[c] && center_val > up[c+1] &&
center_val > up_next[c-1] && center_val > up_next[c] && center_val > up_next[c+1] &&
center_val > down_prev[c-1] && center_val > down_prev[c] && center_val > down_prev[c+1] &&
center_val > down[c-1] && center_val > down[c] && center_val > down[c+1] &&
center_val > down_next[c-1] && center_val > down_next[c] && center_val > down_next[c+1]) {
float sigma = 0.8f * std::pow(1.2f, s);
std::vector<float> hist(36, 0.0f);
float max_hist = compute_orientation_histogram(mags[s], orients[s], cv::Point2f(c, r), sigma, hist);
float peak_thr = 0.8f * max_hist;
for (int b = 0; b < 36; ++b) {
int left = (b == 0) ? 35 : b - 1;
int right = (b == 35) ? 0 : b + 1;
if (hist[b] >= peak_thr && hist[b] >= hist[left] && hist[b] >= hist[right]) {
float bin_offset = 0.5f * (hist[left] - hist[right]) / (hist[left] + hist[right] - 2.0f * hist[b]);
float dominant_angle = (static_cast<float>(b) + bin_offset) * (360.0f / 36.0f);
cv::KeyPoint kp;
kp.pt = cv::Point2f(c, r);
kp.size = sigma * 3.0f;
kp.angle = dominant_angle;
kp.response = center_val;
kp.octave = s;
keypoints.push_back(kp);
}
}
}
}
}
}
return keypoints;
}
The dominant orientation is determined by aggregating gradient directiosn within a spatial window. The histogram approach bins angles and applies Gaussian weighting to reduce boundary effects.
float compute_orientation_histogram(const cv::Mat& mags, const cv::Mat& orients,
cv::Point2f center, float scale, std::vector<float>& hist) {
int radius = cvRound(3.0f * scale);
int rows = mags.rows, cols = mags.cols;
int x_min = std::max(0, cvRound(center.x) - radius);
int x_max = std::min(cols - 1, cvRound(center.x) + radius);
int y_min = std::max(0, cvRound(center.y) - radius);
int y_max = std::min(rows - 1, cvRound(center.y) + radius);
std::vector<float> temp_hist(40, 0.0f);
float inv_2sigma2 = -1.0f / (2.0f * scale * scale);
for (int y = y_min; y <= y_max; ++y) {
const float* m_ptr = mags.ptr<float>(y);
const float* o_ptr = orients.ptr<float>(y);
for (int x = x_min; x <= x_max; ++x) {
float dx = x - center.x;
float dy = y - center.y;
if (dx * dx + dy * dy < radius * radius) {
float dist_weight = std::exp((dx * dx + dy * dy) * inv_2sigma2);
int bin = cvRound(o_ptr[x] / 10.0f);
bin = (bin < 0) ? bin + 36 : (bin >= 36 ? bin - 36 : bin);
temp_hist[bin + 2] += m_ptr[x] * dist_weight;
}
}
}
for (int i = 0; i < 36; ++i) {
hist[i] = (temp_hist[i] * 0.375f) +
(temp_hist[i+1] + temp_hist[i-1]) * 0.25f +
(temp_hist[i+2] + temp_hist[i-2]) * 0.0625f;
}
float max_val = *std::max_element(hist.begin(), hist.end());
return max_val;
}
Descriptor Construction
Once keypoints are localized and oriented, a descriptor vector encodes the local image structure. The GLOH descriptor organizes gradients in a log-polar grid, dividing the region into three concentric rings with varying angular sectors.
#include <algorithm>
#include <numeric>
std::vector<float> extract_gloh_features(const cv::Mat& mags, const cv::Mat& orients,
cv::Point2f center, float scale, float angle) {
float rad = angle * CV_PI / 180.0f;
float cos_a = std::cos(-rad);
float sin_a = std::sin(-rad);
int radius = cvRound(6.0f * scale);
int x_min = std::max(0, cvRound(center.x) - radius);
int x_max = std::min(mags.cols - 1, cvRound(center.x) + radius);
int y_min = std::max(0, cvRound(center.y) - radius);
int y_max = std::min(mags.rows - 1, cvRound(center.y) + radius);
int bin_r = 3;
int bin_c = 8;
int bin_o = 8;
int total_bins = bin_r * bin_c * bin_o;
std::vector<float> descriptor(total_bins, 0.0f);
for (int y = y_min; y <= y_max; ++y) {
const float* m_ptr = mags.ptr<float>(y);
const float* o_ptr = orients.ptr<float>(y);
for (int x = x_min; x <= x_max; ++x) {
float dx = x - center.x;
float dy = y - center.y;
float r_rot = dx * cos_a - dy * sin_a;
float c_rot = dx * sin_a + dy * cos_a;
if (dx*dx + dy*dy < radius*radius) {
float dist = std::sqrt(r_rot*r_rot + c_rot*c_rot);
int r_idx = dist < radius/3.0f ? 0 : (dist < radius*2.0f/3.0f ? 1 : 2);
float ang = std::atan2(c_rot, r_rot) * 180.0f / CV_PI;
ang = (ang < 0) ? ang + 360.0f : ang;
int c_idx = cvRound(ang / 45.0f);
c_idx = (c_idx >= 8) ? 0 : c_idx;
int o_idx = cvRound(o_ptr[x] / 45.0f);
o_idx = (o_idx >= 8) ? 0 : o_idx;
int idx = (r_idx * bin_c + c_idx) * bin_o + o_idx;
descriptor[idx] += m_ptr[x];
}
}
}
float norm = std::sqrt(std::inner_product(descriptor.begin(), descriptor.end(), descriptor.begin(), 0.0f));
std::transform(descriptor.begin(), descriptor.end(), descriptor.begin(),
[norm](float v){ return v / norm; });
std::transform(descriptor.begin(), descriptor.end(), descriptor.begin(),
[](float v){ return std::min(v, 0.2f); });
norm = std::sqrt(std::inner_product(descriptor.begin(), descriptor.end(), descriptor.begin(), 0.0f));
std::transform(descriptor.begin(), descriptor.end(), descriptor.begin(),
[norm](float v){ return v / norm; });
return descriptor;
}
Alternative descriptors like SIFT utilize a 4x4 spatial grid with 8 orientation bins per cell, yielding a 128-dimensional vector. SURF relies on Haar wavelet responses with in sub-regions, while Shape Context matches points based on log-polar distributions of relative coordinates.
Correspondence Establishment
Matching pairs descriptors from different images by measuring similarity. The Nearest Neighbor Distance Ratio (NNDR) compares the closest match to the second-closest to filter ambiguous correspondances.
struct MatchResult {
int idx1, idx2;
float dist, dist2;
};
std::vector<MatchResult> compute_feature_matches(const cv::Mat& des1, const cv::Mat& des2, float ratio_threshold) {
std::vector<MatchResult> matches;
int cols = des1.cols;
for (int i = 0; i < des1.rows; ++i) {
const float* d1 = des1.ptr<float>(i);
float min_d = std::numeric_limits<float>::max();
float sec_d = std::numeric_limits<float>::max();
int idx_min = -1, idx_sec = -1;
for (int j = 0; j < des2.rows; ++j) {
const float* d2 = des2.ptr<float>(j);
float dist_sq = 0.0f;
for (int k = 0; k < cols; ++k) {
float diff = d1[k] - d2[k];
dist_sq += diff * diff;
}
float dist = std::sqrt(dist_sq);
if (dist < min_d) {
sec_d = min_d; idx_sec = idx_min;
min_d = dist; idx_min = j;
} else if (dist < sec_d) {
sec_d = dist; idx_sec = j;
}
}
if (idx_min != -1 && (min_d / sec_d) < ratio_threshold) {
matches.push_back({i, idx_min, min_d, sec_d});
}
}
return matches;
}
For large-scale datasets, brute-force comparison becomes inefficient. Space-partitioning structures like k-d trees accelerate nearest-neighbor searches by recursively dividing the feature space. Alternatively, Locality Sensitive Hashing (LSH) projects high-dimensional vectors into buckets where similar descriptors collide with high probability, enabling sub-linear retrieval time.