Openholo  v4.2
Open Source Digital Holographic Library
Depthmap.py
Go to the documentation of this file.
1 import numpy as np
2 import matplotlib.pyplot as plt
3 from numba import njit
4 from PIL import Image
5 
6 from concurrent.futures import ProcessPoolExecutor
7 import multiprocessing
8 
9 # parameters
10 mm = 1e-3
11 um = mm * mm
12 nm = um * mm
13 
14 
15 @njit(nogil=True, cache=True)
16 def k(wvl):
17  return (np.pi * 2) / wvl
18 
19 
20 @njit(nogil=True)
21 def h_frsn(pixel_pitch_x, pixel_pitch_y, nx, ny, zz, wvl):
22  re = np.zeros((ny, nx))
23  im = np.zeros((ny, nx))
24  for i in np.arange(nx):
25  for j in np.arange(ny):
26  x = (i - nx / 2) * pixel_pitch_x
27  y = (j - ny / 2) * pixel_pitch_y
28  re[j, i] = np.cos((np.pi / (wvl * zz)) * (x * x + y * y))
29  im[j, i] = np.sin((np.pi / (wvl * zz)) * (x * x + y * y))
30  return re + 1j * im
31 
32 
33 @njit(nogil=True)
34 def alphamap(depthmap, n):
35  """ extract alpha map """
36  shape = depthmap.shape
37  ww = shape[1]
38  hh = shape[0]
39  amap = np.zeros((hh, ww))
40  for i in np.arange(ww):
41  for j in np.arange(hh):
42  if depthmap[j, i] == n:
43  amap[j, i] = 1
44  return amap
45 
46 
47 @njit(nogil=True)
48 def refwave(wvl, wr, hr, z, pp, thetaX, thetaY):
49  a = np.zeros((hr, wr))
50  b = np.zeros((hr, wr))
51  for i in np.arange(hr):
52  for j in np.arange(wr):
53  x = (j - wr / 2) * pp
54  y = -(i - hr / 2) * pp
55  a[i, j] = np.cos(k(wvl) * (x * np.sin(thetaX) + y * np.sin(thetaY)))
56  b[i, j] = np.sin(k(wvl) * (x * np.sin(thetaX) + y * np.sin(thetaY)))
57  return (a / z) + 1j * (b / z)
58 
59 
61  """
62  Get Fringe pattern by 2D Depth map image and RGB color image
63 
64  Parameters
65  imgpath : image path
66  f : propagation length
67  angle : phase shift angle
68  Red, Green, Blue : wavelength
69  scale : scaling factor
70 
71  Depth map parameters
72  field_len = 1000e-3
73  near_depth = 800e-3
74  far_depth = 1200e-3
75 
76  Depth quantization (깊이 계조)
77  DQ = 256 # 0 to 255
78 
79  Unit depth
80  UD = -(far_depth - near_depth) / 256
81  http://openholo.org/
82  """
83  def __init__(self, RGBimg, Depthimg, f=1, angleX=0, angleY=0,
84  Red=639*nm, Green=525*nm, Blue=463*nm, SLM_width=3840, multicore=True,
85  SLM_height=2160, scale=0.03, pixel_pitch=3.6*um, zeropadding=3840,
86  field_length=1000e-3, near_depth=800e-3, far_depth=1200e-3, DepthQuantization=256):
87  self.zz = f
88  self.imagein = np.asarray(Image.open(RGBimg))
89  self.depthimg = np.asarray(Image.open(Depthimg))
90  self.depthimg = self.depthimg[:, :, 1]
91  self.thetaX = angleX * (np.pi / 180)
92  self.thetaY = angleY * (np.pi / 180)
93  self.wvl_R = Red
94  self.wvl_G = Green
95  self.wvl_B = Blue
96  self.w = SLM_width
97  self.h = SLM_height
98  self.pp = pixel_pitch
99  self.scale = scale
100  self.UD = -(far_depth - near_depth) / 256
101  self.DQ = DepthQuantization
102  self.nzp = zeropadding
103  self.img_R = np.double(self.resizeimg(self.wvl_R, self.imagein[:, :, 0]))
104  self.img_G = np.double(self.resizeimg(self.wvl_G, self.imagein[:, :, 1]))
105  self.img_B = np.double(self.resizeimg(self.wvl_B, self.imagein[:, :, 2]))
106  if multicore:
107  self.num_cpu = multiprocessing.cpu_count() // 2
108  else:
109  self.num_cpu = 1
110  self.num_point = [i for i in range(self.DQ)]
111  self.num_point = np.array(self.num_point)
112 
113  def resizeimg(self, wvl, img):
114  """RGB 파장에 맞게 원본 이미지를 리사이징 + zero padding"""
115  w_n = int(self.w * (self.wvl_B / wvl))
116  h_n = int(self.h * (self.wvl_B / wvl))
117  img_new = np.zeros((self.h * 2, self.w * 2))
118  im = Image.fromarray(img)
119  im = im.resize((w_n, h_n), Image.BILINEAR) # resize image
120  im = np.asarray(im)
121  img_new[(self.h * 2 - h_n) // 2:(self.h * 2 + h_n) // 2, (self.w * 2 - w_n) // 2:(self.w * 2 + w_n) // 2] = im
122  return img_new
123 
124  def Prop(self, color, image, amap, n):
125  if color == 'green':
126  wvl = self.wvl_G
127  elif color == 'blue':
128  wvl = self.wvl_B
129  else:
130  wvl = self.wvl_R
131  # resize image
132  imgs = image * amap
133  imgs = np.flip(imgs, axis=0)
134  zzz = self.UD * n + self.zz
135  phase = np.random.random((2 * self.h, 2 * self.w)) * 2 * np.pi # random phase
136  ph = np.exp(1j * phase)
137  ph *= imgs
138  ps = self.scale / (2 * self.w)
139  ch2 = ph * h_frsn(ps, ps, self.w + self.w, self.h + self.h, zzz, wvl)
140  CH1 = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(ch2)))
141  result = CH1 * h_frsn(self.pp, self.pp, self.w + self.w, self.h + self.h, zzz, wvl)
142  result = result[self.h // 2: (3 * self.h) // 2, self.w // 2: (3 * self.w) // 2]
143  result *= refwave(wvl, self.w, self.h, zzz, self.pp, self.thetaX, self.thetaY)
144  return result
145 
146  def prop_R(self, n):
147  wvl = self.wvl_R
148  dimg = self.resizeimg(wvl, self.depthimg)
149  amap = alphamap(dimg, n)
150  return self.Prop('red', self.img_R, amap, n)
151 
152  def prop_G(self, n):
153  wvl = self.wvl_G
154  dimg = self.resizeimg(wvl, self.depthimg)
155  amap = alphamap(dimg, n)
156  return self.Prop('green', self.img_G, amap, n)
157 
158  def prop_B(self, n):
159  wvl = self.wvl_B
160  dimg = self.resizeimg(wvl, self.depthimg)
161  amap = alphamap(dimg, n)
162  return self.Prop('blue', self.img_B, amap, n)
163 
164  def parallelCal(self, color):
165  if color == 'green':
166  fun = self.prop_G
167  elif color == 'blue':
168  fun = self.prop_B
169  else:
170  fun = self.prop_R
171  H = np.zeros((self.h, self.w), dtype='complex128')
172  count = np.split(self.num_point, [i * self.num_cpu for i in range(1, self.DQ // self.num_cpu)])
173  for n in count:
174  with ProcessPoolExecutor(self.num_cpu) as ex:
175  cache = [result for result in ex.map(fun, list(n))]
176  cache = np.asarray(cache)
177  print(n, ' depth done')
178  for j in range(len(n)):
179  H += cache[j, :, :]
180  return H
181 
182 
183  def fft(self, f):
184  return np.fft.fftshift(np.fft.fft2(np.fft.fftshift(f)))
185 
186  def ifft(self, f):
187  return np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(f)))
188 
189  def SSB(self, ch):
190  """single side band encoding"""
191  height, width = ch.shape
192  a = np.zeros((height, width), dtype='complex128')
193  CH = self.fft(ch) # fourier transformed image
194  CH = CH[height // 4: (height * 3) // 4, :]
195  a[0:height // 2, :] = CH
196  a[height // 2:, :] = np.conj(CH)
197  return self.ifft(a)
198 
199  def normalize(self, arr, type='angle'):
200  """normalize"""
201  if type == 'phase':
202  arrin = np.copy(np.imag(arr))
203  elif type == 'real':
204  arrin = np.copy(np.real(arr))
205  elif type == 'angle':
206  arrin = np.copy(np.angle(arr))
207  elif type == 'amplitude':
208  arrin = np.copy(np.abs(arr))
209  else:
210  arrin = np.copy(arr)
211  # arrin = np.float(arrin)
212  arrin -= np.min(arrin)
213  # arrin = arrin + np.abs(arrin)
214  arrin = arrin / np.max(arrin)
215  return arrin
216 
217  def getRGBImage(self, R, G, B, fname, type='angle'):
218  """Get RGB image"""
219  h, w = R.shape
220  img = np.zeros((h, w, 3))
221  img[:, :, 0] = self.normalize(R, type)
222  img[:, :, 1] = self.normalize(G, type)
223  img[:, :, 2] = self.normalize(B, type)
224  plt.imsave(fname, img)
225  return img
226 
227  def getMonoImage(self, ch, fname):
228  """Get Single channel image"""
229  im = self.normalize(ch, 'phase')
230  phase = fname + '_IM.bmp'
231  plt.imsave(phase, im, cmap='gray')
232  re = self.normalize(ch, 'real')
233  real = fname + '_RE.bmp'
234  plt.imsave(real, re, cmap='gray')
235  return im, re
236 
237 
def normalize(self, arr, type='angle')
Definition: Depthmap.py:199
def __init__(self, RGBimg, Depthimg, f=1, angleX=0, angleY=0, Red=639 *nm, Green=525 *nm, Blue=463 *nm, SLM_width=3840, multicore=True, SLM_height=2160, scale=0.03, pixel_pitch=3.6 *um, zeropadding=3840, field_length=1000e-3, near_depth=800e-3, far_depth=1200e-3, DepthQuantization=256)
Definition: Depthmap.py:86
def alphamap(depthmap, n)
Definition: Depthmap.py:34
def refwave(wvl, wr, hr, z, pp, thetaX, thetaY)
Definition: Depthmap.py:48
def k(wvl)
Definition: Depthmap.py:16
def Prop(self, color, image, amap, n)
Definition: Depthmap.py:124
def getMonoImage(self, ch, fname)
Definition: Depthmap.py:227
def resizeimg(self, wvl, img)
Definition: Depthmap.py:113
def getRGBImage(self, R, G, B, fname, type='angle')
Definition: Depthmap.py:217
def parallelCal(self, color)
Definition: Depthmap.py:164
def h_frsn(pixel_pitch_x, pixel_pitch_y, nx, ny, zz, wvl)
Definition: Depthmap.py:21