A. Points in Segments

本文介绍了一个算法问题:给定一系列区间,如何找出这些区间未覆盖的所有整数点。通过使用C++实现的示例代码,展示了如何初始化数组、标记区间内的点,并最终找出所有未被区间覆盖的点。

A. Points in Segments

time limit per test

1 second

memory limit per test

256 megabytes

input

standard input

output

standard output

You are given a set of nn segments on the axis OxOx, each segment has integer endpoints between 11 and mm inclusive. Segments may intersect, overlap or even coincide with each other. Each segment is characterized by two integers lili and riri (1≤li≤ri≤m1≤li≤ri≤m) — coordinates of the left and of the right endpoints.

Consider all integer points between 11 and mm inclusive. Your task is to print all such points that don't belong to any segment. The point xxbelongs to the segment [l;r][l;r] if and only if l≤x≤rl≤x≤r.

Input

The first line of the input contains two integers nn and mm (1≤n,m≤1001≤n,m≤100) — the number of segments and the upper bound for coordinates.

The next nn lines contain two integers each lili and riri (1≤li≤ri≤m1≤li≤ri≤m) — the endpoints of the ii-th segment. Segments may intersect, overlap or even coincide with each other. Note, it is possible that li=rili=ri, i.e. a segment can degenerate to a point.

Output

In the first line print one integer kk — the number of points that don't belong to any segment.

In the second line print exactly kk integers in any order — the points that don't belong to any segment. All points you print should be distinct.

If there are no such points at all, print a single integer 00 in the first line and either leave the second line empty or do not print it at all.

Examples

input

Copy

3 5
2 2
1 2
5 5

output

Copy

2
3 4 

input

Copy

1 7
1 7

output

Copy

0

Note

In the first example the point 11 belongs to the second segment, the point 22 belongs to the first and the second segments and the point 55belongs to the third segment. The points 33 and 44 do not belong to any segment.

In the second example all the points from 11 to 77 belong to the first segment.

代码:

#include <iostream>
#include <cstdio>

using namespace std;

int main()
{
    int a[105];
    fill(a,a+105,0);
    int n,m;
    scanf("%d%d",&n,&m);
    for(int i=0;i<n;++i)
    {
        int x,y;
        scanf("%d%d",&x,&y);
        fill(a+x,a+y+1,1);
    }
    int ans=0;
    for(int i=1;i<=m;++i)
    {
        if(a[i]==0) ans++;
    }
    cout<<ans<<endl;
    for(int i=1;i<=m;++i)
    {
        if(a[i]==0) cout<<i<<" ";
    }
    return 0;
}

 

