-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathshrinkwrap.py
160 lines (132 loc) · 6.09 KB
/
shrinkwrap.py
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 2 11:04:38 2017
@author: HonkyT
"""
#from numpy import fft # om man skriver såhär begöver man itne använda np.
# men eftersom jag behöver fler saker från nupy och jag redan skrivit np på alla gör jag inte det.
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from create2Dgaussian import create2Dgaussian #import all functions
from scipy import misc
from math import ceil
def data():
# data = np.genfromtxt("diffractionPattern.txt", dtype='float')
#image = misc.imread('P.png',flatten=True)
#p = misc.imread('star.bmp',flatten=True)
#absF = abs(np.fft.fftshift(np.fft.fft2(p)))
########### Test att använda fft av kristall som inputi shrinkwrap
crystal = np.zeros((201,201), dtype= np.int32)
dx=1
for row in range(60,140,dx):
for col in range(60,140,dx):
crystal[row,col] = 1
crystal_fourier = np.fft.fftshift(np.fft.fft2(crystal))
absF = abs(crystal_fourier)**2
return (absF,crystal)
absF, crystal = data()
#############Test att använda bild som är skjuvad coord
#image = misc.imread('star.bmp',flatten=True)
##image = misc.imread('P.png',flatten=True)
#data = abs(np.fft.fftshift(np.fft.fft2(image)))
#
#theta = 0*np.pi / 180 #rad
#r3 = 1 + 1/np.cos(theta) #antal pixlar som motsvarar 1 pixel i xled i xz systemet
#r1 = 1 + 1/ np.sin(theta)
#
#image_skewed = np.zeros((image.shape[0], ceil((image.shape[1]/np.cos(theta)))+ 50 ))
#
#for i in range(0,image.shape[0]):
# for j in range(0,image.shape[1]):
# xs = ceil(j / (1 + np.tan(theta)**2) + i*np.tan(theta)/ ( 1 + np.tan(theta)**2) )
# #np.disp(xs)
# ys = i
# image_skewed[ys,xs] = image[i,j]
# #image_skewed =
#
## cut ot hte skewed image so that it is ccentered
#image_skewed = image_skewed[:,0:120] # Hamnar rätt om man klipper den rätt
#image_skewed = image_skewed/image_skewed.max()
#fft_image_skewed = abs(np.fft.fftshift(np.fft.fft2(image_skewed)))
#absF = fft_image_skewed
def shrinkwrap(absF):
# the intensity of the diffraction pattern is F^2
absF = pow(absF,0.5)
# parameters
beta = 0.9 # HIO feedback parameter
thresholdMask1 = 0.04; # Initial threshold for mask
thresholdMask2 = 0.2; # threshold for mask
sigma = 3; # initial value of gaussian width
c = absF.shape[0]/2 # center of Gaussian
n = 100 #number of iterations
# initialize vector for error calculation
# err2 = np.zeros(shape=(1,n))
# g represents the unknown object. Create initial guess for g using
# the autocorrelation function of the object. Start by Fourier
# transform the diffraction pattern:
g = np.fft.fftshift(np.fft.ifft2(absF)) #finns några val som symmetric?
# create autocorrelation function of g
#kanske finns än funtion i python som gör det här
g = abs(np.fft.fftshift(np.fft.ifft2(np.fft.fft2(g)*np.conj(np.fft.fft2(g)))))/(g.shape[0]*g.shape[1])
# create a logical matrix that is the Shrinkwrap support
support = g > thresholdMask1*g.max()
# define iteration counter
k = 0
# Start of iterations
while k<n:
# every 20th iteration, update the support
if k%20==0 and k>0: # VARFÖR SKULLE K VA NEG?
# call function create2Dgaussian
gauss2D = create2Dgaussian(sigma,c,absF.shape[0],absF.shape[1])
# calculate the convolution the absolute value of the object wave
# with a gaussian
# finns en djungel av olika laternativ hur man convolverar
support = signal.convolve2d(gauss2D,abs(g),'same')
# Create a logical matrix that is the Shrinkwrap support
support = support > thresholdMask2*support.max() #SÄKERT ATT DEN HITTAR 2D MAXET?
# reduce sigma with 1% every 20th iteration until a limit of 1.5
if sigma >= 1.5:
sigma = sigma*0.99
# STEP 1: Fourier transform of g(x), G(u) :
G = np.fft.fftshift(np.fft.fft2(g))
# STEP 2: make |G| confirm with |F|
Gprime = absF*np.exp(1j*np.angle(G))
# STEP 3: inverse Fourier transform
gprime = np.fft.ifft2(np.fft.ifftshift(Gprime))
# STEP 4: See to that gprime satisfies its constraints:
# create inverse of the support (compuational aid)
support_inv = np.logical_not(support)
# update g(l,j) outside the support
gprime = gprime*support + g*support_inv - beta*gprime*support_inv
# update g'(l,j) inside the support in one of
# the two following ways:
# gprime(l,j) = g(l,j)(in matrix multiplication i.e.:)
#gprime = gprime*support_inv + g*support
#or
# gprime[l,j] = gprime[l,j] (i.e. not at all)
# for l in range(0,absF.shape[0]):
# for j in range(0,absF.shape[1]): ## ska DET VARA RANGE 0 ÄR FÖRSTA INDEXEDT 0?
# if support[l,j] == 0:
# # update g(l,j) outside the support
# gprime[l,j] = g[l,j] - beta * gprime[l,j]
#
# else:
# # update g'(l,j) inside the support in one of
# # the two following ways:
# #gprime(l,j) = g(l,j);
# gprime[l,j] = gprime[l,j]
# overwrite g with result from iteration
g = gprime
# set all negative values of g to zero
g[g<0] = 0
# err2[0][k] = sum(sum( (abs(abs(np.fft.fftshift(np.fft.fft2(((g*np.conj(g))))))) - absF)**2)) / sum(sum(absF**2))
k = k+1
return g
g = shrinkwrap(absF)
plt.figure()
plt.imshow(abs(g), cmap='gray')
plt.title('reconstructed object')
plt.figure()
plt.imshow(abs(crystal), cmap='gray')
plt.title('original object')