BA-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xlwt\n",
    "from scipy import integrate\n",
    "from skimage import measure,data,color\n",
    "\n",
    "\n",
    "def cal_area(src):\n",
    "    src[src > 0] = 1\n",
    "    src_area = np.sum(src)\n",
    "    return src_area\n",
    "\n",
    "def detect_peaks(x, mph=None, mpd=1, threshold=0, edge='rising',\n",
    "                 kpsh=False, valley=False, show=False, ax=None):\n",
    "    \"\"\"Detect peaks in data based on their amplitude and other features.\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : 1D array_like\n",
    "        data.\n",
    "    mph : {None, number}, optional (default = None)\n",
    "        detect peaks that are greater than minimum peak height.\n",
    "    mpd : positive integer, optional (default = 1)\n",
    "        detect peaks that are at least separated by minimum peak distance (in\n",
    "        number of data).\n",
    "    threshold : positive number, optional (default = 0)\n",
    "        detect peaks (valleys) that are greater (smaller) than `threshold`\n",
    "        in relation to their immediate neighbors.\n",
    "    edge : {None, 'rising', 'falling', 'both'}, optional (default = 'rising')\n",
    "        for a flat peak, keep only the rising edge ('rising'), only the\n",
    "        falling edge ('falling'), both edges ('both'), or don't detect a\n",
    "        flat peak (None).\n",
    "    kpsh : bool, optional (default = False)\n",
    "        keep peaks with same height even if they are closer than `mpd`.\n",
    "    valley : bool, optional (default = False)\n",
    "        if True (1), detect valleys (local minima) instead of peaks.\n",
    "    show : bool, optional (default = False)\n",
    "        if True (1), plot data in matplotlib figure.\n",
    "    ax : a matplotlib.axes.Axes instance, optional (default = None).\n",
    "    Returns\n",
    "    -------\n",
    "    ind : 1D array_like\n",
    "        indeces of the peaks in `x`.\n",
    "    Notes\n",
    "    -----\n",
    "    The detection of valleys instead of peaks is performed internally by simply\n",
    "    negating the data: `ind_valleys = detect_peaks(-x)`\n",
    "\n",
    "    The function can handle NaN's\n",
    "    See this IPython Notebook [1]_.\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/DetectPeaks.ipynb\n",
    "    Examples\n",
    "    --------\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    x = np.atleast_1d(x).astype('float64')   # 将输入的数据转换成1维数组，浮点型\n",
    "    if x.size < 3:          # 如果数组尺寸小于3，则返回空，表示没有波峰或波谷\n",
    "        return np.array([], dtype=int)\n",
    "    if valley:      # 表示当求波谷时，先将所有数据取反\n",
    "        x = -x\n",
    "    # find indexes of all peaks\n",
    "    dx = x[1:] - x[:-1]\n",
    "    # handle NaN's\n",
    "    indnan = np.where(np.isnan(x))[0]\n",
    "    if indnan.size:\n",
    "        x[indnan] = np.inf\n",
    "        dx[np.where(np.isnan(dx))[0]] = np.inf\n",
    "    ine, ire, ife = np.array([[], [], []], dtype=int)\n",
    "    if not edge:\n",
    "        ine = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) > 0))[0]\n",
    "    else:\n",
    "        if edge.lower() in ['rising', 'both']:\n",
    "            ire = np.where((np.hstack((dx, 0)) <= 0) & (np.hstack((0, dx)) > 0))[0]\n",
    "        if edge.lower() in ['falling', 'both']:\n",
    "            ife = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) >= 0))[0]\n",
    "    ind = np.unique(np.hstack((ine, ire, ife)))\n",
    "    # handle NaN's\n",
    "    if ind.size and indnan.size:\n",
    "        # NaN's and values close to NaN's cannot be peaks\n",
    "        ind = ind[np.in1d(ind, np.unique(np.hstack((indnan, indnan - 1, indnan + 1))), invert=True)]\n",
    "    # first and last values of x cannot be peaks\n",
    "    if ind.size and ind[0] == 0:\n",
    "        ind = ind[1:]\n",
    "    if ind.size and ind[-1] == x.size - 1:\n",
    "        ind = ind[:-1]\n",
    "    # remove peaks < minimum peak height\n",
    "    if ind.size and mph is not None:\n",
    "        ind = ind[x[ind] >= mph]\n",
    "    # remove peaks - neighbors < threshold\n",
    "    if ind.size and threshold > 0:\n",
    "        dx = np.min(np.vstack([x[ind] - x[ind - 1], x[ind] - x[ind + 1]]), axis=0)\n",
    "        ind = np.delete(ind, np.where(dx < threshold)[0])\n",
    "    # detect small peaks closer than minimum peak distance\n",
    "    if ind.size and mpd > 1:\n",
    "        ind = ind[np.argsort(x[ind])][::-1]  # sort ind by peak height\n",
    "        idel = np.zeros(ind.size, dtype=bool)\n",
    "        for i in range(ind.size):\n",
    "            if not idel[i]:\n",
    "                # keep peaks with the same height if kpsh is True\n",
    "                idel = idel | (ind >= ind[i] - mpd) & (ind <= ind[i] + mpd) \\\n",
    "                       & (x[ind[i]] > x[ind] if kpsh else True)\n",
    "                idel[i] = 0  # Keep current peak\n",
    "        # remove the small peaks and sort back the indexes by their occurrence\n",
    "        ind = np.sort(ind[~idel])\n",
    "\n",
    "    if show:\n",
    "        if indnan.size:\n",
    "            x[indnan] = np.nan\n",
    "        if valley:\n",
    "            x = -x\n",
    "        _plot(x, mph, mpd, threshold, edge, valley, ax, ind)\n",
    "\n",
    "    return ind\n",
    "\n",
    "def rotation2(src):\n",
    "    cnt = get_cnt(src)\n",
    "    rect = cv2.minAreaRect(cnt)\n",
    "    width = int(rect[1][0])\n",
    "    height = int(rect[1][1])\n",
    "    angle = rect[2]\n",
    "\n",
    "    if width < height:  # calculate the angle, prepare for next step\n",
    "        angle = angle - 90\n",
    "\n",
    "    src_pts = cv2.boxPoints(rect) # 使用cv2.boxPoints()可获取该矩形的四个顶点坐标\n",
    "    src_pts = src_pts.astype('float32')\n",
    "    dst_pts = np.array([[0, height],\n",
    "                        [0, 0],\n",
    "                        [width, 0],\n",
    "                        [width, height]], dtype=\"float32\")\n",
    "    matrix = cv2.getPerspectiveTransform(src_pts, dst_pts)   # 返回由源图像中矩形到目标图像矩形变换的矩阵，src：源图像中待测矩形的四点坐标，sdt：目标图像中矩形的四点坐标\n",
    "    warped = cv2.warpPerspective(src, matrix, (width, height)) # 参数为：src，输入图像；M：变换矩阵；dsize：目标图像shape\n",
    "\n",
    "    if angle <= -90:  # rotate for -90 degrees images    # 将图像旋转90度\n",
    "        warped = cv2.transpose(warped)  # 转置图像\n",
    "        warped = cv2.flip(warped, 1)  # anticlockwise 90 degrees, if u want clockwise, turn 0 to 1逆时针旋转90度，如果您想顺时针旋转，则将0旋转至1\n",
    "        warped = np.rot90(warped, k=3) #将矩阵A逆时针旋转（90×k）°以后返回B，k取负数时表示顺时针旋转。\n",
    "    else:\n",
    "        warped = np.rot90(warped, k=1)\n",
    "    return warped\n",
    "\n",
    "def get_cnt(src):\n",
    "    #gray = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)\n",
    "    ret, binary = cv2.threshold(src, 0, 255, cv2.THRESH_BINARY)\n",
    "    if np.sum(binary) == 0:\n",
    "        cnt = np.array([[[0,0]]])\n",
    "    else:\n",
    "        contours, hierarchy = cv2.findContours(binary, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n",
    "        if contours == []:\n",
    "            cnt = np.array([[[0,0]]])\n",
    "        else:\n",
    "            cnt = contours[0]\n",
    "    return cnt\n",
    "\n",
    "def minmax(area):\n",
    "    maxi = max(area)\n",
    "    mini = min(area)\n",
    "    return mini,maxi\n",
    "\n",
    "def fill_contour(img):\n",
    "    contours, _ = cv2.findContours(img, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n",
    "    n = len(contours)  # 轮廓的个数\n",
    "    max_area = 0\n",
    "    for i  in range(n):\n",
    "        if cv2.contourArea(contours[i]) > max_area:\n",
    "            max_area = cv2.contourArea(contours[i])\n",
    "    cv_contours = []\n",
    "    if n != 1:\n",
    "        for contour in contours:\n",
    "            area = cv2.contourArea(contour)\n",
    "            if area < max_area:\n",
    "                cv_contours.append(contour)\n",
    "                x, y, w, h = cv2.boundingRect(contour)\n",
    "                img[y:y + h, x:x + w] = 1\n",
    "            else:\n",
    "                continue\n",
    "    else:\n",
    "        pass\n",
    "    \n",
    "    return img\n",
    "def roi_dect(cav_npy):\n",
    "    frame = cav_npy[:,:,:]\n",
    "    frame  = np.concatenate((frame , frame , frame ), axis=-1)\n",
    "    #frame = cv2.resize(frame, (shape0,shape1))\n",
    "    frame = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)\n",
    "    #plt.imshow(frame)\n",
    "    #print(frame2.shape)\n",
    "    frame = np.array(frame,np.uint8)\n",
    "    frame_detect = rotation2(frame)\n",
    "    #plt.figure()\n",
    "    #plt.imshow(frame_detect)\n",
    "    return frame_detect\n",
    "\n",
    "def lvd_dect(cav_npy):\n",
    "    frame = cav_npy[:, :]\n",
    "    frame  = np.concatenate((frame , frame , frame ), axis=-1)\n",
    "    #print(frame.shape)\n",
    "    #frame = cv2.resize(frame, (shape0,shape1))\n",
    "    frame = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)\n",
    "    #plt.imshow(frame)\n",
    "    #print(frame2.shape)\n",
    "    frame = np.array(frame,np.uint8)\n",
    "    frame_detect = rotation2(frame)\n",
    "    #plt.figure()\n",
    "    #plt.imshow(frame_detect)\n",
    "    return frame_detect\n",
    "\n",
    "def cal_keyframe(srcnpy,i):\n",
    "    #srcnpy = np.load(srcnpy_path)[:, :, :, :, np.newaxis]\n",
    "    #for i in range(srcnpy.shape[0]):\n",
    "    area = []\n",
    "    lvd = []\n",
    "    volume = []\n",
    "    for k in range(srcnpy.shape[1]):\n",
    "        area1 = cal_area(srcnpy[i,k,:,:,:])\n",
    "        area.append(area1)\n",
    "        img = srcnpy[i,k,:,:,:]\n",
    "        img = np.array(img,np.uint8)\n",
    "        img = fill_contour(img)\n",
    "        l1 = int(max(img.sum(axis=0)))\n",
    "        #print(l1)\n",
    "        #lvd1 = lvd_dect(img).shape[0]\n",
    "        lvd1 = l1\n",
    "        lvd.append(lvd1)\n",
    "        volume1 = 2/1000*0.85*(area1**2)/lvd1\n",
    "        volume.append(volume1)\n",
    "    area = np.array(area)\n",
    "    area_avg = np.mean(area)\n",
    "    peak_h = detect_peaks(area,mph=np.mean(area), mpd=len(area)/5)\n",
    "    peak_l = detect_peaks(area, mpd=len(area)/4, valley=True, threshold=-(2*np.mean(area))/3)\n",
    "    minvol,maxvol = minmax(volume)\n",
    "    ef = (maxvol-minvol)/maxvol\n",
    "    peak_list = np.zeros(srcnpy.shape[1])\n",
    "    #print(peak_h)\n",
    "    for i in range(len(peak_h)):\n",
    "        peak_list[peak_h[i]] = 1\n",
    "    for i in range(len(peak_l)):\n",
    "        peak_list[peak_l[i]] = -1\n",
    "    #print(peak_list)\n",
    "    return area,volume,peak_h,peak_l,peak_list,lvd,minvol,maxvol,ef\n",
    "\n",
    "def plot_peaks(x, indexes, algorithm=None, mph=None, mpd=None):\n",
    "    \"\"\"Plot results of the peak dectection.\"\"\"\n",
    "    try:\n",
    "        import matplotlib.pyplot as plt\n",
    "    except ImportError as e:\n",
    "        print('matplotlib is not available.', e)\n",
    "        return\n",
    "    _, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
    "    ax.plot(x, 'b', lw=1)\n",
    "    if indexes.size:\n",
    "        label = 'peak'\n",
    "        label = label + 's' if indexes.size > 1 else label\n",
    "        ax.plot(indexes, x[indexes], '+', mfc=None, mec='r', mew=2, ms=8,\n",
    "                label='%d %s' % (indexes.size, label))\n",
    "        ax.legend(loc='best', framealpha=.5, numpoints=1)\n",
    "    ax.set_xlim(-.02*x.size, x.size*1.02-1)\n",
    "    ymin, ymax = x[np.isfinite(x)].min(), x[np.isfinite(x)].max()\n",
    "    yrange = ymax - ymin if ymax > ymin else 1\n",
    "    # ax.set_ylim(ymin - 0.1*yrange, ymax + 0.1*yrange)\n",
    "#     ax.set_xlabel('frame', fontsize=14)\n",
    "#     ax.set_ylabel('area', fontsize=14)\n",
    "    # ax.set_title('%s (mph=%s, mpd=%s)' % (algorithm, mph, mpd))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "road = '/data/zhangzhenxuan/zhangzhenxuan/zhangzhenxuan'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1990, 30, 128, 128, 1) (1990, 30, 128, 128, 1)\n"
     ]
    }
   ],
   "source": [
    "# path1 = road+'/MTL_LOSS/multi-task/nature_1_pred.npy'\n",
    "# prdnpy = np.load(path1).reshape(1990, 30, 128, 128, 1)\n",
    "path1 = road+'/MTL_LOSS/multi-task/nature_2_pred.npy'\n",
    "prdnpy = np.load(path1).reshape(1990, 30, 128, 128, 1)\n",
    "path2 = '/data/zhangzhenxuan/nature_data/gts_a4c_1.npy'\n",
    "srcnpy = np.load(path2)[:, :, :, :, np.newaxis]\n",
    "print(prdnpy.shape,srcnpy.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "peak_list_all = []\n",
    "ef_pr = []\n",
    "ef_gt = []\n",
    "edv_pr = []\n",
    "edv_gt = []\n",
    "esv_pr = []\n",
    "esv_gt = []\n",
    "area_pr = []\n",
    "area_gt = []\n",
    "vol_pr = []\n",
    "vol_gt = []\n",
    "for i in range(0,1000):\n",
    "    area1,vol1,peak_h1,peak_l1,peak_list1,lvd1,minvol1,maxvol1,ef1 = cal_keyframe(prdnpy,i)\n",
    "    area2,vol2,peak_h2,peak_l2,peak_list2,lvd2,minvol2,maxvol2,ef2 = cal_keyframe(srcnpy,i)\n",
    "    print(ef1,ef2)\n",
    "    print(minvol1,maxvol1,minvol2,maxvol2)\n",
    "    ef_pr.append(ef1)\n",
    "    ef_gt.append(ef2)\n",
    "    edv_pr.append(maxvol1)\n",
    "    edv_gt.append(maxvol2)\n",
    "    esv_pr.append(minvol1)\n",
    "    esv_gt.append(minvol2)\n",
    "    area_pr.append(area1)\n",
    "    area_gt.append(area2)\n",
    "    vol_pr.append(vol1)\n",
    "    vol_gt.append(vol2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vol_pr = np.array(vol_pr)\n",
    "vol_gt = np.array(vol_gt)\n",
    "area_pr = np.array(area_pr)/100\n",
    "area_gt = np.array(area_gt)/100\n",
    "print(vol_pr.shape,vol_gt.shape)\n",
    "print(vol_gt.mean(axis = 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bland(data1,data2, *args, **kwargs):\n",
    "    data1=np.asarray(data1.reshape(-1))\n",
    "    data2=np.asarray(data2.reshape(-1))\n",
    "    #print(data1)\n",
    "    mean=np.mean([data1,data2],axis = 0)\n",
    "    diff=data1-data2\n",
    "#     print(diff)\n",
    "    md=np.mean(diff)\n",
    "    zero = 0\n",
    "    sd=np.std(diff,axis = 0)\n",
    "#     print(sd)\n",
    "    #plt.figure(figsize=(7, 7))\n",
    "    plt.figure(dpi=300,figsize=(7, 7))\n",
    "    # plt.scatter(mean,diff,marker='o',s=22,c='#E5E12F',edgecolor='#6F0707',*args,** kwargs)\n",
    "    plt.scatter(mean,diff,marker='o',s=62,c='#273B66',edgecolor='#23A4E7',*args,** kwargs)\n",
    "    plt.axhline(md,color = 'black',linewidth = 8,linestyle = '--',label='md')\n",
    "    #plt.axhline(zero,color = 'gray',linestyle = '--',label='md')\n",
    "    plt.axhline(md+1.96 * sd,color = 'crimson',linewidth = 8,linestyle = '-.',label='md+1.96*sd')\n",
    "    plt.axhline(md-1.96 * sd,color = 'crimson',linewidth = 8,linestyle = '-.',label='md-1.96*sd')\n",
    "    #plt.legend()\n",
    "    plt.text(x=max(mean), y=(md+1.96 * sd)+sd*0.1,s=str(round(md+1.96 * sd,2)),fontsize=15)\n",
    "    plt.text(x=max(mean), y=md+sd*0.1,s=str(round(md,2)),fontsize=15)\n",
    "    plt.text(x=max(mean), y=(md-1.96 * sd)+sd*0.1,s=str(round(md-1.96 * sd,2)),fontsize=15)\n",
    "    x = np.linspace(min(mean)-20,max(mean)+20,200)\n",
    "    y1 = md+1.96 * sd * np.ones(200)\n",
    "    y2 = md-1.96 * sd * np.ones(200)\n",
    "    #plt.fill_between(x, y1,y2,color = 'crimson', alpha=0.05)\n",
    "    #plt.ylabel('Difference between groundtruth ESV \\n and predicted ESV',fontsize=15)\n",
    "    #plt.xlabel('Average between groundtruth ESV and predicted ESV',fontsize=15)\n",
    "    plt.ylim([-(1.96 * sd)*2,(1.96 * sd)*2])\n",
    "    plt.xlim([min(mean)-(1/5)*max(mean),max(mean)+(1/5)*max(mean)])\n",
    "    #plt.title('Bland-Altman Analyse',fontsize=25)\n",
    "    bwith = 2 #边框宽度设置为2\n",
    "    TK = plt.gca()#获取边框\n",
    "    TK.spines['bottom'].set_linewidth(bwith)#图框下边\n",
    "    TK.spines['left'].set_linewidth(bwith)#图框左边\n",
    "    TK.spines['top'].set_linewidth(bwith)#图框上边\n",
    "    TK.spines['right'].set_linewidth(bwith)#图框右边\n",
    "    plt.savefig('BA-'+indic+'_deepblue.png', format='png', transparent=False, dpi=300)\n",
    "    #plt.show()\n",
    "    print(md+1.96 * sd,md-1.96 * sd,md)\n",
    "indic = 'ef'\n",
    "if indic == 'ef':\n",
    "    bland(np.array(ef_pr),np.array(ef_gt))\n",
    "elif indic == 'esv':\n",
    "    bland(np.array(esv_pr),np.array(esv_gt))\n",
    "elif indic == 'edv':\n",
    "    bland(np.array(edv_pr),np.array(edv_gt))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import gaussian_kde\n",
    "import matplotlib.colors as colors\n",
    "def bland(data1,data2, *args, **kwargs):\n",
    "    data1=np.asarray(data1.reshape(-1))\n",
    "    data2=np.asarray(data2.reshape(-1))\n",
    "    #print(data1)\n",
    "    mean=np.mean([data1,data2],axis = 0)\n",
    "    diff=data1-data2\n",
    "#     print(diff)\n",
    "    md=np.mean(diff)\n",
    "    zero = 0\n",
    "    sd=np.std(diff,axis = 0)\n",
    "#     print(sd)\n",
    "    #plt.figure(figsize=(7, 7))\n",
    "    plt.figure(dpi=300,figsize=(9, 7))\n",
    "    # plt.scatter(mean,diff,marker='o',s=22,c='#E5E12F',edgecolor='#6F0707',*args,** kwargs)\n",
    "    # plt.scatter(mean,diff,marker='o',s=62,c='#273B66',edgecolor='#23A4E7',*args,** kwargs)\n",
    "    xy = np.vstack([mean,diff])\n",
    "    z = gaussian_kde(xy)(xy)\n",
    "    idx = z.argsort()\n",
    "    mean, diff, z = mean[idx], diff[idx], z[idx]\n",
    "    # mean_m, diff_m = np.meshgrid(mean, diff)\n",
    "    # z = [z, z]\n",
    "    scatter = plt.scatter(mean,diff,marker='o',s=62,c=z**(1/8),cmap='PuBu',*args,** kwargs)\n",
    "    # scatter = plt.pcolormesh(mean_m,diff_m, z, cmap ='blues', vmin = z_min, vmax = z_max)\n",
    "    cbar=plt.colorbar(scatter,shrink=1,orientation='vertical',extend='both',pad=0.015,aspect=30)\n",
    "    # plt.pcolormesh(mean,diff, z,norm=colors.SymLogNorm(linthresh=0.4, linscale=1,vmin=0, vmax=8, base=10))\n",
    "    plt.axhline(md,color = 'black',linewidth = 8,linestyle = '--',label='md')\n",
    "    #plt.axhline(zero,color = 'gray',linestyle = '--',label='md')\n",
    "    plt.axhline(md+1.96 * sd,color = 'crimson',linewidth = 8,linestyle = '-.',label='md+1.96*sd')\n",
    "    plt.axhline(md-1.96 * sd,color = 'crimson',linewidth = 8,linestyle = '-.',label='md-1.96*sd')\n",
    "    #plt.legend()\n",
    "    plt.text(x=max(mean), y=(md+1.96 * sd)+sd*0.1,s=str(round(md+1.96 * sd,2)),fontsize=15)\n",
    "    plt.text(x=max(mean), y=md+sd*0.1,s=str(round(md,2)),fontsize=15)\n",
    "    plt.text(x=max(mean), y=(md-1.96 * sd)+sd*0.1,s=str(round(md-1.96 * sd,2)),fontsize=15)\n",
    "    x = np.linspace(min(mean)-20,max(mean)+20,200)\n",
    "    y1 = md+1.96 * sd * np.ones(200)\n",
    "    y2 = md-1.96 * sd * np.ones(200)\n",
    "    #plt.fill_between(x, y1,y2,color = 'crimson', alpha=0.05)\n",
    "    #plt.ylabel('Difference between groundtruth ESV \\n and predicted ESV',fontsize=15)\n",
    "    #plt.xlabel('Average between groundtruth ESV and predicted ESV',fontsize=15)\n",
    "    plt.ylim([-(1.96 * sd)*2,(1.96 * sd)*2])\n",
    "    plt.xlim([min(mean)-(1/5)*max(mean),max(mean)+(1/5)*max(mean)])\n",
    "    #plt.title('Bland-Altman Analyse',fontsize=25)\n",
    "    bwith = 2 #边框宽度设置为2\n",
    "    TK = plt.gca()#获取边框\n",
    "    TK.spines['bottom'].set_linewidth(bwith)#图框下边\n",
    "    TK.spines['left'].set_linewidth(bwith)#图框左边\n",
    "    TK.spines['top'].set_linewidth(bwith)#图框上边\n",
    "    TK.spines['right'].set_linewidth(bwith)#图框右边\n",
    "    plt.savefig('BA_gau-'+indic+'_deepblue.png', format='png', transparent=False, dpi=300)\n",
    "    #plt.show()\n",
    "    print(md+1.96 * sd,md-1.96 * sd,md)\n",
    "indic = 'ef'\n",
    "if indic == 'ef':\n",
    "    bland(np.array(ef_pr),np.array(ef_gt))\n",
    "elif indic == 'esv':\n",
    "    bland(np.array(esv_pr),np.array(esv_gt))\n",
    "elif indic == 'edv':\n",
    "    bland(np.array(edv_pr),np.array(edv_gt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "data_frame=pd.read_excel('./ef.xls')\n",
    "# indic = 'esv'\n",
    "a_ef = data_frame[indic+'gt']\n",
    "b_ef = data_frame[indic+'pr']\n",
    "#c = data_frame['Bedside EF']\n",
    "a_ef = np.array(a_ef)\n",
    "b_ef = np.array(b_ef)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "if indic == 'ef':\n",
    "    a,b = np.array(ef_pr)*100,np.array(ef_gt)*100\n",
    "elif indic == 'esv':\n",
    "    a,b = np.array(esv_pr),np.array(esv_gt)\n",
    "elif indic == 'edv':\n",
    "    a,b = np.array(edv_pr),np.array(edv_gt)\n",
    "plt.figure(dpi=300,figsize=(7, 7))\n",
    "plt.plot(np.arange(0, max(a)+(1/10)*max(a)),np.arange(0, max(a)+(1/10)*max(a)),color = 'k',linewidth = 8,linestyle = '-.')\n",
    "plt.scatter(abs(a),abs(b),marker='o',s=62,c='#273B66',edgecolor='#23A4E7')\n",
    "# plt.scatter(abs(a_ef),abs(b_ef),marker='o',s=18,c='magenta',edgecolor='b')\n",
    "plt.scatter(abs(a_ef),abs(b_ef),marker='o',s=62,c='#273B66',edgecolor='#23A4E7')\n",
    "y_est = a\n",
    "y_err = a.std() * np.sqrt(1/len(a) +(a - a.mean())**2 / np.sum((a - a.mean())**2))\n",
    "y_est = np.array(np.sort(y_est))\n",
    "y_err = np.array(np.sort(y_err))\n",
    "plt.ylim([0,100])\n",
    "plt.xlim([0,100])\n",
    "bwith = 2 #边框宽度设置为2\n",
    "TK = plt.gca()#获取边框\n",
    "TK.spines['bottom'].set_linewidth(bwith)#图框下边\n",
    "TK.spines['left'].set_linewidth(bwith)#图框左边\n",
    "TK.spines['top'].set_linewidth(bwith)#图框上边\n",
    "TK.spines['right'].set_linewidth(bwith)#图框右边\n",
    "plt.savefig('cor'+indic+'_deepblue.png', format='png', transparent=False, dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "indic = 'ef'\n",
    "if indic == 'ef':\n",
    "    a,b = np.array(ef_pr)*100,np.array(ef_gt)*100\n",
    "elif indic == 'esv':\n",
    "    a,b = np.array(esv_pr),np.array(esv_gt)\n",
    "elif indic == 'edv':\n",
    "    a,b = np.array(edv_pr),np.array(edv_gt)\n",
    "plt.figure(dpi=300,figsize=(9, 7))\n",
    "plt.plot(np.arange(0, max(a)+(1/10)*max(a)),np.arange(0, max(a)+(1/10)*max(a)),color = 'k',linewidth = 8,linestyle = '-.')\n",
    "# plt.scatter(abs(a),abs(b),marker='o',s=62,c='#273B66',edgecolor='#23A4E7')\n",
    "# plt.scatter(abs(a_ef),abs(b_ef),marker='o',s=62,c='#273B66',edgecolor='#23A4E7')\n",
    "a = np.concatenate((a,a_ef),axis=0)\n",
    "b = np.concatenate((b,b_ef),axis=0)\n",
    "xy = np.vstack([a,b])\n",
    "z = gaussian_kde(xy)(xy)\n",
    "idx = z.argsort()\n",
    "a,b, z = a[idx], b[idx], z[idx]\n",
    "# mean_m, diff_m = np.meshgrid(mean, diff)\n",
    "# z = [z, z]\n",
    "scatter = plt.scatter(a,b,marker='o',s=62,c=z**(1/8),cmap='PuBu')\n",
    "# scatter = plt.pcolormesh(mean_m,diff_m, z, cmap ='blues', vmin = z_min, vmax = z_max)\n",
    "cbar=plt.colorbar(scatter,shrink=1,orientation='vertical',extend='both',pad=0.015,aspect=30)\n",
    "y_est = a\n",
    "y_err = a.std() * np.sqrt(1/len(a) +(a - a.mean())**2 / np.sum((a - a.mean())**2))\n",
    "y_est = np.array(np.sort(y_est))\n",
    "y_err = np.array(np.sort(y_err))\n",
    "plt.ylim([0,100])\n",
    "plt.xlim([0,100])\n",
    "bwith = 2 #边框宽度设置为2\n",
    "TK = plt.gca()#获取边框\n",
    "TK.spines['bottom'].set_linewidth(bwith)#图框下边\n",
    "TK.spines['left'].set_linewidth(bwith)#图框左边\n",
    "TK.spines['top'].set_linewidth(bwith)#图框上边\n",
    "TK.spines['right'].set_linewidth(bwith)#图框右边\n",
    "#print(y_est.shape,y_err.shape)\n",
    "#plt.plot(y_est,y_est + y_err,color = 'crimson',linewidth = 1,linestyle = '-.')\n",
    "#plt.plot(y_est,y_est - y_err,color = 'crimson',linewidth = 1,linestyle = '-.')\n",
    "#plt.fill_between(y_est, y_est - y_err, y_est + y_err, color = 'crimson', alpha=0.05)\n",
    "#plt.ylabel('Predicted ESV',fontsize=15)\n",
    "#plt.xlabel('GroundTruth ESV',fontsize=15)\n",
    "#plt.title('Correlation Analyse',fontsize=25)\n",
    "plt.savefig('cor_gau'+indic+'_deepblue.png', format='png', transparent=False, dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}