文章目录
前言
最近在研究关于红外图像和自然光影像匹配的问题，参考到了李加元老师的文章：RIFT: Multi-modal image matching based on radiation-variation insensitive feature transform
.文章源码为matlab，以下整理一下python代码的复现思路
一、图像读取后的通道统一
红外图像是单通道，自然管图像为三通道，读取数据后需要将红外影像复制为三通道，使其能够对应到RGB三个通道的特征if img1.ndim == 2:
img = np.zeros([len(img1), len(img1[0]), 3], int)
for i in range(3):
img[..., i] = img1
img0 = img1
img1 = img
if img2.ndim == 2:
img = np.zeros([len(img2), len(img2[0]), 3], int)
for i in range(3):
img[..., i] = img2
img2 = img
二、phasecong函数
2.1 方法介绍
phasecong是一个对比度不变的边缘和角点检测器，计算图像相位一致性的函数。它最开始是matlab的一个函数（文章中用到的为matlab中的phasecong3），幸有Alistair Muldal大佬实现了python的phasecong,可以直接cv来用。
2.2 参数及返回值
phasecong参数和返回值列表与matlab版本的相同，文章模型只需要获取其返回的M和EO：
参数:
------
<Name> <Default> <Description>
img N/A The input image
nscale 5 Number of wavelet scales, try values 3-6
norient 6 Number of filter orientations.
minWaveLength 3 Wavelength of smallest scale filter.
mult 2.1 Scaling factor between successive filters.
sigmaOnf 0.55 Ratio of the standard deviation of the Gaussian
describing the log Gabor filter's transfer function
in the frequency domain to the filter center
frequency.
k 2.0 No. of standard deviations of the noise energy
beyond the mean at which we set the noise threshold
point. You may want to vary this up to a value of
10 or 20 for noisy images.
cutOff 0.5 The fractional measure of frequency spread below
which phase congruency values get penalized.
g 10 Controls the 'sharpness' of the transition in the
sigmoid function used to weight phase congruency
for frequency spread.
noiseMethod -1 Parameter specifies method used to determine
noise statistics.
-1 use median of smallest scale filter responses
-2 use mode of smallest scale filter responses
>=0 use this value as the fixed noise threshold
返回值:
---------
M Maximum moment of phase congruency covariance, which can be used as
a measure of edge strength
m Minimum moment of phase congruency covariance, which can be used as
a measure of corner strength
ori Orientation image, in integer degrees (0-180), positive angles
anti-clockwise.
ft Local weighted mean phase angle at every point in the image. A
value of pi/2 corresponds to a bright line, 0 to a step and -pi/2
to a dark line.
PC A list of phase congruency images (values between 0 and 1), one per
orientation.
EO A list containing the complex-valued convolution results (see
below)
T Calculated noise threshold (can be useful for diagnosing noise
characteristics of images). Once you know this you can then specify
fixed thresholds and save some computation time.
m1, __, __, __, __, eo1, __ = phasecong.phasecong(img=img1, nscale=4, norient=6, minWaveLength=3, mult=1.6,
sigmaOnf=0.75, g=3, k=1)
m2, __, __, __, __, eo2, __ = phasecong.phasecong(img=img2, nscale=4, norient=6, minWaveLength=3, mult=1.6,
sigmaOnf=0.75, g=3, k=1)
三、利用opencv的快速检测提取边缘特征点
使用快速检测算法实现提取点操作
注意：需要先将提取到的特征转为unit8格式，不然会报错
m1, m2 = map(lambda img: (img.astype(np.float) - img.min()) / (img.max() - img.min()), (m1, m2))
cm1 = m1 * 255
cm2 = m2 * 255
fast = cv2.FastFeatureDetector_create(nonmaxSuppression=True, type=cv2.FAST_FEATURE_DETECTOR_TYPE_7_12)
kp1 = fast.detect(np.uint8(cm1), None)
kp2 = fast.detect(np.uint8(cm2), None)
四、复现RIFT_descriptor_no_rotation_invariance函数
matlab可以直接操作矩阵，python里面对矩阵的操作主要依赖于numpy的array，array需要从detect()的返回值keypoint对象中提取坐标（但好像python也有matrix，不太了解就没用）
4.1 复现RIFT_descriptor_no_rotation_invariance函数
def RIFT_descriptor_no_rotation_invariance(img, kps, eo, patch_size, s, o):
KPS = kps.T
(yim, xim, _) = np.shape(img)
CS = np.zeros([yim, xim, o], np.float)
for j in range(o):
for i in range(s):
# 将各个scale的变换结果的幅度相加
CS[..., j] = CS[..., j] + np.abs(np.array(eo[j][i]))
mim = np.argmax(CS, axis=2)
des = np.zeros([36 * o, np.size(KPS, 1)])
kps_to_ignore = np.ones([1, np.size(KPS, 1)], bool)
for k in range(np.size(KPS, 1)):
x = round(KPS[0][k])
y = round(KPS[1][k])
x1 = max(0, x - math.floor(patch_size / 2))
y1 = max(0, y - math.floor(patch_size / 2))
x2 = min(x + math.floor(patch_size / 2), np.size(img, 1))
y2 = min(y + math.floor(patch_size / 2), np.size(img, 0))
if y2 - y1 != patch_size or x2 - x1 != patch_size:
kps_to_ignore[0][i] = 0
continue
patch = mim[y1:y2, x1:x2]
ys, xs = np.size(patch, 0), np.size(patch, 1)
ns = 6;
RIFT_des = np.zeros([ns, ns, o])
for j in range(ns):
for i in range(ns):
clip = patch[round((j) * ys / ns):round((j + 1) * ys / ns),
round((i) * xs / ns): round((i + 1) * xs / ns)]
x, __ = np.histogram(clip.T.flatten(), bins=6, range=(0, o), density=False)
te = RIFT_des[j][i]
RIFT_des[j][i] = x.reshape(1, 1, len(x))
RIFT_des = RIFT_des.T.flatten()
df = np.linalg.norm(RIFT_des)
if df != 0:
RIFT_des = RIFT_des / df
des[:, [k]] = np.expand_dims(RIFT_des, axis=1)
m = repeat(kps_to_ignore, '1 n -> c n', c=2)
v = KPS[m]
KPS_out = v.reshape(2, int(len(v) / 2)).T
w = repeat(kps_to_ignore, '1 n -> c n', c=len(des))
z = des[w]
des_out = z.reshape(len(des), int(len(z) / len(des))).T
des_out = np.float32(des_out) * 100
return KPS_out, des_out
4.2 复现FSC函数
（只保留了’affine’部分）
def FSC(cor1, cor2, change_form, error_t):
(M, N) = np.shape(cor1)
if (change_form == 'similarity'):
n = 2
max_iteration = M * (M - 1) / 2
elif (change_form == 'affine'):
n = 3
max_iteration = M * (M - 1) * (M - 2) / (2 * 3)
elif (change_form == 'perspective'):
n = 4
max_iteration = M * (M - 1) * (M - 2) / (2 * 3)
if (max_iteration > 10000):
iterations = 10000
else:
iterations = max_iteration
most_consensus_number = 0
cor1_new = np.zeros([M, N])
cor2_new = np.zeros([M, N])
for i in range(iterations):
while (True):
a = np.floor(1 + (M - 1) * np.random.rand(1, n)).astype(np.int)[0]
cor11 = cor1[a]
cor22 = cor2[a]
if n == 2 and (a[0] != a[1]) and sum(cor11[0] != cor11[1]) and sum(cor22[0] != cor22[1]):
break
if n == 3 and (a[0] != a[1] and a[0] != a[2] and a[1] != a[2]) and sum(cor11[0] != cor11[1]) and sum(
cor11[0] != cor11[2]) and sum(cor11[1] != cor11[2]) and sum(cor22[0] != cor22[1]) and sum(
cor22[0] != cor22[2]) and sum(cor22[1] != cor22[2]):
break
if n == 4 and (
a[0] != a[1] and a[0] != a[2] and a[0] != a[3] and a[1] != a[2] and a[1] != a[3] and a[2] != a[
3]) and sum(cor11[0] != cor11[1]) and sum(cor11[0] != cor11[2]) and sum(cor11[0] != cor11[3]) and sum(
cor11[1] != cor11[2]) and sum(cor11[1] != cor11[3]) and sum(cor11[2] != cor11[3]) and sum(
cor22[0] != cor11[1]) and sum(cor22[0] != cor22[2]) and sum(cor22[0] != cor22[3]) and sum(
cor22[1] != cor22[2]) and sum(cor22[1] != cor22[3]) and sum(cor22[2] != cor22[3]):
break
parameters, __ = LSM(cor11, cor22, change_form)
solution = np.array([[parameters[0], parameters[1], parameters[4]],
[parameters[2], parameters[3], parameters[5]],
[parameters[6], parameters[7], 1]])
match1_xy = np.ones([3, len(cor1)])
match1_xy[:2] = cor1.T
if change_form == 'affine':
t_match1_xy = solution.dot(match1_xy)
match2_xy = np.ones([3, len(cor1)])
match2_xy[:2] = cor2.T
diff_match2_xy = t_match1_xy - match2_xy
diff_match2_xy = np.sqrt(sum(np.power(diff_match2_xy, 2)))
index_in = np.argwhere(diff_match2_xy < error_t)
consensus_num = len(index_in)
index_in = np.squeeze(index_in)
if consensus_num > most_consensus_number:
most_consensus_number = consensus_num
cor1_new = cor1[index_in]
cor2_new = cor2[index_in]
unil = cor1_new
__, IA = np.unique(unil, return_index=True, axis=0)
IA_new = np.sort(IA)
cor1_new = cor1_new[IA_new]
cor2_new = cor2_new[IA_new]
unil = cor2_new
__, IA = np.unique(unil, return_index=True, axis=0)
IA_new = np.sort(IA)
cor1_new = cor1_new[IA_new]
cor2_new = cor2_new[IA_new]
parameters, rmse = LSM(cor1_new, cor2_new, change_form)
solution = np.array([[parameters[0], parameters[1], parameters[4]],
[parameters[2], parameters[3], parameters[5]],
[parameters[6], parameters[7], 1]])
return solution
4.3 复现LSM
def LSM(match1, match2, change_form):
A = np.zeros([2 * len(match1), 4])
for i in range(len(match1)):
A[2 * i:2 * i + 2] = np.tile(match1[i], (2, 2))
B = np.array([[1, 1, 0, 0], [0, 0, 1, 1]])
B = np.tile(B, (len(match1), 1))
A = A * B
B = np.array([[1, 0], [0, 1]])
B = np.tile(B, (len(match1), 1))
A = np.hstack((A, B))
b = match2.reshape(1, int(len(match2) * len(match2[0]))).T
if change_form == "affine":
Q, R = np.linalg.qr(A)
parameters = np.zeros([8, 1])
parameters[:6] = np.linalg.solve(R, np.dot(Q.T, b))
N = len(match1)
M = np.array([[parameters[0][0], parameters[1][0]], [parameters[2][0], parameters[3][0]]])
match1_test_trans = M.dot(match1.T) + np.tile([parameters[4], parameters[5]], (1, N))
match1_test_trans = match1_test_trans.T
test = match1_test_trans - match2
rmse = math.sqrt(sum(sum(np.power(test, 2))) / N)
return np.squeeze(parameters), rmse
五、复现RIFT_demo
img1 = cv2.imread('./sar-optical/box.png')
img2 = cv2.imread('./sar-optical/box_in_scene.png')
img0 = None
# 统一图像通道数
if img1.ndim == 2:
img = np.zeros([len(img1), len(img1[0]), 3], int)
for i in range(3):
img[..., i] = img1
img0 = img1
img1 = img
if img2.ndim == 2:
img = np.zeros([len(img2), len(img2[0]), 3], int)
for i in range(3):
img[..., i] = img2
img2 = img
# 计算图像相位一致性
m1, __, __, __, __, eo1, __ = phasecong.phasecong(img=img1, nscale=4, norient=6, minWaveLength=3, mult=1.6,
sigmaOnf=0.75, g=3, k=1)
m2, __, __, __, __, eo2, __ = phasecong.phasecong(img=img2, nscale=4, norient=6, minWaveLength=3, mult=1.6,
sigmaOnf=0.75, g=3, k=1)
# 将提取到的特征转为unit8格式
m1, m2 = map(lambda img: (img.astype(np.float) - img.min()) / (img.max() - img.min()), (m1, m2))
cm1 = m1 * 255
cm2 = m2 * 255
fast = cv2.FastFeatureDetector_create(nonmaxSuppression=True, type=cv2.FAST_FEATURE_DETECTOR_TYPE_7_12)
kp1 = fast.detect(np.uint8(cm1), None)
kp2 = fast.detect(np.uint8(cm2), None)
m1_point = kp2point(kp1)
m2_point = kp2point(kp2)
kps1, des1 = RIFT_descriptor_no_rotation_invariance(img1, m1_point, eo1, 96, 4, 6)
kps2, des2 = RIFT_descriptor_no_rotation_invariance(img2, m2_point, eo2, 96, 4, 6)
# 生成匹配，获取匹配点坐标序列，puthon中没有matlab的matchFeatures函数的对应函数，因此用了opencv
bf = cv2.FlannBasedMatcher()
matches = bf.match(des1, des2)
match_point1 = np.zeros([len(matches), 2], int)
match_point2 = np.zeros([len(matches), 2], int)
# 获取匹配点对索引
for m in range(len(matches)):
match_point1[m] = kps1[matches[m].queryIdx]
match_point2[m] = kps2[matches[m].trainIdx]
# 去重
match_point2, IA = np.unique(match_point2, return_index=True, axis=0)
match_point1 = match_point1[IA]
H = FSC(match_point1, match_point2, 'affine', 2)
Y_ = np.ones([3, len(match_point1)])
Y_[:2] = match_point1.T
Y_ = H.dot(Y_)
Y_[0] = Y_[0] / Y_[2]
Y_[1] = Y_[1] / Y_[2]
E = np.sqrt(sum(np.power((Y_[0:2] - match_point2.T), 2)))
inliersIndex = np.squeeze(np.argwhere(E < 3))
clean_point1 = match_point1[inliersIndex]
clean_point2 = match_point2[inliersIndex]
kp1 = [cv2.KeyPoint(float(clean_point1[i][0]), float(clean_point1[i][1]), 1) for i in range(len(clean_point1))]
kp2 = [cv2.KeyPoint(float(clean_point2[i][0]), float(clean_point2[i][1]), 1) for i in range(len(clean_point2))]
matches = [cv2.DMatch(i, i, 1) for i in range(len(clean_point1))]
img3 = cv2.drawMatches(img1, kp1, img2, kp2, matches, None, flags=2)
cv_show('img3', img3)
扫一扫
3896
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
