CTransformer/tem.py at main · chiellini/CTransformer · 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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
"""Demo for calculating cubical complexes.

This example demonstrates how to perform topological operations on
a structured array, such as a grey-scale image.
"""
import nibabel
from scipy import ndimage
from skimage.transform import resize
from matplotlib import pyplot as plt
from gtda.plotting import plot_heatmap
import numpy as np
import matplotlib.pyplot as plt
from torchmetrics.functional import dice,cosine_similarity,mean_squared_log_error,mean_squared_error

from torch_topological.nn import CubicalComplex
from torch_topological.nn import WassersteinDistance
from gudhi.bottleneck import bottleneck_distance
import torch

from tqdm import tqdm
from gtda.images import RadialFiltration,DilationFiltration
from gtda.homology import CubicalPersistence
from gtda.diagrams import Scaler

import torch

if __name__ == '__main__':

    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
    print('Device', device)

    np.random.seed(23)

    volume_seg_nuc_path=r'/home/cimda/zelinli/CellAtlas/DataSource/TrainingandEvaluation/evaluation/SegNuc/200113plc1p2_185_segNuc.nii.gz'
    seg_nuc_volume=nibabel.load(volume_seg_nuc_path).get_fdata()
    seg_nuc_volume_num=len(np.unique(seg_nuc_volume))-1
    prior_PH=np.zeros((seg_nuc_volume_num,2))
    prior_PH[:,1]=2.0

    volume_path=r'/home/cimda/zelinli/CellAtlas/OutputData/OriginalSwinUNETR_validating/SegMemb/200113plc1p2_185_segMemb.nii.gz'
    initvolume=nibabel.load(volume_path).get_fdata()
    volume = resize(initvolume, (256, 384, 128), mode='constant', cval=0, order=0,anti_aliasing=False)
    volume = volume-np.mean(volume)
    volume = volume/np.max(volume)
    # volume=torch.as_tensor(volume,dtype=torch.float,device=device)
    noisess = np.random.normal(0.0, 0.001, volume.shape)
    # X = torch.nn.Parameter(volume, requires_grad=True).to(device)
    #
    # source_gpu = X.clone()
    # source = source_gpu.cpu().detach().numpy()
    # # source_binarized_volume=np.zeros(source.shape)
    # # source_binarized_volume[source>200]=1
    # # source_binarized_volume=torch.as_tensor(source_binarized_volume,dtype=torch.float,device=device)
    # optimizer = torch.optim.Adam([X], lr=1e-2)

    X = torch.autograd.Variable(torch.tensor(volume+noisess).type(torch.float), requires_grad=True)
    source_gpu = torch.tensor(volume, dtype=torch.float, requires_grad=False)
    # y_t = torch.tensor(y, dtype=torch.float, requires_grad=False)
    optimizer = torch.optim.Adam([X], lr=1)


    n_iter = 200
    progress = tqdm(range(n_iter))

    for i in progress:
        volume_cpu_tem=X.cpu().detach().numpy()
        binarized_volume=np.zeros(volume_cpu_tem.shape)
        # random_threshold=np.random.randint(low=204,high=229)
        threshold_std=np.max(volume_cpu_tem)-0.2*(np.max(volume_cpu_tem)-np.min(volume_cpu_tem))
        print('\n',threshold_std)
        binarized_volume[volume_cpu_tem > threshold_std] = 1
        dilation_filtration = DilationFiltration(n_iterations=1)
        volume_filtration = dilation_filtration.fit_transform(binarized_volume[np.newaxis, ...])
        # print('\n',torch.unique(source_gpu,return_counts=True))
        # print('\n',np.unique(binarized_volume,return_counts=True))

        fig, ax = plt.subplots(ncols=3)
        # volume_cpu_tem_showing
        ax[0].imshow(volume_cpu_tem[128,:,:])
        ax[0].set_title('X')

        ax[1].imshow(binarized_volume[128,:,:])
        ax[1].set_title('Binary')

        ax[2].imshow(volume_filtration[0,128,:,:])
        ax[2].set_title('Filtration Binary')
        plt.show()

        mse_loss = (mean_squared_log_error(X, source_gpu)+1e-3).to(device)
        print('MSE LOSSSS', mse_loss)

        cubical_persistence = CubicalPersistence((0,1,2))
        volume_cubical = cubical_persistence.fit_transform(volume_filtration)
        scaler = Scaler()
        volume_cubical_scaled = scaler.fit_transform(volume_cubical)[0]
        # volume_cubical_scaled=volume_cubical_scaled[0]
        # persistence_information = cubical_complex(X)
        # persistence_information = persistence_information[0]
        volume_h_arr = volume_cubical_scaled[:, 2].astype(int)
        # print(volume_cubical_scaled)
        volume_h0_num, volume_h1_num,volume_h2_num = np.unique(volume_h_arr, return_counts=True)[1]
        pred_h2_pd=volume_cubical_scaled[volume_h0_num+volume_h1_num:,:2]
        pred_h2_pd=pred_h2_pd[(pred_h2_pd[:,1]-pred_h2_pd[:,0])>0.5,:]
        h2_bottleneck_dis=bottleneck_distance(pred_h2_pd,prior_PH)
        print('H2 bottle neck distance\n',pred_h2_pd.shape,prior_PH.shape,h2_bottleneck_dis)


        print(abs(volume_h2_num-seg_nuc_volume_num),2-np.mean(np.abs(pred_h2_pd[:,0]-pred_h2_pd[:,1])))

        loss =mse_loss+torch.as_tensor(0.01*np.array(abs(volume_h2_num-seg_nuc_volume_num)+abs(2-np.mean(np.abs(pred_h2_pd[:,0]-pred_h2_pd[:,1])))),dtype=torch.float,
                              device=device)
        # loss.requires_grad=True

        print('LOSS \n',loss)
        # print('X \n',X)

        optimizer.zero_grad()
        loss.backward()
        # print('Gradient of optimizer\n',[x.grad for x in optimizer.param_groups[0]['params']])
        optimizer.step()

        # progress.set_postfix(loss=loss.item())

        # result = X.cpu().detach().numpy()
        #
        # fig, ax = plt.subplots(ncols=3)
        #
        # ax[0].imshow(source[128,:,:])
        # ax[0].set_title('Source')
        #
        # ax[1].imshow(volume_filtration[0,128,:,:])
        # ax[1].set_title('FIltration Source')
        #
        # ax[2].imshow(result[128,:,:])
        # ax[2].set_title('Result')
        #
        # plt.show()