import open3d as o3d import numpy as np import time import random from collections import deque import multiprocessing as mp import sys import traceback import ctypes def region_growing_segmentation(pcd, radius=0.02): """ 区域生长分割算法 :param pcd: 输入点云 :param radius: 邻域搜索半径(米) :return: 分割后的区域列表(每个区域是点索引列表) """ print("开始区域生长分割...") points = np.asarray(pcd.points) n_points = len(points) kdtree = o3d.geometry.KDTreeFlann(pcd) visited = np.zeros(n_points, dtype=bool) segments = [] for idx in range(n_points): if visited[idx]: continue # 创建新区域 segment = [] queue = deque([idx]) visited[idx] = True while queue: current_idx = queue.popleft() segment.append(current_idx) # 搜索半径内的邻域点 [_, idxs, _] = kdtree.search_radius_vector_3d(pcd.points[current_idx], radius) for n_idx in idxs: if not visited[n_idx]: visited[n_idx] = True queue.append(n_idx) segments.append(segment) print(f"分割完成,共 {len(segments)} 个区域") return segments def assign_random_colors(pcd, segments): """ 为每个区域分配随机颜色 :param pcd: 点云对象 :param segments: 分割区域列表 """ colors = np.zeros((len(pcd.points), 3)) for segment in segments: color = [random.random(), random.random(), random.random()] for idx in segment: colors[idx] = color pcd.colors = o3d.utility.Vector3dVector(colors) def process_segment(args): """ 处理单个区域的函数,用于多进程执行 :param args: 包含(区域索引, 区域点索引列表, 共享数组, 点云形状, 法向量半径, 计数半径)的元组 :return: (区域索引, 有效点标记数组, 是否移除整个区域标志, 处理时间, 错误信息) """ seg_idx, segment, shared_array, points_shape, radius_normal, radius_count = args start_time = time.time() error = None try: # 从共享数组重建点云数据 flat_array = np.frombuffer(shared_array.get_obj(), dtype=np.float64) all_points = flat_array.reshape(points_shape) # 创建点云对象 pcd = o3d.geometry.PointCloud() pcd.points = o3d.utility.Vector3dVector(all_points) # 创建KDTree kdtree = o3d.geometry.KDTreeFlann(pcd) # 创建该区域的点云 segment_pcd = pcd.select_by_index(segment) # 下采样 voxel_size = radius_normal / 2 down_pcd = segment_pcd.voxel_down_sample(voxel_size) if down_pcd is None or len(down_pcd.points) < 10: # 点数太少,跳过重建 valid_segment_points = np.ones(len(segment), dtype=bool) remove_entire_segment = len(segment) < 100 end_time = time.time() return seg_idx, segment, valid_segment_points, remove_entire_segment, end_time - start_time, error # 估计法向量 down_pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid( radius=radius_normal*2, max_nn=30)) # 泊松重建 - 使用更稳定的参数 if len(down_pcd.points) > 50: # 确保有足够点进行重建 try: with o3d.utility.VerbosityContextManager(o3d.utility.VerbosityLevel.Debug): mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( down_pcd, depth=8, linear_fit=True, scale=1.1) except Exception as e: mesh = o3d.geometry.TriangleMesh() error = f"泊松重建失败: {str(e)}" else: mesh = o3d.geometry.TriangleMesh() # 计算区域平均法向量A2 segment_normals = np.asarray(down_pcd.normals) A2 = np.mean(segment_normals, axis=0) if np.linalg.norm(A2) > 1e-6: A2 /= np.linalg.norm(A2) # 步骤4: 计算每个点的平均法向量A1并过滤 valid_segment_points = np.ones(len(segment), dtype=bool) for j, orig_idx in enumerate(segment): # 找到该点在下采样后的点云中的对应点 [_, idxs, _] = kdtree.search_knn_vector_3d(pcd.points[orig_idx], 1) down_idx = idxs[0] if idxs else -1 # 如果找不到对应点,跳过 if down_idx == -1: continue # 获取以该点为顶点的所有三角形 triangles = [] if mesh.has_triangles(): for tri in mesh.triangles: if down_idx in tri: triangles.append(tri) if len(triangles) == 0: continue # 计算这些三角形的法向量平均值A1 tri_normals = [] mesh_vertices = np.asarray(mesh.vertices) for tri in triangles: # 三角形三个顶点 pts = [mesh_vertices[i] for i in tri] # 计算三角片法向量 v1 = pts[1] - pts[0] v2 = pts[2] - pts[0] normal = np.cross(v1, v2) norm_val = np.linalg.norm(normal) if norm_val > 1e-6: normal /= norm_val tri_normals.append(normal) if tri_normals: A1 = np.mean(tri_normals, axis=0) norm_val = np.linalg.norm(A1) if norm_val > 1e-6: A1 /= norm_val # 计算A1和A2的夹角(度) if np.linalg.norm(A1) > 1e-6 and np.linalg.norm(A2) > 1e-6: cos_angle = np.clip(np.dot(A1, A2), -1.0, 1.0) angle = np.degrees(np.arccos(cos_angle)) # 如果偏差大于10度,标记为无效点 if angle > 10: valid_segment_points[j] = False # 步骤5: 半径内点云计数过滤 for j, orig_idx in enumerate(segment): if not valid_segment_points[j]: continue [k, idxs, _] = kdtree.search_radius_vector_3d(pcd.points[orig_idx], radius_count) # 排除自身点 neighbor_count = k - 1 if k > 0 else 0 # 如果邻域点少于100个,标记为无效点 if neighbor_count < 100: valid_segment_points[j] = False # 步骤6: 检查是否移除整个区域 remove_entire_segment = np.sum(valid_segment_points) < 100 end_time = time.time() process_time = end_time - start_time return seg_idx, segment, valid_segment_points, remove_entire_segment, process_time, error except Exception as e: # 捕获并记录异常 error = f"区域 {seg_idx} 处理错误: {str(e)}\n{traceback.format_exc()}" valid_segment_points = np.ones(len(segment), dtype=bool) remove_entire_segment = False end_time = time.time() return seg_idx, segment, valid_segment_points, remove_entire_segment, end_time - start_time, error def process_segments(pcd, segments, radius_normal=0.1, radius_count=0.1, max_workers=4): """ 多进程处理点云区域 :param pcd: 点云对象 :param segments: 分割区域列表 :param radius_normal: 法向量计算半径(米) :param radius_count: 点云计数半径(米) :param max_workers: 最大进程数 :return: 处理后的点云 """ print(f"开始多进程处理区域,共 {len(segments)} 个区域,使用 {max_workers} 个进程...") # 准备点云数据 points = np.asarray(pcd.points, dtype=np.float64) # 创建共享数组 shared_array = mp.Array(ctypes.c_double, points.size) # 将点云数据复制到共享数组 shared_np_array = np.frombuffer(shared_array.get_obj(), dtype=np.float64) shared_np_array[:] = points.flatten() # 创建有效点标记数组 valid_points = np.ones(len(pcd.points), dtype=bool) # 准备任务参数 tasks = [] for i, segment in enumerate(segments): tasks.append(( i, segment, shared_array, points.shape, radius_normal, radius_count )) # 用于存储结果的列表 results = [None] * len(segments) # 使用进程池并行处理 with mp.Pool(processes=max_workers) as pool: # 提交所有任务 futures = [pool.apply_async(process_segment, (task,)) for task in tasks] # 处理完成的任务 for i, future in enumerate(futures): try: result = future.get() seg_idx, segment, valid_segment_points, remove_entire_segment, process_time, error = result if error: print(error) # 存储结果 results[seg_idx] = (segment, valid_segment_points, remove_entire_segment, process_time) print(f"区域 {seg_idx+1}/{len(segments)} 处理完成,用时: {process_time:.2f}秒") except Exception as e: print(f"区域 {i} 结果获取失败: {str(e)}") results[i] = (segments[i], np.ones(len(segments[i]), dtype=bool), False, 0) # 处理所有结果 for seg_idx, result in enumerate(results): if result is None: continue segment, valid_segment_points, remove_entire_segment, _ = result if remove_entire_segment: # 移除整个区域 for idx in segment: valid_points[idx] = False print(f"区域 {seg_idx+1} 点数不足,已移除") else: # 更新全局有效点标记 for j, orig_idx in enumerate(segment): if not valid_segment_points[j]: valid_points[orig_idx] = False # 创建有效点的点云 valid_indices = np.where(valid_points)[0] result_pcd = pcd.select_by_index(valid_indices) print(f"处理完成,剩余点数: {len(result_pcd.points)}") return result_pcd def main(): # 步骤1: 读取点云文件 file_path = r"D:\chenle\ms\2.pcd" print(f"加载点云文件: {file_path}") pcd = o3d.io.read_point_cloud(file_path) if not pcd.has_points(): print("错误: 无法读取点云文件") return print(f"原始点云点数: {len(pcd.points)}") # 步骤2: 分区域处理 print("开始区域分割...") segments = region_growing_segmentation(pcd, radius=0.02) assign_random_colors(pcd, segments) # 步骤7: 迭代处理三次 for i in range(3): print(f"\n=== 开始第 {i+1} 次迭代 ===") # 使用CPU核心数一半的进程数,避免内存溢出 worker_count = min(4, max(1, mp.cpu_count() // 2)) pcd = process_segments(pcd, segments, max_workers=worker_count) segments = region_growing_segmentation(pcd) # 重新分割 # 步骤9: 保存结果 output_path = r"D:\chenle\ms\processed.pcd" o3d.io.write_point_cloud(output_path, pcd) print(f"结果已保存至: {output_path}") if __name__ == "__main__": # 设置多进程启动方法 mp.set_start_method('spawn', force=True) main() 出现区域 123 结果获取失败: SynchronizedArray objects should only be shared between processes through inheritance错误,更改,并给出完整解决后的代码
11-19
import open3d as o3d import numpy as np import time import random from collections import deque import multiprocessing as mp import sys import traceback import ctypes import os import tempfile import logging # 配置日志系统 logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s', handlers=[ logging.FileHandler('pointcloud_processing.log'), logging.StreamHandler() ]) logger = logging.getLogger('PointCloudProcessor') def region_growing_segmentation(pcd, radius=0.02): """ 区域生长分割算法 :param pcd: 输入点云 :param radius: 邻域搜索半径(米) :return: 分割后的区域列表(每个区域是点索引列表) """ logger.info("开始区域生长分割...") points = np.asarray(pcd.points) n_points = len(points) kdtree = o3d.geometry.KDTreeFlann(pcd) visited = np.zeros(n_points, dtype=bool) segments = [] for idx in range(n_points): if visited[idx]: continue # 创建新区域 segment = [] queue = deque([idx]) visited[idx] = True while queue: current_idx = queue.popleft() segment.append(current_idx) # 搜索半径内的邻域点 [_, idxs, _] = kdtree.search_radius_vector_3d(pcd.points[current_idx], radius) for n_idx in idxs: if not visited[n_idx]: visited[n_idx] = True queue.append(n_idx) # 只保留足够大的区域 if len(segment) > 50: segments.append(segment) logger.info(f"分割完成,共 {len(segments)} 个区域") return segments def assign_random_colors(pcd, segments): """ 为每个区域分配随机颜色 :param pcd: 点云对象 :param segments: 分割区域列表 """ colors = np.zeros((len(pcd.points), 3)) for segment in segments: color = [random.random(), random.random(), random.random()] for idx in segment: colors[idx] = color pcd.colors = o3d.utility.Vector3dVector(colors) def process_segment(args): """ 处理单个区域的函数,用于多进程执行 :param args: 包含(区域索引, 区域点索引列表, 点云文件路径, 法向量半径, 计数半径)的元组 :return: (区域索引, 有效点标记数组, 是否移除整个区域标志, 处理时间, 错误信息) """ seg_idx, segment, pcd_file_path, radius_normal, radius_count = args start_time = time.time() error = None try: # 从文件加载点云 pcd = o3d.io.read_point_cloud(pcd_file_path) if not pcd.has_points(): error = f"无法加载点云文件: {pcd_file_path}" valid_segment_points = np.ones(len(segment), dtype=bool) remove_entire_segment = False end_time = time.time() return seg_idx, segment, valid_segment_points, remove_entire_segment, end_time - start_time, error # 创建KDTree kdtree = o3d.geometry.KDTreeFlann(pcd) # 创建该区域的点云 segment_pcd = pcd.select_by_index(segment) # 下采样 voxel_size = radius_normal / 2 down_pcd = segment_pcd.voxel_down_sample(voxel_size) if down_pcd is None or len(down_pcd.points) < 10: # 点数太少,跳过重建 valid_segment_points = np.ones(len(segment), dtype=bool) remove_entire_segment = len(segment) < 100 end_time = time.time() return seg_idx, segment, valid_segment_points, remove_entire_segment, end_time - start_time, error # 估计法向量 down_pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid( radius=radius_normal*2, max_nn=30)) # 泊松重建 - 使用更稳定的参数 mesh = o3d.geometry.TriangleMesh() if len(down_pcd.points) > 50: # 确保有足够点进行重建 try: # 使用VerbosityContextManager设置日志级别 mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( down_pcd, depth=8, linear_fit=True, scale=1.1) except Exception as e: error = f"泊松重建失败: {str(e)}" logger.warning(error) else: logger.debug(f"区域 {seg_idx} 点云点数不足({len(down_pcd.points)}),跳过泊松重建") # 计算区域平均法向量A2 segment_normals = np.asarray(down_pcd.normals) A2 = np.mean(segment_normals, axis=0) if np.linalg.norm(A2) > 1e-6: A2 /= np.linalg.norm(A2) else: logger.warning(f"区域 {seg_idx} 平均法向量为零向量") # 步骤4: 计算每个点的平均法向量A1并过滤 valid_segment_points = np.ones(len(segment), dtype=bool) for j, orig_idx in enumerate(segment): # 找到该点在下采样后的点云中的对应点 [_, idxs, _] = kdtree.search_knn_vector_3d(pcd.points[orig_idx], 1) down_idx = idxs[0] if idxs else -1 # 如果找不到对应点,跳过 if down_idx == -1: continue # 获取以该点为顶点的所有三角形 triangles = [] if mesh.has_triangles(): for tri in mesh.triangles: if down_idx in tri: triangles.append(tri) if len(triangles) == 0: continue # 计算这些三角形的法向量平均值A1 tri_normals = [] mesh_vertices = np.asarray(mesh.vertices) for tri in triangles: # 三角形三个顶点 pts = [mesh_vertices[i] for i in tri] # 计算三角片法向量 v1 = pts[1] - pts[0] v2 = pts[2] - pts[0] normal = np.cross(v1, v2) norm_val = np.linalg.norm(normal) if norm_val > 1e-6: normal /= norm_val tri_normals.append(normal) if tri_normals: A1 = np.mean(tri_normals, axis=0) norm_val = np.linalg.norm(A1) if norm_val > 1e-6: A1 /= norm_val else: # 零向量无法计算角度,跳过 continue # 计算A1和A2的夹角(度) if np.linalg.norm(A1) > 1e-6 and np.linalg.norm(A2) > 1e-6: cos_angle = np.clip(np.dot(A1, A2), -1.0, 1.0) angle = np.degrees(np.arccos(cos_angle)) # 如果偏差大于10度,标记为无效点 if angle > 10: valid_segment_points[j] = False logger.debug(f"区域 {seg_idx} 点 {orig_idx} 法向量偏差 {angle:.2f} 度") # 步骤5: 半径内点云计数过滤 for j, orig_idx in enumerate(segment): if not valid_segment_points[j]: continue [k, idxs, _] = kdtree.search_radius_vector_3d(pcd.points[orig_idx], radius_count) # 排除自身点 neighbor_count = k - 1 if k > 0 else 0 # 如果邻域点少于100个,标记为无效点 if neighbor_count < 100: valid_segment_points[j] = False logger.debug(f"区域 {seg_idx} 点 {orig_idx} 邻域点数不足: {neighbor_count}") # 步骤6: 检查是否移除整个区域 valid_count = np.sum(valid_segment_points) remove_entire_segment = valid_count < 100 end_time = time.time() process_time = end_time - start_time logger.info(f"区域 {seg_idx} 处理完成,有效点: {valid_count}/{len(segment)},用时: {process_time:.2f}秒") return seg_idx, segment, valid_segment_points, remove_entire_segment, process_time, error except Exception as e: # 捕获并记录异常 error = f"区域 {seg_idx} 处理错误: {str(e)}\n{traceback.format_exc()}" logger.error(error) valid_segment_points = np.ones(len(segment), dtype=bool) remove_entire_segment = False end_time = time.time() return seg_idx, segment, valid_segment_points, remove_entire_segment, end_time - start_time, error finally: # 清理临时文件 if os.path.exists(pcd_file_path): try: os.remove(pcd_file_path) except Exception as e: logger.error(f"删除临时文件失败: {pcd_file_path}, 错误: {str(e)}") def process_segments(pcd, segments, radius_normal=0.1, radius_count=0.1, max_workers=4): """ 多进程处理点云区域 :param pcd: 点云对象 :param segments: 分割区域列表 :param radius_normal: 法向量计算半径(米) :param radius_count: 点云计数半径(米) :param max_workers: 最大进程数 :return: 处理后的点云 """ logger.info(f"开始多进程处理区域,共 {len(segments)} 个区域,使用 {max_workers} 个进程...") # 创建有效点标记数组 valid_points = np.ones(len(pcd.points), dtype=bool) # 准备任务参数 tasks = [] temp_files = [] # 保存点云到临时文件 temp_dir = tempfile.mkdtemp() pcd_file_path = os.path.join(temp_dir, "temp_pcd.pcd") o3d.io.write_point_cloud(pcd_file_path, pcd) logger.info(f"主点云已保存至: {pcd_file_path}") for i, segment in enumerate(segments): # 创建临时文件副本 temp_file_path = os.path.join(temp_dir, f"temp_pcd_{i}.pcd") try: os.link(pcd_file_path, temp_file_path) except OSError: # 如果硬链接失败(例如跨文件系统),则复制文件 import shutil shutil.copy(pcd_file_path, temp_file_path) temp_files.append(temp_file_path) tasks.append(( i, segment, temp_file_path, radius_normal, radius_count )) logger.debug(f"为区域 {i} 创建临时文件: {temp_file_path}") # 用于存储结果的列表 results = [None] * len(segments) # 使用进程池并行处理 with mp.Pool(processes=max_workers) as pool: # 提交所有任务 futures = [pool.apply_async(process_segment, (task,)) for task in tasks] # 处理完成的任务 for i, future in enumerate(futures): try: result = future.get() seg_idx, segment, valid_segment_points, remove_entire_segment, process_time, error = result if error: logger.error(error) # 存储结果 results[seg_idx] = (segment, valid_segment_points, remove_entire_segment, process_time) logger.info(f"区域 {seg_idx+1}/{len(segments)} 处理完成,用时: {process_time:.2f}秒") except Exception as e: logger.error(f"区域 {i} 结果获取失败: {str(e)}", exc_info=True) results[i] = (segments[i], np.ones(len(segments[i]), dtype=bool), False, 0) # 清理临时文件 for file_path in temp_files: if os.path.exists(file_path): try: os.remove(file_path) except Exception as e: logger.error(f"删除临时文件失败: {file_path}, 错误: {str(e)}") if os.path.exists(pcd_file_path): try: os.remove(pcd_file_path) except Exception as e: logger.error(f"删除主临时文件失败: {pcd_file_path}, 错误: {str(e)}") try: os.rmdir(temp_dir) except Exception as e: logger.error(f"删除临时目录失败: {temp_dir}, 错误: {str(e)}") # 处理所有结果 for seg_idx, result in enumerate(results): if result is None: continue segment, valid_segment_points, remove_entire_segment, _ = result if remove_entire_segment: # 移除整个区域 for idx in segment: valid_points[idx] = False logger.info(f"区域 {seg_idx+1} 点数不足,已移除") else: # 更新全局有效点标记 for j, orig_idx in enumerate(segment): if not valid_segment_points[j]: valid_points[orig_idx] = False # 创建有效点的点云 valid_indices = np.where(valid_points)[0] result_pcd = pcd.select_by_index(valid_indices) logger.info(f"处理完成,剩余点数: {len(result_pcd.points)}") return result_pcd def main(): # 步骤1: 读取点云文件 file_path = r"D:\chenle\ms\2.pcd" logger.info(f"加载点云文件: {file_path}") pcd = o3d.io.read_point_cloud(file_path) if not pcd.has_points(): logger.error("错误: 无法读取点云文件") return logger.info(f"原始点云点数: {len(pcd.points)}") # 步骤2: 分区域处理 logger.info("开始区域分割...") segments = region_growing_segmentation(pcd, radius=0.02) assign_random_colors(pcd, segments) # 步骤7: 迭代处理三次 for i in range(3): logger.info(f"\n=== 开始第 {i+1} 次迭代 ===") # 使用CPU核心数一半的进程数,避免内存溢出 worker_count = min(4, max(1, mp.cpu_count() // 2)) pcd = process_segments(pcd, segments, max_workers=worker_count) # 重新分割前过滤小区域 segments = [seg for seg in region_growing_segmentation(pcd) if len(seg) > 50] logger.info(f"第 {i+1} 次迭代后剩余区域数: {len(segments)}") # 步骤9: 保存结果 output_path = r"D:\chenle\ms\processed.pcd" o3d.io.write_point_cloud(output_path, pcd) logger.info(f"结果已保存至: {output_path}") # 可视化结果 logger.info("显示处理后的点云...") o3d.visualization.draw_geometries([pcd]) if __name__ == "__main__": # 设置多进程启动方法 mp.set_start_method('spawn', force=True) main() 更改为8进程,其他不变
最新发布
11-19
%% ICG特征点检测算法 % 样本周期分割 % === 周期分割优化(同步ECG周期边界)=== % === 信号质量评估 (兼容性修复) === try % 尝试使用原生range函数 signalRange = range(app.icgData_diff); catch % 原生函数不可用时的手动计算 signalRange = max(app.icgData_diff) - min(app.icgData_diff); end % 改进幅度阈值逻辑 % baselineNoise = std(app.icgData_diff(app.icgData_diff < prctile(app.icgData_diff, 10))); % if signalRange < 5 * baselineNoise % 动态阈值 % errordlg(sprintf('信号质量过低 (动态范围: %.4f, 基线噪声: %.4f)', ... % signalRange, baselineNoise), '信号质量警告'); % return; % end % === 动态参数计算 === % 1. 自适应显著度阈值 noiseLevel = std(app.icgData_diff(1:min(100, end))); % 前100点噪声估计 prominenceThresh = max(noiseLevel * 3, 0.1 * signalRange); switch app.DropDown_sele.Value case 'ARM' % 2. 智能最小峰值距离 defaultDistance = 0.5; % 默认值取最小高度 case 'FPGA' defaultDistance = 0.7; % 默认值取最小高度 end adaptiveDistance = 50; % 最小距离 % === 分级峰值检测策略 === [pks, locs] = findpeaks1(app, app.icgData_diff, defaultDistance, adaptiveDistance); % 初级检测失败时启用备用方案 if numel(locs) < 2 % 方案1:降低显著度要求 [pks2, locs2] = findpeaks(... app.icgData_diff, ... 'MinPeakDistance', adaptiveDistance * 0.7, ... 'MinPeakProminence', prominenceThresh * 0.5 ... ); % 方案2:基于局部极大值检测 if numel(locs2) < 2 [~, locs3] = findpeaks(... movmax(app.icgData_diff, 5) == app.icgData_diff, ... 'MinPeakDistance', adaptiveDistance ... ); locs = locs3; else locs = locs2; end end app.icg_C_points = locs; % 新增:保存检测到的ICG C点全局位置 fprintf('成功检测到%d个C峰\n', numel(locs)); % 使用ECG划分的周期数 numCycles = numel(app.ecgsegments); % 计算整体平均心率(用于全局参数) % if numCycles >= 1 % RR_intervals = diff(app.R_points(1:numCycles+1)); % meanRR = mean(RR_intervals); % currentHR = 60 / (meanRR / app.fs); % else % currentHR = 60; % 默认值 % end % 计算周期心率(使用ECG的R峰间隔) % RR_interval = app.R_points(i+1) - app.R_points(i); % app.segments(i).HR = 60 / (RR_interval / app.fs); % app.MinBCDistance = round(0.08 * app.fs); % 最小BC间距80ms (固定值) % app.MaxBCDistance = round(0.16 * currentHR/60 * app.fs); % 基于平均心率计算 % 周期处理 if numCycles < 1 app.segments = []; else app.segments = struct(... 'HR', 0, 'C_point', 0, 'C_value', 0, ... 'C2_point', 0, 'C2_value', 0,... 'B_point', 0, 'B_value', 0, 'X_point', 0, 'X_value', 0, ... 'BX_time', 0, 'startIdx', 0, 'endIdx', 0, ... 'valid', false, 'SV', 0, 'CO', 0, 'CI', 0, 'SVI', 0); app.segments = repmat(app.segments, 1, numCycles); for i = 2:numCycles app.segments(i).valid = false; % 关键修改:直接使用ECG周期边界 % 获取当前周期边界(注意偏移量可能导致越界) segStart = max(1, app.ecgsegments(i).start_idx - 30); % 边界保护 segEnd = min(length(app.icgData_diff), app.ecgsegments(i).end_idx + 30); % 边界保护 % 验证信号段有效性 if segStart >= segEnd || (segEnd - segStart) < 50 fprintf('周期%d: 无效的信号段长度\n', i); continue; % 跳过当前周期 end % 存储周期信息 app.segments(i).startIdx = segStart ; app.segments(i).endIdx = segEnd ; segmentData = app.icgData_diff(segStart:segEnd); segmentData_dz2dt = app.icgData_diff_dz2dt(segStart:segEnd); %% 二阶微分 % 特征点检测 % 特征点检测 try %% 1. 检测C点(一阶微分最大值)并打印 % [c_local, c_value] = detectCpoint(app, segmentData); % app.segments(i).C_point = segStart + c_local - 1; % app.segments(i).C_value = c_value; % fprintf('周期%d: C点检测结果 -> 局部位置: %d, 全局位置: %d, 幅值: %.2fV\n', ... % i, c_local, app.segments(i).C_point, c_value); % 在当前ECG周期内查找预存的ICG C点 idx_in_segment = find(app.icg_C_points >= segStart & app.icg_C_points <= segEnd); if isempty(idx_in_segment) fprintf('周期%d: 未找到C点,跳过特征点检测\n', i); app.segments(i).valid = false; continue; % 无C点时跳过该周期 else % 处理多个C点:取最接近周期中心的点 if numel(idx_in_segment) > 1 segmentCenter = round((segStart + segEnd) / 2); [~, minIdx] = min(abs(app.icg_C_points(idx_in_segment) - segmentCenter)); c_global = app.icg_C_points(idx_in_segment(minIdx)); else c_global = app.icg_C_points(idx_in_segment(1)); end c_local = c_global - segStart + 1; c_value = app.icgData_diff(c_global); app.segments(i).C_point = c_global; app.segments(i).C_value = c_value; fprintf('周期%d: C点使用预检测结果 -> 全局位置: %d, 幅值: %.2fV\n', i, c_global, c_value); end %% 2. 二阶微分处理(C2点检测)并打印 minPeakHeight = 0.1; minPeakDistance = 20; [pks2, locs2] = findpeaks(segmentData_dz2dt, ... 'MinPeakHeight', minPeakHeight, ... 'MinPeakDistance', minPeakDistance); if isempty(locs2) [c2_value, c2_local] = max(segmentData_dz2dt); locs2 = c2_local; pks2 = c2_value; end [~, idx] = min(abs(locs2 - c_local)); app.segments(i).C2_point = segStart + locs2(idx) - 1; app.segments(i).C2_value = pks2(idx); fprintf('周期%d: C2点检测结果 -> 局部位置: %d, 全局位置: %d, 幅值: %.2fV\n', ... i, locs2(idx), app.segments(i).C2_point, pks2(idx)); %% 3. 检测B点并打印 % 检查C点有效性 if isempty(app.segments(i).C_point) || isnan(app.segments(i).C_value) fprintf('周期%d: C点无效,跳过B点检测\n', i); app.segments(i).B_point = NaN; app.segments(i).B_value = NaN; continue; end % 检查C2点有效性 if ~isfield(app.segments(i), 'C2_point') || isempty(locs2) || idx > numel(locs2) fprintf('周期%d: C2点无效,使用单参数模式\n', i); [b_local, b_value] = detectBpoint_simple(app, segmentData, app.segments(i).C_point - segStart + 1, app.segments(i).C_value); else % 标准调用(带C2点) [b_local, b_value] = detectBpoint(app, segmentData, c_local, c_value, locs2(idx) ); % 确保局部位置 end app.segments(i).B_point = segStart + b_local - 1; app.segments(i).B_value = b_value; fprintf('周期%d: B点检测结果 -> 局部位置: %d, 全局位置: %d, 幅值: %.2fV\n', i, b_local, app.segments(i).B_point, b_value); %% 4. 检测X点并打印 [x_local, x_value] = detectXpoint(app, segmentData, c_local, app.segments(i).B_point, segStart, app.fs); app.segments(i).X_point = segStart + x_local - 1; app.segments(i).X_value = x_value; fprintf('周期%d: X点检测结果 -> 局部位置: %d, 全局位置: %d, 幅值: %.2fV\n', i, x_local, app.segments(i).X_point, x_value); %% 计算并打印BX时间 if app.segments(i).B_point < app.segments(i).X_point app.segments(i).BX_time = (app.segments(i).X_point - app.segments(i).B_point) / app.fs; fprintf('周期%d: BX时间 = %.3fs\n', i, app.segments(i).BX_time); else app.segments(i).BX_time = NaN; fprintf('周期%d: B点位置超过X点,BX时间无效\n', i); end app.segments(i).valid = true; fprintf('周期%d特征点检测成功!\n', i); catch ME fprintf('周期%d特征点检测失败: %s\n', i, ME.message); app.segments(i).valid = false; end end end % === 结果输出 === validCount = sum([app.segments.valid]); fprintf('成功处理%d/%d个周期:\n', validCount, numel(app.segments)); app.segments = app.segments; % 在计算segments的回调函数中执行 % 调用一次show_plotButtonPushed来绘制所有特征点 app.show_plotButtonPushed(); % 初始化累加变量 total_C_value = 0; total_BX_time = 0; total_C2_value = 0; valid_count = 0; for i = 1:numel(app.segments) if app.segments(i).valid % 打印当前周期信息 fprintf('周期%d: B@%d(%.2fV), C@%d(%.2fV), C2@%d(%.2fV), X@%d(%.2fV), BX=%.3fs\n',... i, app.segments(i).B_point, app.segments(i).B_value,... app.segments(i).C_point, app.segments(i).C_value,... app.segments(i).C2_point, app.segments(i).C2_value,... app.segments(i).X_point, app.segments(i).X_value,... app.segments(i).BX_time); % 累加有效周期的值 total_C_value = total_C_value + app.segments(i).C_value; total_BX_time = total_BX_time + app.segments(i).BX_time; total_C2_value = total_C2_value + app.segments(i).C2_value; valid_count = valid_count + 1; end end % 计算平均值(避免除零错误) if valid_count > 0 app.avg_C_value = total_C_value / valid_count; app.avg_BX_time = total_BX_time / valid_count; app.avg_C2_value = total_C2_value / valid_count; fprintf('\n==== 统计结果 ====\n'); fprintf('有效周期数: %d\n', valid_count); fprintf('平均C值: %.2fV\n', app.avg_C_value); fprintf('平均C2值: %.2fV\n', app.avg_C2_value); fprintf('平均BX时间: %.3fs\n', app.avg_BX_time); else fprintf('\n警告: 未发现有效周期数据\n'); end 根据代码生成特征点识别与特征参数提取的处理流程
10-29
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值