numpy-np.ix_函数的使用

Numpy中np.ix_函数用法解析
博客围绕Numpy的np.ix_函数展开,给出官方文档链接,虽表示对用法和操作代码不太懂,但通过参考链接分析得出,该函数第一个参数返回行索引,第二个返回列索引等。还举例说明了二维和三维array使用该函数获取元素的情况。

官方文档:https://numpy.org/doc/stable/reference/generated/numpy.ix_.html#numpy.ix_,给出的用法是:

numpy.ix_(*args)
'''
Construct an open mesh from multiple sequences.

This function takes N 1-D sequences and returns N outputs with N dimensions each, such that the shape is 1 in all but one dimension and the dimension with the non-unit shape value cycles through all N dimensions.

Using ix_ one can quickly construct index arrays that will index the cross product. a[np.ix_([1,3],[2,5])] returns the array [[a[1,2] a[1,5]], [a[3,2] a[3,5]]].
'''

看不太懂,操作代码也是不懂:

a = np.arange(10).reshape(2, 5)
print(a)
ixgrid = np.ix_([0, 1], [2, 4])
print("-" * 10)
print(ixgrid)
print("-" * 10)
print(a[ixgrid])
'''
输出:
[[0 1 2 3 4]
 [5 6 7 8 9]]
----------
(array([[0],
       [1]]), array([[2, 4]]))
----------
[[2 4]
 [7 9]]
'''

参考链接中给出的:

'''
Its main use is to create an open mesh to select specific indices from an array (specific sub-array). 
An easy example to understand it is: Say you have an array of shape (5,5), 
and you would like to select the sub-array that is constructed from selecting 
rows 1 and 3 and columns 0 and 3. 
You can use np.ix_ to create such (index) mesh to be able to select sub-array as follows:
'''
a = np.arange(5*5).reshape(5,5)
print(a)
print("-" * 10)
sub_indices = np.ix_([1,3],[0,3])
print(sub_indices)
print("-" * 10)
print(a[sub_indices])
'''
输出:
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]
 [15 16 17 18 19]
 [20 21 22 23 24]]
----------
(array([[1],
       [3]]), array([[0, 3]]))
----------
[[ 5  8]
 [15 18]]
'''

通过上面的操作,得到的是'a'中行索引为1、3和列索引为0、3的元素:即'a[1, 0], a[1, 3], a[3, 0], a[3, 3]'。

col 0    col 3
   |        |
   v        v
[[ 0  1  2  3  4]   
 [ 5  6  7  8  9]   <- row 1
 [10 11 12 13 14]
 [15 16 17 18 19]   <- row 3
 [20 21 22 23 24]]

从np.ix_返回的结果看,函数第一个参数返回‘行的索引’,第二个参数返回‘列的索引’,第三个参数返回‘深度的索引’,以此类推。

当索引三维array时:

b = np.arange(5*5*5).reshape(5, 5, 5)
print(b)
print("-" * 10)
sub_indices = np.ix_([1,3],[0,3],[1,3])
print(sub_indices)
print("-" * 10)
print(b[sub_indices])
'''
输出:
[[[  0   1   2   3   4]
  [  5   6   7   8   9]
  [ 10  11  12  13  14]
  [ 15  16  17  18  19]
  [ 20  21  22  23  24]]

 [[ 25  26  27  28  29]
  [ 30  31  32  33  34]
  [ 35  36  37  38  39]
  [ 40  41  42  43  44]
  [ 45  46  47  48  49]]

 [[ 50  51  52  53  54]
  [ 55  56  57  58  59]
  [ 60  61  62  63  64]
  [ 65  66  67  68  69]
  [ 70  71  72  73  74]]

 [[ 75  76  77  78  79]
  [ 80  81  82  83  84]
  [ 85  86  87  88  89]
  [ 90  91  92  93  94]
  [ 95  96  97  98  99]]

 [[100 101 102 103 104]
  [105 106 107 108 109]
  [110 111 112 113 114]
  [115 116 117 118 119]
  [120 121 122 123 124]]]
----------
(array([[[1]],

       [[3]]]), array([[[0],
        [3]]]), array([[[1, 3]]]))
----------
[[[26 28]
  [41 43]]

 [[76 78]
  [91 93]]]
'''

