cv/main.ipynb

{
 "cells": [
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-08T15:04:41.951180Z",
     "start_time": "2024-05-08T15:04:41.479489Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ],
   "outputs": [],
   "execution_count": 1
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-08T16:52:58.149750Z",
     "start_time": "2024-05-08T16:52:58.140930Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def fit_contour(contour):\n",
    "    # 轮廓逼近\n",
    "    epsilon = 0.01 * cv2.arcLength(contour, True)\n",
    "    approx = cv2.approxPolyDP(contour, epsilon, True)\n",
    "    return approx\n",
    "\n",
    "\n",
    "def cal_curvature(contour):\n",
    "    # 计算曲率\n",
    "    curvatures = []\n",
    "    length = len(contour)\n",
    "    for i in range(1, length):\n",
    "        prev = contour[(i - 1) % length].ravel()  # 将二维数组转换为一维数组\n",
    "        curr = contour[i].ravel()\n",
    "        next = contour[(i + 1) % length].ravel()\n",
    "        # 使用向量角度来近似曲率\n",
    "        vector1 = prev - curr\n",
    "        vector2 = next - curr\n",
    "        dot_product = np.dot(vector1, vector2)\n",
    "        norm_product = np.linalg.norm(vector1) * np.linalg.norm(vector2)\n",
    "        if norm_product == 0:\n",
    "\n",
    "            angle = 0  # 避免除以零的情况\n",
    "        else:\n",
    "            cosine_angle = np.clip(dot_product / norm_product, -1.0, 1.0)  # 避免计算错误导致的数值问题\n",
    "            angle = np.arccos(cosine_angle)\n",
    "        curvatures.append(angle)\n",
    "    length = len(curvatures)\n",
    "    first_order_differences = [abs(curvatures[i] - curvatures[(i + 1) % length]) for i in range(length)]\n",
    "    # length = len(first_order_differences)\n",
    "    # first_order_differences = [abs(first_order_differences[i] - first_order_differences[(i + 1) % length]) for i in range(length)]\n",
    "    return first_order_differences\n",
    "\n",
    "\n",
    "def fourier_trans(contour):\n",
    "    # 傅里叶描述子\n",
    "    contour_array = contour[:, 0, :]\n",
    "    contour_complex = np.empty(contour_array.shape[:-1], dtype=complex)\n",
    "    contour_complex.real = contour_array[:, 0]\n",
    "    contour_complex.imag = contour_array[:, 1]\n",
    "    fourier_result = np.fft.fft(contour_complex)\n",
    "    return fourier_result\n",
    "\n",
    "\n",
    "# 读取图片\n",
    "def show_contours(path_image, contour1=None):\n",
    "    image = cv2.imread(path_image)\n",
    "\n",
    "    # 将图片转换为RGB格式\n",
    "    image_rgb = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n",
    "\n",
    "    # 增强饱和度\n",
    "    hsv = cv2.cvtColor(image_rgb, cv2.COLOR_RGB2HSV)\n",
    "    hsv[..., 1] = hsv[..., 1] * 6  # 增强饱和度\n",
    "    enhanced_image_rgb = cv2.cvtColor(hsv, cv2.COLOR_HSV2RGB)\n",
    "\n",
    "    # 进行Prewitt边缘检测\n",
    "    gray_image = cv2.cvtColor(enhanced_image_rgb, cv2.COLOR_RGB2GRAY)\n",
    "    prewittx = cv2.Sobel(gray_image, cv2.CV_64F, 1, 0, ksize=3)\n",
    "    prewitty = cv2.Sobel(gray_image, cv2.CV_64F, 0, 1, ksize=3)\n",
    "    prewitt_edge = cv2.magnitude(prewittx, prewitty)\n",
    "\n",
    "    _, thresholded_image = cv2.threshold(prewitt_edge, 50, 255, cv2.THRESH_BINARY)\n",
    "\n",
    "    # 寻找轮廓\n",
    "    contours, hierarchy = cv2.findContours(thresholded_image.astype(np.uint8), cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)\n",
    "\n",
    "    # 过滤面积小于阈值的轮廓\n",
    "    filtered_contours = []\n",
    "    # 计算和显示每个轮廓的最小斜率\n",
    "    c = 0\n",
    "    for contour in contours:\n",
    "        rect = cv2.minAreaRect(contour)  # 获取最小外接矩形\n",
    "        box = cv2.boxPoints(rect)  # 计算矩形的四个角点\n",
    "        box = np.intp(box)  # 角点坐标整数化  # 角点坐标整数化\n",
    "\n",
    "        # 求边长\n",
    "        widths = [np.linalg.norm(box[i] - box[(i + 1) % 4]) for i in range(4)]\n",
    "        sorted_widths = sorted(widths)\n",
    "\n",
    "        # 宽和高\n",
    "        width, height = sorted_widths[0], sorted_widths[1]\n",
    "\n",
    "        # 角度调整\n",
    "        angle = rect[2]\n",
    "        # if width < height:\n",
    "        #     angle += 90\n",
    "        area = width * height\n",
    "        # 计算斜率\n",
    "        slope = np.tan(np.radians(angle))\n",
    "        # if angle < 45 and area > 5000:\n",
    "        if area <= 5000:\n",
    "            continue\n",
    "        \n",
    "        # 计算曲率\n",
    "        curvatures = cal_curvature(contour)\n",
    "\n",
    "        # 计算标准差和平均值\n",
    "        mean_change = np.mean(curvatures)\n",
    "        std_change = np.std(curvatures)\n",
    "        threshold = mean_change + 1 * std_change  # 设定阈值为平均值加两倍标准差\n",
    "\n",
    "        # 寻找异常值\n",
    "        outliers = [change for change in curvatures if change > threshold]\n",
    "\n",
    "        if len(outliers) < 20:\n",
    "            continue\n",
    "\n",
    "        # curvatures = pd.Series(curvatures)\n",
    "        # std = curvatures.describe()['std']\n",
    "        # if std < 0.15:\n",
    "        #     continue\n",
    "        filtered_contours.append(contour)\n",
    "        c += 1\n",
    "\n",
    "    # 将二值图像转换为RGB格式，以便能在其上绘制彩色轮廓\n",
    "    contoured_image = cv2.cvtColor(thresholded_image.astype(np.uint8), cv2.COLOR_GRAY2RGB)\n",
    "\n",
    "    # 使用红色来绘制轮廓\n",
    "    if contour1 is None:\n",
    "        cv2.drawContours(contoured_image, filtered_contours, -1, (255, 0, 0), 2)\n",
    "    else:\n",
    "        cv2.drawContours(contoured_image, [contour1], -1, (255, 0, 0), 2)\n",
    "\n",
    "    # 使用matplotlib显示图像\n",
    "    plt.figure(figsize=(8, 8))\n",
    "    plt.imshow(contoured_image)\n",
    "    plt.title('Filtered Image with Red Contours and Slope Info')\n",
    "    plt.axis('off')\n",
    "    plt.show()\n",
    "    return filtered_contours"
   ],
   "outputs": [],
   "execution_count": 183
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-08T16:57:20.434078Z",
     "start_time": "2024-05-08T16:57:20.137351Z"
    }
   },
   "cell_type": "code",
   "source": [
    "\n",
    "tmp = show_contours(\"00.jpg\")"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAKSCAYAAACdh68YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3xN5x/A8c8duTd7SSKyJSISm9pi701RtVWNVlvVodtu1ShtVbWqKFpao0btnZghQhorG5HI3jv3Pr8/0tyfLIIgOO/X67xenPvcc557cs75nvNMmRBCIJFIJBLJUyZ/2hmQSCQSiQSkgCSRSCSSKkIKSBKJRCKpEqSAJJFIJJIqQQpIEolEIqkSpIAkkUgkkipBCkgSiUQiqRKkgCSRSCSSKkEKSBKJRCKpEio1IEVGRiKTyVi7dq1u3axZs5DJZJW5m8fqWctvVdShQwc6dOhQ4bT16tV7vBmqJC4uLowdO/ZpZ0NSSdauXYtMJiMyMrLStimTyZg1a1alba+qKSgoYPr06Tg6OiKXyxkwYEClbv+BAlLRH7Cs5eOPP67wdr766iu2b9/+oHmtUsaOHYuxsfHTzsYzITo6mlmzZnHx4sVK37aLi0ux89DIyIjmzZuzbt26St/Xg8jJyWHp0qW0aNECMzMz9PX1qV27Nm+99RbBwcGPdd/Pw/VV1Zw4cYKePXtib2+Pvr4+Tk5O9O3blz/++ONpZ+2BHDt2DJlMxpYtWx7q+6tXr2bRokUMHjyY3377jWnTplVq/pQP86U5c+ZQs2bNYuvq1auHs7Mz2dnZ6Onp3fP7X331FYMHD6706CqpGg4cOFDs/9HR0cyePRsXFxcaNWpU6ftr1KgR77//PgAxMTGsWrWKMWPGkJuby4QJEyp9f/eTkJBAjx498Pf3p0+fPgwfPhxjY2OuX7/Opk2bWLlyJXl5eY9t/9L1Vbk2b97MK6+8QqNGjZg6dSoWFhZERETg4+PDL7/8wvDhw592Fp+YI0eOYG9vz9KlSx/L9h8qIPXs2ZOXXnqpzM/09fUfKUMPKycnB5VKhVwuVYs9bSqV6onuz97enpEjR+r+P3bsWFxdXVm6dOlTCUhjx44lICCALVu28PLLLxf7bO7cuXz22WdPPE9PihCCnJwcDAwMnnZWKs2sWbPw8vLizJkzpc7tuLi4p5SrpyMuLg5zc/PHtv3HXodUkkwmIzMzk99++01XzHJ3ufzt27d57bXXqF69Omq1mrp167J69epi2yh67dy0aROff/459vb2GBoakpaWBsDZs2fp0aMHZmZmGBoa0r59e06ePFkqLydOnKBZs2bo6+vj5ubGzz///Ei/38XFhT59+nDs2DFeeuklDAwMqF+/PseOHQNg27Zt1K9fH319fZo2bUpAQECx7wcGBupupvr6+tja2vLaa6+RmJhYal9F+7g77+XVf23YsIGmTZtiYGCApaUlw4YN49atW/f8LYGBgchkMnbu3Klb5+/vj0wmo0mTJsXS9uzZkxYtWuj+f3cd0rFjx2jWrBkA48aN0/3NS54jV65coWPHjhgaGmJvb8/ChQvvmb97sba2pk6dOoSFhRVbr9Vq+fbbb6lbty76+vpUr16dSZMmkZycXCydEIJ58+bh4OCAoaEhHTt25PLlyxXa99mzZ9m9ezfjx48vFYwA1Go1ixcvLrbuyJEjeHt7Y2RkhLm5Of379+fq1avF0hT9bUNDQxk7dizm5uaYmZkxbtw4srKydOnud30FBATQs2dPTE1NMTY2pnPnzpw5c6bMfZVUVp1L0Tm/f/9+3TlfdB0dPHiQtm3bYm5ujrGxMR4eHnz66af3PYZr1qyhU6dO2NjYoFar8fLyYsWKFaXSFe37xIkTNG/eHH19fVxdXcssrr18+TKdOnXCwMAABwcH5s2bh1arvW9eAMLCwmjWrFmZD1o2Njb3/X5FjnnRsfXx8WHSpElUq1YNU1NTRo8eXer8BNi7d6/unDExMaF3794VPkdLqsi5VXRvP3r0KJcvX9adW0X3tszMTN5//30cHR1Rq9V4eHiwePFiHnQyiYd6Q0pNTSUhIaHYOisrqwp9d/369bz++us0b96ciRMnAuDm5gZAbGwsLVu2RCaT8dZbb2Ftbc3evXsZP348aWlpvPvuu8W2NXfuXFQqFR988AG5ubmoVCqOHDlCz549adq0KTNnzkQul+tOcF9fX5o3bw7Av//+S7du3bC2tmbWrFkUFBQwc+ZMqlev/jCHRCc0NJThw4czadIkRo4cyeLFi+nbty8//fQTn376KW+++SYA8+fPZ+jQoVy/fl33Vnfw4EHCw8MZN24ctra2XL58mZUrV3L58mXOnDmju0kEBATQo0cPatSowezZs9FoNMyZMwdra+tS+fnyyy/54osvGDp0KK+//jrx8fEsW7aMdu3aERAQUO7TTr169TA3N8fHx4d+/foB4Ovri1wu59KlS6SlpWFqaopWq+XUqVO6v2VJnp6ezJkzhxkzZjBx4kS8vb0BaN26tS5NcnIyPXr0YNCgQQwdOpQtW7bw0UcfUb9+fXr27PnAf4OCggKioqKwsLAotn7SpEmsXbuWcePG8c477xAREcEPP/xAQEAAJ0+e1BU1z5gxg3nz5tGrVy969erFhQsX6NatW4WK2YoC+KhRoyqU10OHDtGzZ09cXV2ZNWsW2dnZLFu2jDZt2nDhwgVcXFyKpR86dCg1a9Zk/vz5XLhwgVWrVmFjY8OCBQuAe19fly9fxtvbG1NTU6ZPn46enh4///wzHTp04Pjx48UeKh7E9evXefXVV5k0aRITJkzAw8ODy5cv06dPHxo0aMCcOXNQq9WEhoaW+WBY0ooVK6hbty79+vVDqVSya9cu3nzzTbRaLVOmTCmWNjQ0lMGDBzN+/HjGjBnD6tWrGTt2LE2bNqVu3boA3Llzh44dO1JQUMDHH3+MkZERK1eurPBbnLOzM4cPHyYqKgoHB4cHOjYPeszfeustzM3NmTVrFtevX2fFihXcuHFD9xAOhX/jMWPG0L17dxYsWEBWVhYrVqygbdu2BAQElDpnKupe55a1tTXr16/nyy+/JCMjg/nz5wOF17cQgn79+nH06FHGjx9Po0aN2L9/Px9++CG3b99+sOI98QDWrFkjgDIXIYSIiIgQgFizZo3uOzNnzhQld2NkZCTGjBlTavvjx48XNWrUEAkJCcXWDxs2TJiZmYmsrCwhhBBHjx4VgHB1ddWtE0IIrVYr3N3dRffu3YVWq9Wtz8rKEjVr1hRdu3bVrRswYIDQ19cXN27c0K27cuWKUCgUpfJbljFjxggjI6Ni65ydnQUgTp06pVu3f/9+AQgDA4Ni+/r5558FII4ePVosnyVt3LhRAMLHx0e3rm/fvsLQ0FDcvn1bty4kJEQolcpieY+MjBQKhUJ8+eWXxbb577//CqVSWWp9Sb179xbNmzfX/X/QoEFi0KBBQqFQiL179wohhLhw4YIAxI4dO3Tp2rdvL9q3b6/7/7lz50qdF3enBcS6det063Jzc4Wtra14+eWX75k/IQqPebdu3UR8fLyIj48X//77rxg1apQAxJQpU3TpfH19BSB+//33Yt/ft29fsfVxcXFCpVKJ3r17FzuHPv30UwGUed7ebeDAgQIQycnJ9827EEI0atRI2NjYiMTERN26S5cuCblcLkaPHq1bV3Qdvfbaa6X2V61atWLryru+BgwYIFQqlQgLC9Oti46OFiYmJqJdu3al9lVS0fUfERGhW1d0zu/bt69Y2qVLlwpAxMfH3/sAlKGs66B79+7C1dW12Lqifd99bcTFxQm1Wi3ef/993bp3331XAOLs2bPF0pmZmZX6PWX59ddfBSBUKpXo2LGj+OKLL4Svr6/QaDSl0gJi5syZuv9X9JgXHdumTZuKvLw83fqFCxcWu77S09OFubm5mDBhQrH93rlzR5iZmZVaX1LRvXPz5s26dQ9ybrVv317UrVu32Lrt27cLQMybN6/Y+sGDBwuZTCZCQ0Pvmae7PVSR3fLlyzl48GCx5VEJIdi6dSt9+/ZFCEFCQoJu6d69O6mpqVy4cKH
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 188
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-08T16:53:47.763789Z",
     "start_time": "2024-05-08T16:53:47.324568Z"
    }
   },
   "cell_type": "code",
   "source": [
    "i = 8\n",
    "tmp = show_contours(\"10.jpg\", tmp[i])\n",
    "curvatures = cal_curvature(tmp[i])\n",
    "curvatures = pd.Series(curvatures)\n",
    "curvatures.describe()"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAKSCAYAAACdh68YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydd1hUx/f/33d7Ydml9w4iRUVFkChi7y32mKgxJtEYY2JJYjT2XmKJRqOJxqixdxMrVizYsKIoIh2k9wUWds/vD397w9JBLPl89/U88yh3770zd+7cOTNnzpzDEBFBjx49evToectw3nYB9OjRo0ePHkAvkPTo0aNHzzuCXiDp0aNHj553Ar1A0qNHjx497wR6gaRHjx49et4J9AJJjx49evS8E+gFkh49evToeSfQCyQ9evTo0fNOoBdIevTo0aPnnaBBBVJMTAwYhsHWrVvZY3PmzAHDMA2ZzWvlv1bed5H27dujffv2tT7X29v79RaogXB0dMTHH3/8touhp4HYunUrGIZBTExMg92TYRjMmTOnwe73rlFaWorvvvsOdnZ24HA46N+/f4Pev04CSfsCK0vTpk2r9X0WLVqEw4cP17Ws7xQff/wxDAwM3nYx/hMkJSVhzpw5uHv3boPf29HRUacdSqVS+Pn5Ydu2bQ2eV10oKirCqlWr4O/vD7lcDpFIhEaNGmHChAl4+vTpa837f+H7ete4fPkyevToARsbG4hEItjb26NPnz7YuXPn2y5anbhw4QIYhsH+/fvrdf2WLVuwfPlyDBo0CH/++ScmTZrUoOXj1eeiefPmwcnJSeeYt7c3HBwcUFhYCD6fX+31ixYtwqBBgxpcuup5Nzh9+rTO30lJSZg7dy4cHR3h4+PT4Pn5+PhgypQpAIDk5GT8/vvvGDVqFIqLi/HZZ581eH41kZ6eju7du+P27dvo3bs3hg8fDgMDAzx58gS7d+/Gpk2boFKpXlv++u+rYdm3bx+GDh0KHx8ffP311zAyMkJ0dDQuXbqE3377DcOHD3/bRXxjnDt3DjY2Nli1atVruX+9BFKPHj3g6+tb6W8ikeiVClRfioqKIBAIwOHol8XeNgKB4I3mZ2Njg48++oj9++OPP4azszNWrVr1VgTSxx9/jDt37mD//v0YOHCgzm/z58/HjBkz3niZ3hREhKKiIojF4rddlAZjzpw58PT0RGhoaIW2nZqa+pZK9XZITU2FQqF4bfd/7WtI5WEYBgUFBfjzzz9ZNUtZvXxiYiI++eQTWFhYQCgUwsvLC1u2bNG5h3bauXv3bvz444+wsbGBRCJBbm4uAOD69evo3r075HI5JBIJgoKCcOXKlQpluXz5Mlq1agWRSAQXFxds3LjxlZ7f0dERvXv3xoULF+Dr6wuxWIwmTZrgwoULAICDBw+iSZMmEIlEaNmyJe7cuaNz/f3799nOVCQSwdLSEp988gkyMjIq5KXNo2zZq1r/2rFjB1q2bAmxWAxjY2MMGzYM8fHx1T7L/fv3wTAMjh49yh67ffs2GIZBixYtdM7t0aMH/P392b/LriFduHABrVq1AgCMHj2afefl28ijR4/QoUMHSCQS2NjYYNmyZdWWrzrMzMzQuHFjREVF6RzXaDRYvXo1vLy8IBKJYGFhgbFjxyIrK0vnPCLCggULYGtrC4lEgg4dOiA8PLxWeV+/fh3//PMPxowZU0EYAYBQKMSKFSt0jp07dw6BgYGQSqVQKBTo168fHj9+rHOO9t0+e/YMH3/8MRQKBeRyOUaPHg2lUsmeV9P3defOHfTo0QOGhoYwMDBAp06dEBoaWmle5alszUXb5k+dOsW2ee13dObMGbRt2xYKhQIGBgZwd3fH9OnTa6zDP/74Ax07doS5uTmEQiE8PT2xYcOGCudp8758+TL8/PwgEong7Oxcqbo2PDwcHTt2hFgshq2tLRYsWACNRlNjWQAgKioKrVq1qnSgZW5uXuP1talzbd1eunQJY8eOhYmJCQwNDTFy5MgK7RMATpw4wbYZmUyGXr161bqNlqc2bUvbt58/fx7h4eFs29L2bQUFBZgyZQrs7OwgFArh7u6OFStWoK7BJOo1Q8rJyUF6errOMVNT01pdu337dnz66afw8/PD559/DgBwcXEBAKSkpKB169ZgGAYTJkyAmZkZTpw4gTFjxiA3NxfffPONzr3mz58PgUCAqVOnori4GAKBAOfOnUOPHj3QsmVLzJ49GxwOh23gISEh8PPzAwA8ePAAXbt2hZmZGebMmYPS0lLMnj0bFhYW9akSlmfPnmH48OEYO3YsPvroI6xYsQJ9+vTBr7/+iunTp2P8+PEAgMWLF2PIkCF48uQJO6s7c+YMnj9/jtGjR8PS0hLh4eHYtGkTwsPDERoaynYSd+7cQffu3WFlZYW5c+dCrVZj3rx5MDMzq1CehQsXYubMmRgyZAg+/fRTpKWlYe3atWjXrh3u3LlT5WjH29sbCoUCly5dQt++fQEAISEh4HA4uHfvHnJzc2FoaAiNRoOrV6+y77I8Hh4emDdvHmbNmoXPP/8cgYGBAID33nuPPScrKwvdu3fHgAEDMGTIEOzfvx/ff/89mjRpgh49etT5HZSWliIhIQFGRkY6x8eOHYutW7di9OjRmDhxIqKjo7Fu3TrcuXMHV65cYVXNs2bNwoIFC9CzZ0/07NkTYWFh6Nq1a63UbFoBPmLEiFqVNTg4GD169ICzszPmzJmDwsJCrF27Fm3atEFYWBgcHR11zh8yZAicnJywePFihIWF4ffff4e5uTmWLl0KoPrvKzw8HIGBgTA0NMR3330HPp+PjRs3on379rh48aLOoKIuPHnyBB988AHGjh2Lzz77DO7u7ggPD0fv3r3RtGlTzJs3D0KhEM+ePat0YFieDRs2wMvLC3379gWPx8OxY8cwfvx4aDQafPnllzrnPnv2DIMGDcKYMWMwatQobNmyBR9//DFatmwJLy8vAMCLFy/QoUMHlJaWYtq0aZBKpdi0aVOtZ3EODg44e/YsEhISYGtrW6e6qWudT5gwAQqFAnPmzMGTJ0+wYcMGxMbGsoNw4OU7HjVqFLp164alS5dCqVRiw4YNaNu2Le7cuVOhzdSW6tqWmZkZtm/fjoULFyI/Px+LFy8G8PL7JiL07dsX58+fx5gxY+Dj44NTp07h22+/RWJiYt3Ue1QH/vjjDwJQaSIiio6OJgD0xx9/sNfMnj2bymcjlUpp1KhRFe4/ZswYsrKyovT0dJ3jw4YNI7lcTkqlkoiIzp8/TwDI2dmZPUZEpNFoyM3Njbp160YajYY9rlQqycnJibp06cIe69+/P4lEIoqNjWWPPXr0iLhcboXyVsaoUaNIKpXqHHNwcCAAdPXqVfbYqVOnCACJxWKdvDZu3EgA6Pz58zrlLM+uXbsIAF26dIk91qdPH5JIJJSYmMgei4yMJB6Pp1P2mJgY4nK5tHDhQp17PnjwgHg8XoXj5enVqxf5+fmxfw8YMIAGDBhAXC6XTpw4QUREYWFhBICOHDnCnhcUFERBQUHs3zdv3qzQLsqeC4C2bdvGHisuLiZLS0saOHBgteUjelnnXbt2pbS0NEpLS6MHDx7QiBEjCAB9+eWX7HkhISEEgP766y+d60+ePKlzPDU1lQQCAfXq1UunDU2fPp0AVNpuy/L+++8TAMrKyqqx7EREPj4+ZG5uThkZGeyxe/fuEYfDoZEjR7LHtN/RJ598UiE/ExMTnWNVfV/9+/cngUBAUVFR7LGkpCSSyWTUrl27CnmVR/v9R0dHs8e0bf7kyZM6565atYoAUFpaWvUVUAmVfQfdunUjZ2dnnWPavMt+G6mpqSQUCmnKlCnssW+++YYA0PXr13XOk8vlFZ6nMjZv3kwASCAQUIcOHWjmzJkUEhJCarW6wrkAaPbs2ezfta1zbd22bNmSVCoVe3zZsmU631deXh4pFAr67LPPdPJ98eIFyeXyCsfLo+079+3bxx6rS9sKCgoiLy8vnWOHDx8mALRgwQKd44MGDSKGYejZs2fVlqks9VLZ/fLLLzhz5oxOelWICAcOHECfPn1AREhPT2dTt27
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "count    457.000000\n",
       "mean       0.103116\n",
       "std        0.320896\n",
       "min        0.000000\n",
       "25%        0.000000\n",
       "50%        0.000000\n",
       "75%        0.000000\n",
       "max        2.356194\n",
       "dtype: float64"
      ]
     },
     "execution_count": 187,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 187
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-08T16:50:32.975145Z",
     "start_time": "2024-05-08T16:50:32.971255Z"
    }
   },
   "cell_type": "code",
   "source": "curvatures",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0      0.0\n",
       "1      0.0\n",
       "2      0.0\n",
       "3      0.0\n",
       "4      0.0\n",
       "      ... \n",
       "352    0.0\n",
       "353    0.0\n",
       "354    0.0\n",
       "355    0.0\n",
       "356    0.0\n",
       "Length: 357, dtype: float64"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 174
  },
  {
   "metadata": {},
   "cell_type": "code",
   "outputs": [],
   "execution_count": null,
   "source": "[change for change in curvature_changes if change > threshold]"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-08T16:22:44.242067Z",
     "start_time": "2024-05-08T16:22:44.157432Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 假设 curvatures 是已经计算好的曲率列表\n",
    "length = len(curvatures)\n",
    "curvature_changes = [abs(curvatures[i] - curvatures[(i + 1) % length]) for i in range(length)]\n",
    "\n",
    "# 计算标准差和平均值\n",
    "mean_change = np.mean(curvature_changes)\n",
    "std_change = np.std(curvature_changes)\n",
    "threshold = mean_change + 1 * std_change  # 设定阈值为平均值加两倍标准差\n",
    "\n",
    "# 寻找异常值\n",
    "outliers = [i for i, change in enumerate(curvature_changes) if change > threshold]\n",
    "\n",
    "# 可视化结果\n",
    "print(len(curvatures))\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(curvatures, label='Curvature')\n",
    "plt.scatter(outliers, [curvatures[i] for i in outliers], color='red', label='Outliers')  # 标记异常点\n",
    "plt.legend()\n",
    "plt.title('Curvature and Outliers')\n",
    "plt.xlabel('Index')\n",
    "plt.ylabel('Curvature')\n",
    "plt.show()\n"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "63\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAIjCAYAAAA0vUuxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABopklEQVR4nO3deXiTVfr/8U+SNqFAaQHZKYssIqLIgAsyKjgsojIiKiD8hsVd64royLigM2rRUccNdNQR9CsI4gDuKIOyqCCyuSKyFItaQFFb1jZNnt8fJSGhabM9abb367pyDUmepKfNeE7uc+5zH4thGIYAAAAAANWyxrsBAAAAAJDoCJwAAAAAIAgCJwAAAAAIgsAJAAAAAIIgcAIAAACAIAicAAAAACAIAicAAAAACILACQAAAACCIHACAAAAgCAInAAAMJHFYtE999wT72YEtG3bNlksFs2YMcP72D333COLxRK/RgFAkiBwAoAEsmXLFl111VU6+uijVadOHTVo0EB9+vTR448/rgMHDsS7eTV65513EjZgSFS7d+/WrbfeqmOOOUZ16tRRo0aNNGjQIL311ltRve+sWbP02GOPmdNIAIAkAicASBhvv/22jj/+eL366qsaMmSInnzySRUUFKhNmza69dZbdeONN8a7iTV65513dO+998a7GUlj48aN6t69u5544gn169dPTz31lP72t79p165dGjJkiG699daI3zucwOnOO+9M+KAcABJBRrwbAACQCgsLNXLkSLVt21YffPCBWrRo4X0uPz9fmzdv1ttvv23Kz9q3b5/q1atnynvVhmRrbyicTqcuuugi/fbbb1q2bJlOOeUU73M333yzRo8erYcffli9evXSiBEjYtqWjIwMZWSY93Vg//79qlu3rmnvBwCJghUnAEgADz30kPbu3av//Oc/fkGTR8eOHb0rToH2qXgcub/Gs3/lm2++0ahRo9SwYUP98Y9/1MMPPyyLxaLvv/++yntMmjRJdrtdv/32myRp+fLluvjii9WmTRs5HA7l5eXp5ptv9lulGDdunKZOneptg+cmSUuWLJHFYtGSJUv8fk6g32PcuHGqX7++tmzZonPOOUfZ2dkaPXq0JMntduuxxx7Tcccdpzp16qhZs2a66qqrvO2syRdffKFx48Z5UyCbN2+uSy+9VLt37/a7zvP32rx5s8aNG6fc3Fzl5ORo/Pjx2r9/v9+1ZWVluvnmm9WkSRNlZ2frz3/+s3744YegbZGk//73v/rqq690++23+wVNkmSz2fTvf/9bubm5fp/ljBkzZLFYtG3bNr/rj/z79u3bV2+//ba+//577+fQrl27attS3R6nl19+WT179lRWVpYaNWqkkSNHavv27X7X9O3bV926ddOaNWt0xhlnqG7duvrb3/4mSVq9erUGDRqko446SllZWWrfvr0uvfTSkP4+AJCIWHECgATw5ptv6uijj9Zpp50Wk/e/+OKL1alTJz3wwAMyDEPnnXeebrvtNr366qtVUsJeffVVDRw4UA0bNpQkzZ07V/v379c111yjxo0ba9WqVXryySf1ww8/aO7cuZKkq666Sj/99JMWLVqk//u//4uqrRUVFRo0aJA3wPOsXlx11VWaMWOGxo8frxtuuEGFhYV66qmntG7dOn388cfKzMys9j0XLVqkrVu3avz48WrevLm+/vprPfvss/r666+1cuXKKoHD8OHD1b59exUUFGjt2rV6/vnn1bRpUz344IPeay6//HK9/PLLGjVqlE477TR98MEHOvfcc0P6Hd98801J0pgxYwI+n5OTo/PPP18vvviiNm/erI4dO4b0vpJ0xx13qKSkRD/88IP+9a9/SZLq168f8usl6f7779ddd92l4cOH6/LLL9fPP/+sJ598UmeccYbWrVun3Nxc77W7d+/W4MGDNXLkSP2///f/1KxZM+3atUsDBw5UkyZNdPvttys3N1fbtm3TvHnzwmoHACQUAwAQVyUlJYYk4/zzzw/p+sLCQkOSMX369CrPSTImT57svT958mRDknHJJZdUubZ3795Gz549/R5btWqVIcl46aWXvI/t37+/ymsLCgoMi8VifP/9997H8vPzjUDDyocffmhIMj788MOgv8fYsWMNScbtt9/ud+3y5csNScbMmTP9Hl+4cGHAx48U6Hd45ZVXDEnGsmXLvI95/l6XXnqp37UXXHCB0bhxY+/99evXG5KMa6+91u+6UaNGVfkMAjnxxBONnJycGq959NFHDUnGG2+8YRiGYUyfPt2QZBQWFvpdF+jve+655xpt27at8p6B/uae39lj27Zths1mM+6//36/13755ZdGRkaG3+NnnnmmIcl45pln/K6dP3++Icn47LPPavwdASCZkKoHAHFWWloqScrOzo7Zz7j66qurPDZixAitWbNGW7Zs8T42Z84cORwOnX/++d7HsrKyvP/et2+ffvnlF5122mkyDEPr1q2LSXuvueYav/tz585VTk6OBgwYoF9++cV769mzp+rXr68PP/ywxvfz/R0OHjyoX375Raeeeqokae3atVWuP/Lvdfrpp2v37t3ez+qdd96RJN1www1+1910000h/X579uwJ+nl7nvf8zNoyb948ud1uDR8+3O9v3bx5c3Xq1KnK39rhcGj8+PF+j3lWpN566y05nc7aajoAxBSBEwDEWYMGDSRVfpmOlfbt21d57OKLL5bVatWcOXMkSYZhaO7cuRo8eLC3TZJUVFSkcePGqVGjRqpfv76aNGmiM888U5JUUlJielszMjLUunVrv8c2bdqkkpISNW3aVE2aNPG77d27V7t27arxPX/99VfdeOONatasmbKystSkSRPv3yTQ79CmTRu/+560Rc9+qu+//15Wq1UdOnTwu+6YY44J6XfMzs4O+nl7no9lQB3Ipk2bZBiGOnXqVOVvvWHDhip/61atWslut/s9duaZZ+rCCy/Uvffeq6OOOkrnn3++pk+frrKystr8VQDAVOxxAoA4a9CggVq2bKmvvvoqpOurO6zU5XJV+xrfFRePli1b6vTTT9err76qv/3tb1q5cqWKior89vG4XC4NGDBAv/76q/7617+qS5cuqlevnn788UeNGzdObrfb9PY6HA5Zrf7zem63W02bNtXMmTMDvqZJkyY1tmH48OH65JNPdOutt+rEE09U/fr15Xa7dfbZZwf8HWw2W8D3MQyjxp8TqmOPPVbr169XUVFRlSDN44svvpAkde3aVVJkn3sk3G63LBaL3n333YB/hyP3SwX6/5bFYtFrr72mlStX6s0339R7772nSy+9VI888ohWrlwZ9p4rAEgEBE4AkADOO+88Pfvss1qxYoV69+5d47We1Y/ff//d7/FAFfKCGTFihK699lpt3LhRc+bMUd26dTVkyBDv819++aW+++47vfjii36FDBYtWlTlvar7Ym9Gezt06KD//e9/6tOnT8Av6jX57bfftHjxYt177726++67vY9v2rQprPfx1bZtW7ndbm3ZssVvlWnjxo0hvf68887TK6+8opdeekl33nlnledLS0v1+uuvq0uXLt7CEOH8Hav7LELRoUMHGYah9u3bq3PnzhG/jySdeuqpOvXUU3X//fdr1qxZGj16tGbPnq3LL788qvcFgHggVQ8AEsBtt92mevXq6fLLL9fOnTurPL9lyxY9/vjjkipXqI466igtW7bM75pp06aF/XMvvPBC2Ww2vfLKK5o7d67OO+88vzOTPCsOvisthmF42+LL87ojv9i3bdtWNpstqvYOHz5cLpdL//jHP6o8V1FRUeVn+gr0O0gK+YDYQAYPHixJeuKJJyJ6z4suukhdu3bVlClTtHr1ar/n3G63rrnmGv3222+aPHmy93FPWqDv39HlcunZZ5+t8v716tWLOI1y2LBhstlsuvfee6v8zQzDqFLCPZDffvutymtPPPFESSJdD0DSYsUJABJAhw4dNGvWLI0YMULHHnusxowZo27duqm8vFyffPKJ5s6dq3Hjxnmvv/zyyzVlyhRdfvnl6tWrl5YtW6bvvvsu7J/btGlT9evXT48++qj27Nl
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 98
  },
  {
   "metadata": {},
   "cell_type": "code",
   "outputs": [],
   "execution_count": null,
   "source": ""
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}