{"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":"Dali code is from @tivfrvqhs5, see\nhttps://www.kaggle.com/competitions/rsna-breast-cancer-detection/discussion/371534\n\nHow to use:\n1. make a dictionary of { machine_id: TransferSyntaxUID } using the function make_transfer_syntax_uid()\n\n2. in your dataloader, iterate over each patient id. For (each patient id, machine_id) pair, load all images using function read_image_with_dicomsdl() or read_image_with_dali(), according to dictionary TransferSyntaxUID[machine_id]","metadata":{}},{"cell_type":"code","source":"# todo\n# ! pip download\n! pip install --extra-index-url https://developer.download.nvidia.com/compute/redist/nightly --upgrade nvidia-dali-nightly-cuda110\n! pip install dicomsdl\n! pip install pylibjpeg\n! pip install python-gdcm\n\nprint('install ok')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-12-12T15:24:07.454975Z","iopub.execute_input":"2022-12-12T15:24:07.455344Z","iopub.status.idle":"2022-12-12T15:24:47.373467Z","shell.execute_reply.started":"2022-12-12T15:24:07.455311Z","shell.execute_reply":"2022-12-12T15:24:47.370508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import dicomsdl\ndef __dataset__to_numpy_image(self, index=0):\n    info = self.getPixelDataInfo()\n    dtype = info['dtype']\n    if info['SamplesPerPixel'] != 1:\n        raise RuntimeError('SamplesPerPixel != 1')\n    else:\n        shape = [info['Rows'], info['Cols']]\n    outarr = np.empty(shape, dtype=dtype)\n    self.copyFrameData(index, outarr)\n    return outarr\ndicomsdl._dicomsdl.DataSet.to_numpy_image = __dataset__to_numpy_image\n\n\nimport nvidia.dali as dali\n#https://github.com/NVIDIA/DALI/issues/2330\n#https://medium.datadriveninvestor.com/gpu-accelerated-data-loading-with-dali-part-2-pipelines-and-data-loaders-99f51548e8a6\n#https://github.com/NVIDIA/DALI/blob/main/dali/python/nvidia/dali/plugin/pytorch.py\n#https://zhuanlan.zhihu.com/p/518240063\n\nclass ExternalInputIterator(object):\n    def __init__(self,jpeg_stream):\n        self.jpeg_stream = jpeg_stream\n    def __iter__(self):\n        return self\n    def __next__(self):\n        return self.jpeg_stream\n\n@dali.pipeline_def\ndef j2k_decode_pipeline(eii, is_resize=False):\n    jpeg  = dali.fn.external_source(source=eii, dtype=dali.types.UINT8)\n    image = dali.fn.experimental.decoders.image(jpeg, device='mixed', output_type=dali.types.ANY_DATA, dtype=dali.types.UINT16)\n    if is_resize:\n        image = dali.fn.resize(image, size=[image_size, image_size])\n    return image\n\n##################################################################################\nimport pydicom\nimport pandas as pd\nimport numpy as np\nimport cv2\n\n\nfrom timeit import default_timer as timer\n#from tqdm.notebook import tqdm\nfrom tqdm import tqdm\nfrom joblib import Parallel, delayed\nfrom glob import glob\n\n\ndef time_to_str(t, mode='min'):\n    if mode=='min':\n        t  = int(t)/60\n        hr = t//60\n        min = t%60\n        return '%2d hr %02d min'%(hr,min)\n\n    elif mode=='sec':\n        t   = int(t)\n        min = t//60\n        sec = t%60\n        return '%2d min %02d sec'%(min,sec)\n\n    else:\n        raise NotImplementedError\n\n\nprint('import ok!')","metadata":{"execution":{"iopub.status.busy":"2022-12-12T15:24:47.376585Z","iopub.execute_input":"2022-12-12T15:24:47.377017Z","iopub.status.idle":"2022-12-12T15:24:47.391986Z","shell.execute_reply.started":"2022-12-12T15:24:47.376956Z","shell.execute_reply":"2022-12-12T15:24:47.390571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nimage_size = 1024\n\nmode = 'local-dicom'   #\n\nif 'local' in mode:\n    csv_file = '/kaggle/input/rsna-breast-cancer-detection/train.csv'\n    #csv_file  = '/home/titanx/hengck/share1/kaggle/2022/rsna-breast-mammography/data/valid_df.fold0.csv'\n    dcm_dir   = '/kaggle/input/rsna-breast-cancer-detection/train_images'\n \nif 'submit' in mode:\n    csv_file = '/kaggle/input/rsna-breast-cancer-detection/test.csv'\n    dcm_dir  = '/kaggle/input/rsna-breast-cancer-detection/test_images'\n\n\ndef make_transfer_syntax_uid(df):\n    machine_id_to_transfer = {}\n    machine_id = df.machine_id.unique()\n    for i in machine_id:\n        d = df[df.machine_id==i].iloc[0]\n        f = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n        dicom = pydicom.dcmread(f)\n        machine_id_to_transfer[i]= dicom.file_meta.TransferSyntaxUID\n    return machine_id_to_transfer\n'''\n1.2.840.10008.1.2.4.70 = JPEG Lossless, Nonhierarchical, First- Order Prediction (Processes 14)\n1.2.840.10008.1.2.4.90 = JPEG 2000 Image Compression (Lossless Only)\n'''\n\n\ntest_df = pd.read_csv(csv_file)\nmachine_id_to_transfer = make_transfer_syntax_uid(test_df)\n\ntest_df.loc[:, 'i'] = np.arange(len(test_df))\ntest_df.loc[:,'TransferSyntaxUID']=test_df.machine_id.map(machine_id_to_transfer)\nif 'local' in mode:\n    test_df.loc[:, 'prediction_id'] = test_df.patient_id.astype(str) + '_' + test_df.laterality\n    if 'dicom' in mode:\n        # if 1:\n        test_id = [\n            #  1.2.840.10008.1.2.4.70  count(14)\n            65, 127, 152, 272, 282, 308, 477, 505, 2989, 3542, 7780, 9014, 11094, 11937,\n            #  1.2.840.10008.1.2.4.90  count(26)\n            30, 36, 90, 111, 122, 158, 204, 289, 299, 399, 425, 454, 826, 1703, 1759, 2346, 3021, 4340, 4824, 5059, 5769, 6654, 6658, 7053, 7493, 14292\n        ]\n        test_df = test_df[test_df.patient_id.isin(test_id)].reset_index(drop=True)\nprint(test_df.shape)\nprint(test_df)\n\n\n### read dicom\ndef norm_data(m):\n    max=m.max()\n    min=m.min()\n    m = (m-min)/(max-min+1e-6) #float32\n    m = m.astype(np.float16)\n    return m\n\n\ndef read_image_with_dicomsdl_non_parallel(df): \n    image = []\n    for t,d in df.iterrows():\n        dcm_file = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n        dicom = dicomsdl.open(dcm_file)\n        m = dicom.to_numpy_image()\n        m = norm_data(m)\n        if dicom.PhotometricInterpretation== 'MONOCHROME1':\n            m = 1-m\n        image.append(m)\n\n    image = np.stack(image)\n    return image\n\n\ndef parallel_process(f):\n    dicom = dicomsdl.open(f)\n    m = dicom.to_numpy_image()\n    m = norm_data(m)\n    if dicom.PhotometricInterpretation== 'MONOCHROME1':\n        m = 1-m\n    return m\n\ndef read_image_with_dicomsdl(df):\n\n    #https://stackoverflow.com/questions/56659294/does-joblib-parallel-keep-the-original-order-of-data-passed\n    image =  Parallel(n_jobs=2, backend='multiprocessing')(\n        delayed(parallel_process)(f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm')\n        for t,d in df.iterrows()\n    )\n    image = np.stack(image)\n    return image\n\n\n\n\ndef read_image_with_dali(df):\n\n    jpeg_stream = []\n    invert = []\n    for t,d in df.iterrows():\n        dcm_file = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n        dicom = pydicom.dcmread(dcm_file)\n        s = dicom.PixelData\n        offset = s.find(b'\\x00\\x00\\x00\\x0C')  # <---- the jpeg2000 header info we're looking for\n        jpeg_stream.append(\n            np.array(bytearray(s[offset:]), np.uint8)\n        )\n        invert.append(\n            dicom.PhotometricInterpretation == 'MONOCHROME1'\n        )\n\n    pipe = j2k_decode_pipeline(ExternalInputIterator(jpeg_stream), batch_size=len(jpeg_stream), num_threads=2, device_id=0, debug=True)\n    pipe.build()\n    pipe_out = pipe.run()\n\n    #todo : use feed_ndarray() /nvidia/dali/plugin/pytorch.py to directly convert to gpu pytorch tensor\n    pipe_out = pipe_out[0].as_cpu().as_array().squeeze(-1)\n    #print(pipe_out.shape)\n\n    image = []\n    for i in range(len(df)):\n        m = pipe_out[i]\n        m = norm_data(m)\n        if invert[i]:\n            m = 1-m\n        image.append(m)\n\n    image = np.stack(image)\n    return image","metadata":{"execution":{"iopub.status.busy":"2022-12-12T15:24:47.39384Z","iopub.execute_input":"2022-12-12T15:24:47.394297Z","iopub.status.idle":"2022-12-12T15:24:47.617195Z","shell.execute_reply.started":"2022-12-12T15:24:47.394258Z","shell.execute_reply":"2022-12-12T15:24:47.613257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n##debug function:\nif 1: \n    test_df = test_df[test_df.TransferSyntaxUID == '1.2.840.10008.1.2.4.90'].reset_index(drop=True)\n    patient_id = sorted(test_df.patient_id.unique())\n    print('num_image', len(test_df))\n    print('num_patient', len(patient_id))\n    print('')\n\n    start_timer = timer()\n    image0 = []\n    for i in patient_id:\n        df = test_df[test_df.patient_id==i].reset_index(drop=True)\n        image0.append(read_image_with_dicomsdl(df))\n    print('read_image_with_dicomsdl',time_to_str(timer() - start_timer, 'sec'))\n   \n\n    start_timer = timer()\n    image1 = []\n    for i in patient_id:\n        df = test_df[test_df.patient_id == i].reset_index(drop=True)\n        image1.append(read_image_with_dali(df))\n    print('read_image_with_dali', time_to_str(timer() - start_timer, 'sec'))\n    print('')\n","metadata":{"execution":{"iopub.status.busy":"2022-12-12T15:24:47.619659Z","iopub.execute_input":"2022-12-12T15:24:47.622613Z","iopub.status.idle":"2022-12-12T15:27:00.281916Z","shell.execute_reply.started":"2022-12-12T15:24:47.62254Z","shell.execute_reply":"2022-12-12T15:27:00.280686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\n#check difference\nfor b,i in enumerate(patient_id):\n\n    m0 = image0[b]\n    m1 = image1[b]\n\n    #image = read_image_from_dicom_070(df)\n    \n    diff = m0.astype(np.float32) - m1.astype(np.float32)\n    print(i, m0.shape)\n    print('diff', diff.min(),diff.max(),diff.mean())\n    print('')\n    \n    if b<=4:\n        plt.figure(figsize=(15,20))\n        plt.imshow(np.hstack([m for m in m0]))\n\n        plt.figure(figsize=(15,20))\n        plt.imshow(np.hstack([m for m in m1]))\n\n        plt.show()\n\n  ","metadata":{"execution":{"iopub.status.busy":"2022-12-12T15:27:00.283688Z","iopub.execute_input":"2022-12-12T15:27:00.284835Z","iopub.status.idle":"2022-12-12T15:28:21.78415Z","shell.execute_reply.started":"2022-12-12T15:27:00.28479Z","shell.execute_reply":"2022-12-12T15:28:21.783039Z"},"trusted":true},"execution_count":null,"outputs":[]}]}