cva_cd.py
1 #!/usr/bin/env python3 2 # -*- coding: utf-8 -*- 3 """ 4 Created on Thu Mar 10 09:40:56 2022 5 6 @author: aleoikon 7 """ 8 9 import sys 10 sys.path.append('/home/dvalsamis/Documents/projects/Change_detection_SSL_Siamese') 11 import time 12 13 14 from skimage.filters import threshold_otsu, threshold_triangle 15 from tensorflow import keras 16 from architectures.branch import branches_nopool, branch_cva, branch_cva_aspp,two_branch_cva_with_aspp,two_branch_cva_with_aspp_fmaps 17 #from tests import change_detection_noup, change_detection_noup_1x1convs 18 from architectures.similarity_detection import pretext_task_one_nopool 19 from tensorflow.keras.optimizers import Adam 20 from tensorflow.keras.utils import plot_model 21 import pandas as pd 22 import tensorflow as tf 23 import numpy as np 24 import os 25 import random 26 from numpy import expand_dims 27 from skimage.morphology import remove_small_objects 28 import matplotlib 29 matplotlib.use('TkAgg') 30 import matplotlib.pyplot as plt 31 from utils.layer_select import feature_selector_cva, feature_selector_cva_aspp,two_feature_selector_cva_aspp 32 from utils.log_params import log_params_cva 33 from architectures.conv_classifier import conv_classifier, conv_classifier_two,conv_classifier_two_with_aspp 34 from keras.models import Sequential 35 from skimage.filters import gaussian 36 from utils.my_metrics import recall, accuracy, specificity, precision, f_measure, get_roc 37 import uuid 38 from datetime import datetime 39 os.environ["CUDA_VISIBLE_DEVICES"]="1" 40 41 42 #Euclidean Distance 43 def calculate_distancemap(f1, f2): 44 """ 45 calcualtes pixelwise euclidean distance between images with multiple imput channels 46 47 Parameters 48 ---------- 49 f1 : np.ndarray of shape (N,M,D) 50 image 1 with the channels in the third dimension 51 f2 : np.ndarray of shape (N,M,D) 52 image 2 with the channels in the third dimension 53 54 Returns 55 ------- 56 np.ndarray of shape(N,M) 57 pixelwise euclidean distance between image 1 and image 2 58 59 """ 60 dist_per_fmap= [(f2[i,:,:]-f1[i,:,:])**2 for i in range(f1.shape[0])] 61 62 return np.sqrt(sum(dist_per_fmap)) 63 64 def create_rgb_onera(x,channel): 65 if channel == 'red': 66 r = x[:,:,2] 67 r = np.expand_dims(r, axis=2) 68 return r 69 if channel == 'green': 70 g = x[:,:,1] 71 g = np.expand_dims(g, axis=2) 72 return g 73 if channel == 'blue': 74 b = x[:,:,0] 75 b = np.expand_dims(b, axis=2) 76 return b 77 if channel == 'rgb': 78 r = x[:,:,2] 79 g = x[:,:,1] 80 b = x[:,:,0] 81 rgb = np.dstack((r,g,b)) 82 return(rgb) 83 if channel == 'rgbvnir': 84 r = x[:,:,2] 85 g = x[:,:,1] 86 b = x[:,:,0] 87 vnir = x[:,:,3] 88 rgbvnir = np.stack((r,g,b,vnir),axis=2).astype('float') 89 return(rgbvnir) 90 else: 91 return x 92 print("NOT CORRECT CHANNELS") 93 94 def normalize(x): 95 img =((x - x.mean()) / x.std()) 96 return img 97 98 def scaleMinMax(x): 99 return ((x - np.nanpercentile(x,2)) / (np.nanpercentile(x,98) - np.nanpercentile(x,2))) 100 101 102 def create_rgb(x, channels): 103 if channels == 'red': 104 r = x[:,:,0] 105 r = scaleMinMax(r) 106 return r 107 if channels == 'green': 108 g = x[:,:,0] 109 g = scaleMinMax(g) 110 return g 111 if channels == 'blue': 112 b = x[:,:,0] 113 b = scaleMinMax(b) 114 return b 115 if channels == 'rgb': 116 r = x[:,:,2] 117 r = scaleMinMax(r) 118 g = x[:,:,1] 119 g = scaleMinMax(g) 120 b = x[:,:,0] 121 b = scaleMinMax(b) 122 rgb = np.dstack((r,g,b)) 123 return(rgb) 124 125 126 127 #---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 128 129 depth = 2 130 dropout = 0.1 131 decay = 0.0001 132 NORM = True 133 ImageSize = 96 134 n_ch = 3 135 channel = 'rgb' 136 137 # Models ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 138 source_model = conv_classifier_two_with_aspp(depth, dropout, decay, ImageSize, ImageSize, n_ch) 139 mtype = 'conv_classifier_two' 140 141 cd_model_path = '/home/dvalsamis/Documents/projects/Change_detection_SSL_Siamese/saved_models/CD_Simple_Onera_0a8d.h5' 142 cd_model_name = 'CD_Simple_CBMI_0a8d' 143 model_id = cd_model_name.split('_')[-1].split('.')[0] 144 145 146 source_model.load_weights(cd_model_path) 147 148 149 branch_model = two_branch_cva_with_aspp(dropout, decay, depth, ImageSize, ImageSize, n_ch) 150 branch_model = two_feature_selector_cva_aspp(depth, source_model, branch_model) 151 #plot_model(branch_model, to_file='/home/dvalsamis/Documents/projects/Change_detection_SSL_Siamese/graphs_cva/'+model_id+'_model_plot.png', show_shapes=True, show_layer_names=True) 152 153 154 # Data ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 155 156 onera_train_target = '/home/dvalsamis/Documents/data/CBMI/CBMI_0.3/CBMI_0.3/NPY_dataset/aug_train_data/' 157 onera_test_target = '/home/dvalsamis/Documents/data/CBMI/CBMI_0.3/CBMI_0.3/NPY_dataset/aug_test_data/' 158 159 160 onera_test_df = pd.read_csv(onera_test_target + "dataset_test.csv") 161 162 dsize = len(onera_test_df) # Total size of the dataset 163 trsize = len(onera_train_target) # Size of the training set 164 tesize = len(onera_test_target) # Size of the testing set 165 166 167 X1 = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize,n_ch)) 168 X2 = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize,n_ch)) 169 170 nonNorm1 = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize,n_ch)) 171 nonNorm2 = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize,n_ch)) 172 173 y = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize)) 174 y_pred = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize)) 175 y_pred_otsu = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize)) 176 y_pred_triangle = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize)) 177 178 179 180 for i in range(len(onera_test_df)): 181 img1 = np.load(onera_test_target+onera_test_df['pair1'][i]) 182 img2 = np.load(onera_test_target+onera_test_df['pair2'][i]) 183 img1 = create_rgb_onera(img1, channel) 184 img2 = create_rgb_onera(img2, channel) 185 186 nonNorm1[i] = img1 187 nonNorm2[i] = img2 188 189 if NORM: 190 X1[i] = normalize(img1) 191 X2[i] = normalize(img2) 192 else: 193 X1[i] = img1 194 X2[i] = img2 195 y[i] = np.load(onera_test_target+onera_test_df['change_mask'][i]) 196 197 198 branch_model.summary() 199 200 201 # Setting Predictions---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 202 203 # INITIAL code 204 205 # FILTER = False 206 207 # for index in range(len(X1)): 208 # img1 = expand_dims(X1[index], axis=0) 209 # img2 = expand_dims(X2[index], axis=0) 210 211 # feature_maps_left = branch_model.predict(img1) 212 # feature_maps_right = branch_model.predict(img2) 213 214 # left = np.ndarray(shape=(32,ImageSize,ImageSize)) 215 # right = np.ndarray(shape=(32,ImageSize,ImageSize)) 216 # for i in range(left.shape[0]): 217 # left[i] = feature_maps_left[0,:,:,i] 218 # right[i] = feature_maps_right[0,:,:,i] 219 220 # distmap = calculate_distancemap(left, right) 221 # if FILTER == True: 222 # distmap = gaussian(distmap, sigma=1, mode='constant', cval=0.0) 223 224 # #otsu 225 # binary_otsu = distmap > threshold_otsu(distmap) 226 # binary_otsu = remove_small_objects(binary_otsu,min_size=100) 227 228 # #triange 229 # binary_triangle = distmap > threshold_triangle(distmap) 230 # binary_triangle = remove_small_objects(binary_triangle,min_size=100) 231 232 # y_pred_otsu[index]=binary_otsu 233 # y_pred_triangle[index]=binary_triangle 234 235 # y_pred_conv = source_model.predict([X1,X2]) 236 # y_pred_conv = np.argmax(y_pred_conv, axis=3) 237 238 239 #New method with ASPP WORKING one Branch -------------------------------------------------------------------------------------------------------------------------------------------- 240 # FILTER = False 241 242 # for index in range(len(X1)): 243 # img1 = np.expand_dims(X1[index], axis=0) 244 # img2 = np.expand_dims(X2[index], axis=0) 245 246 # feature_maps_left = branch_model.predict(img1) 247 # feature_maps_right = branch_model.predict(img2) 248 249 # # Reduce the dimensionality from (96, 96, 32) to (96, 96) by averaging across channels 250 # feature_maps_left = np.mean(feature_maps_left, axis=-1) 251 # feature_maps_right = np.mean(feature_maps_right, axis=-1) 252 253 # # Now, you have a single-channel (96, 96) distance map 254 # distmap = np.abs(feature_maps_left - feature_maps_right) # absolute difference between averaged maps 255 256 # if FILTER: 257 # distmap = gaussian(distmap, sigma=1, mode='constant', cval=0.0) 258 259 # # Apply Otsu and Triangle thresholding 260 # binary_otsu = distmap > threshold_otsu(distmap) 261 # binary_otsu = remove_small_objects(binary_otsu, min_size=100) 262 263 # binary_triangle = distmap > threshold_triangle(distmap) 264 # binary_triangle = remove_small_objects(binary_triangle, min_size=100) 265 266 # y_pred_otsu[index] = binary_otsu 267 # y_pred_triangle[index] = binary_triangle 268 269 # y_pred_conv = source_model.predict([X1,X2]) 270 # y_pred_conv = np.argmax(y_pred_conv, axis=3) 271 272 #---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 273 274 # Final two Branch Model 275 276 for index in range(len(X1)): 277 img1 = np.expand_dims(X1[index], axis=0) 278 img2 = np.expand_dims(X2[index], axis=0) 279 280 # Get the output from the model, which is (1, 96, 96, 32) 281 combined_features = branch_model.predict([img1, img2]) 282 283 # Calculate the mean across the channels to get a (96, 96) map 284 mean_feature_map = np.mean(combined_features, axis=-1)[0] 285 286 287 # Otsu's thresholding 288 binary_otsu = mean_feature_map > threshold_otsu(mean_feature_map) 289 binary_otsu = remove_small_objects(binary_otsu, min_size=100) 290 291 # Triangle thresholding 292 binary_triangle = mean_feature_map > threshold_triangle(mean_feature_map) 293 binary_triangle = remove_small_objects(binary_triangle, min_size=100) 294 295 y_pred_otsu[index] = binary_otsu 296 y_pred_triangle[index] = binary_triangle 297 298 y_pred_conv = source_model.predict([X1,X2]) 299 y_pred_conv = np.argmax(y_pred_conv, axis=3) 300 301 print('Recall',recall(y,y_pred_conv)) 302 print('Specificity',specificity(y,y_pred_conv)) 303 print('Precision',precision(y,y_pred_conv)) 304 print('F1',f_measure(y,y_pred_conv)) 305 print('Accuracy',accuracy(y,y_pred_conv)) 306 307 308 309 310 311 #---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 312 313 pos = random.randint(0, len(y)) 314 print(pos) 315 print(len(y)) 316 #pos = 157 317 #pos = 712 318 319 fig, ax = plt.subplots(3, 3, figsize=(10,10),constrained_layout=True) 320 321 ax[0,0].imshow(create_rgb(X1[pos], channel)) 322 ax[0,0].set_title('Left Image', fontsize=20) 323 ax[0,0].axis('off') 324 325 ax[0,1].imshow(y[pos], cmap='gray') 326 ax[0,1].set_title('Ground Truth', fontsize=20) 327 ax[0,1].axis('off') 328 329 ax[0,2].imshow(create_rgb(X2[pos], channel)) 330 ax[0,2].set_title('Right Image', fontsize=20) 331 ax[0,2].axis('off') 332 333 ax[1,0].imshow(create_rgb(X1[pos], channel)) 334 ax[1,0].set_title('Left Image', fontsize=20) 335 ax[1,0].axis('off') 336 337 ax[1,1].imshow(y_pred_otsu[pos], cmap='gray') 338 ax[1,1].set_title('CVA(with Otsu)', fontsize=20) 339 ax[1,1].axis('off') 340 341 ax[1,2].imshow(create_rgb(X2[pos], channel)) 342 ax[1,2].set_title('Right Image', fontsize=20) 343 ax[1,2].axis('off') 344 345 ax[2,0].imshow(create_rgb(X1[pos], channel)) 346 ax[2,0].set_title('Left Image', fontsize=20) 347 ax[2,0].axis('off') 348 349 ax[2,1].imshow(y_pred_triangle[pos], cmap='gray') 350 ax[2,1].set_title('CVA(with Triagle)', fontsize=20) 351 ax[2,1].axis('off') 352 353 ax[2,2].imshow(create_rgb(X2[pos], channel)) 354 ax[2,2].set_title('Right Image', fontsize=20) 355 ax[2,2].axis('off') 356 357 358 359 recall(y,y_pred_otsu) 360 get_roc(y,y_pred_otsu) 361 362 363 log_params_cva('CBMI_Test', model_id, cd_model_name, mtype, depth, dropout, decay, ImageSize, n_ch, channel, NORM) 364 365 #####Metrics##### 366 data_dict = {'Pretext':'Task 1', 367 'Model ID': model_id, 368 'Downstream':'CVA+Otsu(min size = 100)', 369 'Sensitivity/Recall':recall(y,y_pred_otsu), 370 'Specificity':specificity(y,y_pred_otsu), 371 'Precision':precision(y,y_pred_otsu), 372 'F1':f_measure(y,y_pred_otsu), 373 'Accuracy':accuracy(y,y_pred_otsu), 374 'Set':'CBMI Test', 375 'ImageSize':ImageSize, 376 'Norm':NORM, 377 'Pretext Model':'-', 378 'CD model':cd_model_name} 379 380 # Make data frame of above data 381 df = pd.DataFrame(data_dict, index=[0]) 382 383 results_path = '/home/dvalsamis/Documents/projects/Change_detection_SSL_Siamese/logs/cva_results_1.csv' 384 385 # append data frame to CSV file 386 #df.to_csv(results_path, mode='a', index=False, header=False) 387 388 # Append data to the CSV file 389 with open(results_path, 'a') as f: 390 # If the file is empty, write the headers 391 if f.tell() == 0: 392 pd.DataFrame([data_dict.keys()]).to_csv(f, header=False, index=False) 393 # Append data 394 pd.DataFrame([data_dict.values()]).to_csv(f, header=False, index=False) 395 396 # print message 397 print("Data logged successfully.") 398 399 data_dict = {'Pretext':'Task 1', 400 'Model ID': model_id, 401 'Downstream':'CVA+Triangle(min size = 100)', 402 'Sensitivity/Recall':recall(y,y_pred_triangle), 403 'Specificity':specificity(y,y_pred_triangle), 404 'Precision':precision(y,y_pred_triangle), 405 'F1':f_measure(y,y_pred_triangle), 406 'Accuracy':accuracy(y,y_pred_triangle), 407 'Set':'CBMI Test', 408 'ImageSize':ImageSize, 409 'Norm':NORM, 410 'Pretext Model':'-', 411 'CD model':cd_model_name} 412 413 # Make data frame of above data 414 df = pd.DataFrame(data_dict, index=[0]) 415 # append data frame to CSV file 416 417 # Append data to the CSV file 418 with open(results_path, 'a') as f: 419 # If the file is empty, write the headers 420 if f.tell() == 0: 421 pd.DataFrame([data_dict.keys()]).to_csv(f, header=False, index=False) 422 # Append data 423 pd.DataFrame([data_dict.values()]).to_csv(f, header=False, index=False) 424 # print message 425 print("Data logged successfully.") 426 427 428 cd_results_df = pd.read_csv(results_path) 429 430 ################fusion#################### 431 from architectures.fusion_maria import fusion 432 433 y_pred_fusion = np.ndarray(shape=(len(onera_test_df),ImageSize,ImageSize)) 434 435 for param in range(0, 13): 436 for cm_pos in range(len(y_pred_fusion)): 437 fused_mask = fusion(y_pred_triangle[cm_pos],y_pred_otsu[cm_pos],y_pred_conv[cm_pos], param) 438 y_pred_fusion[cm_pos] = fused_mask 439 440 data_dict = {'Pretext':'Task 1', 441 'Model ID': model_id, 442 'Downstream':'Fusion', 443 'Sensitivity/Recall':recall(y,y_pred_fusion), 444 'Specificity':specificity(y,y_pred_fusion), 445 'Precision':precision(y,y_pred_fusion), 446 'F1':f_measure(y,y_pred_fusion), 447 'Accuracy':accuracy(y,y_pred_fusion), 448 'Set':'CBMI Test', 449 'ImageSize':ImageSize, 450 'Norm':"param="+str(param), 451 'Pretext Model':'-', 452 'CD model':cd_model_name} 453 454 # Make data frame of above data 455 df = pd.DataFrame(data_dict, index=[0]) 456 # Append data to the CSV file 457 with open(results_path, 'a') as f: 458 # If the file is empty, write the headers 459 if f.tell() == 0: 460 pd.DataFrame([data_dict.keys()]).to_csv(f, header=False, index=False) 461 # Append data 462 pd.DataFrame([data_dict.values()]).to_csv(f, header=False, index=False) 463 # print message 464 print("Data logged successfully.", param) 465 466 #cd_results_df = pd.read_csv('cd_results.csv') 467 468 #plot fusion 469 cm_pos = 0 470 for cm_pos in range(len(y_pred_fusion)): 471 fused_mask = fusion(y_pred_triangle[cm_pos],y_pred_otsu[cm_pos],y_pred_conv[cm_pos], 10) 472 y_pred_fusion[cm_pos] = fused_mask 473 474 pos = random.randint(0, len(y)) 475 print(pos) 476 print(len(y)) 477 #pos = 157 478 #pos = 712 479 #pos = 577 480 fig, ax = plt.subplots(5, 1, figsize=(10,40),constrained_layout=True) 481 482 ax[0].imshow(y[pos], cmap='gray') 483 ax[0].set_title('Ground Truth', fontsize=30) 484 ax[0].axis('off') 485 486 ax[1].imshow(y_pred_conv[pos], cmap='gray') 487 ax[1].set_title('Conv classifier', fontsize=30) 488 ax[1].axis('off') 489 490 ax[2].imshow(y_pred_otsu[pos], cmap='gray') 491 ax[2].set_title('CVA(with Otsu)', fontsize=30) 492 ax[2].axis('off') 493 494 ax[3].imshow(y_pred_triangle[pos], cmap='gray') 495 ax[3].set_title('CVA(with Triagle)', fontsize=30) 496 ax[3].axis('off') 497 498 ax[4].imshow(y_pred_fusion[pos], cmap='gray') 499 ax[4].set_title('Fusion', fontsize=30) 500 ax[4].axis('off') 501 502 fig, ax = plt.subplots(1, 2, figsize=(10,40),constrained_layout=True) 503 ax[0].imshow(create_rgb(X1[pos], 'rgb')) 504 ax[0].set_title('Left Image', fontsize=20) 505 ax[0].axis('off') 506 507 508 509 ax[1].imshow(create_rgb(X2[pos], 'rgb')) 510 ax[1].set_title('Right Image', fontsize=20) 511 ax[1].axis('off') 512 513 514 515 516 print("End of first excecution") 517