把下面的代码从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)