基于Opencv的矩阵三角函数封装 cvAtan2Mat cvSinMat cvAtan2Mat功能等于matlab里的atan2 sin cos...

本文介绍了在Opencv中如何封装缺失的三角函数操作,包括正弦(sin)、余弦(cos)及反正切(atan2)函数,这些函数能够直接对矩阵进行操作,极大提升了图像处理任务中的效率。

Opencv中没有提供对矩阵进行操作的三角函数,我封装了三个,其他三角函数可以仿照此进行封装。

头文件中要包含

#include <cmath>


CvMat* cvAtan2Mat(CvMat *a, CvMat *b)
{
int rows = a->rows;
int cols = a->cols;
CvMat *out = cvCreateMat(rows, cols, a->type);
for(int i=0; i<rows; i++)
{
float* ptra = ( float*)(a->data.ptr+i*a->step);
float* ptrb = ( float*)(b->data.ptr+i*b->step);
float* ptrout = ( float*)(out->data.ptr+i*out->step);
for(int j=0; j<cols; j++)
{
*ptrout = atan2(*ptra,*ptrb);
ptra++;
ptrb++;
ptrout++;
}
}
return out;
}


CvMat* cvSinMat(CvMat *a)
{
int rows = a->rows;
int cols = a->cols;
CvMat *out = cvCreateMat(rows, cols, a->type);
for(int i=0; i<rows; i++)
{
float* ptra = ( float*)(a->data.ptr+i*a->step);
float* ptrout = ( float*)(out->data.ptr+i*out->step);
for(int j=0; j<cols; j++)
{
*ptrout = sin(*ptra);
ptra++;
ptrout++;
}
}
return out;
}


CvMat* cvCosMat(CvMat *a)
{
int rows = a->rows;
int cols = a->cols;
CvMat *out = cvCreateMat(rows, cols, a->type);
for(int i=0; i<rows; i++)
{
float* ptra = ( float*)(a->data.ptr+i*a->step);
float* ptrout = ( float*)(out->data.ptr+i*out->step);
for(int j=0; j<cols; j++)
{
*ptrout = cos(*ptra);
ptra++;
ptrout++;
}
}
return out;
}

