-
Notifications
You must be signed in to change notification settings - Fork 0
/
2D-Algo
179 lines (118 loc) · 4.27 KB
/
2D-Algo
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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
# coding: utf-8
#2-D signals
import numpy as np
import matplotlib.pyplot as plt
from cmath import polar
from PIL import Image
# in signal_noisy, we put the padded noisy image (we put the inital image into a bigger one to be able to shift it)
image_file = Image.open("/Users/victorstorchan/Desktop/RA/ra/apple.png") # open colour image
image_signal = image_file.convert('L') # convert image to black and white
signal=np.asarray(image_signal)
(a,b)=signal.shape
randommat=np.zeros((400,400),dtype=complex)
for i in range(400):
for j in range(400):
randommat[i][j]+=complex(np.random.normal(0,1,1))
signal_noisy=np.zeros((400,400),dtype=complex)
for i in range(a):
for j in range(b):
signal_noisy[10+i][10+j]=signal[i][j]
for i in range(400):
for j in range(400):
signal_noisy[i][j]+=randommat[i][j]
im3=signal_noisy.astype(np.uint8)
img=Image.fromarray(im3,'L')
img.save('/Users/victorstorchan/Desktop/RA/ra/apple_true_sig_nooise1.png')
#(m,n) is the size of the noisy image, vect_of_shift is the vector of shift,
# in the files 2Dsignal_to_noise_* we will let this vector vary between (0,40)
m=signal_noisy.shape[0]
n=signal_noisy.shape[1]
vect_of_shift=[(0,0),(2,0),(2,5),(4,3)]
len_shift=len(vect_of_shift)
#signal with shift:
shifted_signals=[]
shifted_signals_1=[]
shifted_signals_2=[]
for s in range(len_shift):
y=np.zeros((n,m),dtype=complex)
y1=np.zeros((n,m),dtype=complex)
y2=np.zeros((n,m),dtype=complex)
for k in range(m):
for l in range(n):
if (l<10+vect_of_shift[s][0] or l>315+vect_of_shift[s][0]) and (k<10+vect_of_shift[s][1] or k>324+vect_of_shift[s][1]):
y[k][l]=randommat[k][l]
y1[k][l]=randommat[k][l]*np.exp(2J*np.pi*k/m)
y2[k][l]=randommat[k][l]*np.exp(2J*np.pi*l/n)
else:
y[k][l]=signal_noisy[k-vect_of_shift[s][0]][l-vect_of_shift[s][1]]
y1[k][l]=signal_noisy[k-vect_of_shift[s][0]][l-vect_of_shift[s][1]]*np.exp(2J*np.pi*k/m)
y2[k][l]=signal_noisy[k-vect_of_shift[s][0]][l-vect_of_shift[s][1]]*np.exp(2J*np.pi*l/n)
shifted_signals.append(y)
shifted_signals_1.append(y1)
shifted_signals_2.append(y2)
A=[]
A_1=[]
A_2=[]
for i in range(len(shifted_signals)):
A.append(np.matrix(np.fft.fft2(shifted_signals[i])))
A_1.append(np.matrix(np.fft.fft2(shifted_signals_1[i])).conjugate())
A_2.append(np.matrix(np.fft.fft2(shifted_signals_2[i])).conjugate())
G1=[]
G2=[]
for s in range(len_shift):
G1.append(np.multiply(A[s],A_1[s]))
G2.append(np.multiply(A[s],A_2[s]))
A_mat1=np.zeros((len_shift,n**2),dtype=complex)
for s in range(len_shift):
for i in range(n):
for k in range(n):
A_mat1[s][400*i+k]=G1[s][i,k]
A_mat2=np.zeros((len_shift,n**2),dtype=complex)
for s in range(len_shift):
for i in range(n):
for k in range(n):
A_mat2[s][400*i+k]=G2[s][i,k]
A_mat1_mat=np.matrix(A_mat1)
A_mat2_mat=np.matrix(A_mat2)
A_mat1_transpose=A_mat1_mat.getH()
A_mat2_transpose=A_mat2_mat.getH()
A_prod1=A_mat1_mat*A_mat1_transpose
A_prod2=A_mat2_mat*A_mat2_transpose
A_final1=A_prod1/A_prod1[0,0]
A_final2=A_prod2/A_prod2[0,0]
#perform the SVD
(V1,sigma1,V_star1)=np.linalg.svd(A_final1)
(V2,sigma2,V_star2)=np.linalg.svd(A_final2)
v1=V_star1[0].getH()
v2=V_star2[0].getH()
#the shifts are recovered:
output1=np.zeros(len_shift)
for i in range(len_shift):
output1[i]=-n*polar(-v1[i,0])[1]/(2*np.pi).real
output1
#the shifts are recovered:
output2=np.zeros(len_shift)
for i in range(len_shift):
output2[i]=-n*polar(-v2[i,0])[1]/(2*np.pi).real
output2
recover_A=[]
for i in range(len_shift):
M=np.matrix(np.zeros((m,n),dtype=complex))
for l in range(m):
for k in range(n):
M[l,k]=A[i][l,k]/np.exp(-2J*np.pi*(l*output1[i]+k*output2[i])/n)
recover_A.append(M)
A_final=[]
for i in range(len_shift):
A_final.append(np.fft.ifft2(recover_A[i]))
recov_signal=np.zeros((m,n))
for i in range(m):
for j in range(n):
k=0
for s in range(len_shift):
k+=A_final[s][i][j].real
recov_signal[i][j]=k/len_shift
recov_signal1=recov_signal.astype(np.uint8)
img_recov = Image.fromarray(recov_signal1,'L')
img_recov
img_recov.save('/Users/victorstorchan/Desktop/RA/ra/apple_recov_noise1.png')