我的代码长这样：import numpy as np
from scipy import ndimage, spatial
import cv2
import os
import matplotlib.pyplot as plt
ROOT_DIR = ‘/Users/mac/Desktop/Problem_Set_2-3’
IMGDIR = os.path.join(ROOT_DIR, ‘Problem2Images’)
OUTPUT_DIR = os.path.join(ROOT_DIR, ‘outputs’)
os.makedirs(OUTPUT_DIR, exist_ok=True)
def gradient_x(img):
if len(img.shape) == 3:
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
else:
gray = img
gray = ndimage.gaussian_filter(gray.astype(np.float32), sigma=1)
grad_x = ndimage.sobel(gray, axis=1, mode=‘reflect’)
return grad_x
def gradient_y(img):
if len(img.shape) == 3:
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
else:
gray = img
gray = ndimage.gaussian_filter(gray.astype(np.float32), sigma=1)
grad_y = ndimage.sobel(gray, axis=0, mode=‘reflect’)
return grad_y
def harris_response(img, alpha, win_size):
if img is None:
raise ValueError(“输入图像为空”)
if len(img.shape) == 3:
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
else:
gray = img
grad_x = cv2.Sobel(gray, cv2.CV_32F, 1, 0, ksize=3)
grad_y = cv2.Sobel(gray, cv2.CV_32F, 0, 1, ksize=3)
Ixx = grad_x * grad_x
Ixy = grad_x * grad_y
Iyy = grad_y * grad_y
sigma = win_size / 4
A = cv2.GaussianBlur(Ixx, (win_size, win_size), sigma)
B = cv2.GaussianBlur(Ixy, (win_size, win_size), sigma)
C = cv2.GaussianBlur(Iyy, (win_size, win_size), sigma)
det_M = A * C - B * B
trace_M = A + C
R = det_M - alpha * (trace_M ** 2)
return R
def corner_selection(R, thresh, min_dist):
size = min_dist * 2 + 1
R_max = ndimage.maximum_filter(R, size=size)
mask = (R == R_max) & (R > thresh * 0.08)
corners_y, corners_x = np.where(mask)
corner_responses = R[corners_y, corners_x]
sorted_indices = np.argsort(corner_responses)[::-1]
max_corners = min(2000, len(sorted_indices))
selected_indices = sorted_indices[:max_corners]
corners_x = corners_x[selected_indices]
corners_y = corners_y[selected_indices]
pix = list(zip(corners_x, corners_y))
return pix
def histogram_of_gradients(img, pix):
if img is None:
raise ValueError(“输入图像为空”)
if len(img.shape) == 3:
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
else:
gray = img
grad_x = cv2.Sobel(gray, cv2.CV_32F, 1, 0, ksize=3)
grad_y = cv2.Sobel(gray, cv2.CV_32F, 0, 1, ksize=3)
grad_mag = np.sqrt(grad_x2 + grad_y2)
grad_ori = np.arctan2(grad_y, grad_x) * 180 / np.pi
grad_ori = np.mod(grad_ori, 360)
cell_size = 8
block_size = 2
num_bins = 9
bin_width = 360 / num_bins
features = []
for point in pix:
x, y = point
patch_size = cell_size * block_size
if (x < patch_size//2 or x >= gray.shape[1] - patch_size//2 or
y < patch_size//2 or y >= gray.shape[0] - patch_size//2):
features.append(np.zeros(block_size * block_size * num_bins))
continue
keypoint_feature = []
for i in range(block_size):
for j in range(block_size):
start_x = x - patch_size//2 + j * cell_size
start_y = y - patch_size//2 + i * cell_size
cell_hist = np.zeros(num_bins)
for cy in range(cell_size):
for cx in range(cell_size):
px = start_x + cx
py = start_y + cy
magnitude = grad_mag[py, px]
orientation = grad_ori[py, px]
bin_idx = int(orientation / bin_width)
bin_idx = bin_idx % num_bins
cell_hist[bin_idx] += magnitude
keypoint_feature.extend(cell_hist)
keypoint_feature = np.array(keypoint_feature)
norm = np.linalg.norm(keypoint_feature)
if norm > 0.1:
keypoint_feature = keypoint_feature / norm
features.append(keypoint_feature)
features = np.array(features)
return features
def feature_matching(img_1, img_2):
R1 = harris_response(img_1, 0.04, 9)
R2 = harris_response(img_2, 0.04, 9)
cor1 = corner_selection(R1, 0.0005 * np.max(R1) if np.max(R1) > 0 else 0, 8)
cor2 = corner_selection(R2, 0.0005 * np.max(R2) if np.max(R2) > 0 else 0, 8)
if len(cor1) < 4 or len(cor2) < 4:
return [], []
fea1 = histogram_of_gradients(img_1, cor1)
fea2 = histogram_of_gradients(img_2, cor2)
dis = spatial.distance.cdist(fea1, fea2, metric=‘euclidean’)
matches_1to2 = []
matches_2to1 = []
threshold = 0.75
for i in range(len(cor1)):
if len(dis[i]) < 2:
continue
sorted_idx = np.argsort(dis[i])
if dis[i][sorted_idx[0]] < threshold * dis[i][sorted_idx[1]]:
matches_1to2.append((i, sorted_idx[0]))
for j in range(len(cor2)):
if len(dis[:, j]) < 2:
continue
sorted_idx = np.argsort(dis[:, j])
if dis[sorted_idx[0], j] < threshold * dis[sorted_idx[1], j]:
matches_2to1.append((sorted_idx[0], j))
consistent = set(matches_1to2) & set(matches_2to1) pixels_1, pixels_2 = [], [] for i, j in consistent: pixels_1.append(cor1[i]) pixels_2.append(cor2[j]) if len(pixels_1) < 6: all_possible = [] for i in range(len(cor1)): if len(dis[i]) == 0: continue j = np.argmin(dis[i]) all_possible.append((i, j, dis[i][j])) all_possible.sort(key=lambda x: x[2]) for i, j, _ in all_possible: if (i, j) not in consistent: pixels_1.append(cor1[i]) pixels_2.append(cor2[j]) if len(pixels_1) >= 6: break return pixels_1, pixels_2
def test_matching():
img_1 = cv2.imread(os.path.join(IMGDIR, ‘1_1.jpg’))
img_2 = cv2.imread(os.path.join(IMGDIR, ‘1_2.jpg’))
if img_1 is None or img_2 is None:
print(“⚠️ 测试图像加载失败，请检查路径”)
return
img_gray_1 = cv2.cvtColor(img_1, cv2.COLOR_BGR2GRAY)
img_gray_2 = cv2.cvtColor(img_2, cv2.COLOR_BGR2GRAY)
pixels_1, pixels_2 = feature_matching(img_1, img_2)
H_1, W_1 = img_gray_1.shape
H_2, W_2 = img_gray_2.shape
img = np.zeros((max(H_1, H_2), W_1 + W_2, 3))
img[:H_1, :W_1, (2, 1, 0)] = img_1 / 255
img[:H_2, W_1:, (2, 1, 0)] = img_2 / 255
plt.figure(figsize=(20, 10), dpi=300)
plt.imshow(img)
N = len(pixels_1)
for i in range(N):
x1, y1 = pixels_1[i]
x2, y2 = pixels_2[i]
plt.plot([x1, x2+W_1], [y1, y2])
plt.savefig(os.path.join(OUTPUT_DIR, ‘test_matching.jpg’), bbox_inches=‘tight’)
plt.close()
def compute_homography(pixels_1, pixels_2):
if len(pixels_1) < 4:
raise ValueError(“Need at least 4 point pairs to compute homography”)
A = []
for (x1, y1), (x2, y2) in zip(pixels_1, pixels_2):
A.append([x1, y1, 1, 0, 0, 0, -x2x1, -x2y1, -x2])
A.append([0, 0, 0, x1, y1, 1, -y2x1, -y2y1, -y2])
A = np.array(A)
U, S, Vt = np.linalg.svd(A)
H = Vt[-1, :].reshape(3, 3)
homo_matrix = H / H[2, 2]
return homo_matrix
def align_pair(pixels_1, pixels_2):
if len(pixels_1) < 4:
raise ValueError(“Need at least 4 point pairs for RANSAC”)
est_homo = None best_inliers = 0 max_iterations = 3000 threshold = 4.0 for _ in range(max_iterations): indices = np.random.choice(len(pixels_1), 4, replace=False) sample_1 = [pixels_1[i] for i in indices] sample_2 = [pixels_2[i] for i in indices] try: H = compute_homography(sample_1, sample_2) inlier_count = 0 for (x1, y1), (x2, y2) in zip(pixels_1, pixels_2): p1 = np.array([x1, y1, 1]) p2_est = H @ p1 p2_est = p2_est / p2_est[2] error = np.sqrt((p2_est[0] - x2)**2 + (p2_est[1] - y2)**2) if error < threshold: inlier_count += 1 if inlier_count > best_inliers: best_inliers = inlier_count est_homo = H except np.linalg.LinAlgError: continue if est_homo is None: est_homo = compute_homography(pixels_1, pixels_2) return est_homo
def stitch_blend(img_1, img_2, est_homo):
if img_1 is None or img_2 is None:
raise ValueError(“输入图像为空”)
h1, w1, d1 = img_1.shape
h2, w2, d2 = img_2.shape
img_1_uint8 = img_1.copy()
img_2_uint8 = img_2.copy()
corners_img1 = np.array([[0,0,1], [w1,0,1], [0,h1,1], [w1,h1,1]], dtype=np.float32)
corners_img1_proj = (est_homo @ corners_img1.T).T
corners_img1_proj = corners_img1_proj[:, :2] / corners_img1_proj[:, 2:]
all_x = np.concatenate([[0, w2], corners_img1_proj[:, 0]])
all_y = np.concatenate([[0, h2], corners_img1_proj[:, 1]])
x_min = int(np.floor(all_x.min()))
x_max = int(np.ceil(all_x.max()))
y_min = int(np.floor(all_y.min()))
y_max = int(np.ceil(all_y.max()))
canvas_h = y_max - y_min + 1 canvas_w = x_max - x_min + 1 est_img1 = np.zeros((canvas_h, canvas_w, 3), dtype=np.uint8) est_img2 = np.zeros_like(est_img1) x_range = np.arange(x_min, x_max + 1, 1, dtype=np.float64) y_range = np.arange(y_min, y_max + 1, 1, dtype=np.float64) x, y = np.meshgrid(x_range, y_range) try: homo_inv = np.linalg.inv(est_homo) except np.linalg.LinAlgError: homo_inv = np.linalg.pinv(est_homo) ones = np.ones_like(x, dtype=np.float64) coords = np.stack([x, y, ones], axis=-1) coords_proj = coords @ homo_inv.T coords_proj = coords_proj / coords_proj[..., 2:3] trans_x = coords_proj[..., 0].astype(np.float32) trans_y = coords_proj[..., 1].astype(np.float32) est_img1 = cv2.remap( img_1_uint8, trans_x, trans_y, interpolation=cv2.INTER_CUBIC, borderMode=cv2.BORDER_CONSTANT, borderValue=(0, 0, 0) ) start_y = max(0, -y_min) end_y = start_y + h2 start_x = max(0, -x_min) end_x = start_x + w2 end_y = min(end_y, canvas_h) end_x = min(end_x, canvas_w) img2_crop_h = end_y - start_y img2_crop_w = end_x - start_x if img2_crop_h > 0 and img2_crop_w > 0: est_img2[start_y:end_y, start_x:end_x] = img_2_uint8[:img2_crop_h, :img2_crop_w] alpha1 = cv2.remap( np.ones((h1, w1), dtype=np.float32), trans_x, trans_y, cv2.INTER_CUBIC, borderMode=cv2.BORDER_CONSTANT, borderValue=0 ) alpha2 = np.zeros_like(alpha1) alpha2[start_y:end_y, start_x:end_x] = 1.0 alpha_sum = alpha1 + alpha2 alpha_sum[alpha_sum == 0] = 1.0 alpha1_norm = alpha1 / alpha_sum alpha2_norm = alpha2 / alpha_sum alpha1_3ch = np.stack([alpha1_norm, alpha1_norm, alpha1_norm], axis=-1) alpha2_3ch = np.stack([alpha2_norm, alpha2_norm, alpha2_norm], axis=-1) est_img = (est_img1.astype(np.float32) * alpha1_3ch + est_img2.astype(np.float32) * alpha2_3ch) est_img = np.clip(est_img, 0, 255).astype(np.uint8) return est_img
def generate_panorama(ordered_img_seq):
length = len(ordered_img_seq)
mid = length // 2
principle_img = ordered_img_seq[mid].copy()
for j in range(mid + 1, length): try: pixels1, pixels2 = feature_matching(ordered_img_seq[j], principle_img) if len(pixels1) < 4: print(f"匹配点不足，跳过图像 {j}") continue homo_matrix = align_pair(pixels1, pixels2) principle_img = stitch_blend(ordered_img_seq[j], principle_img, homo_matrix) principle_img = np.uint8(principle_img) except Exception as e: print(f"右侧拼接失败 {j}: {e}") for i in range(mid - 1, -1, -1): try: pixels1, pixels2 = feature_matching(ordered_img_seq[i], principle_img) if len(pixels1) < 4: print(f"匹配点不足，跳过图像 {i}") continue homo_matrix = align_pair(pixels1, pixels2) principle_img = stitch_blend(ordered_img_seq[i], principle_img, homo_matrix) principle_img = np.uint8(principle_img) except Exception as e: print(f"左侧拼接失败 {i}: {e}") return principle_img
test_matching()
blend_pairs = [
(os.path.join(IMGDIR, “1_1.jpg”), os.path.join(IMGDIR, “1_2.jpg”)),
(os.path.join(IMGDIR, “2_1.jpg”), os.path.join(IMGDIR, “2_2.jpg”)),
(os.path.join(IMGDIR, “3_1.jpg”), os.path.join(IMGDIR, “3_2.jpg”))
]
for idx, (path1, path2) in enumerate(blend_pairs, 1):
img1 = cv2.imread(path1)
img2 = cv2.imread(path2)
if img1 is None:
print(f"无法加载图像1: {path1}“)
continue
if img2 is None:
print(f"无法加载图像2: {path2}”)
continue
try:
pixels1, pixels2 = feature_matching(img1, img2)
homo = align_pair(pixels1, pixels2)
blended_img = stitch_blend(img1, img2, homo)
save_path = os.path.join(OUTPUT_DIR, f"blend_{idx}.png")
cv2.imwrite(save_path, blended_img)
print(f"成功保存: {save_path}“)
except Exception as e:
print(f"混合组{idx}失败：{e}”)
panorama_tasks = [
{
“name”: “grail”,
“files”: [f"grail{i:02d}.jpg" for i in range(5)],
“output”: “panorama_1.png”
},
{
“name”: “library”,
“files”: [f"{i}.jpg" for i in range(10, 15)],
“output”: “panorama_2.png”
},
{
“name”: “parrington”,
“files”: [f"prtn{i:02d}.jpg" for i in range(5)],
“output”: “panorama_3.png”
},
{
“name”: “Xue-Mountain-Entrance”,
“files”: [f"DSC_017{i}.jpg" for i in range(4, 9)],
“output”: “panorama_4.png”
}
]
for task in panorama_tasks:
img_list = []
for file in task[“files”]:
img_path = os.path.join(IMGDIR, “panoramas”, task[“name”], file)
img = cv2.imread(img_path)
if img is None:
print(f"无法加载图像: {img_path}“)
continue
img_list.append(img)
if len(img_list) < 2:
print(f"图像数量不足，跳过任务: {task[‘name’]}”)
continue
try:
pano = generate_panorama(img_list)
if pano is not None:
save_path = os.path.join(OUTPUT_DIR, task[“output”])
cv2.imwrite(save_path, pano)
print(f"成功保存全景图: {save_path}“)
else:
print(f"全景图生成失败: {task[‘name’]}”)
except Exception as e:
print(f"景任务失败 {task[‘name’]}: {e}")
对的吗？
Mac'sgraybox样式实现
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
