{"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":"[Nvidia Dali ](http://https://docs.nvidia.com/deeplearning/dali/user-guide/docs/index.html) is a GPU based library for very fast data loading and preprocessing.  It contains GPU based image decoders which can be used for fast and parallel decoding of jpeg2000 images.  This notebook contains a minimal example of extracting the jpeg2000 encoded images conatained in a dicom container and decoding them on GPU.","metadata":{}},{"cell_type":"markdown","source":"![](https://arcwiki.rs.gsu.edu/nvidia_dali_pipeline.png)","metadata":{}},{"cell_type":"markdown","source":"**Note:** Need to use a nightly build of DALI because UINT16 support was recently added and is not in the main prod wheel","metadata":{}},{"cell_type":"code","source":"!pip install -qU python-gdcm pydicom pylibjpeg\n!pip install --extra-index-url https://developer.download.nvidia.com/compute/redist/nightly --upgrade nvidia-dali-nightly-cuda110","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:28:37.469425Z","iopub.execute_input":"2022-12-10T16:28:37.469835Z","iopub.status.idle":"2022-12-10T16:29:16.552221Z","shell.execute_reply.started":"2022-12-10T16:28:37.46975Z","shell.execute_reply":"2022-12-10T16:29:16.551061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport glob, os\nimport pydicom\nfrom pydicom.filebase import DicomBytesIO\nfrom tqdm.notebook import tqdm\nfrom joblib import Parallel, delayed\n\nfrom nvidia.dali import pipeline_def\nimport nvidia.dali.fn as fn\nimport nvidia.dali.types as types\nfrom nvidia.dali.types import DALIDataType","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:16.555183Z","iopub.execute_input":"2022-12-10T16:29:16.555622Z","iopub.status.idle":"2022-12-10T16:29:17.033147Z","shell.execute_reply.started":"2022-12-10T16:29:16.555534Z","shell.execute_reply":"2022-12-10T16:29:17.032038Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#here's the magic of hacking the jpeg2000 encoded bitstream.  Function saves jp2 encoded image files contained within dicom which will be later decoded by DALI \ndef convert_dicom_to_j2k(file):\n    patient = file.split('/')[-2]\n    image = file.split('/')[-1][:-4]\n    dcmfile = pydicom.dcmread(f'../input/rsna-breast-cancer-detection/train_images/{file}')\n    if dcmfile.file_meta.TransferSyntaxUID=='1.2.840.10008.1.2.4.90':\n        with open(f'../input/rsna-breast-cancer-detection/train_images/{file}', 'rb') as fp:\n            raw = DicomBytesIO(fp.read())\n            ds = pydicom.dcmread(raw)\n        offset = ds.PixelData.find(b\"\\x00\\x00\\x00\\x0C\")  #<---- the jpeg2000 header info we're looking for\n        hackedbitstream = bytearray()\n        hackedbitstream.extend(ds.PixelData[offset:])\n        with open(f\"../working/{patient}_{image}.jp2\", \"wb\") as binary_file:\n            binary_file.write(hackedbitstream)","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:17.034675Z","iopub.execute_input":"2022-12-10T16:29:17.035044Z","iopub.status.idle":"2022-12-10T16:29:17.044348Z","shell.execute_reply.started":"2022-12-10T16:29:17.035008Z","shell.execute_reply":"2022-12-10T16:29:17.042377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv('/kaggle/input/dali-rsna-helpers/dicom_jpg_formats.csv')\ndel df['label']\ndf['uid'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:32:07.523548Z","iopub.execute_input":"2022-12-10T16:32:07.523951Z","iopub.status.idle":"2022-12-10T16:32:07.576034Z","shell.execute_reply.started":"2022-12-10T16:32:07.523894Z","shell.execute_reply":"2022-12-10T16:32:07.574918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1.2.840.10008.1.2.4.70  ----> JPEG Lossless, Nonhierarchical, First- Order Prediction (Processes 14)\n\n1.2.840.10008.1.2.4.90  ----> JPEG 2000 Image Compression (Lossless Only) \n\n**For now this solution applies only to the JPEG 2000 Image Compression standard**\n","metadata":{}},{"cell_type":"code","source":"allj2kdicoms = df[df['uid']=='1.2.840.10008.1.2.4.90']['dicom'].tolist()","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:17.164162Z","iopub.execute_input":"2022-12-10T16:29:17.165038Z","iopub.status.idle":"2022-12-10T16:29:17.178511Z","shell.execute_reply.started":"2022-12-10T16:29:17.165002Z","shell.execute_reply":"2022-12-10T16:29:17.176368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's get our initial baseline of time to decode a single image using pydicom (CPU based decoding).","metadata":{}},{"cell_type":"code","source":"def pydicom_benchmark(file):\n    dicom = pydicom.dcmread(f'../input/rsna-breast-cancer-detection/train_images/{file}')\n    img = dicom.pixel_array","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:17.179861Z","iopub.execute_input":"2022-12-10T16:29:17.18046Z","iopub.status.idle":"2022-12-10T16:29:17.186531Z","shell.execute_reply.started":"2022-12-10T16:29:17.180427Z","shell.execute_reply":"2022-12-10T16:29:17.185428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_ = Parallel(n_jobs=2)(delayed(pydicom_benchmark)(uid) for uid in tqdm(allj2kdicoms[:32]))","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:17.188187Z","iopub.execute_input":"2022-12-10T16:29:17.188514Z","iopub.status.idle":"2022-12-10T16:29:54.53177Z","shell.execute_reply.started":"2022-12-10T16:29:17.188483Z","shell.execute_reply":"2022-12-10T16:29:54.530645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now to the GPU based vesion. \n\nStep 1: convert dicoms to .jp2 files (jpeg2000 encoded images)","metadata":{}},{"cell_type":"code","source":"_ = Parallel(n_jobs=2)(delayed(convert_dicom_to_j2k)(uid) for uid in tqdm(allj2kdicoms[:32]))","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:54.533314Z","iopub.execute_input":"2022-12-10T16:29:54.533704Z","iopub.status.idle":"2022-12-10T16:29:54.919207Z","shell.execute_reply.started":"2022-12-10T16:29:54.53366Z","shell.execute_reply":"2022-12-10T16:29:54.918094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Step 2: create [DALI pipeline](https://docs.nvidia.com/deeplearning/dali/user-guide/docs/pipeline.html).  \n\nAll we're doing in the pipeline is reading and decoding the j2k files and returining the uint16 array (equivalent to pydicom.pixel_array function)","metadata":{}},{"cell_type":"code","source":"j2kfiles = [f'../working/{thing.split(\"/\")[-2]}_{thing.split(\"/\")[-1][:-4]}.jp2' for thing in allj2kdicoms[:32]]","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:54.920659Z","iopub.execute_input":"2022-12-10T16:29:54.921072Z","iopub.status.idle":"2022-12-10T16:29:54.930099Z","shell.execute_reply.started":"2022-12-10T16:29:54.921038Z","shell.execute_reply":"2022-12-10T16:29:54.929072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"@pipeline_def\ndef j2k_decode_pipeline():\n    jpegs, _ = fn.readers.file(files = j2kfiles)\n    images = fn.experimental.decoders.image(jpegs, device='mixed', output_type=types.ANY_DATA, dtype=DALIDataType.UINT16)\n    return images","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:54.934762Z","iopub.execute_input":"2022-12-10T16:29:54.9351Z","iopub.status.idle":"2022-12-10T16:29:54.944132Z","shell.execute_reply.started":"2022-12-10T16:29:54.935075Z","shell.execute_reply":"2022-12-10T16:29:54.943098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"max_batch_size = 32\npipe = j2k_decode_pipeline(batch_size=max_batch_size, num_threads=2, device_id=0, debug=True)\npipe.build()","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:54.945557Z","iopub.execute_input":"2022-12-10T16:29:54.946002Z","iopub.status.idle":"2022-12-10T16:29:55.156984Z","shell.execute_reply.started":"2022-12-10T16:29:54.945966Z","shell.execute_reply":"2022-12-10T16:29:55.155964Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\npipe_out = pipe.run()","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:55.160832Z","iopub.execute_input":"2022-12-10T16:29:55.161174Z","iopub.status.idle":"2022-12-10T16:29:57.111366Z","shell.execute_reply.started":"2022-12-10T16:29:55.161145Z","shell.execute_reply":"2022-12-10T16:29:57.110525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Wow nice speedup, that's the time to decode **32 images** AND the data is sitting on the GPU already.  How much memory are we using?","metadata":{}},{"cell_type":"code","source":"!nvidia-smi","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:57.112843Z","iopub.execute_input":"2022-12-10T16:29:57.113258Z","iopub.status.idle":"2022-12-10T16:29:58.95448Z","shell.execute_reply.started":"2022-12-10T16:29:57.11322Z","shell.execute_reply":"2022-12-10T16:29:58.953363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's test and make sure the array we get from pydicom is the same as what DALI gives us after decode","metadata":{}},{"cell_type":"code","source":"pydicom_image = pydicom.dcmread(f'../input/rsna-breast-cancer-detection/train_images/{allj2kdicoms[30]}').pixel_array","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:58.956528Z","iopub.execute_input":"2022-12-10T16:29:58.956974Z","iopub.status.idle":"2022-12-10T16:29:59.552274Z","shell.execute_reply.started":"2022-12-10T16:29:58.95693Z","shell.execute_reply":"2022-12-10T16:29:59.55121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pydicom_image","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:59.553939Z","iopub.execute_input":"2022-12-10T16:29:59.554328Z","iopub.status.idle":"2022-12-10T16:29:59.564576Z","shell.execute_reply.started":"2022-12-10T16:29:59.554288Z","shell.execute_reply":"2022-12-10T16:29:59.563436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"testfile = [f'../working/{thing.split(\"/\")[-2]}_{thing.split(\"/\")[-1][:-4]}.jp2' for thing in [allj2kdicoms[30]]]\n\n\n@pipeline_def\ndef j2k_decode_pipeline():\n    jpegs, _ = fn.readers.file(files = testfile)\n    images = fn.experimental.decoders.image(jpegs, device='mixed', output_type=types.ANY_DATA, dtype=DALIDataType.UINT16)\n    return images\n\nmax_batch_size = 1\npipe = j2k_decode_pipeline(batch_size=max_batch_size, num_threads=2, device_id=0, debug=True)\npipe.build()\n\npipe_out = pipe.run()\n\ndali_image = pipe_out[0].as_cpu().as_array()[0,:,:,0]","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:59.566587Z","iopub.execute_input":"2022-12-10T16:29:59.567532Z","iopub.status.idle":"2022-12-10T16:29:59.694468Z","shell.execute_reply.started":"2022-12-10T16:29:59.567493Z","shell.execute_reply":"2022-12-10T16:29:59.693485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dali_image","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:59.695717Z","iopub.execute_input":"2022-12-10T16:29:59.698051Z","iopub.status.idle":"2022-12-10T16:29:59.706574Z","shell.execute_reply.started":"2022-12-10T16:29:59.698023Z","shell.execute_reply":"2022-12-10T16:29:59.705904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.all(dali_image==pydicom_image)","metadata":{"execution":{"iopub.status.busy":"2022-12-10T16:29:59.708163Z","iopub.execute_input":"2022-12-10T16:29:59.709051Z","iopub.status.idle":"2022-12-10T16:29:59.721509Z","shell.execute_reply.started":"2022-12-10T16:29:59.709014Z","shell.execute_reply":"2022-12-10T16:29:59.720648Z"},"trusted":true},"execution_count":null,"outputs":[]}]}