OCTA_SEG/depth_map_generation.py at main · xiaobisheng/OCTA_SEG · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
import os
import cv2
import numpy as np
import open3d as o3d

root_path = './Data/OCTA_6mm/Projection Maps'
def generate_depth_maps(data, i, threshold=10):
    data = data.transpose(2, 1, 0)

    # OCTA projection map
    img = np.average(data, axis=1) / 255
    cv2.imshow('img', img * 40)
    cv2.waitKey(200)

    # DCCM
    output_dccm = np.sum(data > threshold, axis=1)
    output_dccm = output_dccm / data.shape[1]
    cv2.imshow('DCCM', output_dccm)
    cv2.waitKey(2000)
    # save_path = root_path + '/OCTA_dccm_%d/%d.bmp' % (threshold, i)
    # cv2.imwrite(save_path, output_dccm)

    # DCIM
    data_ = data.transpose(0, 2, 1)
    indices = np.tile(np.arange(data_.shape[2]), (data_.shape[0], data_.shape[1], 1))  # generate the channel index matrix, (W, H, C)
    valid_indices = np.where(data_ > threshold, indices, 0)
    output_dcim = np.sum(valid_indices, axis=2)
    num = np.sum(data_ > threshold, axis=2) + 0.01
    output_dcim = output_dcim / num
    output_dcim = output_dcim / data.shape[1]
    cv2.imshow('DCIM', output_dcim)
    cv2.waitKey(20000)
    # save_path = root_path + '/OCTA_dcim_%d/%d.bmp' % (threshold, i)
    # cv2.imwrite(save_path, output_dcim)

def adaptativeDCIM_square(data, i):
    """
    This is the code to compute DCIM without parameters.
    """
    data = data.transpose(2, 1, 0)
    data = data.transpose(1, 0, 2)

    # OCTA projection map
    img = np.average(data, axis=0) / 255
    cv2.imshow('img', img * 10)
    cv2.waitKey(200)

    # Create a depth index array [0, 1, ..., C-1], shape (C, 1, 1)
    data = data.astype(np.int32)
    data2 = data ** 2
    depth_indices = np.arange(data2.shape[0], dtype=np.float32).reshape(-1, 1, 1)
    # Numerator: weighted sum of depth indices by intensity
    numerator = np.sum(data2 * depth_indices, axis=0)
    # Denominator: total intensity at each (i, j)
    denominator = np.sum(data2, axis=0)
    # Avoid division by zero -> result=0 when denominator=0
    dcim = np.divide(
        numerator,
        denominator,
        out=np.zeros_like(numerator, dtype=np.float32),
        where=denominator > 0
    )
    dcim = dcim / 128.0

    # dcim = dcim * 255
    # dcim = np.array(dcim, dtype=np.uint8)

    cv2.imshow('dcim', dcim)
    cv2.waitKey(100)

    save_path = root_path + '/OCTA_DCIM2/'
    if not os.path.exists(save_path):
        os.makedirs(save_path)
    save_file = save_path + '%d.bmp' % (i)
    cv2.imwrite(save_file, dcim)

    return dcim

def visualize_3d(data):
    data = data[:, :, :] * 3
    W, H, C = data.shape
    threshold = 25  # setting a threshold for ignoring low intensity voxel
    coords = np.argwhere(data > threshold)  # only choose voxel > threshold
    colors = data[coords[:, 0], coords[:, 1], coords[:, 2]]
    colors = np.tile(colors[:, None], (1, 3)) / 255.0  # normalization
    # create Open3D point cloud object
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(coords)  # setting coordinates
    pcd.colors = o3d.utility.Vector3dVector(colors)  # setting color (gray)
    # visualization
    o3d.visualization.draw_geometries([pcd])

if __name__ == '__main__':
    for i in range(10001, 10300, 1):
        octa_data = np.load(root_path + '/OCTA_npy/%d.npy' % (i))
        visualize_3d(octa_data)
        generate_depth_maps(octa_data, i)