demo.py 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116
  1. import os
  2. import numpy as np
  3. import cv2
  4. import matplotlib.pyplot as plt
  5. from pyfeats import *
  6. #%% Pre-processing
  7. def padImage(im, pad=2, value=0):
  8. TDLU=[1, 1, 1, 1] #top, down, left, right pad
  9. out = im.copy()
  10. for _ in range(pad):
  11. out = cv2.copyMakeBorder(out, TDLU[0], TDLU[1], TDLU[2], TDLU[3],\
  12. cv2.BORDER_CONSTANT, None, value)
  13. return out
  14. # Load image
  15. path = os.getcwd() + '/data/'
  16. image_name ='ultrasound.bmp'
  17. image = cv2.imread(path + image_name, cv2.IMREAD_GRAYSCALE)
  18. plt.imshow(image, cmap='gray')
  19. plt.title('Initial Image')
  20. plt.show()
  21. # Load mask
  22. points_name = 'points.out'
  23. points = np.loadtxt(path + points_name, delimiter=',')
  24. points = np.array(points, np.int32)
  25. mask = cv2.fillPoly(np.zeros(image.shape,np.double), [points.reshape((-1,1,2))], color=1).astype('i')
  26. perimeter = cv2.polylines(np.zeros(image.shape), [points.reshape((-1,1,2))],isClosed=True,color=1,thickness=1).astype('i')
  27. plt.imshow(mask, cmap='gray')
  28. plt.title('Initial Mask')
  29. plt.show()
  30. # Fill image with zeros outside ROI
  31. image = np.multiply(image.astype(np.double), mask).astype(np.uint8)
  32. max_point_x = max(points[:,0])
  33. max_point_y = max(points[:,1])
  34. min_point_x = min(points[:,0])
  35. min_point_y = min(points[:,1])
  36. a,b,c,d = min_point_y, (max_point_y+1), min_point_x, (max_point_x+1)
  37. image, mask, perimeter = image[a:b,c:d], mask[a:b,c:d], perimeter[a:b,c:d]
  38. image, mask, perimeter = padImage(image, 2, 0), padImage(mask, 2, 0), padImage(perimeter, 2, 0) # ALWAYS PAD IMAGE SO AS SHAPE PARAMETERS WORK
  39. plt.imshow(image, cmap='gray')
  40. plt.title('Image cropped to ROI')
  41. plt.show()
  42. plt.imshow(mask, cmap='gray')
  43. plt.title('Mask of image cropped to ROI')
  44. plt.show()
  45. plt.imshow(perimeter, cmap='gray')
  46. plt.title('Perimeter of image cropped to ROI')
  47. plt.show()
  48. #%% A1. Texture features
  49. features = {}
  50. features['A_FOS'] = fos(image, mask)
  51. features['A_GLCM'] = glcm_features(image, ignore_zeros=True)
  52. features['A_GLDS'] = glds_features(image, mask, Dx=[0,1,1,1], Dy=[1,1,0,-1])
  53. features['A_NGTDM'] = ngtdm_features(image, mask, d=1)
  54. features['A_SFM'] = sfm_features(image, mask, Lr=4, Lc=4)
  55. features['A_LTE'] = lte_measures(image, mask, l=7)
  56. features['A_FDTA'] = fdta(image, mask, s=3)
  57. features['A_GLRLM'] = glrlm_features(image, mask, Ng=256)
  58. features['A_FPS'] = fps(image, mask)
  59. features['A_Shape_Parameters'] = shape_parameters(image, mask, perimeter, pixels_per_mm2=1)
  60. features['A_HOS'] = hos_features(image, th=[135,140])
  61. features['A_LBP'] = lbp_features(image, image, P=[8,16,24], R=[1,2,3])
  62. features['A_GLSZM'] = glszm_features(image, mask)
  63. plot_sinogram(image, image_name)
  64. #%% B. Morphological features
  65. features['B_Morphological_Grayscale_pdf'], features['B_Morphological_Grayscale_cdf'] = grayscale_morphology_features(image, N=30)
  66. features['B_Morphological_Binary_L_pdf'], features['B_Morphological_Binary_M_pdf'], features['B_Morphological_Binary_H_pdf'], features['B_Morphological_Binary_L_cdf'], \
  67. features['B_Morphological_Binary_M_cdf'], features['B_Morphological_Binary_H_cdf'] = multilevel_binary_morphology_features(image, mask, N=30, thresholds=[25,50])
  68. plot_pdf_cdf(features['B_Morphological_Grayscale_pdf'], features['B_Morphological_Grayscale_cdf'], image_name)
  69. plot_pdfs_cdfs(features['B_Morphological_Binary_L_pdf'], features['B_Morphological_Binary_M_pdf'], features['B_Morphological_Binary_H_pdf'], features['B_Morphological_Binary_L_cdf'], features['B_Morphological_Binary_M_cdf'], features['B_Morphological_Binary_H_cdf'])
  70. #%% C. Histogram Based features
  71. features['C_Histogram'] = histogram(image, mask, bins=32)
  72. features['C_MultiregionHistogram'] = multiregion_histogram(image, mask, bins=32, num_eros=3, square_size=3)
  73. features['C_Correlogram'] = correlogram(image, mask, bins_digitize=32, bins_hist=32, flatten=True)
  74. plot_histogram(image, mask, bins=32, name=image_name)
  75. plot_correlogram(image, mask, bins_digitize=32, bins_hist=32, name=image_name)
  76. #%% D. Multi-Scale features
  77. features['D_DWT'] = dwt_features(image, mask, wavelet='bior3.3', levels=3)
  78. features['D_SWT'] = swt_features(image, mask, wavelet='bior3.3', levels=3)
  79. features['D_WP'] = wp_features(image, mask, wavelet='coif1', maxlevel=3)
  80. features['D_GT'] = gt_features(image, mask)
  81. features['D_AMFM'] = amfm_features(image)
  82. #%% E. Other
  83. features['E_HOG'] = hog_features(image, ppc=8, cpb=3)
  84. features['E_HuMoments'] = hu_moments(image)
  85. features['E_TAS'] = tas_features(image)
  86. features['E_ZernikesMoments'] = zernikes_moments(image, radius=9)
  87. #%% Print
  88. for x, y in features.items():
  89. if x.startswith("B") | x.startswith("C") | x.startswith("E"):
  90. continue
  91. print('-' * 50)
  92. print(x)
  93. print('-' * 50)
  94. if len(y[1]) == 1:
  95. print(y[1][0], '=', y[0])
  96. else:
  97. for i in range(len(y[1])):
  98. print(y[1][i], '=', y[0][i])