-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathgenerate_spectra.py
155 lines (122 loc) · 4.94 KB
/
generate_spectra.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
#
# Copyright (c) 2023 Image Processing Research Group of University Federico II of Naples ('GRIP-UNINA').
# All rights reserved.
# This work should only be used for nonprofit purposes.
#
# By downloading and/or using any of these files, you implicitly agree to all the
# terms of the license, as specified in the document LICENSE.txt
# (included in this package) and online at
# http://www.grip.unina.it/download/LICENSE_OPEN.txt
#
import matplotlib.pyplot as plt
from tqdm import tqdm
import numpy as np
from PIL import Image
import glob
import os
import argparse
import cv2
import random
def imread(filename):
return np.asarray(Image.open(filename).convert('RGB'))/256.0
def rescale_img(img, siz):
h, w = img.shape[:2]
m = min(w, h)
if m != siz:
dim = (siz*w//m, siz*h//m)
# resize image
if siz < m:
img = cv2.resize(img, dim, interpolation=cv2.INTER_AREA)
else:
img = cv2.resize(img, dim, interpolation=cv2.INTER_LINEAR)
h, w = img.shape[:2]
assert min(w, h) == siz
py = (h - siz)//2
px = (w - siz)//2
return img[py:(py+siz), px:(px+siz)]
def get_fft2(x):
x = np.float64(x)
x = x - np.mean(x, (-3, -2, -1), keepdims=True)
x = x/np.sqrt(np.mean(np.abs(x**2), (-3, -2, -1), keepdims=True))
x = np.fft.fft2(x, axes=(-3, -2), norm='ortho')
x = np.abs(x)**2
return x
def get_spectrum(power_spec, q_step=None):
power_spec = np.mean(power_spec, -1)
power_spec = power_spec / power_spec.size
H, W = power_spec.shape
h, w = np.meshgrid(np.fft.fftfreq(H), np.fft.fftfreq(W), indexing='ij')
r = np.sqrt(h**2 + w**2)
if q_step is None:
q_step = 1.0/min(H, W)
r_quant = np.round(r/q_step)
freq = np.sort(np.unique(r_quant))
y = np.asarray([np.sum(power_spec[r_quant == f]) for f in freq])
return y, q_step*freq
def get_spectrum_angular(power_spec, num = 16):
power_spec = np.mean(power_spec, -1)
power_spec = power_spec / power_spec.size
H, W = power_spec.shape
h, w = np.meshgrid(np.fft.fftfreq(H), np.fft.fftfreq(W), indexing='ij')
r = np.sqrt(h**2 + w**2)
angular = np.round(num * np.arctan2(h, w) / np.pi) % num
ang_freq = np.sort(np.unique(angular))
y = np.asarray([np.sum(power_spec[(angular==f) & (r>0.1)]) for f in ang_freq])
return y, ang_freq/num
def get_spectra(files_path, output_dir, output_code):
print("Starting generation of spectra")
print(files_path)
filenames = glob.glob(files_path + "/*")
print(len(filenames))
random.shuffle(filenames)
siz = 256
filenames = filenames[:1000]
print("Starting to generate fingerprints")
img_fft2 = [get_fft2(rescale_img(imread(_), siz)) for _ in tqdm(filenames)]
freq = get_spectrum(img_fft2[0])[1]
ang_freq = np.pi*get_spectrum_angular(img_fft2[0])[1]
spectra = [get_spectrum(_)[0] for _ in tqdm(img_fft2)]
ang_spectra = [get_spectrum_angular(_)[0] for _ in tqdm(img_fft2)]
spectra_mean = np.mean(spectra, 0)
ang_spectra_mean = np.mean(ang_spectra, 0)
spectra_var = np.var(spectra, 0)
ang_spectra_var = np.var(ang_spectra, 0)
figures_output_dir = os.path.join(output_dir, output_code)
os.makedirs(figures_output_dir, exist_ok=True)
dict_out = dict()
dict_out['freq'] = freq
dict_out['spectra_mean'] = spectra_mean
dict_out['spectra_var'] = spectra_var
dict_out['ang_freq'] = ang_freq
dict_out['ang_spectra_mean'] = ang_spectra_mean
dict_out['ang_spectra_var'] = ang_spectra_var
np.savez(figures_output_dir+'/spectra.npz', **dict_out)
# save figures
fig = plt.figure(figsize=(6, 5))
plt.plot(freq, spectra_mean, linewidth=2)
plt.xlabel('$freq$', fontsize=10)
plt.yticks(fontsize=10)
plt.xticks(fontsize=10)
plt.xlim([0.2, 0.5])
plt.ylim([0.0, 0.0012])
plt.grid()
fig.savefig(figures_output_dir+'/spectra.png',
bbox_inches='tight', pad_inches=0.0)
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': 'polar'}, figsize=(6, 6))
ang_spectra_mean = np.concatenate((ang_spectra_mean, ang_spectra_mean, ang_spectra_mean[...,:1]),-1)
ang_freq = np.concatenate((ang_freq, np.pi+ang_freq, 2*np.pi+ang_freq[:1]),0)
ax.plot(ang_freq, ang_spectra_mean, linewidth=2)
ax.set_yticks(ax.get_yticks(), list())
ax.grid('on')
fig.savefig(figures_output_dir+'/ang_spectra.png',
bbox_inches='tight', pad_inches=0.0)
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--files_path", type=str,
help="The path where the images are stored")
parser.add_argument("--out_dir", type=str,
help="The path where to save the images")
parser.add_argument("--out_name", type=str,
help="The name of the folder in which to save the images and the numpy arrays")
args = vars(parser.parse_args())
get_spectra(args['files_path'], args['out_dir'], args['out_name'])