{"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":"#from __future__ import print_function\nimport argparse\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim\nfrom torchvision import datasets, transforms\nfrom torch.autograd import Variable\nimport cv2\nimport sys\nimport numpy as np\nfrom skimage import segmentation\nimport torch.nn.init\n\nimport pydicom\nimport matplotlib.pyplot as plt\nfrom skimage.transform import resize\nfrom skimage import measure\nfrom scipy.fft import fft\nfrom scipy.optimize import curve_fit","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-06-24T13:46:13.416149Z","iopub.execute_input":"2021-06-24T13:46:13.416462Z","iopub.status.idle":"2021-06-24T13:46:13.423744Z","shell.execute_reply.started":"2021-06-24T13:46:13.416432Z","shell.execute_reply":"2021-06-24T13:46:13.422686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"use_cuda = torch.cuda.is_available()\nuse_cuda","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:02:44.455137Z","iopub.execute_input":"2021-06-24T13:02:44.455479Z","iopub.status.idle":"2021-06-24T13:02:44.464909Z","shell.execute_reply.started":"2021-06-24T13:02:44.455443Z","shell.execute_reply":"2021-06-24T13:02:44.464195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\nparser = argparse.ArgumentParser(description='PyTorch Unsupervised Segmentation')\nparser.add_argument('--nChannel', metavar='N', default=100, type=int, \n                    help='number of channels')\nparser.add_argument('--maxIter', metavar='T', default=1000, type=int, \n                    help='number of maximum iterations')\nparser.add_argument('--minLabels', metavar='minL', default=3, type=int, \n                    help='minimum number of labels')\nparser.add_argument('--lr', metavar='LR', default=0.1, type=float, \n                    help='learning rate')\nparser.add_argument('--nConv', metavar='M', default=2, type=int, \n                    help='number of convolutional layers')\nparser.add_argument('--num_superpixels', metavar='K', default=10000, type=int, \n                    help='number of superpixels')\nparser.add_argument('--compactness', metavar='C', default=100, type=float, \n                    help='compactness of superpixels')\nparser.add_argument('--visualize', metavar='1 or 0', default=1, type=int, \n                    help='visualization flag')\nparser.add_argument('--input', metavar='FILENAME',\n                    help='input image file name', required=True)\nargs = parser.parse_args()\n\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CNN model\nnChannel = 100\nnConv = 2\n\nclass MyNet(nn.Module):\n    def __init__(self, input_dim):\n        super(MyNet, self).__init__()\n        self.conv1 = nn.Conv2d(input_dim, nChannel, kernel_size=3, stride=1, padding=1 )\n        self.bn1 = nn.BatchNorm2d(nChannel)\n        self.conv2 = nn.ModuleList()\n        self.bn2 = nn.ModuleList()\n        for i in range(nConv-1):\n            self.conv2.append( nn.Conv2d(nChannel, nChannel, kernel_size=3, stride=1, padding=1 ) )\n            self.bn2.append( nn.BatchNorm2d(nChannel) )\n        self.conv3 = nn.Conv2d(nChannel, nChannel, kernel_size=1, stride=1, padding=0 )\n        self.bn3 = nn.BatchNorm2d(nChannel)\n\n    def forward(self, x):\n        x = self.conv1(x)\n        x = F.relu( x )\n        x = self.bn1(x)\n        for i in range(nConv-1):\n            x = self.conv2[i](x)\n            x = F.relu( x )\n            x = self.bn2[i](x)\n        x = self.conv3(x)\n        x = self.bn3(x)\n        return x\n","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:02:53.550255Z","iopub.execute_input":"2021-06-24T13:02:53.550583Z","iopub.status.idle":"2021-06-24T13:02:53.559301Z","shell.execute_reply.started":"2021-06-24T13:02:53.550551Z","shell.execute_reply":"2021-06-24T13:02:53.558505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# load image\nfilepath = '/kaggle/input/siim-covid19-detection/train/000c9c05fd14/e555410bd2cd/51759b5579bc.dcm'\nim = resize(pydicom.dcmread(filepath).pixel_array, (512, 512))\nim = im - im.min()\nim = im.astype(float) / im.max()\ndata = torch.from_numpy(im)\nif use_cuda:\n    data = data.cuda()\ndata = Variable(data)\n\nplt.figure(figsize=(10, 10))\nplt.imshow(im, cmap=plt.cm.bone)","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:02:56.787022Z","iopub.execute_input":"2021-06-24T13:02:56.787340Z","iopub.status.idle":"2021-06-24T13:03:02.475912Z","shell.execute_reply.started":"2021-06-24T13:02:56.787311Z","shell.execute_reply":"2021-06-24T13:03:02.475096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# segmentation\n\nmin_sizes = [50, 100, 200]\n\nfig, axs = plt.subplots(2, 2, figsize=(20, 20))\n\nfor i in range(len(min_sizes)):\n    labels = segmentation.felzenszwalb(im, scale=1, sigma=0.8, min_size=min_sizes[i])\n\n    axs[i // 2, i % 2].imshow(segmentation.mark_boundaries(im, labels))\n    \n    \naxs[1, 1].imshow(im, cmap=plt.cm.bone)","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:15:07.744283Z","iopub.execute_input":"2021-06-24T13:15:07.744604Z","iopub.status.idle":"2021-06-24T13:15:09.666514Z","shell.execute_reply.started":"2021-06-24T13:15:07.744574Z","shell.execute_reply":"2021-06-24T13:15:09.665716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labels = segmentation.felzenszwalb(im, scale=1, sigma=0.8, min_size=100)\nlabels = labels.reshape(im.shape[0]*im.shape[1])\nu_labels = np.unique(labels)\nl_inds = []\nfor i in range(len(u_labels)):\n    l_inds.append( np.where( labels == u_labels[ i ] )[ 0 ] )\n    \nplt.figure(figsize=(10, 10))\nplt.imshow(segmentation.mark_boundaries(im, labels.reshape(im.shape)))","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:37:37.404154Z","iopub.execute_input":"2021-06-24T13:37:37.404468Z","iopub.status.idle":"2021-06-24T13:37:38.132632Z","shell.execute_reply.started":"2021-06-24T13:37:37.404439Z","shell.execute_reply":"2021-06-24T13:37:38.131703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train\nmodel = MyNet(1)\nif use_cuda:\n    model.cuda()\nmodel.train()\nloss_fn = torch.nn.CrossEntropyLoss()\n\nlr = 0.1\nmaxIter = 1000\nminLabels = 5\n\noptimizer = optim.SGD(model.parameters(), lr=lr, momentum=0.9)\nfor batch_idx in range(maxIter):\n    # forwarding\n    optimizer.zero_grad()\n    output = model( data.reshape((1, 1, im.shape[0], im.shape[1])).float() )[ 0 ]\n    output = output.permute( 1, 2, 0 ).contiguous().view( -1, nChannel )\n    ignore, target = torch.max( output, 1 )\n    im_target = target.data.cpu().numpy()\n    nLabels = len(np.unique(im_target))\n\n    # superpixel refinement\n    # TODO: use Torch Variable instead of numpy for faster calculation\n    for i in range(len(l_inds)):\n        labels_per_sp = im_target[ l_inds[ i ] ]\n        u_labels_per_sp = np.unique( labels_per_sp )\n        hist = np.zeros( len(u_labels_per_sp) )\n        for j in range(len(hist)):\n            hist[ j ] = len( np.where( labels_per_sp == u_labels_per_sp[ j ] )[ 0 ] )\n        im_target[ l_inds[ i ] ] = u_labels_per_sp[ np.argmax( hist ) ]\n    target = torch.from_numpy( im_target )\n    if use_cuda:\n        target = target.cuda()\n    target = Variable( target )\n    loss = loss_fn(output, target)\n    loss.backward()\n    optimizer.step()\n\n    #print (batch_idx, '/', args.maxIter, ':', nLabels, loss.data[0])\n    print (batch_idx, '/', maxIter, ':', nLabels, loss.item())\n\n    if nLabels <= minLabels:\n        print (\"nLabels\", nLabels, \"reached minLabels\", minLabels, \".\")\n        break\n","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:38:19.128702Z","iopub.execute_input":"2021-06-24T13:38:19.129042Z","iopub.status.idle":"2021-06-24T13:38:21.861793Z","shell.execute_reply.started":"2021-06-24T13:38:19.129014Z","shell.execute_reply":"2021-06-24T13:38:21.860601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output = model( data.reshape((1, 1, im.shape[0], im.shape[1])).float() )[0]\noutput = output.permute( 1, 2, 0 ).contiguous().view( -1, nChannel )\nignore, target = torch.max( output, 1 )\nim_target = target.data.cpu().numpy().reshape(im.shape[0], im.shape[1])\nprint(im_target.shape)\n\nfig, axs = plt.subplots(1, 2, figsize=(20, 10))\naxs[0].imshow(im, cmap=plt.cm.bone)\naxs[1].imshow(im_target)","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:38:23.723445Z","iopub.execute_input":"2021-06-24T13:38:23.723791Z","iopub.status.idle":"2021-06-24T13:38:24.167358Z","shell.execute_reply.started":"2021-06-24T13:38:23.723756Z","shell.execute_reply":"2021-06-24T13:38:24.166649Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im_target","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:40:54.351649Z","iopub.execute_input":"2021-06-24T13:40:54.351979Z","iopub.status.idle":"2021-06-24T13:40:54.357267Z","shell.execute_reply.started":"2021-06-24T13:40:54.351950Z","shell.execute_reply":"2021-06-24T13:40:54.356435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(1, 5, figsize=(20, 10))\nfor i in range(5):\n    axs[i].plot(range(512), im[50 + 100 * i])","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:41:44.191656Z","iopub.execute_input":"2021-06-24T13:41:44.191984Z","iopub.status.idle":"2021-06-24T13:41:44.717451Z","shell.execute_reply.started":"2021-06-24T13:41:44.191954Z","shell.execute_reply":"2021-06-24T13:41:44.716538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xs = np.array(range(512))\nplt.plot(xs, 0.5 + np.sin(xs * 0.03 - a) * 0.3)\nplt.plot(xs, im[150])","metadata":{"execution":{"iopub.status.busy":"2021-06-24T13:54:03.314523Z","iopub.execute_input":"2021-06-24T13:54:03.314883Z","iopub.status.idle":"2021-06-24T13:54:03.616663Z","shell.execute_reply.started":"2021-06-24T13:54:03.314850Z","shell.execute_reply":"2021-06-24T13:54:03.615727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(1, 5, figsize=(20, 10))\n\ndef f(x, a, b):\n    return 0.5 + np.sin(x * 0.03 - a) * b\n\nfor i in range(5):\n    popt, pcov = curve_fit(f, xs, im[50 + 100 * i], p0=[0, 0.3])\n    ys = f(xs, *popt)\n    print(pcov, np.max(np.abs(ys - im[50 + 100 * i])))\n    axs[i].plot(range(512), im[50 + 100 * i])\n    axs[i].plot(xs, ys)","metadata":{"execution":{"iopub.status.busy":"2021-06-24T14:01:32.504289Z","iopub.execute_input":"2021-06-24T14:01:32.504946Z","iopub.status.idle":"2021-06-24T14:01:33.251821Z","shell.execute_reply.started":"2021-06-24T14:01:32.504900Z","shell.execute_reply":"2021-06-24T14:01:33.250878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"i = 2\npopt, pcov = curve_fit(f, xs, im[50 + 100 * i], p0=[0, 0.3])\nys = f(xs, *popt)\n#print(pcov, np.max(np.abs(ys - im[50 + 100 * i])))\naxs[i].plot(range(512), im[50 + 100 * i])\naxs[i].plot(xs, ys)\n\nprint(np.argsort(ys)[:2])\n\nindexes = components[250, np.argsort(ys)[:2]]\nplt.figure(figsize=(10, 10))\nmask = (components == indexes[0]) | (components == indexes[1])\nplt.imshow(mask)","metadata":{"execution":{"iopub.status.busy":"2021-06-24T14:40:06.895545Z","iopub.execute_input":"2021-06-24T14:40:06.895920Z","iopub.status.idle":"2021-06-24T14:40:07.084314Z","shell.execute_reply.started":"2021-06-24T14:40:06.895888Z","shell.execute_reply":"2021-06-24T14:40:07.083408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CNN model\nnChannel = 100\nnConv = 2\n\nclass MyNet(nn.Module):\n    def __init__(self, input_dim):\n        super(MyNet, self).__init__()\n        self.conv1 = nn.Conv2d(input_dim, nChannel, kernel_size=3, stride=1, padding=1 )\n        self.bn1 = nn.BatchNorm2d(nChannel)\n        self.conv2 = nn.ModuleList()\n        self.bn2 = nn.ModuleList()\n        for i in range(nConv-1):\n            self.conv2.append( nn.Conv2d(nChannel, nChannel, kernel_size=3, stride=1, padding=1 ) )\n            self.bn2.append( nn.BatchNorm2d(nChannel) )\n        self.conv3 = nn.Conv2d(nChannel, nChannel, kernel_size=1, stride=1, padding=0 )\n        self.bn3 = nn.BatchNorm2d(nChannel)\n\n    def forward(self, x):\n        x = self.conv1(x)\n        x = F.relu( x )\n        x = self.bn1(x)\n        for i in range(nConv-1):\n            x = self.conv2[i](x)\n            x = F.relu( x )\n            x = self.bn2[i](x)\n        x = self.conv3(x)\n        x = self.bn3(x)\n        return x\n\n\ndef detect_lungs(filepath):\n    # load image\n    #print(filepath)\n    im = resize(pydicom.dcmread(filepath).pixel_array, (256, 256))\n    im = im - im.min()\n    im = im.astype(float) / im.max()\n    data = torch.from_numpy(im)\n    if use_cuda:\n        data = data.cuda()\n    data = Variable(data)\n    \n    # Segmentation\n    labels = segmentation.felzenszwalb(im, scale=1, sigma=0.8, min_size=50)\n    labels = labels.reshape(im.shape[0]*im.shape[1])\n    u_labels = np.unique(labels)\n    l_inds = []\n    for i in range(len(u_labels)):\n        l_inds.append( np.where( labels == u_labels[ i ] )[ 0 ] )\n\n    # train\n    model = MyNet(1)\n    if use_cuda:\n        model.cuda()\n    model.train()\n    loss_fn = torch.nn.CrossEntropyLoss()\n\n    lr = 0.1\n    maxIter = 100\n    minLabels = 10\n\n    optimizer = optim.SGD(model.parameters(), lr=lr, momentum=0.9)\n    for batch_idx in range(maxIter):\n        # forwarding\n        optimizer.zero_grad()\n        output = model( data.reshape((1, 1, im.shape[0], im.shape[1])).float() )[ 0 ]\n        output = output.permute( 1, 2, 0 ).contiguous().view( -1, nChannel )\n        ignore, target = torch.max( output, 1 )\n        im_target = target.data.cpu().numpy()\n        nLabels = len(np.unique(im_target))\n\n        # superpixel refinement\n        # TODO: use Torch Variable instead of numpy for faster calculation\n        for i in range(len(l_inds)):\n            labels_per_sp = im_target[ l_inds[ i ] ]\n            u_labels_per_sp = np.unique( labels_per_sp )\n            hist = np.zeros( len(u_labels_per_sp) )\n            for j in range(len(hist)):\n                hist[ j ] = len( np.where( labels_per_sp == u_labels_per_sp[ j ] )[ 0 ] )\n            im_target[ l_inds[ i ] ] = u_labels_per_sp[ np.argmax( hist ) ]\n        target = torch.from_numpy( im_target )\n        if use_cuda:\n            target = target.cuda()\n        target = Variable( target )\n        loss = loss_fn(output, target)\n        loss.backward()\n        optimizer.step()\n\n        #print (batch_idx, '/', args.maxIter, ':', nLabels, loss.data[0])\n        #print (batch_idx, '/', maxIter, ':', nLabels, loss.item())\n\n        if nLabels <= minLabels:\n            #print (\"nLabels\", nLabels, \"reached minLabels\", minLabels, \".\")\n            break\n\n    # predict\n    output = model( data.reshape((1, 1, im.shape[0], im.shape[1])).float() )[0]\n    output = output.permute( 1, 2, 0 ).contiguous().view( -1, nChannel )\n    ignore, target = torch.max( output, 1 )\n    components = target.data.cpu().numpy().reshape(im.shape[0], im.shape[1])\n\n    # Find a horizontal line with lungs\n    def f(x, a, b):\n        return 0.5 + np.sin(x * 0.03 - a) * b\n\n    xs = np.array(range(im.shape[1]))\n    best_line = 0\n    best_score = np.inf\n    for i in range(7):\n        line = im.shape[1] // 8 * (i + 1)\n        popt, pcov = curve_fit(f, xs, im[line], p0=[0, 0.3])\n        ys = f(xs, *popt)\n        score = np.max(np.abs(ys - im[line]))\n        if score < best_score:\n            best_line = line\n            best_score = score\n    \n    # Find the lungs on the best line\n    popt, pcov = curve_fit(f, xs, im[best_line], p0=[0, 0.3])\n    ys = f(xs, *popt)\n\n    indexes = components[best_line, np.argsort(ys)[:2]]\n    mask = (components == indexes[0]) | (components == indexes[1])\n    return mask","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:08:07.076665Z","iopub.execute_input":"2021-06-24T15:08:07.076999Z","iopub.status.idle":"2021-06-24T15:08:07.104487Z","shell.execute_reply.started":"2021-06-24T15:08:07.076970Z","shell.execute_reply":"2021-06-24T15:08:07.103647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#filepath = '/kaggle/input/siim-covid19-detection/train/000c9c05fd14/e555410bd2cd/51759b5579bc.dcm'\nfilepath = '/kaggle/input/siim-covid19-detection/train/ff0879eb20ed/d8a644cc4f93/000c3a3f293f.dcm'\nmask = detect_lungs(filepath)\n\nim = resize(pydicom.dcmread(filepath).pixel_array, (512, 512))\nim = im - im.min()\nim = im.astype(float) / im.max()\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 10))\naxs[0].imshow(im, cmap=plt.cm.bone)\naxs[1].imshow(mask)","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:08:07.383194Z","iopub.execute_input":"2021-06-24T15:08:07.383513Z","iopub.status.idle":"2021-06-24T15:08:08.906352Z","shell.execute_reply.started":"2021-06-24T15:08:07.383485Z","shell.execute_reply":"2021-06-24T15:08:08.905494Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f = '/kaggle/input/siim-covid19-detection/train/9d514ce429a7/22897cd1daa0/0012ff7358bc.dcm'\nmask = detect_lungs(f)\n\nim = resize(pydicom.dcmread(f).pixel_array, (256, 256))\nim = im - im.min()\nim = im.astype(float) / im.max()\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 10))\naxs[0].imshow(im, cmap=plt.cm.bone)\naxs[1].imshow(mask)","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:08:52.474406Z","iopub.execute_input":"2021-06-24T15:08:52.474741Z","iopub.status.idle":"2021-06-24T15:08:54.134743Z","shell.execute_reply.started":"2021-06-24T15:08:52.474709Z","shell.execute_reply":"2021-06-24T15:08:54.133947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labels = segmentation.felzenszwalb(im, scale=1, sigma=0.8, min_size=1000)\nlabels = labels.reshape(im.shape[0]*im.shape[1])\nu_labels = np.unique(labels)\nl_inds = []\nfor i in range(len(u_labels)):\n    l_inds.append( np.where( labels == u_labels[ i ] )[ 0 ] )\n    \nplt.figure(figsize=(10, 10))\nplt.imshow(segmentation.mark_boundaries(im, labels.reshape(im.shape)))","metadata":{"execution":{"iopub.status.busy":"2021-06-24T15:08:37.967958Z","iopub.execute_input":"2021-06-24T15:08:37.968286Z","iopub.status.idle":"2021-06-24T15:08:38.582608Z","shell.execute_reply.started":"2021-06-24T15:08:37.968254Z","shell.execute_reply":"2021-06-24T15:08:38.581705Z"},"trusted":true},"execution_count":null,"outputs":[]}]}