{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"}],"dockerImageVersionId":30761,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"from scipy.linalg import norm\nimport numpy as np\nfrom math import atan\nfrom math import sin, cos\nimport cv2\nimport pydicom\n\n\nimport matplotlib\nimport matplotlib.pyplot as plt\n\nprint('import ok!')","metadata":{"execution":{"iopub.status.busy":"2024-08-23T02:25:48.098016Z","iopub.execute_input":"2024-08-23T02:25:48.098607Z","iopub.status.idle":"2024-08-23T02:25:48.107178Z","shell.execute_reply.started":"2024-08-23T02:25:48.098534Z","shell.execute_reply":"2024-08-23T02:25:48.105544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# shape alignment code from\n# https://medium.com/@olga_kravchenko/generalized-procrustes-analysis-with-python-numpy-c571e8e8a421\n\nREFERENCE_H = 512\nREFERENCE_W = 512\n\nlevel_name_to_color={\n    'L1/L2':[255,0,0],\n    'L2/L3':[0,255,0],\n    'L3/L4':[0,0,255],\n    'L4/L5':[255,255,0],\n    'L5/S1':[0,255,255],\n}\nlevel_color= [\n    v for k,v in level_name_to_color.items()\n]\nlevel_color1=[\n    [128,0,0],\n    [0,128,0],\n    [0,0,128],\n    [128,128,0],\n    [0,128,128],\n]\n\n\n\ndef draw_shape(\n    shape,  #x,y,x,y,x,y ... format\n    image=None, \n    is_shift_center=True, \n    is_line=True, line_thickness=1,\n    is_circle=True, circle_radius=8,\n    level_color=level_color\n):\n    \n    if image is None:\n        H,W = REFERENCE_H,REFERENCE_W\n        image = np.zeros((H,W,3))\n    else:\n        H,W = image.shape[:2]\n\n    point = shape.reshape(5,2,2) #assume 5x2=10 points\n    if is_shift_center:\n        point = point + np.array([W//2,H//2]).reshape(1,1,2)\n\n    point = (np.round(point)).astype(np.int32)\n    for i in range(5):\n        color = level_color[i]\n        px0,py0 = point[i,0]\n        px1,py1 = point[i,1]\n        if is_circle:\n            cv2.circle(image, (px0, py0), circle_radius, color, -1, cv2.LINE_AA)\n            cv2.circle(image, (px1, py1), circle_radius, color, -1, cv2.LINE_AA)\n        \n        if is_line:\n            cv2.line(image, (px0, py0), (px1, py1), color, line_thickness, cv2.LINE_AA) \n            if i!=0:\n                qx0,qy0 = point[i-1,0]\n                qx1,qy1 = point[i-1,1]\n                cv2.line(image, (px0, py0), (qx0, qy0), color, line_thickness, cv2.LINE_AA)\n                cv2.line(image, (px1, py1), (qx1, qy1), color, line_thickness, cv2.LINE_AA)\n    return image\n\n\n\ndef get_rotation_scale(reference_shape, shape, num_point=10):\n    \n    reference_shape = reference_shape.reshape(-1) \n    shape = shape.reshape(-1)\n\n    a = np.dot(shape, reference_shape) / norm(reference_shape) ** 2\n\n    # separate x and y for the sake of convenience\n    ref_x = reference_shape[::2]\n    ref_y = reference_shape[1::2]\n\n    x = shape[::2]\n    y = shape[1::2]\n\n    b = np.sum(x * ref_y - ref_x * y) / norm(reference_shape) ** 2\n\n    scale = np.sqrt(a ** 2 + b ** 2)\n    theta = atan(b / max(a, 10 ** -10))  # avoid dividing by 0\n\n    return scale, theta\n\n\ndef get_rotation_matrix(theta):\n    return np.array([[cos(theta), -sin(theta)], [sin(theta), cos(theta)]])\n\n\ndef do_zero_mean(shape):\n    shape = shape.reshape(-1,2)\n    shape = shape-shape.mean(0,keepdims=True)\n    shape = shape.reshape(-1)\n    return shape\n\n\ndef do_scale(shape, scale):\n    return shape / scale\n\n\ndef do_rotate(shape, theta): \n    mat = get_rotation_matrix(theta) \n    shape = shape.reshape((-1, 2)).T \n    rotated_shape = np.dot(mat, shape)\n    rotated_shape = rotated_shape.T.reshape(-1)\n    return rotated_shape\n\n\ndef do_align_shape(shape, reference_shape, is_rotate=True):\n\n    reference_shape = np.copy(reference_shape)\n    shape = np.copy(shape)\n\n    # get scale and rotation\n    scale, theta = get_rotation_scale(reference_shape, shape)\n \n    aligned_shape = shape / scale\n    if is_rotate:\n        aligned_shape = do_rotate(aligned_shape, theta)\n\n    return aligned_shape\n\n\n#---helper ---\ndef convert_to_8bit(x):\n    lower, upper = np.percentile(x, (1, 99))\n    x = np.clip(x, lower, upper)\n    x = x - np.min(x)\n    x = x / np.max(x) \n    return (x * 255).astype(\"uint8\")\n\ndef read_dicom_as_image(dicom_file):\n    \n    d = pydicom.dcmread(dicom_file)\n    image = d.pixel_array.astype(\"float32\")\n    image = convert_to_8bit(image)\n    return image\n\nclass dotdict(dict):\n    __setattr__ = dict.__setitem__\n    __delattr__ = dict.__delitem__\n\n    def __getattr__(self, name):\n        try:\n            return self[name]\n        except KeyError:\n            raise AttributeError(name)\n            \nprint('helper ok!')","metadata":{"execution":{"iopub.status.busy":"2024-08-23T02:25:48.109472Z","iopub.execute_input":"2024-08-23T02:25:48.109996Z","iopub.status.idle":"2024-08-23T02:25:48.150568Z","shell.execute_reply.started":"2024-08-23T02:25:48.109950Z","shell.execute_reply":"2024-08-23T02:25:48.149194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#some dummy data from\n# https://www.kaggle.com/competitions/rsna-2024-lumbar-spine-degenerative-classification/discussion/528653\ndummy_data=[\n    dotdict(\n        study_id = 1217477368,\n        series_id = 370109,  \n        instance_number =[9],\n        point=np.array(\n            [[200.        , 134.        ],\n             [288.86107085, 144.5533225 ],\n             [188.        , 204.        ],\n             [279.55700992, 211.74931814],\n             [174.        , 272.        ],\n             [270.76984115, 274.810176  ],\n             [172.        , 344.        ],\n             [277.48944077, 324.43183411],\n             [200.        , 408.        ],\n             [285.2428247 , 375.08727706]]\n        ),\n    ),\n    dotdict(\n        study_id = 2155667219,\n        series_id = 2142829,  \n        instance_number =[8],\n        point=np.array(\n            [[ 81.        ,  77.        ],\n             [137.31788928,  90.70113075],\n             [ 72.        , 117.        ],\n             [128.2193431 , 118.82390938],\n             [ 74.        , 160.        ],\n             [126.97863245, 152.3231017 ],\n             [ 82.        , 197.        ],\n             [136.07717837, 175.89660749],\n             [ 95.        , 229.        ],\n             [149.31142733, 203.1922455 ]]\n        ),\n    ),\n    dotdict(\n        study_id = 425970461,\n        series_id = 8693307,  \n        instance_number =[10],\n        point=np.array(\n            [[135.       ,   46.25      ],\n             [180.7934336,   67.4145008 ],\n             [120.       ,   80.        ],\n             [166.7852256,   98.49521216],\n             [105.        , 120.        ],\n             [155.84131328, 133.95348832],\n             [ 95.        , 166.25      ],\n             [150.5882352 , 164.15868672],\n             [105.        , 217.5       ],\n             [159.34336544, 192.6128592 ]]\n        ),\n    ),\n    dotdict(\n        study_id = 376723024,\n        series_id = 9334563,  \n        instance_number =[8],\n        point=np.array(\n            [[273.        , 226.3125    ],\n             [410.00000026, 244.99999962],\n             [264.        , 342.796875  ],\n             [398.00000026, 334.99999978],\n             [261.        , 459.28125   ],\n             [415.76475955, 429.25688076],\n             [288.        , 575.765625  ],\n             [435.09713203, 510.534375  ],\n             [348.        , 692.25      ],\n             [465.71588198, 573.103125  ]]\n        ),\n    ),\n]\n\n\ndata_shape = []\ndata_image = []\nfor i in range(len(dummy_data)):\n    d = dummy_data[i]\n\n    image = [] \n    for n in d.instance_number: \n        dicom_file =f'/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images/{d.study_id}/{d.series_id}/{n}.dcm'\n        m = read_dicom_as_image(dicom_file)\n        m = cv2.cvtColor(m, cv2.COLOR_GRAY2RGB)\n        image.append(m)\n\n    image = np.stack(image).mean(0)\n    image = image.astype(np.uint8)\n    shape = d.point.reshape(-1,2)\n    \n    #for debug and illustrative purpose, we mark the image\n    image =  draw_shape(\n        shape,  #x,y,x,y,x,y ... format\n        image=image, \n        is_shift_center=False, \n        is_line=False, line_thickness=1,\n        is_circle=True, circle_radius=8,\n        level_color=level_color\n    )\n    data_image.append(image)\n    data_shape.append(shape)\n    if i==0:\n        plt.imshow(image, cmap='gray')\n        plt.show()\n\ndata_shape = np.stack(data_shape)\nnum_shape=N=len(data_shape)\nnum_point=P=10","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-08-23T02:25:48.152277Z","iopub.execute_input":"2024-08-23T02:25:48.152804Z","iopub.status.idle":"2024-08-23T02:25:48.590181Z","shell.execute_reply.started":"2024-08-23T02:25:48.152747Z","shell.execute_reply":"2024-08-23T02:25:48.588924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. make reference shape","metadata":{}},{"cell_type":"code","source":"# make mean shape (reference shape)\nmean = data_shape.mean(0).reshape(-1)\nzero_mean = do_zero_mean(mean)\n\nprint('mean', mean.shape)\nprint('zero_mean', zero_mean.shape)\n\naligned_shape = []\nfor i in range(N):\n    #print(i)\n    p = data_shape[i]\n    p = do_zero_mean(p)\n    p_align = do_align_shape(p, zero_mean)\n    aligned_shape.append(p_align)\n\naligned_shape = np.stack(aligned_shape)\nprint('aligned_shape', aligned_shape.shape)\nzero_mean = aligned_shape.mean(0)\nzero_mean = do_zero_mean(zero_mean)\n\nprint('zero_mean shape\\n', zero_mean)\nmean_overlay =  draw_shape(\n    zero_mean,  #x,y,x,y,x,y ... format\n    image=None, \n    is_shift_center=True, \n    is_line=True, line_thickness=2,\n    is_circle=True, circle_radius=8,\n    level_color=level_color\n)\nplt.imshow(mean_overlay, cmap='gray')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-23T02:25:48.592765Z","iopub.execute_input":"2024-08-23T02:25:48.593192Z","iopub.status.idle":"2024-08-23T02:25:48.870026Z","shell.execute_reply.started":"2024-08-23T02:25:48.593149Z","shell.execute_reply":"2024-08-23T02:25:48.868754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. DEMO: make normalised data\n- resize to e.g. 512x512\n- align kaggle data to reference shape\n- you just need to compute the affine parameters (scale, rotation, shift). You can leave the image and points untouched.","metadata":{}},{"cell_type":"code","source":"#we usually normalised training data to fixed size. here we use reference size 512x512\n\nimage_512=[]\nshape_512=[]\n\nfor i in range(N):\n    #print(i)\n    p = data_shape[i].reshape(-1)\n    m = data_image[i]\n \n    H,W = m.shape[:2]\n    s = 512/max(H,W)\n     \n    p_512 = p*s \n    m_512 = cv2.resize(m, dsize=None,fx=s,fy=s)\n    h,w = m_512.shape[:2]\n    m_512 = np.pad(m_512,[[0,512-h],[0,512-w], [0,0]],mode='constant',constant_values=0)\n    \n    image_512.append(m_512)\n    shape_512.append(p_512)\n    \n# compute affine parameters\nmean_512 = (zero_mean.reshape(-1,2) + [[256,256]]).reshape(-1)\noverlay_512_overlay = mean_overlay\n\naffine =[]\nfor i in range(N):\n    p_512 = shape_512[i].copy()\n    \n    #get kaggle to mean   \n    #print(p_512.shape, mean_512.shape)\n    mat, inlier = cv2.estimateAffinePartial2D(p_512.reshape(-1,2), mean_512.reshape(-1,2))\n\n    #nomalised shape\n    p_align = np.concatenate([p_512.reshape(-1,2), np.ones((10,1))], axis=1)@mat.T\n    #recover parameter to make normalised shape\n    mat, inlier = cv2.estimateAffinePartial2D(p_512.reshape(-1,2),p_align.reshape(-1,2))\n    scale = np.sqrt(mat[0,0]**2+mat[0,1]**2)\n    theta = np.arctan2(mat[1,0,],mat[0,0])/np.pi*180\n    shift = mat[:,2].tolist()\n    affine.append((scale,theta, shift))\n    \n#--------------------------------\n#let's visualise and debug\nshape_aligned=[]\nimage_aligned=[]\nfor i in range(N):\n    m_512 = image_512[i].copy()\n    p_512 = shape_512[i].copy()\n    \n    scale, theta, shift = affine[i]\n    theta = theta/180*np.pi\n    mat = np.array([\n        scale*cos(theta), -scale*sin(theta), shift[0],\n        scale*sin(theta),  scale*cos(theta), shift[1],\n    ]).reshape(2,3)\n    \n    p_align = np.concatenate([p_512.reshape(-1,2), np.ones((10,1))], axis=1)@mat.T\n    p_align = p_align.reshape(-1)\n    m_align = cv2.warpAffine(m_512, mat, (512, 512))\n    m_align =  draw_shape(\n        p_align,  #x,y,x,y,x,y ... format\n        image=m_align, \n        is_shift_center=False, \n        is_line=True, line_thickness=2,\n        is_circle=False, circle_radius=8,\n        level_color=[[255,0,0]]*5\n    )\n\n    shape_aligned.append(p_align)\n    image_aligned.append(m_align)\n    \n     \n    \nimage_aligned = np.stack(image_aligned)\nmean_image_aligned = image_aligned.mean(0).astype(np.uint8)\nplt.imshow(mean_image_aligned, cmap='gray')\nplt.show()\n\n\n    \nshape_aligned = np.stack(shape_aligned)\nmean_shape_aligned = shape_aligned.mean(0)\n\nprint('mean_shape_aligned\\n', mean_shape_aligned)\nprint('mean_512\\n', mean_512)","metadata":{"execution":{"iopub.status.busy":"2024-08-23T02:25:48.871698Z","iopub.execute_input":"2024-08-23T02:25:48.872093Z","iopub.status.idle":"2024-08-23T02:25:49.263076Z","shell.execute_reply.started":"2024-08-23T02:25:48.872051Z","shell.execute_reply":"2024-08-23T02:25:49.261758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. DEMO: how to use affine paramters for augmentation\n- generate parameters: random scale, rotation, shift\n- combine random parameters with affine parameters\n- apply new parameters to image and shape","metadata":{}},{"cell_type":"code","source":"for i in range(N):\n    m_512 = image_512[i].copy()\n    p_512 = shape_512[i].copy()\n    \n    scale0, theta0, shift0 = affine[i]\n       \n    shape_augment=[]\n    image_augment=[]\n    \n    # 1 orginal + 1 normalised + 5 augmentation\n    for j in range(7):\n        if j==0:\n            mat = np.array([1,0,0,0,1,0], dtype=np.float32).reshape(2,3)\n        else:\n            if j==1:\n                scale,theta, shift= scale0,theta0, shift0\n            else:\n                dscale = np.random.uniform(-0.25,0.25)\n                dtheta = np.random.uniform(-20,20)\n                dshift = np.random.uniform(-30,30,2).tolist()\n                #print(dscale,dtheta,dshift)\n                scale = scale0 + dscale\n                theta = theta0 + dtheta\n                shift = [s+ds for s,ds in zip(shift0,dshift)]\n\n            theta = theta/180*np.pi\n            mat = np.array([\n                scale*cos(theta), -scale*sin(theta), shift[0],\n                scale*sin(theta),  scale*cos(theta), shift[1],\n            ]).reshape(2,3)\n\n        p_augment = np.concatenate([p_512.reshape(-1,2), np.ones((10,1))], axis=1)@mat.T\n        p_augment = p_augment.reshape(-1)\n        m_augment = cv2.warpAffine(m_512, mat, (512, 512))\n         \n        m_augment =  draw_shape(\n            p_augment,  #x,y,x,y,x,y ... format\n            image=m_augment, \n            is_shift_center=False, \n            is_line=j!=0, line_thickness=3,\n            is_circle=j!=0, circle_radius=20,\n            level_color=level_color\n        ) \n        shape_augment.append(p_augment)\n        image_augment.append(m_augment)\n\n    image_augment = np.hstack(image_augment) \n    plt.figure(figsize=(15,5))\n    plt.imshow(image_augment, cmap='gray')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-23T02:25:49.264493Z","iopub.execute_input":"2024-08-23T02:25:49.264917Z","iopub.status.idle":"2024-08-23T02:25:51.566889Z","shell.execute_reply.started":"2024-08-23T02:25:49.264874Z","shell.execute_reply":"2024-08-23T02:25:51.565658Z"},"trusted":true},"execution_count":null,"outputs":[]}]}