通过上面的操作,得到的是'b'中'b[1, 0, 1], b[1, 0, 3], b[1, 3, 1], b[1, 3, 3],b[3, 0, 1], b[3, 0, 3], b[3, 3, 1], b[3, 3, 3]'。

 

参考:

https://stackoverflow.com/questions/62505046/what-does-numpy-ix-function-do-and-what-is-the-output-used-for

 

import cv2 import numpy as np import math class CDomPho: def init(self): “”“初始化摄影测量参数”“” self.m_f = self.m_Pxz = self.m_k1 = self.m_k2 = self.m_p1 = self.m_p2 = 0.0 self.m_X0 = self.m_Y0 = 0.0 self.m_cols = self.m_rows = 0 self.m_aopX = self.m_aopY = self.m_aopZ = self.m_aopP = self.m_aopO = self.m_aopK = 0.0 self.m_rotM = np.zeros(9) # 旋转矩阵 self.m_ref = np.zeros(6) # 参考参数 def POK2M(self, p, o, k, r): """将俯仰角(p)、翻滚角(o)、偏航角(k)转换为旋转矩阵""" cosp = math.cos(p) cosw = math.cos(o) cosk = math.cos(k) sinp = math.sin(p) sinw = math.sin(o) sink = math.sin(k) # 计算旋转矩阵元素 r[0] = cosp * cosk - sinp * sinw * sink r[1] = -cosp * sink - sinp * sinw * cosk r[2] = -sinp * cosw r[3] = cosw * sink r[4] = cosw * cosk r[5] = -sinw r[6] = sinp * cosk + cosp * sinw * sink r[7] = -sinp * sink + cosp * sinw * cosk r[8] = cosp * cosw def SetPara(self, aopX, aopY, aopZ, aopP, aopO, aopK, f, pxz, cols, rows, k1, k2, p1, p2, X0, Y0): """设置摄影测量参数""" self.m_rotM = np.zeros(9) self.POK2M(aopP, aopO, aopK, self.m_rotM) # 计算旋转矩阵 # 存储参数 self.m_aopX, self.m_aopY, self.m_aopZ = aopX, aopY, aopZ self.m_aopP, self.m_aopO, self.m_aopK = aopP, aopO, aopK self.m_f, self.m_Pxz = f, pxz self.m_k1, self.m_k2, self.m_p1, self.m_p2 = k1, k2, p1, p2 self.m_X0, self.m_Y0 = X0, Y0 self.m_cols, self.m_rows = cols, rows # 设置参考参数 self.m_ref[0] = -(cols / 2) * pxz self.m_ref[1] = pxz self.m_ref[2] = 0 self.m_ref[3] = -(rows / 2) * pxz self.m_ref[4] = 0 self.m_ref[5] = -pxz def Scan2Pho(self, scX, scY): """从扫描坐标转换到摄影坐标""" phoX = self.m_ref[0] + self.m_ref[1] * scX - self.m_ref[2] * scY phoY = self.m_ref[3] + self.m_ref[4] * scX - self.m_ref[5] * scY # 应用畸变校正 if abs(self.m_k1) > 1e-12: x, y = phoX, phoY r2 = x * x + y * y r4 = r2 * r2 dx = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y dy = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y phoX -= dx phoY -= dy return phoX, phoY def Pho2Scan(self, phoX, phoY): """从摄影坐标转换到扫描坐标""" # 应用畸变校正 if abs(self.m_k1) > 1e-12: x, y = phoX, phoY r2 = x * x + y * y r4 = r2 * r2 dx0 = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y dy0 = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y x = phoX + dx0 y = phoY + dy0 r2 = x * x + y * y r4 = r2 * r2 dx = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y dy = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y phoX += dx phoY += dy phoX -= self.m_ref[0] phoY -= self.m_ref[3] tt = self.m_ref[1] * self.m_ref[5] - self.m_ref[2] * self.m_ref[4] scX = (self.m_ref[5] * phoX - self.m_ref[2] * phoY) / tt scY = -(-self.m_ref[4] * phoX + self.m_ref[1] * phoY) / tt return scX, scY def PhoZ2Grd(self, px, py, gz): """从摄影坐标和高程转换到地面坐标""" A = self.m_rotM[0] * px + self.m_rotM[1] * py - self.m_rotM[2] * self.m_f B = self.m_rotM[3] * px + self.m_rotM[4] * py - self.m_rotM[5] * self.m_f C = self.m_rotM[6] * px + self.m_rotM[7] * py - self.m_rotM[8] * self.m_f if abs\(C\) > 1e-12: N = (gz - self.m_aopZ) / C else: N = 0.0 gx = self.m_aopX + A * N gy = self.m_aopY + B * N return gx, gy def Grd2Pho_(self, gx, gy, gz): """从地面坐标转换到摄影坐标""" dx = gx - self.m_aopX dy = gy - self.m_aopY dz = gz - self.m_aopZ denom = self.m_rotM[2] * dx + self.m_rotM[5] * dy + self.m_rotM[8] * dz fm = -self.m_f / denom if abs(denom) > 1e-12 else 0 px = (self.m_rotM[0] * dx + self.m_rotM[3] * dy + self.m_rotM[6] * dz) * fm py = (self.m_rotM[1] * dx + self.m_rotM[4] * dy + self.m_rotM[7] * dz) * fm return px, py def ImgZ2Grd(self, xyz): """从图像坐标和高程转换到地面坐标""" results = [] for i in range(0, len(xyz), 3): scX, scY, Z = xyz[i], xyz[i + 1], xyz[i + 2] phoX, phoY = self.Scan2Pho(scX, scY) gx, gy = self.PhoZ2Grd(phoX, phoY, Z) results.extend([gx, gy, Z]) return results def Grd2Img(self, gx, gy, gz): """从地面坐标转换到图像坐标""" px, py = self.Grd2Pho_(gx, gy, gz) scX, scY = self.Pho2Scan(px, py) return scX, scY 全局变量用于图像交互 scale = 1.0 angle = 0.0 offset_x = 0 offset_y = 0 last_mouse_pos = (0, 0) dragging = False def mouse_callback(event, x, y, flags, param): “”“鼠标回调函数,处理图像交互”“” global scale, angle, offset_x, offset_y, dragging, last_mouse_pos if event == cv2.EVENT_MBUTTONDOWN: # 鼠标中键按下 dragging = True last_mouse_pos = (x, y) elif event == cv2.EVENT_MOUSEMOVE and dragging: # 拖动图像 dx = x - last_mouse_pos[0] dy = y - last_mouse_pos[1] offset_x += dx offset_y += dy last_mouse_pos = (x, y) update_display() elif event == cv2.EVENT_MBUTTONUP: # 鼠标中键释放 dragging = False elif event == cv2.EVENT_MOUSEWHEEL: # 鼠标滚轮缩放 if flags > 0: scale *= 1.1 # 放大 else: scale /= 1.1 # 缩小 update_display() elif event == cv2.EVENT_RBUTTONDOWN: # 鼠标右键旋转 angle += 10 if angle >= 360: angle -= 360 update_display() def update_display(): “”“根据当前变换参数更新显示图像”“” global domImg, scale, angle, offset_x, offset_y height, width = domImg.shape[:2] center = (width // 2, height // 2) # 创建变换矩阵 M_scale = cv2.getRotationMatrix2D(center, 0, scale) M_rot = cv2.getRotationMatrix2D(center, angle, 1) M_translate = np.float32([[1, 0, offset_x], [0, 1, offset_y]]) # 组合变换:先缩放,再旋转,最后平移 M_combined = np.vstack([M_translate, [0, 0, 1]]) @ np.vstack([M_rot, [0, 0, 1]]) @ np.vstack([M_scale, [0, 0, 1]]) M_final = M_combined[:2, :] # 转换为2x3矩阵用于warpAffine # 应用变换并显示 transformed_img = cv2.warpAffine(domImg, M_final, (width, height)) cv2.imshow('Orthophoto', transformed_img) if name == “main”: # 读取摄影测量参数 try: with open(r’D:\shuiwei\para.txt’) as file: file.readline() # 跳过标题行 params = list(map(float, file.readline().split())) except Exception as e: print(f"读取参数文件时出错: {e}") exit(1) # 解析参数 aopX, aopY, aopZ, aopP, aopO, aopK = params[0:6] f, Pxz, k1, k2, p1, p2, X0, Y0 = params[6:14] # 读取原始图像 m_img = cv2.imread(r'D:\shuiwei\8.jpg') if m_img is None: raise FileNotFoundError("无法读取图像文件") rows, cols = m_img.shape[:2] # 初始化摄影测量模型 dompho = CDomPho() dompho.SetPara(aopX, aopY, aopZ, aopP, aopO, aopK, f, Pxz, cols, rows, k1, k2, p1, p2, X0, Y0) # 定义图像四个角点及其高程 corPt = [0.0, 0.0, 64.3, cols, 0, 64.3, 0, rows, 64.3, cols, rows, 64.3] # 计算地面坐标范围 grd_coords = dompho.ImgZ2Grd(corPt) gxs = grd_coords[0::3] gys = grd_coords[1::3] mingx, maxgx = min(gxs), max(gxs) mingy, maxgy = min(gys), max(gys) dx = maxgx - mingx dy = maxgy - mingy # 创建正射影像 resolution = 0.5 # 地面分辨率(米/像素) domCols = int(dx / resolution) + 1 domRows = int(dy / resolution) + 1 domImg = np.zeros((domRows, domCols, 3), dtype=np.uint8) h = 64.3 # 高程 gy = mingy # 起始Y坐标 total_pixels = domRows * domCols processed = 0 # 生成正射影像(修改了镜像问题) for r in range(domRows): gx = mingx for c in range(domCols): # 计算图像坐标 ix, iy = dompho.Grd2Img(gx, gy, h) i = int(ix + 0.5) # 列坐标 j = int(iy + 0.5) # 行坐标 # 修正:直接使用j坐标,避免垂直翻转 if 0 <= i < cols and 0 <= j < rows: domImg[r, c] = m_img[j, i] # 关键修改:使用j而不是rows-1-j gx += resolution processed += 1 # 显示进度 if processed % 10000 == 0: percent = processed / total_pixels * 100 print(f"处理进度: {percent:.2f}%") gy += resolution # 创建交互窗口 cv2.namedWindow('Orthophoto', cv2.WINDOW_NORMAL) cv2.setMouseCallback('Orthophoto', mouse_callback) update_display() print("操作说明:") print("- 鼠标中键拖动:平移图像") print("- 鼠标滚轮:缩放图像") print("- 鼠标右键:旋转图像") print("- ESC键:退出程序") # 主循环 while True: key = cv2.waitKey(1) if key == 27: # ESC键退出 break cv2.destroyAllWindows() 再正射校正之前,帮我添加一个选择roi感兴趣区域的代码,这样我可以直接选择感兴趣区域作为我的校正的部分,其他部分就不考虑了
06-27
import cv2 import numpy as np import math from tqdm import tqdm # 添加进度条库 class CDomPho: def __init__(self): """初始化摄影测量参数""" self.m_f = self.m_Pxz = self.m_k1 = self.m_k2 = self.m_p1 = self.m_p2 = 0.0 self.m_X0 = self.m_Y0 = 0.0 self.m_cols = self.m_rows = 0 self.m_aopX = self.m_aopY = self.m_aopZ = self.m_aopP = self.m_aopO = self.m_aopK = 0.0 self.m_rotM = np.zeros(9) # 旋转矩阵 self.m_ref = np.zeros(6) # 参考参数 def POK2M(self, p, o, k, r): """将俯仰角(p)、翻滚角(o)、偏航角(k)转换为旋转矩阵""" cosp = math.cos(p) cosw = math.cos(o) cosk = math.cos(k) sinp = math.sin(p) sinw = math.sin(o) sink = math.sin(k) # 计算旋转矩阵元素 r[0] = cosp * cosk - sinp * sinw * sink r[1] = -cosp * sink - sinp * sinw * cosk r[2] = -sinp * cosw r[3] = cosw * sink r[4] = cosw * cosk r[5] = -sinw r[6] = sinp * cosk + cosp * sinw * sink r[7] = -sinp * sink + cosp * sinw * cosk r[8] = cosp * cosw def SetPara(self, aopX, aopY, aopZ, aopP, aopO, aopK, f, pxz, cols, rows, k1, k2, p1, p2, X0, Y0): """设置摄影测量参数""" self.m_rotM = np.zeros(9) self.POK2M(aopP, aopO, aopK, self.m_rotM) # 计算旋转矩阵 # 存储参数 self.m_aopX, self.m_aopY, self.m_aopZ = aopX, aopY, aopZ self.m_aopP, self.m_aopO, self.m_aopK = aopP, aopO, aopK self.m_f, self.m_Pxz = f, pxz self.m_k1, self.m_k2, self.m_p1, self.m_p2 = k1, k2, p1, p2 self.m_X0, self.m_Y0 = X0, Y0 self.m_cols, self.m_rows = cols, rows # 设置参考参数 self.m_ref[0] = -(cols / 2) * pxz self.m_ref[1] = pxz self.m_ref[2] = 0 self.m_ref[3] = -(rows / 2) * pxz self.m_ref[4] = 0 self.m_ref[5] = -pxz def Scan2Pho(self, scX, scY): """从扫描坐标转换到摄影坐标""" phoX = self.m_ref[0] + self.m_ref[1] * scX - self.m_ref[2] * scY phoY = self.m_ref[3] + self.m_ref[4] * scX - self.m_ref[5] * scY # 应用畸变校正 if abs(self.m_k1) > 1e-12: x, y = phoX, phoY r2 = x * x + y * y r4 = r2 * r2 dx = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y dy = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y phoX -= dx phoY -= dy return phoX, phoY def Pho2Scan(self, phoX, phoY): """从摄影坐标转换到扫描坐标""" # 应用畸变校正 if abs(self.m_k1) > 1e-12: x, y = phoX, phoY r2 = x * x + y * y r4 = r2 * r2 dx0 = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y dy0 = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y x = phoX + dx0 y = phoY + dy0 r2 = x * x + y * y r4 = r2 * r2 dx = x * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p1 * (r2 + 2 * x * x) + 2 * self.m_p2 * x * y dy = y * (self.m_k1 * r2 + self.m_k2 * r4) + self.m_p2 * (r2 + 2 * y * y) + 2 * self.m_p1 * x * y phoX += dx phoY += dy phoX -= self.m_ref[0] phoY -= self.m_ref[3] tt = self.m_ref[1] * self.m_ref[5] - self.m_ref[2] * self.m_ref[4] scX = (self.m_ref[5] * phoX - self.m_ref[2] * phoY) / tt scY = -(-self.m_ref[4] * phoX + self.m_ref[1] * phoY) / tt return scX, scY def PhoZ2Grd(self, px, py, gz): """从摄影坐标和高程转换到地面坐标""" A = self.m_rotM[0] * px + self.m_rotM[1] * py - self.m_rotM[2] * self.m_f B = self.m_rotM[3] * px + self.m_rotM[4] * py - self.m_rotM[5] * self.m_f C = self.m_rotM[6] * px + self.m_rotM[7] * py - self.m_rotM[8] * self.m_f if abs(C) > 1e-12: N = (gz - self.m_aopZ) / C else: N = 0.0 gx = self.m_aopX + A * N gy = self.m_aopY + B * N return gx, gy def Grd2Pho_(self, gx, gy, gz): """从地面坐标转换到摄影坐标 - 修复数组处理问题""" # 确保输入是标量 if isinstance(gx, np.ndarray) or isinstance(gy, np.ndarray) or isinstance(gz, np.ndarray): # 处理数组输入 px = np.zeros_like(gx) py = np.zeros_like(gy) for i in range(len(gx)): dx = gx[i] - self.m_aopX dy = gy[i] - self.m_aopY dz = gz[i] - self.m_aopZ denom = self.m_rotM[2] * dx + self.m_rotM[5] * dy + self.m_rotM[8] * dz if abs(denom) > 1e-12: fm = -self.m_f / denom else: fm = 0.0 px[i] = (self.m_rotM[0] * dx + self.m_rotM[3] * dy + self.m_rotM[6] * dz) * fm py[i] = (self.m_rotM[1] * dx + self.m_rotM[4] * dy + self.m_rotM[7] * dz) * fm return px, py else: # 处理标量输入 dx = gx - self.m_aopX dy = gy - self.m_aopY dz = gz - self.m_aopZ denom = self.m_rotM[2] * dx + self.m_rotM[5] * dy + self.m_rotM[8] * dz if abs(denom) > 1e-12: fm = -self.m_f / denom else: fm = 0.0 px = (self.m_rotM[0] * dx + self.m_rotM[3] * dy + self.m_rotM[6] * dz) * fm py = (self.m_rotM[1] * dx + self.m_rotM[4] * dy + self.m_rotM[7] * dz) * fm return px, py def ImgZ2Grd(self, xyz): """从图像坐标和高程转换到地面坐标""" results = [] for i in range(0, len(xyz), 3): scX, scY, Z = xyz[i], xyz[i + 1], xyz[i + 2] phoX, phoY = self.Scan2Pho(scX, scY) gx, gy = self.PhoZ2Grd(phoX, phoY, Z) results.extend([gx, gy, Z]) return results def Grd2Img(self, gx, gy, gz): """从地面坐标转换到图像坐标 - 修复数组处理问题""" # 确保输入是标量 if isinstance(gx, np.ndarray) or isinstance(gy, np.ndarray) or isinstance(gz, np.ndarray): # 处理数组输入 scX = np.zeros_like(gx) scY = np.zeros_like(gy) for i in range(len(gx)): px, py = self.Grd2Pho_(gx[i], gy[i], gz[i]) scX[i], scY[i] = self.Pho2Scan(px, py) return scX, scY else: # 处理标量输入 px, py = self.Grd2Pho_(gx, gy, gz) scX, scY = self.Pho2Scan(px, py) return scX, scY # 全局变量用于图像交互 scale = 1.0 angle = 0.0 offset_x = 0 offset_y = 0 last_mouse_pos = (0, 0) dragging = False domImg = None # 全局正射影像 def mouse_callback(event, x, y, flags, param): """鼠标回调函数,处理图像交互""" global scale, angle, offset_x, offset_y, dragging, last_mouse_pos if event == cv2.EVENT_MBUTTONDOWN: # 鼠标中键按下 dragging = True last_mouse_pos = (x, y) elif event == cv2.EVENT_MOUSEMOVE and dragging: # 拖动图像 dx = x - last_mouse_pos[0] dy = y - last_mouse_pos[1] offset_x += dx offset_y += dy last_mouse_pos = (x, y) update_display() elif event == cv2.EVENT_MBUTTONUP: # 鼠标中键释放 dragging = False elif event == cv2.EVENT_MOUSEWHEEL: # 鼠标滚轮缩放 if flags > 0: scale *= 1.1 # 放大 else: scale /= 1.1 # 缩小 update_display() elif event == cv2.EVENT_RBUTTONDOWN: # 鼠标右键旋转 angle += 10 if angle >= 360: angle -= 360 update_display() elif event == cv2.EVENT_LBUTTONDOWN: # 左键点击显示坐标 print(f"点击位置: x={x}, y={y}") def update_display(): """根据当前变换参数更新显示图像""" global domImg, scale, angle, offset_x, offset_y if domImg is None: return height, width = domImg.shape[:2] center = (width // 2, height // 2) # 创建变换矩阵 M_rot = cv2.getRotationMatrix2D(center, angle, scale) M_translate = np.float32([[1, 0, offset_x], [0, 1, offset_y]]) # 组合变换:先旋转缩放,再平移 M_combined = np.vstack([M_translate, [0, 0, 1]]) @ np.vstack([M_rot, [0, 0, 1]]) M_final = M_combined[:2, :] # 转换为2x3矩阵用于warpAffine # 应用变换并显示 transformed_img = cv2.warpAffine(domImg, M_final, (width, height)) # 添加坐标信息 cv2.putText(transformed_img, f"Scale: {scale:.2f} | Angle: {angle}°", (10, 30), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 0), 2) cv2.putText(transformed_img, f"Offset: ({offset_x}, {offset_y})", (10, 60), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 0), 2) cv2.imshow('Orthophoto', transformed_img) def select_roi_custom(image): """自定义ROI选择函数,添加了辅助线""" print("请用鼠标选择感兴趣区域,按空格键或回车确认,按ESC取消") clone = image.copy() roi_selected = False roi_box = (0, 0, 0, 0) drawing = False ix, iy = -1, -1 def mouse_callback(event, x, y, flags, param): nonlocal ix, iy, drawing, clone, roi_box, roi_selected if event == cv2.EVENT_LBUTTONDOWN: drawing = True ix, iy = x, y clone = image.copy() elif event == cv2.EVENT_MOUSEMOVE: if drawing: clone = image.copy() # 绘制矩形 cv2.rectangle(clone, (ix, iy), (x, y), (0, 255, 0), 2) # 绘制中心十字线 center_x = (ix + x) // 2 center_y = (iy + y) // 2 cv2.line(clone, (center_x, 0), (center_x, clone.shape[0]), (0, 255, 255), 1) cv2.line(clone, (0, center_y), (clone.shape[1], center_y), (0, 255, 255), 1) # 显示尺寸信息 width = abs(ix - x) height = abs(iy - y) cv2.putText(clone, f"{width}x{height}", (x + 10, y + 20), cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 255, 255), 1) cv2.imshow('Select ROI', clone) elif event == cv2.EVENT_LBUTTONUP: drawing = False roi_box = (min(ix, x), min(iy, y), abs(ix - x), abs(iy - y)) clone = image.copy() cv2.rectangle(clone, (roi_box[0], roi_box[1]), (roi_box[0] + roi_box[2], roi_box[1] + roi_box[3]), (0, 255, 0), 2) cv2.imshow('Select ROI', clone) roi_selected = True cv2.namedWindow('Select ROI', cv2.WINDOW_NORMAL) cv2.setMouseCallback('Select ROI', mouse_callback) cv2.imshow('Select ROI', image) while True: key = cv2.waitKey(1) & 0xFF if key == 13 or key == 32: # 回车或空格确认 break elif key == 27: # ESC取消 roi_box = (0, 0, 0, 0) break cv2.destroyWindow('Select ROI') if roi_selected: x, y, w, h = roi_box cropped = image[y:y + h, x:x + w] print(f"已选择ROI: x={x}, y={y}, width={w}, height={h}") return cropped, (x, y, w, h) else: print("未选择ROI,将使用整张图像") return image, (0, 0, image.shape[1], image.shape[0]) def generate_dom(dompho, m_img, mingx, maxgx, mingy, maxgy, resolution, h=64.3): """生成正射影像的优化函数""" dx = maxgx - mingx dy = maxgy - mingy domCols = int(dx / resolution) + 1 domRows = int(dy / resolution) + 1 # 创建空图像 domImg = np.zeros((domRows, domCols, 3), dtype=np.uint8) # 获取ROI尺寸 rows, cols = m_img.shape[:2] print(f"生成正射影像: {domCols}x{domRows} 像素 (分辨率: {resolution}米/像素)") # 使用进度条 progress_bar = tqdm(total=domRows, desc="生成DOM") for r in range(domRows): gy = mingy + r * resolution for c in range(domCols): gx = mingx + c * resolution # 计算图像坐标 ix, iy = dompho.Grd2Img(gx, gy, h) i = int(ix + 0.5) # 列坐标 j = int(iy + 0.5) # 行坐标 # 确保坐标在图像范围内 if 0 <= i < cols and 0 <= j < rows: domImg[r, c] = m_img[j, i] # 使用ROI内的图像数据 progress_bar.update(1) progress_bar.close() return domImg if __name__ == "__main__": # 读取摄影测量参数 try: with open(r'D:\shuiwei\para.txt') as file: file.readline() # 跳过标题行 params = list(map(float, file.readline().split())) except Exception as e: print(f"读取参数文件时出错: {e}") exit(1) # 解析参数 aopX, aopY, aopZ, aopP, aopO, aopK = params[0:6] f, Pxz, k1, k2, p1, p2, X0, Y0 = params[6:14] # 读取原始图像 original_img = cv2.imread(r'D:\shuiwei\8.jpg') if original_img is None: raise FileNotFoundError("无法读取图像文件") # 使用自定义的ROI选择函数 m_img, roi_coords = select_roi_custom(original_img) x0, y0, w, h = roi_coords # 调整像主点偏移 - 关键修正! X0_adjusted = X0 - x0 Y0_adjusted = Y0 - y0 # 获取ROI后的图像尺寸 rows, cols = m_img.shape[:2] # 初始化摄影测量模型 dompho = CDomPho() dompho.SetPara(aopX, aopY, aopZ, aopP, aopO, aopK, f, Pxz, cols, rows, k1, k2, p1, p2, X0_adjusted, Y0_adjusted) # 使用调整后的像主点 # 定义图像四个角点及其高程(基于ROI) corPt = [ 0.0, 0.0, 64.3, # ROI左上角 cols, 0, 64.3, # ROI右上角 0, rows, 64.3, # ROI左下角 cols, rows, 64.3 # ROI右下角 ] # 计算地面坐标范围 grd_coords = dompho.ImgZ2Grd(corPt) gxs = grd_coords[0::3] gys = grd_coords[1::3] mingx, maxgx = min(gxs), max(gxs) mingy, maxgy = min(gys), max(gys) # 创建正射影像 resolution = 3 # 地面分辨率(米/像素) print("开始生成正射影像...") # 使用优化函数生成DOM domImg = generate_dom(dompho, m_img, mingx, maxgx, mingy, maxgy, resolution) # 创建交互窗口 cv2.namedWindow('Orthophoto', cv2.WINDOW_NORMAL) cv2.setMouseCallback('Orthophoto', mouse_callback) update_display() print("操作说明:") print("- 鼠标中键拖动:平移图像") print("- 鼠标滚轮:缩放图像") print("- 鼠标右键:旋转图像") print("- 鼠标左键:显示点击位置坐标") print("- ESC键:退出程序") print("- S键:保存当前视图") # 主循环 while True: key = cv2.waitKey(1) if key == 27: # ESC键退出 break elif key == ord('s'): # 添加保存功能 cv2.imwrite(r'D:\shuiwei\current_view1.jpg', domImg) print("当前视图已保存") cv2.destroyAllWindows()检查一下这个代码,几个坐标系换来换去很容易出错,尤其是opencv图像坐标系的坐标原点在图像左上角,而扫描坐标系坐标原点在左下角
最新发布
07-23
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值