{"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"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71447,"databundleVersionId":8208918,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"#%%\nimport pandas as pd\nimport seaborn as sns\nimport glob,pydicom,os,sys\nimport numpy as np\ntra_df = pd.read_csv(\"/home/share/dataset/CTage/train.csv\")\nprint(tra_df.columns)\n# %%\n\n\nprint(tra_df.shape)\n\nsns.histplot(tra_df[\"Age\"],bins=100)\n\n# %%\nex_tra = pd.read_csv(\"/home/share/dataset/CTage/test.csv\")\nprint(ex_tra.shape)\n\ntra_id = tra_df[\"StudyID\"].unique()\nex_id = ex_tra[\"StudyID\"].unique()\nex_id = [str(i).zfill(6) for i in ex_id]\n\ntest = pd.read_csv(\"/home/share/dataset/CTage/submission2.csv\")\ntest_id = test[\"StudyID\"].unique()\ntest_id = [str(i).zfill(6) for i in test_id]\n\n\nprint(len(tra_id),len(ex_id),len(test_id))\n\nprint(set(tra_id)&set(ex_id))\nprint(set(tra_id)&set(test_id))\nprint(tra_id[0],test_id[0])\n\n\nsns.histplot(ex_tra[\"Age\"],bins=100)\n\n\n# %%\ntra_df.head()\n# %%\n\n# %%\nimage_size_seg = (128, 128, 128)\n\n# %%\ndef window_image(img, window_center,window_width, intercept, slope, rescale=True):\n    '''\n    This fucntion came from this notebook https://www.kaggle.com/code/redwankarimsony/ct-scans-dicom-files-windowing-explained\n    If you want to understand more about windowing the referenced notebook is a good read.\n    '''\n    img = (img*slope +intercept) #for translation adjustments given in the dicom file. \n    img_min = window_center - window_width//2 #minimum HU level\n    img_max = window_center + window_width//2 #maximum HU level\n    img[img<img_min] = img_min #set img_min for all HU levels less than minimum HU level\n    img[img>img_max] = img_max #set img_max for all HU levels higher than maximum HU level\n    if rescale: \n        img = (img - img_min) / (img_max - img_min)*255.0 \n    return img\n\ndef normalize_image(image):\n    \"\"\"\n    Normalize image to the range [0, 1].\n    \"\"\"\n    image = image - np.min(image)\n    return image / np.max(image)\ndef load_dicom(dcm, window_center=None, window_width=None):\n    \"\"\"\n    Process a DICOM file and save it as a PNG file.\n    \"\"\"    \n    try:\n        image = dcm.pixel_array\n    except:\n        return\n    if window_center is not None and window_width is not None:\n        \n        image = window_image(image, window_center, window_width, dcm.RescaleIntercept, dcm.RescaleSlope)\n\n    if dcm.PhotometricInterpretation == \"MONOCHROME1\":\n        image = np.invert(image)\n    \n    normalized_image = normalize_image(image)\n    return normalized_image\ndef load_dicom_line_par(path):\n\n    t_paths = sorted(glob.glob(os.path.join(path, \"*\"))[1:], key=lambda x: int(x.split('/')[-1].split(\".\")[0]))\n\n    #indices = np.quantile(list(range(n_scans)), np.linspace(0., 1., image_size_seg[2])).round().astype(int)\n    #t_paths = [t_paths[i] for i in indices]\n    dicoms = [pydicom.dcmread(d, force= True) for d in t_paths]\n    sereis_instance_uid = [d.SeriesInstanceUID for d in dicoms]\n    print(set(sereis_instance_uid))\n        \n    images = {i:[] for i in set(sereis_instance_uid)}\n    z_pos = {i:[] for i in set(sereis_instance_uid)}\n    for dcm,s_id in zip(dicoms,sereis_instance_uid):\n        d = load_dicom(dcm,40,80)\n        #d = load_dicom(dcm,500,2000)\n        if d is None:\n            continue\n        images[s_id].append(d)\n        z_pos[s_id].append(float(dcm.ImagePositionPatient[-1]))\n        \n    #z_pos = [float(d.ImagePositionPatient[-1]) for d in dicoms]\n    for s_id in set(sereis_instance_uid):\n        images[s_id] = np.stack(images[s_id], -1)[:,:,np.argsort(z_pos[s_id])]\n        z_pos[s_id] = np.sort(z_pos[s_id])\n    return images,z_pos\n\n\ndef save_3stack(imgs):\n    new_img = []\n    for i in range(imgs.shape[-1]-2):\n        img = imgs[:,:,i:i+3]\n        print(img.shape)\n       \n\nbase_dir = \"/home/share/dataset/CTage/dataset_jpr_train/dataset_jpr_train\"\noutput_base_dir = \"/home/share/dataset/CTage/window_png/\"  # Update this path as necessary\nos.makedirs(output_base_dir,exist_ok=True)\noutput_z_pos = output_base_dir+\"z_pos\"\noutput_3d = output_base_dir+\"3d\"\nos.makedirs(output_3d,exist_ok=True)\n\nos.makedirs(output_z_pos,exist_ok=True)\n\nfor root_dir in [\"3\"]:\n    root_path = os.path.join(base_dir, root_dir)\n    for accession_number in os.listdir(root_path)[::-1]: \n        if len(glob.glob(f\"{output_3d}/{accession_number}*\"))>0:\n            continue\n        accession_path = os.path.join(root_path, accession_number)\n        img_3d,z_pos = load_dicom_line_par(accession_path)\n        for i in img_3d.keys():\n            #np.save(f\"{output_z_pos}/{accession_number}____{i}.npy\",z_pos[i])\n            np.save(f\"{output_3d}/{accession_number}____{i}.npy\",img_3d[i])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport numpy as np\nimport glob,cv2,multiprocessing\nfrom PIL import Image\nimport matplotlib.pyplot as plt\n\n\"\"\"\n3dで保存したwindow済みの画像を25dに変換する\n\n\n\"\"\"\noutput_base_dir = \"/home/share/dataset/CTage/window_png/\"  # Update this path as necessary\noutput_z_pos = output_base_dir+\"z_pos\"\noutput_3d = output_base_dir+\"3d\"\nimgs = glob.glob(f\"{output_3d}/*\")\n\ndef func(img):\n    x = 255-cv2.morphologyEx(cv2.cvtColor(img, cv2.COLOR_BGR2HSV)[:,:,1], cv2.MORPH_CLOSE, np.ones((5,5),np.uint8))\n    retval, im_bw = cv2.threshold(x, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)\n    contours, hierarchy = cv2.findContours(255-im_bw, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n    max_rect = 0\n    xywh = [0,0,1,1]\n    for i in range(len(contours)):\n        \n        x, y, w, h = cv2.boundingRect(contours[i])\n        rect = w*h\n        if max_rect < rect:\n            max_rect=rect\n            xywh = [x, y, w, h]\n    size_ = xywh[2]*xywh[3]/(img.shape[0]*img.shape[1])\n\n        \n    return xywh,size_\n\noutput_25d = output_base_dir+\"25d\"\nos.makedirs(output_25d,exist_ok=True)\n\ndef save_3stack(path):\n    imgs = np.load(path)\n    p = path.split(\"/\")[-1].replace(\".npy\",\"\")#studyid____seriesid\n    new_img = []\n    \n    os.makedirs(f\"{output_25d}/{p}\",exist_ok=True)\n    for i in range(imgs.shape[-1]-2):\n        img = (imgs[:,:,i:i+3]*255.).astype(np.uint8)\n        cv2.imwrite(f\"{output_25d}/{p}/{p}_{i}.png\",img)\n    return \n\noutput_25d_2 = output_base_dir+\"25d_2\"\nos.makedirs(output_25d_2,exist_ok=True)\ndef save_3stack_2(path):\n    imgs = np.load(path)\n    p = path.split(\"/\")[-1].replace(\".npy\",\"\")#studyid____seriesid\n    new_img = []\n    \n    os.makedirs(f\"{output_25d_2}/{p}\",exist_ok=True)\n    for i in range(imgs.shape[-1]-4):\n        img = (imgs[:,:,[i,i+2,i+4]]*255.).astype(np.uint8)\n        cv2.imwrite(f\"{output_25d_2}/{p}/{p}_{i}.png\",img)\n    return \n\noutput_25d_3 = output_base_dir+\"25d_3\"\nos.makedirs(output_25d_3,exist_ok=True)\ndef save_5stack_3(path):\n    imgs = np.load(path)\n    p = path.split(\"/\")[-1].replace(\".npy\",\"\")#studyid____seriesid\n    new_img = []\n    \n    os.makedirs(f\"{output_25d_3}/{p}\",exist_ok=True)\n    for i in range(imgs.shape[-1]-10):\n        img = (imgs[:,:,[i,i+2,i+4,i+6,i+10]]*255.).astype(np.uint8)\n        np.save(f\"{output_25d_3}/{p}/{p}_{i}.npy\",img)\n    return \n\noutput_25d_4 = output_base_dir+\"25d_4\"\nos.makedirs(output_25d_4,exist_ok=True)\n\ndef save_7stack_int2(path):\n    imgs = np.load(path)\n    p = path.split(\"/\")[-1].replace(\".npy\",\"\")#studyid____seriesid\n    new_img = []\n    \n    os.makedirs(f\"{output_25d_4}/{p}\",exist_ok=True)\n    os.makedirs(f\"{output_25d_4}/{p}\".replace(\"25d_4\",\"25d_4_xywh\"),exist_ok=True)\n    for i in range(imgs.shape[-1]-14):\n        img = (imgs[:,:,[i,i+2,i+4,i+6,i+10,i+12,i+14]]*255.).astype(np.uint8)\n        np.save(f\"{output_25d_4}/{p}/{p}_{i}.npy\",img)\n        \n        xywh,size_ = func(img[:,:,[3,4,5]])\n        xywh.append(size_)\n        path_npy = f\"{output_25d_4}/{p}/{p}_{i}.npy\".replace(\"25d_4\",\"25d_4_xywh\")\n        #print(path_npy,xywh)\n        np.save(path_npy,np.array(xywh))\n        #print(img.shape)\n    #exit()\n        \n    return \n    \nimport time\nwith multiprocessing.Pool(10) as p:\n    p.map(save_7stack_int2,imgs)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#%%\nimport numpy as np\nimport cv2,glob,os,multiprocessing\nimport matplotlib.pyplot as plt\n\nimg_paths = glob.glob(\"/home/share/dataset/CTage/window_png/25d_2/*\")\n\ndef func(img):\n    x = 255-cv2.morphologyEx(cv2.cvtColor(img, cv2.COLOR_BGR2HSV)[:,:,1], cv2.MORPH_CLOSE, np.ones((5,5),np.uint8))\n    retval, im_bw = cv2.threshold(x, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)\n    contours, hierarchy = cv2.findContours(255-im_bw, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n    max_rect = 0\n    xywh = [0,0,1,1]\n    for i in range(len(contours)):\n        \n        x, y, w, h = cv2.boundingRect(contours[i])\n        rect = w*h\n        if max_rect < rect:\n            max_rect=rect\n            xywh = [x, y, w, h]\n    size_ = xywh[2]*xywh[3]/(img.shape[0]*img.shape[1])\n\n        \n    return xywh,size_\n\n\ndef crop(paths_):\n    os.makedirs(paths_.replace(\"25d_2\",\"25d_2_xywh\"),exist_ok=True)\n    \n    paths_ =  glob.glob(f\"{paths_}/*\")\n    for p in paths_:\n\n        img  = cv2.imread(p)\n\n        xywh,size_ = func(img)\n        xywh.append(size_)\n        path_npy = p.replace(\"25d_2\",\"25d_2_xywh\").replace(\".png\",\".npy\")\n        #print(path_npy,xywh)\n        np.save(path_npy,np.array(xywh))\n        \ndef crop_5(paths_):\n    # 5ch stack\n    os.makedirs(paths_.replace(\"25d_3\",\"25d_3_xywh\"),exist_ok=True)\n    \n    paths_ =  glob.glob(f\"{paths_}/*\")\n    print(paths_)\n    for p in paths_:\n\n        img  = np.load(p)[:,:,[1,2,3]]\n\n        xywh,size_ = func(img)\n        xywh.append(size_)\n        path_npy = p.replace(\"25d_3\",\"25d_3_xywh\").replace(\".npy\",\".npy\")\n        #print(path_npy,xywh)\n        np.save(path_npy,np.array(xywh))\n        \ndef crop_7(paths_):\n    # 7ch stack\n    os.makedirs(paths_.replace(\"25d_4\",\"25d_4_xywh\"),exist_ok=True)\n    \n    paths_ =  glob.glob(f\"{paths_}/*\")\n    for p in paths_:\n\n        img  = np.load(p)[:,:,[3,4,5]]\n\n        xywh,size_ = func(img)\n        xywh.append(size_)\n        path_npy = p.replace(\"25d_4\",\"25d_4_xywh\").replace(\".npy\",\".npy\")\n        #print(path_npy,xywh)\n        np.save(path_npy,np.array(xywh))\n        \nwith multiprocessing.Pool(4) as p:\n    p.map(crop,img_paths[:])\n\n","metadata":{},"execution_count":null,"outputs":[]}]}