{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install /kaggle/input/mypackages/timm-0.3.2-py3-none-any.whl\n#!pip install /kaggle/input/mypackages/PyYAML-5.4.1-cp37-cp37m-manylinux1_x86_64.whl\n!pip install /kaggle/input/mypackages/Pillow-7.2.0-cp37-cp37m-manylinux1_x86_64.whl\n!pip install /kaggle/input/mypackages/thop-0.0.31.post2005241907-py3-none-any.whl","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nsys.path.insert(0, \"../input/yolov5\")\nsys.path.insert(1, \"../input/mypackages\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')\nimport argparse\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport numpy as np\nimport pandas as pd\nimport pathlib, sys, os, random, time\nimport numba, cv2, gc\nimport glob\nfrom PIL import Image\nfrom sklearn.model_selection import KFold\nimport matplotlib.pyplot as plt\nimport segmentation_models_pytorch as smp\nfrom tqdm import tqdm\nimport albumentations as A\nimport rasterio\nfrom rasterio.windows import Window\nimport torch\nimport torch.nn as nn\nfrom tqdm.notebook import tqdm\nimport torch.nn.functional as F\nimport torch.utils.data as D\nimport torchvision\nfrom torchvision import transforms as T\nfrom models.experimental import attempt_load\nfrom utils.datasets import LoadStreams, LoadImages\nfrom utils.general import check_img_size, check_requirements, check_imshow, non_max_suppression, apply_classifier, \\\n    scale_coords, xyxy2xywh, strip_optimizer, set_logging, increment_path\nfrom utils.plots import plot_one_box\nfrom utils.torch_utils import select_device, load_classifier, time_synchronized","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_PATH = '../input/hubmap-kidney-segmentation/'\np = pathlib.Path(DATA_PATH)\nDEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'\nimport segmentation_models_pytorch as smp\nprint(DEVICE)\nos.environ['CUDA_VISIBLE_DEVICES'] = '0'\nEPOCHES = 50\nBATCH_SIZE = 64\nWINDOW = 1024\nMIN_OVERLAP = 100\nNEW_SIZE = 512\nidentity = rasterio.Affine(1, 0, 0, 0, 1, 0)\nsubm = {}\nstart_time = time.time()\nconf_thres=0.25\niou_thres=0.5\nclasses=False\nagnostic_nms =False","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def rle_numba(pixels):\n    size = len(pixels)\n    points = []\n    if pixels[0] == 1: points.append(0)\n    flag = True\n    for i in range(1, size):\n        if pixels[i] != pixels[i - 1]:\n            if flag:\n                points.append(i + 1)\n                flag = False\n            else:\n                points.append(i + 1 - points[-1])\n                flag = True\n    if pixels[-1] == 1: points.append(size - points[-1] + 1)\n    return points\ndef rle_numba_encode(image):\n    pixels = image.flatten(order='F')\n    points = rle_numba(pixels)\n    return ' '.join(str(x) for x in points)\ndef make_grid(shape, window=256, min_overlap=32):\n    \"\"\"\n        Return Array of size (N,4), where N - number of tiles,\n        2nd axis represente slices: x1,x2,y1,y2\n    \"\"\"\n    x, y = shape\n    nx = x // (window - min_overlap) + 1\n    x1 = np.linspace(0, x, num=nx, endpoint=False, dtype=np.int64)\n    x1[-1] = x - window\n    x2 = (x1 + window).clip(0, x)\n    ny = y // (window - min_overlap) + 1\n    y1 = np.linspace(0, y, num=ny, endpoint=False, dtype=np.int64)\n    y1[-1] = y - window\n    y2 = (y1 + window).clip(0, y)\n    slices = np.zeros((nx, ny, 4), dtype=np.int64)\n\n    for i in range(nx):\n        for j in range(ny):\n            slices[i, j] = x1[i], x2[i], y1[j], y2[j]\n    return slices.reshape(nx * ny, 4)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trfm_512 = T.Compose([\n    T.ToPILImage(),\n    T.Resize(512),\n    T.ToTensor(),\n    T.Normalize([0.64254668, 0.49338829, 0.68244302],\n                [0.1927673, 0.25079806, 0.17813688]),\n])\n# trfm_384 = T.Compose([\n#     T.ToPILImage(),\n#     T.Resize(384),\n#     T.ToTensor(),\n#     T.Normalize([0.64254668, 0.49338829, 0.68244302],\n#                 [0.1927673, 0.25079806, 0.17813688]),\n# ])\nyolo_trfm = T.Compose([\n    T.ToPILImage(),\n    T.Resize(NEW_SIZE),\n    T.ToTensor(),\n    # T.Normalize([0.64254668, 0.49338829, 0.68244302],\n    #             [0.1927673,  0.25079806, 0.17813688]),\n])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"segmodel_512_1 = smp.Unet(\n    encoder_name=\"timm-efficientnet-b4\",  # choose encoder, e.g. mobilenet_v2 or efficientnet-b7  se_resnext50_32x4d\n    encoder_weights=None,  # use `imagenet` pretreined weights for encoder initialization\n    in_channels=3,  # model input channels (1 for grayscale images, 3 for RGB, etc.)\n    classes=1,  # model output channels (number of classes in your dataset)\n)\nsegmodel_512_1 = nn.DataParallel(segmodel_512_1)\nsegmodel_512_1.load_state_dict(torch.load(\"../input/hubmapmodel/out915_0.932_Unet_b4_v0428_fold1_epoch4_vdice0.949331.pth\"))\nsegmodel_512_1.to(DEVICE);\nsegmodel_512_1.eval()\n\n# segmodel_512_2 = smp.UnetPlusPlus(\n#     encoder_name=\"timm-efficientnet-b2\",  # choose encoder, e.g. mobilenet_v2 or efficientnet-b7  se_resnext50_32x4d\n#     encoder_weights=None,  # use `imagenet` pretreined weights for encoder initialization\n#     in_channels=3,  # model input channels (1 for grayscale images, 3 for RGB, etc.)\n#     classes=1,  # model output channels (number of classes in your dataset)\n# )\n# segmodel_512_2 = nn.DataParallel(segmodel_512_2)\n# segmodel_512_2.load_state_dict(torch.load(\"../input/hubmapmodel/out890_0.932_UNetPlusPluse_b2_v0429_fold0_epoch20_vdice0.949854.pth\"))\n# segmodel_512_2.to(DEVICE);\n# segmodel_512_2.eval()\n\n\n# segmodel_512_3 = smp.Unet(\n#     encoder_name=\"timm-efficientnet-b2\",  # choose encoder, e.g. mobilenet_v2 or efficientnet-b7  se_resnext50_32x4d\n#     encoder_weights=None,  # use `imagenet` pretreined weights for encoder initialization\n#     in_channels=3,  # model input channels (1 for grayscale images, 3 for RGB, etc.)\n#     classes=1,  # model output channels (number of classes in your dataset)\n# )\n# segmodel_512_3 = nn.DataParallel(segmodel_512_3)\n# segmodel_512_3.load_state_dict(torch.load(\"../input/hubmapmodel/out898_0.932_Unet_b2_v0426_fold4_epoch8_vdice0.955260.pth\"))\n# segmodel_512_3.to(DEVICE);\n# segmodel_512_3.eval()\n\nsegmodel_512_4 = smp.MAnet(\n    encoder_name=\"timm-efficientnet-b4\",  # choose encoder, e.g. mobilenet_v2 or efficientnet-b7  se_resnext50_32x4d\n    encoder_weights=None,  # use `imagenet` pretreined weights for encoder initialization\n    in_channels=3,  # model input channels (1 for grayscale images, 3 for RGB, etc.)\n    classes=1,  # model output channels (number of classes in your dataset)\n)\nsegmodel_512_4 = nn.DataParallel(segmodel_512_4)\nsegmodel_512_4.load_state_dict(torch.load(\"../input/hubmapmodel/out918_0.931_MAnet_b4_v0502_fold1_epoch4_vdice0.954013.pth\"))\nsegmodel_512_4.to(DEVICE);\nsegmodel_512_4.eval()\n\n\n# segmodel_384_1 = smp.Unet(\n#     encoder_name=\"timm-efficientnet-b3\",  # choose encoder, e.g. mobilenet_v2 or efficientnet-b7  se_resnext50_32x4d\n#     encoder_weights=None,  # use `imagenet` pretreined weights for encoder initialization\n#     in_channels=3,  # model input channels (1 for grayscale images, 3 for RGB, etc.)\n#     classes=1,  # model output channels (number of classes in your dataset)\n# )\n# segmodel_384_1.load_state_dict(torch.load(\"../input/hubmapnewmodel/timm_b3_all_384_fold_1_15.pth\"))\n# segmodel_384_1.to(DEVICE);\n# segmodel_384_1.eval()\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nsubm = {}\n''' yolo模型定义 '''\nweights = '/kaggle/input/yolov5/runs/train/exp37/weights/best.pt'\nimgsz = 512\n# Initialize\ndevice = select_device('0')\nhalf = device.type != 'cpu'  # half precision only supported on CUDA\n# Load model\nmodel_yolo = attempt_load(weights, map_location=device)  # load FP32 model\nstride = int(model_yolo.stride.max())  # model stride\nimgsz = check_img_size(imgsz, s=stride)  # check img_size\nif half:\n    model_yolo.half()  # to FP16\n# Run inference\nif device.type != 'cpu':\n    model_yolo(torch.zeros(1, 3, imgsz, imgsz).to(device).type_as(next(model_yolo.parameters())))  # run once\n\nprint('1')\n'''seg&yolo'''\nfor filei, filename in enumerate(p.glob('test/*.tiff')):\n    dataset = rasterio.open(filename.as_posix(), transform = identity)\n    slices = make_grid(dataset.shape, window=WINDOW, min_overlap=MIN_OVERLAP)\n    print(filename)\n    print(filei)\n    check = 0\n    yolotime = 0\n    change_seg=0\n#     aotu = 0\n\n\n    preds = np.zeros(dataset.shape, dtype=np.uint8)\n\n    for (x1, x2, y1, y2) in tqdm(slices):\n        '''读入一张切片'''\n        save_yolo = 0\n#         print(check)\n#         check +=1\n        if dataset.count == 3:\n            oimage = dataset.read([1, 2, 3],\n                                 window=Window.from_slices((x1, x2), (y1, y2)))\n            oimage = np.moveaxis(oimage, 0, -1)\n        else:\n            subdatasets = dataset.subdatasets\n            #     if len(subdatasets) > 0:\n            #  image = np.zeros(WINDOW, WINDOW, 3), dtype=np.uint8)\n            oimage = np.zeros((1024, 1024, 3), dtype=np.uint8)\n            for m, subdataset in enumerate(subdatasets, 0):\n                with rasterio.open(subdataset) as layer:\n                    oimage[:, :, m] = layer.read(1, window=Window.from_slices((x1, x2), (y1, y2)))\n\n        img_save = oimage\n\n        '''YOLO输入图片'''\n        image_yolo = yolo_trfm(oimage)\n        '''分割图片'''\n        image_512 = trfm_512(oimage)\n#         image_384 = trfm_384(oimage)\n\n        '''seg分割'''\n        with torch.no_grad():\n            image_512 = image_512.to(DEVICE)[None]\n#             image_384 = image_384.to(DEVICE)[None]\n\n            score_512_1 = segmodel_512_1(image_512)[0][0]\n#             score_512_2 = segmodel_512_2(image_512)[0][0]\n#             score_512_3 = segmodel_512_3(image_512)[0][0]\n            score_512_4 = segmodel_512_4(image_512)[0][0]    \n#             score_384_1 = segmodel_384_1(image_384)[0][0]\n\n\n            score_512_1 = score_512_1\n            score_512_1_sigmoid = score_512_1.sigmoid().cpu().numpy()\n            score_512_1_sigmoid = cv2.resize(score_512_1_sigmoid, (WINDOW, WINDOW))\n            \n#             score_512_2 = score_512_2\n#             score_512_2_sigmoid = score_512_2.sigmoid().cpu().numpy()\n#             score_512_2_sigmoid = cv2.resize(score_512_2_sigmoid, (WINDOW, WINDOW))\n            \n#             score_512_3 = score_512_3\n#             score_512_3_sigmoid = score_512_3.sigmoid().cpu().numpy()\n#             score_512_3_sigmoid = cv2.resize(score_512_3_sigmoid, (WINDOW, WINDOW))\n            \n            score_512_4 = score_512_4\n            score_512_4_sigmoid = score_512_4.sigmoid().cpu().numpy()\n            score_512_4_sigmoid = cv2.resize(score_512_4_sigmoid, (WINDOW, WINDOW))\n            \n#             score_384_1 = score_384_1\n#             score_384_1_sigmoid = score_384_1.sigmoid().cpu().numpy()\n#             score_384_1_sigmoid = cv2.resize(score_384_1_sigmoid, (WINDOW, WINDOW))\n\n\n    \n            score_mean = (score_512_1_sigmoid  + score_512_4_sigmoid ) / 2.0\n#             score_mean = (score_512_1_sigmoid + score_512_2_sigmoid +\n#                           score_512_3_sigmoid + score_512_4_sigmoid ) / 4.0\n\n            \n            '''填充空洞'''\n            score_sigmoid = ( score_mean > 0.4).astype(np.uint8)\n            if (score_sigmoid >0.5 ).mean() > 0.002:\n                contours, hierarchy = cv2.findContours(score_sigmoid, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)#检测轮廓\n                image_fill = np.zeros((1024,1024), dtype=np.uint8)\n                cv2.fillPoly(score_sigmoid, contours, 1)\n                area_ao = np.sum(np.greater(score_sigmoid, 0))#凹填充的面积计算\n                for c in range(len(contours)):  # 凸包检测并填充区域 \n                    points = cv2.convexHull(contours[c])\n                    cv2.fillConvexPoly(image_fill, points, 1)  \n                area_tu = np.sum(np.greater(image_fill, 0))  #凸填充面积计算\n                area_rate = float(area_ao / area_tu)\n            else :\n                area_rate = 1\n            preds[x1:x2, y1:y2] = (score_sigmoid > 0.5).astype(np.uint8)#凹填充的结果或者是没填充\n            # 取出单张切片，等下用来和检测结果比较\n            one_seg_pred = preds[x1:x2, y1:y2]\n\n        '''单张切片分割结果有标注才跑检测。也可以反过来先检测有检测结果再分割，不过我估计效果很差'''\n        if area_rate >0.6 and area_rate <0.79:\n            yolotime +=1\n            \n            image_fill = (image_fill > 0.5).astype(np.uint8)   #凸填充的结果\n            # new_seg_pred = one_seg_pred\n            img_yolo = image_yolo.to(device)[None]\n            img_yolo = img_yolo.half()\n            # 没有batch_size的话则在最前面添加一个轴\n            if img_yolo.ndimension() == 3:\n                img_yolo = img_yolo.unsqueeze(0)\n            # Inference\n            \"\"\"\n                    前向传播 返回pred的shape是(1, num_boxes, 5+num_class)\n                    h,w为传入网络图片的长和宽，注意dataset在检测时使用了矩形推理，所以这里h不一定等于w\n                    num_boxes = h/32 * w/32 + h/16 * w/16 + h/8 * w/8\n                    pred[..., 0:4]为预测框坐标\n                    预测框坐标为xywh(中心点+宽长)格式\n                    pred[..., 4]为objectness置信度\n                    pred[..., 5:-1]为分类结果\n            \"\"\"\n            augment =True\n            pred_yolo = model_yolo(img_yolo, augment=augment)[0]\n            # 进行NMS\n            \"\"\"\n            pred:前向传播的输出\n            conf_thres:置信度阈值\n            iou_thres:iou阈值\n            classes:是否只保留特定的类别\n            agnostic:进行nms是否也去除不同类别之间的框\n            经过nms之后，预测框格式：xywh-->xyxy(左上角右下角)\n            pred是一个列表list[torch.tensor]，长度为batch_size\n            \"\"\"\n            pred_yolo = non_max_suppression(pred_yolo, conf_thres, iou_thres, classes=classes,\n                                            agnostic=agnostic_nms)\n\n            # 依次导出每一个检测框\n            for i1, det in enumerate(pred_yolo):\n                # 判断是否有检测框。若无检测框，则不变\n                if len(det):\n                    # 检测框从512*512 到1024*1024\n                    det[:, :4] = scale_coords(img_yolo.shape[2:], det[:, :4], oimage.shape).round()\n                    box_x1 = int(det[0][0])\n                    box_x2 = int(det[0][2])\n                    box_y1 = int(det[0][1])\n                    box_y2 = int(det[0][3])\n                    \n                    #  这里可以检查x1,x2,y1,y2的位置，检测之前已经凸填充了\n\n                    # 画出检测框，\n                    box = image_fill[ box_y1 :box_y2,box_x1:box_x2] \n                    # 计算检测框内的分割面积\n                    seg_area=box.sum()\n                    # 检测框面积\n                    box_area = (box_x2-box_x1)*(box_y2-box_y1)\n                    # 分割面积和监测面积的比值\n                    area_rate = float(seg_area / box_area)\n                    # 比值小于某值时，填充\n                    if box_area > 5000 and area_rate<0.5 and area_rate>0.05:\n                        change_seg +=1\n \n                        # 填充方式:椭圆，\n                        x = box.shape[0]\n                        y = box.shape[1]\n\n                        #调用cv2的包，先绘制rgb的椭圆，再转成单通道\n                        fill = np.zeros((x,y,3),np.uint8)\n\n                        fill_circle = cv2.circle(img=fill, center=(int(y / 2), int(x / 2)),\n                                                 radius=min(int(x / 3), int(y / 3)),\n                                                 color=(1, 1, 1), thickness=-1)\n                        fill_one = cv2.cvtColor(fill_circle, cv2.COLOR_RGB2GRAY)\n\n                        new_pred = fill_one+box\n                        new_pred = (new_pred > 0.5).astype(np.uint8)\n\n                        contours, hierarchy = cv2.findContours(new_pred, cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_NONE)  # 检测轮廓\n                        for c in range(len(contours)):  # 凸包检测并填充区域\n                            points = cv2.convexHull(contours[c])\n                            cv2.fillConvexPoly(new_pred, points, 1)  # 填充轮廓\n\n                        image_fill[box_y1:box_y2, box_x1:box_x2] = (new_pred > 0.5).astype(np.uint8)\n\n\n                # 将新填充结果保存到最大的切片preds\n            preds[x1:x2, y1:y2] = image_fill  #有可能没有目标直接执行这一行，有可能只是凸填充的结果\n       \n\n    print(filename.stem,'yolo time:', yolotime )\n    print(filename.stem,'change seg time:', change_seg )\n\n        # break\n\n    #生成csv\n    subm[filei] = {'id': filename.stem, 'predicted': rle_numba_encode(preds)}\n\n    del preds\n    gc.collect();\n    \nprint('time:', (time.time() - start_time) // 60, 'min', (time.time() - start_time) % 60, 's')\nsubmission = pd.DataFrame.from_dict(subm, orient='index')\nsubmission.to_csv('./submission.csv', index=False)\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}