把下面的代码从python转成matlab代码 import os import cv2 import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import curve_fit # from sklearn.cluster import DBSCAN # import networkx as nx import warnings warnings.filterwarnings("ignore", category=RuntimeWarning) DRILLHOLE_COORDS = { '1#孔': (500, 2000, 0), '2#孔': (1500, 2000, 0), '3#孔': (2500, 2000, 0), '4#孔': (500, 1000, 0), '5#孔': (1500, 1000, 0), '6#孔': (2500, 1000, 0), } DRILLHOLE_DIAMETER = 94.25 def sine_model(x_angle, R, P_angle, phase_angle, C): return R * np.sin(x_angle / P_angle * 2 * np.pi + phase_angle) + C def process_image(image_path, depth_start, depth_end): try: img_array = np.fromfile(image_path, dtype=np.uint8) img = cv2.imdecode(img_array, cv2.IMREAD_GRAYSCALE) if img is None: raise ValueError("图像解码失败") except Exception as e: print(f"无法读取图像: {image_path}, 错误: {e}") return [] h, w = img.shape num_fractures = np.random.randint(0, 3) fractures = [] for _ in range(num_fractures): fractures.append({ 'R': np.random.uniform(5, 50), 'C': depth_start + np.random.uniform(0.1, 0.9) * (depth_end - depth_start), 'phase_deg': np.random.uniform(0, 360), 'source_image': os.path.basename(image_path), 'jrc': np.random.randint(5, 15) }) return fractures def calculate_3d_properties(fracture_info, drillhole_id): R = fracture_info['R'] C = fracture_info['C'] phase_deg = fracture_info['phase_deg'] hole_coords = DRILLHOLE_COORDS[drillhole_id] # 将所有单位统一为米 dip_rad = np.arctan(2 * R / (DRILLHOLE_DIAMETER)) # R和直径单位需一致 dip_deg = np.degrees(dip_rad) phase_rad = np.radians(phase_deg) dip_direction_deg = (3 * np.pi / 2 - phase_rad) * 180 / np.pi dip_direction_deg = dip_direction_deg % 360 # C是深度,假设孔口z=0, 向下为正 intersection_pt = np.array([hole_coords[0], hole_coords[1], C]) dip_dir_rad = np.radians(dip_direction_deg) # 法向量(x-东, y-北, z-下) normal_vector = np.array([ np.sin(dip_rad) * np.sin(dip_dir_rad), np.sin(dip_rad) * np.cos(dip_dir_rad), np.cos(dip_rad) ]) fracture_info.update({ 'drillhole_id': drillhole_id, 'dip': dip_deg, 'dip_direction': dip_direction_deg, 'intersection_pt': intersection_pt, # 这就是我们模型中的 C_i 'normal_vector': normal_vector, # 这就是我们模型中的 n_i # 'jrc' 属性已在process_image中添加 }) return fracture_info MODEL_PARAMS = { 'weights': {'dist': 0.5, 'orient': 0.4, 'jrc': 0.1}, 'sigma_dist': 1000, # 空间距离评分的尺度参数 (mm), 对应钻孔间距 'jrc_range': 20, 'connectivity_threshold': 0.6, # 判断为连通的概率阈值 'fracture_radius': 400 # 可视化半径 } def calculate_probabilistic_connectivity(all_fractures): print("\n--- Calculating probabilistic fracture connectivity ---") adj_pairs = [( '1#孔', '2#孔'), ('2#孔', '3#孔'), ('4#孔', '5#孔'), ('5#孔', '6#孔'), ( '1#孔', '4#孔'), ('2#孔', '5#孔'), ('3#孔', '6#孔')] connections = [] for id1, id2 in adj_pairs: fracs1 = [f for f in all_fractures if f['drillhole_id'] == id1] fracs2 = [f for f in all_fractures if f['drillhole_id'] == id2] for f1 in fracs1: for f2 in fracs2: # 1. 计算空间距离评分 S_dist dist = np.linalg.norm(f1['intersection_pt'] - f2['intersection_pt']) s_dist = np.exp(-(dist / MODEL_PARAMS['sigma_dist'])**2) # 2. 计算产状相似性评分 S_orient s_orient = np.abs(np.dot(f1['normal_vector'], f2['normal_vector'])) # 3. 计算JRC相似性评分 S_jrc s_jrc = max(0, 1 - (np.abs(f1['jrc'] - f2['jrc']) / MODEL_PARAMS['jrc_range'])) # 4. 加权计算总概率 w = MODEL_PARAMS['weights'] p_conn = w['dist']*s_dist + w['orient']*s_orient + w['jrc']*s_jrc if p_conn > MODEL_PARAMS['connectivity_threshold']: connections.append({ 'p': p_conn, 'c1': f1['intersection_pt'], 'c2': f2['intersection_pt'] }) print(f"Found {len(connections)} high-probability connections (P > {MODEL_PARAMS['connectivity_threshold']})") return connections def optimize_borehole_placement(num_new_boreholes=3): print("\n--- Optimizing supplementary borehole placement ---") boreholes_df = pd.DataFrame([ {'id': k, 'x_start': v[0], 'y_start': v[1], 'z_start': v[2]} for k, v in DRILLHOLE_COORDS.items() ]) print(" - 优先级 1 补充钻孔位置: (1500, 1500)") print(" - 优先级 2 补充钻孔位置: (1500, 0)") print(" - 优先级 3 补充钻孔位置: (3500, 1500)") pass def plot_disc(ax, center, normal, radius): pass def visualize_results(all_fractures, connections): print("\n--- Visualizing 3D probabilistic fracture network ---") fig = plt.figure(figsize=(12, 10)) ax = fig.add_subplot(111, projection='3d') # 绘制钻孔 (使用真实坐标) for hole_id, coords in DRILLHOLE_COORDS.items(): # 假设钻孔深度为7000mm ax.plot([coords[0], coords[0]], [coords[1], coords[1]], [coords[2], coords[2] + 7000], lw=2, color='black', label=f'Drillhole {hole_id}' if hole_id=='1#孔' else "") ax.text(coords[0], coords[1], coords[2], f' {hole_id}', color='red') for fracture in all_fractures: p = fracture['intersection_pt'] ax.scatter(p[0], p[1], p[2], color='blue', s=40) # plot_disc(ax, p, fracture['normal_vector'], MODEL_PARAMS['fracture_radius']) for conn in connections: c1, c2 = conn['c1'], conn['c2'] ax.plot([c1[0], c2[0]], [c1[1], c2[1]], [c1[2], c2[2]], color='red', ls='--', lw=1.5) ax.set_xlabel('East (mm)') ax.set_ylabel('North (mm)') ax.set_zlabel('Depth (mm)') ax.set_title('3D Probabilistic Fracture Network Model') if ax.get_zlim()[1] > ax.get_zlim()[0]: ax.invert_zaxis() ax.legend() plt.show() def run_analysis(root_path): all_fractures = [] drillhole_folders = { '1#孔': '1#孔', '2#孔': '2#孔', '3#孔': '3#孔', '4#孔': '4#孔', '5#孔': '5#孔', '6#孔': '6#孔' } for drillhole_id, folder_name in drillhole_folders.items(): drillhole_path = os.path.join(root_path, folder_name) if not os.path.isdir(drillhole_path): print(f"Warning: Directory not found for {drillhole_id}") continue print(f"--- Processing Drillhole: {drillhole_id} ---") for filename in sorted(os.listdir(drillhole_path)): if filename.lower().endswith(('.png', '.jpg', '.jpeg')): try: parts = os.path.splitext(filename)[0].replace('m', '').split('-') depth_start, depth_end = float(parts[0]) * 1000, float(parts[1]) * 1000 except (ValueError, IndexError): continue image_path = os.path.join(drillhole_path, filename) fractures_2d = process_image(image_path, depth_start, depth_end) for f in fractures_2d: f_3d = calculate_3d_properties(f, drillhole_id) all_fractures.append(f_3d) if not all_fractures: print("No fractures detected. Exiting.") return print(f"\nTotal fractures extracted from images: {len(all_fractures)}") prob_connections = calculate_probabilistic_connectivity(all_fractures) visualize_results(all_fractures, prob_connections) optimize_borehole_placement() if __name__ == '__main__': print("Starting integrated 3D fracture network analysis...") attachment_4_path = 'D:\hw杯数模题目\C题\附件4' if not os.path.exists(attachment_4_path): print(f"Error: The path '{attachment_4_path}' does not exist.") else: run_analysis(attachment_4_path)
09-25
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值