Blame view

python/muve-align.py 4.95 KB
a296b833   David Mayerich   added the muve al...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
  # -*- coding: utf-8 -*-
  """
  Created on Fri May  4 17:15:07 2018
  
  @author: david
  """
  
  import imagestack
  import numpy
  import matplotlib.pyplot as plt
  import cv2
  import progressbar
  import scipy.ndimage
  import sys
2119be83   David Mayerich   updated libraries...
15
  import os
a296b833   David Mayerich   added the muve al...
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
  
  
  def press(event):
      global i
      global ax
      global fig
      if event.key == 'z':
          i = i - 1
          if i < 0:
              i = 0
      if event.key == 'x':
          i = i + 1
          if i == S.shape[0]:
              i = S.shape[0]
      if event.key == "up":
          for n in range(i, S.shape[0]):
              warps[n][1, 2] = warps[n][1, 2] + 1
      if event.key == "down":
          for n in range(i, S.shape[0]):
              warps[n][1, 2] = warps[n][1, 2] - 1
      if event.key == "left":
          for n in range(i, S.shape[0]):
              warps[n][0, 2] = warps[n][0, 2] + 1
      if event.key == "right":
          for n in range(i, S.shape[0]):
              warps[n][0, 2] = warps[n][0, 2] - 1
      aligned = cv2.warpAffine(S[i, :, :, :], warps[i], (S.shape[2], S.shape[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP);
      ax.cla()
      ax.imshow(aligned.astype(numpy.uint8))
      ax.set_title('Slice ' + str(i))
      fig.canvas.draw() 
      sys.stdout.flush
  
  #apply the alignment adjustments to I based on the list of warp matrices in t
  def apply_alignment(I, t):
      A = numpy.zeros(I.shape, dtype=numpy.uint8)
      for i in range(0, I.shape[0]):
          A[i, :, :, :] = cv2.warpAffine(I[i, :, :, :], t[i], (I.shape[2], I.shape[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP);
      
      return A
  
  #get the transform that aligns image B to image A
  def align(A, B, max_power=5):
      # Define the motion model
      warp_mode = cv2.MOTION_TRANSLATION
               
      # Specify the number of iterations.
      number_of_iterations = 5000;
       
      # Specify the threshold of the increment
      # in the correlation coefficient between two iterations
      termination_eps = 1e-10;
       
      # Define termination criteria
      criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations,  termination_eps)
      
      #attempt to fit the image, increasing the blur as necessary
      i = 0
      while i <= max_power:    
          #blur the image used for fitting
          im1_blur = scipy.ndimage.filters.gaussian_filter(A, (2 ** i, 2 ** i), mode='constant')
          im2_blur = scipy.ndimage.filters.gaussian_filter(B, (2 ** i, 2 ** i), mode='constant')
           
          # Run the ECC algorithm. The results are stored in warp_matrix.
          # Define 2x3 matrix and initialize the matrix to identity
          warp_matrix = numpy.eye(2, 3, dtype=numpy.float32)
          try:
              (cc, warp_matrix) = cv2.findTransformECC (im1_blur,im2_blur,warp_matrix, warp_mode, criteria)
          except:
              #print("Error aligning at p = " + str(i))
              i = i + 1
          else:
              #print("Successful alignment at p = " + str(i))
              break
      #enforce the fact that the x-axis is already aligned
      #warp_matrix[0, 2] = 0
      return warp_matrix
      
      #if i > 0:
      #    (cc, warp_matrix) = cv2.findTransformECC (A,B,warp_matrix, warp_mode, criteria)
          #warp_matrix[0, 2] = 0
      #    return warp_matrix
  
59b0e5f5   David Mayerich   updated envi file...
99
100
  fmask = "D:/Dropbox/todo/lab-presentations/jack/TinkParaffinLung0.005S/*.png"
  out_dir = "D:/Dropbox/todo/lab-presentations/jack/TinkParaffinLung0.005S/aligned"
a296b833   David Mayerich   added the muve al...
101
102
103
104
105
106
  
  if not os.path.isdir(out_dir):
      os.mkdir(out_dir)
  
  #read the image stack for alignment
  print("Loading image stack...")
2119be83   David Mayerich   updated libraries...
107
108
  S = imagestack.load(fmask, dtype=numpy.uint8)
  
a296b833   David Mayerich   added the muve al...
109
  #convert to grayscale
2119be83   David Mayerich   updated libraries...
110
  G = imagestack.rgb2gray(S.astype(numpy.float32))
a296b833   David Mayerich   added the muve al...
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
  
  # Define the motion model
  warp_mode = cv2.MOTION_TRANSLATION
  
  # Define the motion model
  warp_mode = cv2.MOTION_TRANSLATION
   
  # Specify the number of iterations.
  number_of_iterations = 5000
   
  # Specify the threshold of the increment
  # in the correlation coefficient between two iterations
  termination_eps = 1e-10
   
  # Define termination criteria
  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations,  termination_eps)
  
  print("\n\nCalculating alignment transformations...")
  bar = progressbar.ProgressBar(max_value=G.shape[0])
  #for each image in the stack, calculate a transformation matrix to align it to the previous image
  warps = [numpy.eye(2, 3, dtype=numpy.float32)]
  for i in range(1, G.shape[0]):
      warps.append(align(G[i-1, :, :], G[i, :, :]))
      bar.update(i+1)
  
  
  print("\n\nApplying alignment transformation...")
  # Define 2x3 matrix and initialize the matrix to identity
  for i in range(1, G.shape[0]):
      warps[i][0, 2] = warps[i][0, 2] + warps[i-1][0, 2]
      warps[i][1, 2] = warps[i][1, 2] + warps[i-1][1, 2]
  
  
  
  i = 0
  aligned = cv2.warpAffine(S[i, :, :, :], warps[i], (S.shape[2], S.shape[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)
  
  fig, ax = plt.subplots()
  
  fig.canvas.mpl_connect('key_press_event', press)
  
  ax.imshow(aligned.astype(numpy.uint8))
  ax.set_title('Slice ' + str(i))
  ax.set_xlabel("Press 'z' and 'x' to change slices and arrow keys to align")
  plt.show()
  
  I = apply_alignment(S, warps)
59b0e5f5   David Mayerich   updated envi file...
158
  imagestack.save(I, out_dir + "/aligned_", "bmp")