{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":94689,"databundleVersionId":11605086,"sourceType":"competition"}],"dockerImageVersionId":31012,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Foraminifera classification\n\nThe task of this competition is to classify volumetric scans of foraminifera into 14 types, plus a 15th 'unknown' class, which includes scans of foram fragments, sand, or other similar objects. Each image is 128-by-128-by-128 voxels.\n\nAcknowledgement: Some elements in this notebook are taken from https://www.kaggle.com/code/muhammadaashirirshad/multiview-2d-cnn-f1-score-0-51.\n\nReferences:\n- [Forams Classification 2025](https://www.kaggle.com/competitions/forams-classification-2025)\n- [Foraminifera on Wikipedia](https://en.wikipedia.org/wiki/Foraminifera)","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nimport seaborn as sns\nfrom glob import glob\nimport PIL\nfrom tifffile import imread\nfrom tqdm import tqdm\n\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import StandardScaler, FunctionTransformer\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis\nfrom sklearn.model_selection import cross_val_score, cross_val_predict\nfrom sklearn.metrics import f1_score, confusion_matrix, log_loss, accuracy_score\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:34:38.033007Z","iopub.execute_input":"2025-04-26T07:34:38.033307Z","iopub.status.idle":"2025-04-26T07:34:40.047301Z","shell.execute_reply.started":"2025-04-26T07:34:38.033283Z","shell.execute_reply":"2025-04-26T07:34:40.046202Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Visualization\n\nWe show two samples per class for all 14 classes:","metadata":{}},{"cell_type":"code","source":"vis_list = sorted(glob('/kaggle/input/forams-classification-2025/visualizations/visualizations/labelled/*.jpg'))\nfor j in [0, 5]:\n    vis_list_subset = vis_list[j::15]\n    plt.figure(figsize=(20, 4))\n    for i, filename in enumerate(vis_list_subset):\n        # print(filename)\n        plt.subplot(1, len(vis_list_subset), i+1)\n        img = np.asarray(PIL.Image.open(filename))\n        plt.imshow(img)\n        plt.title(str(i))\n        plt.axis('off')\n    plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:34:40.049114Z","iopub.execute_input":"2025-04-26T07:34:40.050389Z","iopub.status.idle":"2025-04-26T07:34:42.820097Z","shell.execute_reply.started":"2025-04-26T07:34:40.050355Z","shell.execute_reply":"2025-04-26T07:34:42.819107Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Now we show some cross sections (click on the output to open):","metadata":{}},{"cell_type":"code","source":"def visualize_volume(volume_path, title=None, c=64):\n    # Load the volume\n    volume = imread(volume_path)\n    \n    # Get slices in each dimension\n    slice_x = volume[c, :, :]\n    slice_y = volume[:, c, :]\n    slice_z = volume[:, :, c]\n    \n    # Create a figure with subplots\n    fig, axes = plt.subplots(1, 3, figsize=(25, 8))\n    \n    # Plot slices\n    axes[0].imshow(slice_x, cmap='gray', vmin=0, vmax=255)\n    axes[0].set_title(f'X-Slice ({c})')\n    \n    axes[1].imshow(slice_y, cmap='gray', vmin=0, vmax=255)\n    axes[1].set_title(f'Y-Slice ({c})')\n    \n    axes[2].imshow(slice_z, cmap='gray', vmin=0, vmax=255)\n    axes[2].set_title(f'Z-Slice ({c})')\n    \n    plt.tight_layout()\n    if title:\n        plt.suptitle(title, y=1.05, fontsize=24)\n    plt.show()\n\nvol_list = sorted(glob('/kaggle/input/forams-classification-2025/volumes/volumes/labelled/*.tif'))\n\nfor i, filename in enumerate(vol_list[::15]):\n    visualize_volume(filename, title=f\"Class {i}\")\n    for c in [40, 20]:\n        visualize_volume(filename, c=c)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:34:42.821338Z","iopub.execute_input":"2025-04-26T07:34:42.822370Z","iopub.status.idle":"2025-04-26T07:35:15.882330Z","shell.execute_reply.started":"2025-04-26T07:34:42.822332Z","shell.execute_reply":"2025-04-26T07:35:15.881225Z"},"_kg_hide-input":true,"_kg_hide-output":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Simple feature engineering and model","metadata":{}},{"cell_type":"code","source":"# Read the 210 true labels\nlabels = pd.read_csv(f'/kaggle/input/forams-classification-2025/labelled.csv', index_col='id')['label'] # 210 rows, id and label\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:35:15.884525Z","iopub.execute_input":"2025-04-26T07:35:15.885095Z","iopub.status.idle":"2025-04-26T07:35:15.918526Z","shell.execute_reply.started":"2025-04-26T07:35:15.885066Z","shell.execute_reply":"2025-04-26T07:35:15.917438Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Squared distance from center for all 128*128*128 pixels\nr0 = np.square(np.arange(128) - 63.5)\nr = r0.reshape(-1, 1, 1) + r0.reshape(1, -1, 1) + r0.reshape(1, 1, -1)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:35:15.919594Z","iopub.execute_input":"2025-04-26T07:35:15.919950Z","iopub.status.idle":"2025-04-26T07:35:15.938669Z","shell.execute_reply.started":"2025-04-26T07:35:15.919921Z","shell.execute_reply":"2025-04-26T07:35:15.937763Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Feature engineering\n\ndef make_feature_array(vol_list):\n    \"\"\"Compute the features\n    \n    Parameter\n    vol_list: list of filenames with 128*128*128 tiff files\n\n    Return value\n    X: array of shape (n_samples, n_features)\n    \"\"\"\n    X = np.full((len(vol_list), 14), np.nan)\n    for i, filename in enumerate(tqdm(vol_list)):\n        volume = imread(filename)\n        X[i] = np.array([volume.mean(),\n                         volume.ravel()[volume.ravel() > 0].mean(),\n                         volume.max(axis=0).mean(),\n                         volume.max(axis=1).mean(),\n                         volume.max(axis=2).mean(),\n                         volume.max(axis=1).max(axis=0).mean(),\n                         volume.max(axis=2).max(axis=1).mean(),\n                         volume.max(axis=2).max(axis=0).mean(),\n                         (r * volume).mean(),\n                         (volume != 0).mean(),\n                         np.abs(np.diff(volume, axis=0)).mean(),\n                         np.abs(np.diff(volume, axis=1)).mean(),\n                         np.abs(np.diff(volume, axis=2)).mean(),\n                         float(filename[filename.index('_sc_') + 4 : -4].replace('_', '.')),\n                        ])\n    return X\n\nvol_list = sorted(glob('/kaggle/input/forams-classification-2025/volumes/volumes/labelled/*.tif'))\nX = make_feature_array(vol_list)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:35:15.939556Z","iopub.execute_input":"2025-04-26T07:35:15.939799Z","iopub.status.idle":"2025-04-26T07:35:32.670198Z","shell.execute_reply.started":"2025-04-26T07:35:15.939779Z","shell.execute_reply":"2025-04-26T07:35:32.669352Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Linear Discriminant Analysis\n\nmodel = LinearDiscriminantAnalysis()\n\noof = cross_val_predict(model, X, labels)\nprint(f\"# F1:      {f1_score(labels, oof, average='macro'):.3f}\", end='   ')\nprint(f\"Accuracy:  {accuracy_score(labels, oof):.3f}\")\n\noof_prob = cross_val_predict(model, X, labels, method='predict_proba')\nprint(f\"# Logloss: {log_loss(labels, oof_prob):.3f}\")\n\nsns.heatmap(confusion_matrix(labels, oof), annot=True)\nplt.show()\n# F1:      0.528   Acc:     0.529\n# Logloss: 1.985","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:35:32.671407Z","iopub.execute_input":"2025-04-26T07:35:32.671773Z","iopub.status.idle":"2025-04-26T07:35:33.224098Z","shell.execute_reply.started":"2025-04-26T07:35:32.671739Z","shell.execute_reply":"2025-04-26T07:35:33.223231Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"%%time\n# Fit the model to the full dataset\ndisplay(model)\nmodel.fit(X, labels)\n\n# Compute the test predictions\nvol_list = sorted(glob('/kaggle/input/forams-classification-2025/volumes/volumes/unlabelled/*.tif'))\n# vol_list = np.random.default_rng().choice(vol_list, size=2000, replace=False) # subset to save time\nX_test = make_feature_array(vol_list)\ny_pred = model.predict(X_test)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:35:33.224999Z","iopub.execute_input":"2025-04-26T07:35:33.225280Z","iopub.status.idle":"2025-04-26T07:38:22.816196Z","shell.execute_reply.started":"2025-04-26T07:35:33.225240Z","shell.execute_reply":"2025-04-26T07:38:22.815155Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"x, y = np.unique(y_pred, return_counts=True)\nplt.title('Predicted classes')\nplt.bar(x, y)\nplt.xticks(x)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:40:08.476244Z","iopub.execute_input":"2025-04-26T07:40:08.476615Z","iopub.status.idle":"2025-04-26T07:40:08.665663Z","shell.execute_reply.started":"2025-04-26T07:40:08.476588Z","shell.execute_reply":"2025-04-26T07:40:08.664740Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"unlabelled_index = pd.read_csv(f'/kaggle/input/forams-classification-2025/unlabelled.csv', index_col='id').index # 18216 rows, id and label (label is NaN)\npd.Series(y_pred, index=unlabelled_index[:len(y_pred)], name='label').to_csv('submission.csv')\n!head submission.csv","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-26T07:38:22.817176Z","iopub.execute_input":"2025-04-26T07:38:22.817413Z","iopub.status.idle":"2025-04-26T07:38:22.976902Z","shell.execute_reply.started":"2025-04-26T07:38:22.817394Z","shell.execute_reply":"2025-04-26T07:38:22.975769Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# To do\n\n- Find more and better features\n- Add class 14 ('unknown')\n- Convert model to semi-supervised learning","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}