{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":61446,"databundleVersionId":6962461,"sourceType":"competition"}],"dockerImageVersionId":30579,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"\n\n\n### XGBoost is all you need? as the title implies, the below exercise attempts to use xgboost to classify vessel/not-vessel at the pixel level.\n\n<img src=\"https://i.imgflip.com/85tx5h.jpg\" width=420>\n\n\n\n### Preface\n#### Long long time ago, before the era of GPUs, your grand parents had to fine tune and extract features from image patches using matlab and spend 2 whole years to convert the algo to C++ libraries to match the features coming out from Matlab, only to realize their algo have been replaced by conv-nets 1 year after deployment to production. ;)\n\n<img src='https://media.npr.org/assets/img/2023/05/26/honest-work-meme-cb0f0fb2227fb84b77b3c9a851ac09b095ab74d8-s1100-c50.jpg' width=420>\n\n","metadata":{}},{"cell_type":"code","source":"\nfrom xgboost import XGBClassifier\n","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:46:16.815951Z","iopub.execute_input":"2023-11-13T05:46:16.816449Z","iopub.status.idle":"2023-11-13T05:46:16.822129Z","shell.execute_reply.started":"2023-11-13T05:46:16.816411Z","shell.execute_reply":"2023-11-13T05:46:16.820924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n%matplotlib inline\n\nimport random\nimport time\nimport os\nimport pathlib\nimport pandas as pd\nimport numpy as np\nimport skimage\nimport tifffile as tiff\n\nfrom sklearn.mixture import GaussianMixture\nfrom scipy import ndimage as ndi\nimport skimage\nimport sklearn\nfrom sklearn.model_selection import train_test_split\n\ndef prepare():\n\n    base_path = '/kaggle/input/blood-vessel-segmentation/train'\n\n    for x in os.listdir(base_path):\n        folder_path = os.path.join(base_path,x)\n        print(folder_path)\n        print(os.listdir(folder_path))\n\n    mask_path_list = [str(x) for x in pathlib.Path(base_path).rglob('*.tif') if '/labels/' in str(x)]\n\n    mylist = []\n    for mask_path in mask_path_list:\n        if 'kidney_3_sparse' in mask_path:\n            # we skip it since dense is provided.\n            continue\n        image_path = mask_path.replace(\"/labels/\",\"/images/\")\n        if 'kidney_3_dense' in mask_path:\n            image_path = image_path.replace(\"kidney_3_dense\",\"kidney_3_sparse\")\n        assert(os.path.exists(image_path))\n        assert(os.path.exists(mask_path))\n        myitem=dict(\n            image_path=image_path,\n            mask_path=mask_path\n        )\n        mylist.append(myitem)\n\n    df = pd.DataFrame(mylist)\n    train_df=df.sample(frac=0.01,random_state=420) # tiny!\n    df=df.drop(train_df.index)\n    val_df=df.sample(frac=0.01,random_state=420) # tiny!\n    test_df=df.drop(val_df.index)\n    \n    os.makedirs('csvs',exist_ok=True)\n    \n    train_df=train_df.reset_index()\n    train_df.to_csv('csvs/train.csv',index=False)\n\n    val_df=val_df.reset_index()\n    val_df.to_csv('csvs/val.csv',index=False)\n    \n    test_df=test_df.reset_index()\n    test_df.to_csv('csvs/test.csv',index=False)\n\n    print(train_df.shape,val_df.shape,test_df.shape)\n\ndef readimage(mypath):\n    image = tiff.imread(mypath)\n    assert(len(image.shape)==2)\n    image = image.astype(np.float64)\n    return image\n\ndef readmask(mypath):\n    mask = tiff.imread(mypath)\n    assert(len(mask.shape)==2)\n    mask = mask > 0\n    mask = mask.astype(np.float64)\n    return mask\n\ndef normalize_image(image):\n    image = (image-np.mean(image))/np.std(image)\n    return image\n\n# https://www.kaggle.com/code/paulorzp/run-length-encode-and-decode/script\n# ref.: https://www.kaggle.com/stainsby/fast-tested-rle\ndef rle_encode(img):\n    '''\n    img: numpy array, 1 - mask, 0 - background\n    Returns run length as string formated\n    '''\n    pixels = img.flatten()\n    pixels = np.concatenate([[0], pixels, [0]])\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 1\n    runs[1::2] -= runs[::2]\n    mystr = ' '.join(str(x) for x in runs)\n    if mystr == \"\":\n        mystr = \"1 0\"\n    return mystr\n\ndef rle_decode(mask_rle, shape):\n    '''\n    mask_rle: run-length as string formated (start length)\n    shape: (height,width) of array to return \n    Returns numpy array, 1 - mask, 0 - background\n\n    '''\n    s = mask_rle.split()\n    starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]\n    starts -= 1\n    ends = starts + lengths\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    for lo, hi in zip(starts, ends):\n        img[lo:hi] = 1\n\n    return img.reshape(shape)\n","metadata":{"execution":{"iopub.status.busy":"2023-11-13T04:36:02.216847Z","iopub.execute_input":"2023-11-13T04:36:02.217275Z","iopub.status.idle":"2023-11-13T04:36:02.244821Z","shell.execute_reply.started":"2023-11-13T04:36:02.217242Z","shell.execute_reply":"2023-11-13T04:36:02.243756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prepare()","metadata":{"execution":{"iopub.status.busy":"2023-11-13T03:52:45.381983Z","iopub.execute_input":"2023-11-13T03:52:45.382361Z","iopub.status.idle":"2023-11-13T03:53:07.379255Z","shell.execute_reply.started":"2023-11-13T03:52:45.382333Z","shell.execute_reply":"2023-11-13T03:53:07.377911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"! ls -lha csvs","metadata":{"execution":{"iopub.status.busy":"2023-11-13T03:53:07.381537Z","iopub.execute_input":"2023-11-13T03:53:07.382300Z","iopub.status.idle":"2023-11-13T03:53:08.405533Z","shell.execute_reply.started":"2023-11-13T03:53:07.382255Z","shell.execute_reply":"2023-11-13T03:53:08.404482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n\n<img src=\"https://pictures.abebooks.com/isbn/9780130085191-us.jpg\" width=240>\n<br>\nAre you old enough to remember the above book (version 1 cover!)? ;)\n","metadata":{}},{"cell_type":"code","source":"#\n# sample and visualize data\n#\n\ndf = pd.read_csv('csvs/train.csv')\nfor n,row in df.iterrows():\n    \n    t0 = time.time() \n    basename = os.path.basename(row.image_path)\n    image = readimage(row.image_path)\n    image = normalize_image(image)\n    mask = readmask(row.mask_path)\n    t1 = time.time() \n    \n    t2 = time.time()\n    # speed up processing time?\n    image_mod = skimage.transform.resize(image, (512,512), order=0)\n    mask_mod = skimage.transform.resize(mask, (512,512), order=0)\n    t3 = time.time()\n    if np.sum(mask) > 0:\n        vessel_mean = np.mean(image_mod[mask_mod>0])\n    else:\n        vessel_mean = np.nan\n    \n    X = np.expand_dims(image_mod.ravel(),axis=-1)\n    t4 = time.time()\n    gmm = GaussianMixture(n_components=2, random_state=0).fit(X)\n    t5 = time.time()\n     \n    print(t1-t0,t3-t2,t5-t4)\n    mu0 =gmm.means_[0][0]\n    var0 = gmm.covariances_[0][0][0]\n    mu1 =gmm.means_[1][0]\n    var1 = gmm.covariances_[1][0][0]\n\n    blocks=skimage.util.view_as_blocks(image_mod,(32,32))\n    background_values = []\n    for idx_x,idx_y in [[0,0],[0,-1],[-1,0],[-1,-1]]:\n        background_values.append(np.mean(blocks[:,:,idx_x,idx_y]))\n    background_mean = np.mean(background_values)\n    \n    plt.figure(n)\n    plt.subplot(212)\n    n_bins = 100\n    plt.hist(image.ravel(), bins=n_bins)\n    plt.axvline(vessel_mean,color='red')\n    plt.axvline(mu1,color='green')\n    plt.axvline(mu1-var1,color='blue',linestyle='--')\n    plt.axvline(mu0,color='green')\n    plt.axvline(mu0-var0,color='blue',linestyle='--')\n    plt.axvline(background_mean,color='black')\n    plt.title('vessel(red),mu(green),mu-sd(blue),bkgd(black)')\n    \n    plt.grid(True)\n    plt.subplot(221)\n    plt.imshow(image,cmap='gray')\n    plt.colorbar()\n    plt.subplot(222)\n    plt.imshow(mask,cmap='gray')\n    plt.title(f'mean {vessel_mean:1.1f}')\n    plt.colorbar()\n    # png_path = f'mytmp/eda-{n}.png'\n    #plt.savefig(png_path)\n    #plt.close()\n    if n>10:\n        break","metadata":{"execution":{"iopub.status.busy":"2023-11-13T03:53:08.407735Z","iopub.execute_input":"2023-11-13T03:53:08.408150Z","iopub.status.idle":"2023-11-13T03:53:40.212621Z","shell.execute_reply.started":"2023-11-13T03:53:08.408110Z","shell.execute_reply":"2023-11-13T03:53:40.211459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"prelminary thoughts formed from above samples:\n\n+ image intensity distribution likely is bimodal or unimodal.\n+ looks like image intensity distribution is bimodal at lower resolution (in histology terms: \"wholemount\" versus \"micrographs\"?).\n+ regardless of resolution, vessel structures are darker than its surroundings.\n+ shape wise, on the 2d images, vessel shape is circular or tubular due to image orientation.\n+ given the small sample size, most likely we can assume the flavors in train-set will not be representative of those in testset / real-world data.\n+ for old times sake, why not give old school image processing and analysis / ml techniques a try first ? ;) \n","metadata":{}},{"cell_type":"code","source":"\ndef process(image):\n    # prepare filter bank kernels\n    feat_list = [image]\n    for sigma in (1, 3, 10):\n        gauss_image = skimage.filters.gaussian(image, sigma=sigma)\n        feat_list.append(gauss_image)\n\n    for sigma in (1, 3):\n        for theta in range(4):\n            theta = theta / 4. * np.pi\n            for frequency in (0.05, 0.25):\n                gabor_real,_=skimage.filters.gabor(\n                    image, frequency, theta=theta, bandwidth=1,\n                    sigma_x=sigma,sigma_y=sigma,offset=0, mode='reflect', cval=0)\n                feat_list.append(gabor_real)\n\n    #image = skimage.transform.resize(image, (512,512), order=0)\n    X = np.expand_dims(image.ravel(),axis=-1)\n    gmm = GaussianMixture(n_components=2, random_state=0).fit(X)\n    mu0 =gmm.means_[0][0]\n    var0 = gmm.covariances_[0][0][0]\n    mu1 =gmm.means_[1][0]\n    var1 = gmm.covariances_[1][0][0]\n    \n    blocks=skimage.util.view_as_blocks(image_mod,(32,32))\n    background_values = []\n    for idx_x,idx_y in [[0,0],[0,-1],[-1,0],[-1,-1]]:\n        background_values.append(np.mean(blocks[:,:,idx_x,idx_y]))\n    background_mean = np.mean(background_values)\n    for val in [mu0,var0,mu1,var1,background_mean]:\n        tmp = val*np.ones_like(image)\n        feat_list.append(tmp)\n\n    new_shape = (len(feat_list),np.prod(image.shape))\n    feat = np.reshape(np.array(feat_list), new_shape).T\n    return feat\n","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:25:27.681085Z","iopub.execute_input":"2023-11-13T05:25:27.681943Z","iopub.status.idle":"2023-11-13T05:25:27.693384Z","shell.execute_reply.started":"2023-11-13T05:25:27.681905Z","shell.execute_reply":"2023-11-13T05:25:27.692388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nylist = []\nxlist = []\ndf = pd.read_csv('csvs/train.csv')\nprint(df.shape)\nwidth,height = 1024,1024\n\nfor n,row in df.iterrows():\n    image = readimage(row.image_path)\n    image = skimage.transform.resize(image, (width,height), order=0)\n    image = normalize_image(image)\n    feat = process(image)\n\n    mask = readmask(row.mask_path)\n    mask = skimage.transform.resize(mask, (width,height), order=0)\n    mask = mask>0\n    mask = np.expand_dims(mask.ravel(),axis=-1)    \n    negative_idx,_ = np.where(mask==0)\n    positive_idx,_ = np.where(mask==1)\n    random.shuffle(negative_idx)\n    random.shuffle(positive_idx)\n    print(len(positive_idx[:1000]),len(negative_idx[:3000]))\n    xlist.append(feat[positive_idx[:1000],:])\n    ylist.append(mask[positive_idx[:1000],:])\n    xlist.append(feat[negative_idx[:3000],:])\n    ylist.append(mask[negative_idx[:3000],:])\n    \n    print(n,len(df))\n    #\n    # too slow! performance wise, traditional feature extraction + ml\n    # will be hard to compete with GPU-algos (ofcourse you can figure out \n    # how to pipe the feature compute to GPU... but if you have \"big data\" \n    # why not just opt for gradient-descent/DL with GPUs )\n    #\n    if n > 10:\n        break\n","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:07:08.558570Z","iopub.execute_input":"2023-11-13T05:07:08.558992Z","iopub.status.idle":"2023-11-13T05:09:22.909958Z","shell.execute_reply.started":"2023-11-13T05:07:08.558957Z","shell.execute_reply":"2023-11-13T05:09:22.908838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Xarr = np.concatenate(xlist,axis=0)\nYarr = np.concatenate(ylist,axis=0)\nprint(Xarr.shape,Yarr.shape)","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:14:32.710145Z","iopub.execute_input":"2023-11-13T05:14:32.710628Z","iopub.status.idle":"2023-11-13T05:14:32.720113Z","shell.execute_reply.started":"2023-11-13T05:14:32.710592Z","shell.execute_reply":"2023-11-13T05:14:32.719154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(Xarr,Yarr, test_size=.2)\nbst = XGBClassifier(n_estimators=2, max_depth=2, learning_rate=1, objective='binary:logistic')\n\nbst.fit(X_train, y_train)\nbst.save_model('xgb_model.json')","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:14:33.415239Z","iopub.execute_input":"2023-11-13T05:14:33.415646Z","iopub.status.idle":"2023-11-13T05:14:33.520838Z","shell.execute_reply.started":"2023-11-13T05:14:33.415604Z","shell.execute_reply":"2023-11-13T05:14:33.519948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# bst = XGBClassifier()\n# bst.load_model('xgb_model.json')","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:14:33.924712Z","iopub.execute_input":"2023-11-13T05:14:33.925877Z","iopub.status.idle":"2023-11-13T05:14:33.933472Z","shell.execute_reply.started":"2023-11-13T05:14:33.925831Z","shell.execute_reply":"2023-11-13T05:14:33.932557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# make predictions\ny_pred_train = bst.predict(X_train)\nprint(np.unique(y_pred_train))\ny_pred_test = bst.predict(X_test)\nprint(np.unique(y_pred_test))\n# huh\nacc = sklearn.metrics.accuracy_score(y_train, y_pred_train)\nprint('acc train',acc)\nacc = sklearn.metrics.accuracy_score(y_test, y_pred_test)\nprint('acc \"test\"',acc)\ncm = sklearn.metrics.confusion_matrix(y_test, y_pred_test)\nprint(cm)","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:14:34.376795Z","iopub.execute_input":"2023-11-13T05:14:34.377514Z","iopub.status.idle":"2023-11-13T05:14:34.403931Z","shell.execute_reply.started":"2023-11-13T05:14:34.377478Z","shell.execute_reply":"2023-11-13T05:14:34.402888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def inference(tif_path):\n    image = readimage(tif_path)\n    original_shape = image.shape\n    image = normalize_image(image)\n    image = skimage.transform.resize(image, (width,height), order=0)\n    feat = process(image)\n    y_hat = bst.predict(feat)\n    output_mask = np.reshape(y_hat,image.shape)\n    output_mask = skimage.transform.resize(output_mask, original_shape, order=0)\n    output_mask = output_mask > 0\n    return output_mask","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:14:37.280872Z","iopub.execute_input":"2023-11-13T05:14:37.281877Z","iopub.status.idle":"2023-11-13T05:14:37.287357Z","shell.execute_reply.started":"2023-11-13T05:14:37.281838Z","shell.execute_reply":"2023-11-13T05:14:37.286233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef nosubmit():\n    df = pd.read_csv('csvs/test.csv')\n    mylist = []\n    for n,row in df.iterrows():\n        tif_path = row.image_path\n        mask_path = row.mask_path\n        dataset_id = os.path.basename(os.path.dirname(os.path.dirname(tif_path)))\n        slice_id = os.path.basename(tif_path).replace(\".tif\",\"\")\n        case_id = f'{dataset_id}_{slice_id}'\n        output_mask = inference(tif_path)\n        \n        image = readimage(tif_path)\n        mask = readmask(mask_path)\n        plt.figure(n)\n        plt.subplot(131)\n        plt.imshow(image)\n        plt.subplot(132)\n        plt.imshow(mask)\n        plt.subplot(133)\n        plt.imshow(output_mask)\n        rle_str = rle_encode(output_mask)\n        print(len(rle_str),case_id)\n        if n > 10:\n            break\n\nnosubmit()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:14:38.013634Z","iopub.execute_input":"2023-11-13T05:14:38.014997Z","iopub.status.idle":"2023-11-13T05:16:45.335405Z","shell.execute_reply.started":"2023-11-13T05:14:38.014949Z","shell.execute_reply":"2023-11-13T05:16:45.334217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\ndef submit():\n    test_folder = \"/kaggle/input/blood-vessel-segmentation/test\"\n    tif_list = sorted([str(x) for x in pathlib.Path(test_folder).rglob(\"*.tif\")])\n    mylist = []\n    for n,tif_path in enumerate(tif_list):\n        dataset_id = os.path.basename(os.path.dirname(os.path.dirname(tif_path)))\n        slice_id = os.path.basename(tif_path).replace(\".tif\",\"\")\n        case_id = f'{dataset_id}_{slice_id}'\n        print(case_id)\n        output_mask = inference(tif_path)\n        #plt.figure(n)\n        #plt.imshow(output_mask)\n        rle_str = rle_encode(output_mask)\n        print(len(rle_str),tif_path)\n        \n        myitem={\n            \"id\":case_id,\n            \"rle\":rle_str\n        }\n        mylist.append(myitem)\n\n    df = pd.DataFrame(mylist)\n    df.to_csv(\"submission.csv\",index=False)\n\nsubmit()","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:17:10.220655Z","iopub.execute_input":"2023-11-13T05:17:10.221411Z","iopub.status.idle":"2023-11-13T05:18:09.264809Z","shell.execute_reply.started":"2023-11-13T05:17:10.221373Z","shell.execute_reply":"2023-11-13T05:18:09.263730Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv('submission.csv')\nfor n,row in df.iterrows():\n    print(row.id,len(row.rle))","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:19:32.518916Z","iopub.execute_input":"2023-11-13T05:19:32.519352Z","iopub.status.idle":"2023-11-13T05:19:32.543571Z","shell.execute_reply.started":"2023-11-13T05:19:32.519321Z","shell.execute_reply":"2023-11-13T05:19:32.541978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#\n#\n# Thanks for scrolling this way down. That means you have really liked this post, remember to give it a thumbs up!\n# And by no means am I advocating xgboost for computer vision! conv/u/swin-net is the way to go.\n#\n##\n#\n#    I didn’t choose the XGBoost life.\n#    The XGBoost life chose me.\n#\n#                    - Bojan Tunguz\n#                      https://twitter.com/tunguz/status/1723785270590611903\n#\n#","metadata":{"execution":{"iopub.status.busy":"2023-11-13T05:48:14.930548Z","iopub.execute_input":"2023-11-13T05:48:14.931123Z","iopub.status.idle":"2023-11-13T05:48:14.936147Z","shell.execute_reply.started":"2023-11-13T05:48:14.931082Z","shell.execute_reply":"2023-11-13T05:48:14.934819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}