{"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":"markdown","source":"# Read Some `RSNA 2022 Cervical Spine Fracture Detection` Competetion Data is Slow\n\nAs title, after I did the [data overview](https://www.kaggle.com/code/shihhsuanchen/cspine-overview/data), I found that reading pixel array of some study is very slow. In this notebook, I realize that PixelData of dicom file with header starts with `\\xfe\\xff`, which is likely to be `utf-16` encoding, has reading speed (number of slices per seconds) two times lower than other studies.","metadata":{}},{"cell_type":"code","source":"!pip install -qU ../input/for-pydicom/python_gdcm-3.0.14-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl ../input/for-pydicom/pylibjpeg-1.4.0-py3-none-any.whl --find-links frozen_packages --no-index","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:27:39.636008Z","iopub.execute_input":"2022-08-25T05:27:39.636493Z","iopub.status.idle":"2022-08-25T05:27:51.31669Z","shell.execute_reply.started":"2022-08-25T05:27:39.636459Z","shell.execute_reply":"2022-08-25T05:27:51.315305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport time\nimport nibabel as nib\nimport pandas as pd\nimport numpy as np\nimport torch\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\n\nfrom pydicom import dcmread","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:27:51.318724Z","iopub.execute_input":"2022-08-25T05:27:51.319089Z","iopub.status.idle":"2022-08-25T05:27:51.326191Z","shell.execute_reply.started":"2022-08-25T05:27:51.319055Z","shell.execute_reply":"2022-08-25T05:27:51.325178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"basedir = '../input/rsna-2022-cervical-spine-fracture-detection'\nimage_dir = os.path.join(basedir, 'train_images')\nmask_dir = os.path.join(basedir, 'segmentations')","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:27:51.328121Z","iopub.execute_input":"2022-08-25T05:27:51.328478Z","iopub.status.idle":"2022-08-25T05:27:51.339061Z","shell.execute_reply.started":"2022-08-25T05:27:51.328447Z","shell.execute_reply":"2022-08-25T05:27:51.33748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_volume = pd.read_csv('../input/cspine-overview/df_volume.csv').set_index('sid')\ndf_volume","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:27:51.34165Z","iopub.execute_input":"2022-08-25T05:27:51.342162Z","iopub.status.idle":"2022-08-25T05:27:51.390305Z","shell.execute_reply.started":"2022-08-25T05:27:51.342119Z","shell.execute_reply":"2022-08-25T05:27:51.389137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_volume[(df_volume['read_array_cost'] >12.5) & (abs(df_volume['num_slices'] - 500)<10) ]","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:27:51.393481Z","iopub.execute_input":"2022-08-25T05:27:51.393843Z","iopub.status.idle":"2022-08-25T05:27:51.4205Z","shell.execute_reply.started":"2022-08-25T05:27:51.393814Z","shell.execute_reply":"2022-08-25T05:27:51.418719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_volume[(df_volume['read_array_cost'] <7.5) & (abs(df_volume['num_slices'] - 500)<10) ]","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:27:51.422199Z","iopub.execute_input":"2022-08-25T05:27:51.422568Z","iopub.status.idle":"2022-08-25T05:27:51.446402Z","shell.execute_reply.started":"2022-08-25T05:27:51.422536Z","shell.execute_reply":"2022-08-25T05:27:51.445572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nseries_info = dict()\nfor uid in tqdm(df_volume.index):\n    d = os.path.join(image_dir, uid)\n    idxs = [int(f.split('.')[0]) for f in os.listdir(d)]\n    first = min(idxs)\n    obj = dcmread(os.path.join(d, f'{first}.dcm'))\n    series_info[uid] = {\n        'pixel_representation': int(obj.PixelRepresentation),\n        'bits_stored': int(obj.BitsStored),\n        'patient_id': str(obj.PatientID),\n        'rescale_intercept': float(obj.RescaleIntercept),\n        'header': obj.PixelData[:3],\n    }\ndf_volume['pixel_representation'] = [series_info[uid]['pixel_representation'] for uid in df_volume.index]\ndf_volume['bits_stored'] = [series_info[uid]['bits_stored'] for uid in df_volume.index]\ndf_volume['patient_id'] = [series_info[uid]['patient_id'] for uid in df_volume.index]\ndf_volume['rescale_intercept'] = [series_info[uid]['rescale_intercept'] for uid in df_volume.index]\ndf_volume['header'] = [series_info[uid]['header'] for uid in df_volume.index]","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:27:51.447895Z","iopub.execute_input":"2022-08-25T05:27:51.448197Z","iopub.status.idle":"2022-08-25T05:29:46.583207Z","shell.execute_reply.started":"2022-08-25T05:27:51.44817Z","shell.execute_reply":"2022-08-25T05:29:46.581799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def header_map(v):\n    head_list = [\n        b'0\\xf80',\n        b'\\xfe\\xff\\x00',\n        b'\\x00\\xfc\\x00',\n        b'\\x00\\xf8\\x00',\n    ]\n    if v in head_list:\n        return head_list.index(v)\n    else:\n        return -1\n    \ndef is_utf16(v):\n    return int(v.startswith(b'\\xfe\\xff') or v.startswith(b'\\xff\\xfe'))\n\ndf_volume['slow'] = df_volume['read_array_cost']/df_volume['num_slices'] >= 10./550\ndf_volume['header_tag'] = df_volume['header'].apply(header_map)\ndf_volume['is_utf16'] = df_volume['header'].apply(is_utf16)\n\ndf_volume.to_csv('df_volume.csv')","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:29:46.585507Z","iopub.execute_input":"2022-08-25T05:29:46.585988Z","iopub.status.idle":"2022-08-25T05:29:46.628516Z","shell.execute_reply.started":"2022-08-25T05:29:46.585945Z","shell.execute_reply":"2022-08-25T05:29:46.627251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_volume['patient_id'].value_counts().value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:29:46.630018Z","iopub.execute_input":"2022-08-25T05:29:46.630331Z","iopub.status.idle":"2022-08-25T05:29:46.646727Z","shell.execute_reply.started":"2022-08-25T05:29:46.630303Z","shell.execute_reply":"2022-08-25T05:29:46.645185Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Ncat = 5\ncolors = np.linspace(0, 1, Ncat)\ncmap = dict(enumerate(colors))\n\ntar = 'is_utf16'\nvmin, vmax = min(df_volume[tar]), max(df_volume[tar])\ncolor = df_volume[tar].apply(lambda x: cmap[int((x-vmin)/(vmax-vmin)*(Ncat-1))])\n    \nplt.scatter(df_volume['num_slices'], df_volume['read_array_cost'], marker='.', c=color)\nplt.xlabel('num_slices')\nplt.ylabel('read volume total cost (s)')\ncbar = plt.colorbar()\n\nplt.scatter(*df_volume.loc['1.2.826.0.1.3680043.19388', ['num_slices', 'read_array_cost']], marker='x', s=85)\nplt.scatter(*df_volume.loc['1.2.826.0.1.3680043.25770', ['num_slices', 'read_array_cost']], marker='x', s=85)\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:29:46.65017Z","iopub.execute_input":"2022-08-25T05:29:46.650825Z","iopub.status.idle":"2022-08-25T05:29:47.027589Z","shell.execute_reply.started":"2022-08-25T05:29:46.650781Z","shell.execute_reply":"2022-08-25T05:29:47.02579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vc = df_volume['header'].value_counts()\nvc","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:29:47.029513Z","iopub.execute_input":"2022-08-25T05:29:47.029967Z","iopub.status.idle":"2022-08-25T05:29:47.042575Z","shell.execute_reply.started":"2022-08-25T05:29:47.029922Z","shell.execute_reply":"2022-08-25T05:29:47.040941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\n\ndf_retype = df_volume.copy()\nfor c, typ in df_retype.dtypes.iteritems():\n    if str(typ) == 'bool':\n        df_retype[c] = df_retype[c].astype('int64')\n        \ndf_retype.dtypes\nsns.pairplot(df_retype[['slice_spacing', 'rescale_intercept', 'header_tag', 'is_utf16', 'slow']], hue='slow', diag_kind=\"hist\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:29:47.044278Z","iopub.execute_input":"2022-08-25T05:29:47.044766Z","iopub.status.idle":"2022-08-25T05:29:53.55465Z","shell.execute_reply.started":"2022-08-25T05:29:47.044721Z","shell.execute_reply":"2022-08-25T05:29:53.553197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = df_volume[df_volume['slow']]['read_array_cost'].to_list()\nY = df_volume[df_volume['slow']]['num_slices'].to_list()\nfit_slow = np.polyfit(X, Y, 1)\nprint('[slow] speed:', fit_slow[0], 'intercept:', fit_slow[1])\n\nX = df_volume[~df_volume['slow']]['read_array_cost'].to_list()\nY = df_volume[~df_volume['slow']]['num_slices'].to_list()\nfit_fast = np.polyfit(X, Y, 1)\nprint('[fast] speed:', fit_fast[0], 'intercept:', fit_fast[1])","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:29:53.556248Z","iopub.execute_input":"2022-08-25T05:29:53.55662Z","iopub.status.idle":"2022-08-25T05:29:53.575015Z","shell.execute_reply.started":"2022-08-25T05:29:53.556586Z","shell.execute_reply":"2022-08-25T05:29:53.573445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(10,5))\nNcat = 5\ncolors = np.linspace(0, 1, Ncat)\ncmap = dict(enumerate(colors))\n\ntar = 'is_utf16'\nvmin, vmax = min(df_volume[tar]), max(df_volume[tar])\ncolor = df_volume[tar].apply(lambda x: cmap[int((x-vmin)/(vmax-vmin)*(Ncat-1))])\n    \nplt.scatter(df_volume['num_slices'], df_volume['read_array_cost'], marker='.', c=color)\nplt.xlabel('num_slices')\nplt.ylabel('read volume total cost (s)')\ncbar = plt.colorbar(ticks=[0,1], label='is_utf16')\n\nplt.scatter(*df_volume.loc['1.2.826.0.1.3680043.19388', ['num_slices', 'read_array_cost']], marker='x', s=85, label='slow example: 1.2.826.0.1.3680043.19388')\nplt.scatter(*df_volume.loc['1.2.826.0.1.3680043.25770', ['num_slices', 'read_array_cost']], marker='x', s=85, label='fast example: 1.2.826.0.1.3680043.25770')\n\nx = [0,17.5]\nplt.plot(np.poly1d(fit_slow)(x), x, label='fit speed: slow (~%.0f slices/seconds)' % fit_slow[0])\nplt.plot(np.poly1d(fit_fast)(x), x, label='fit speed: other (~%.0f slices/seconds)'% fit_fast[0])\nplt.legend()\nplt.title('Read Volume Speed: slow when PixelData in UTF-16 encoding')\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T05:29:53.580347Z","iopub.execute_input":"2022-08-25T05:29:53.581207Z","iopub.status.idle":"2022-08-25T05:29:53.999859Z","shell.execute_reply.started":"2022-08-25T05:29:53.58116Z","shell.execute_reply":"2022-08-25T05:29:53.998664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}