-
Notifications
You must be signed in to change notification settings - Fork 0
/
m_specfun.py
1700 lines (1560 loc) · 70.3 KB
/
m_specfun.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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -------------------------------------------------------------------
# m_specfun functions for m_spec
# Author: Martin Dubs, 2020
# -------------------------------------------------------------------
import configparser
import ctypes
import logging
import os
import os.path as path
import platform
import subprocess
import time
import warnings
from datetime import datetime, date
import PySimpleGUI as sg
import numpy as np
from astropy.io import fits
from astropy.time import Time
from scipy import optimize, interpolate
from scipy.ndimage import map_coordinates
from skimage import img_as_float
from skimage import transform as tf
from skimage import io as ios
from PIL import Image
import PIL
import io
import base64
if platform.system() == 'Windows':
ctypes.windll.user32.SetProcessDPIAware() # Set unit of GUI to pixels
version = '0.9.21'
# today = date.today()
logfile = 'm_spec' + date.today().strftime("%y%m%d") + '.log'
# turn off other loggers
for handler in logging.root.handlers[:]:
logging.root.removeHandler(handler)
logging.basicConfig(filename=logfile, format='%(asctime)s %(message)s', level=logging.INFO)
# -------------------------------------------------------------------
# initialize dictionaries for configuration
parv = ['1' for x in range(10)] + ['' for x in range(10, 15)]
parkey = ['f_lam0', 'f_scalxy', 'b_fitxy', 'i_imx', 'i_imy', 'f_f0', 'f_pix', 'f_grat', 'f_rotdeg', 'i_binning',
's_comment', 's_infile', 's_outfil', 's_linelist', 'b_sqrt']
par_dict = dict(list(zip(parkey, parv)))
resv = [0.0 for x in range(7)]
reskey = ['scalxy', 'x00', 'y00', 'rot', 'disp0', 'a3', 'a5']
res_dict = dict(list(zip(reskey, resv)))
fitsv = ['' for x in range(7)]
fitskey = ['DATE-OBS', 'OBSERVER', 'VERSION', 'INSTRUME', 'TELESCOP', 'M_STATIO', 'COMMENT']
fits_dict = dict(list(zip(fitskey, fitsv)))
bc_enabled = ('white', 'green') # button_color
bc_disabled = (None, 'darkblue')
# default values for setup, if no m_set.ini
debug = False
fit_report = False
win2ima = False
zoom = 1.0
wsize = (1060, 660)
wloc = (50, 20)
(xoff_calc, yoff_calc) = (355, 50)
(xoff_setup, yoff_setup) = (250, 120)
wloc_calc = (wloc[0] + xoff_calc, wloc[1] + yoff_calc)
wloc_setup = (wloc[0] + xoff_setup, wloc[1] + yoff_setup)
(max_width, max_height) = (700, 500)
opt_comment = ''
pngdir = ''
png_name = 'tmp/m_'
outpath = 'out'
mdist = 'mdist'
colorflag = False
bob_doubler = False
plot_w = 1000
plot_h = 500
i_min = -0.5
i_max = 5
graph_size = 2000
show_images = True
optkey = ['zoom', 'win_width', 'win_height', 'win_x', 'win_y', 'calc_off_x',
'calc_off_y', 'setup_off_x', 'setup_off_y', 'debug', 'fit-report',
'scale_win2ima', 'comment', 'png_name', 'outpath', 'mdist', 'colorflag', 'bob',
'plot_w', 'plot_h', 'i_min', 'i_max', 'graph_size', 'show_images']
optvar = [zoom, wsize[0], wsize[1], wloc[0], wloc[1], xoff_calc, yoff_calc,
xoff_setup, yoff_setup, debug, fit_report, win2ima, opt_comment, png_name,
outpath, mdist, colorflag, bob_doubler, plot_w, plot_h, i_min, i_max, graph_size, show_images]
opt_dict = dict(list(zip(optkey, optvar)))
# -------------------------------------------------------------------
def read_configuration(conf, par_dict, res_dict, opt_dict):
"""
read configuration file for m_calib and m-spec
:param conf: filename of configuration with extension .ini
:param par_dict: parameters for m_calib, partly used in m-spec
:param res_dict: results of m_calib, used for distortion
:param opt_dict: options, used in m_calib and m_spec
:return:
partext: multiline text of configuration
updated values of par_dict, res_dict, fits_dict, opt_dict
# from readconf in m_config3.py
"""
partext = ''
pngdir = ''
if path.exists(conf):
config = configparser.ConfigParser()
config.read(conf)
for section in config.sections():
partext += f'[{section}]\n'
for key in config[section]:
partext += f'- [{key}] = {config[section][key]}\n'
for key in config['Lasercal'].keys():
k = key[0]
if k == 'b':
if config['Lasercal'][key] == '0':
par_dict[key] = False
else:
par_dict[key] = True
elif k == 'f':
par_dict[key] = float(config['Lasercal'][key])
elif k == 'i':
par_dict[key] = int(config['Lasercal'][key])
elif k == 's':
par_dict[key] = config['Lasercal'][key]
else:
print('unknown key in readconf: ', key)
if 'Calib' in config.sections():
for key in config['Calib']:
res_dict[key] = float(config['Calib'][key])
if 'Fits' in config.sections():
for key in config['Fits']:
fits_dict[key.upper()] = config['Fits'][key]
if 'Options' in config.sections():
for key in config['Options'].keys():
if key in (
'win_width', 'win_height', 'win_x', 'win_y', 'calc_off_x', 'calc_off_y', 'setup_off_x',
'setup_off_y', 'graph_size'):
opt_dict[key] = int(config['Options'][key])
elif key in ('debug', 'fit-report', 'scale_win2ima', 'scale_ima2win',
'colorflag', 'bob', 'show_images'):
opt_dict[key] = bool(int(config['Options'][key]))
elif key in ('zoom', 'i_min', 'i_max'):
opt_dict[key] = float(config['Options'][key])
else:
if key == 'pngdir':
pngdir = config['Options'][key]
else:
opt_dict[key] = config['Options'][key]
opt_dict['png_name'] = path.join(pngdir, opt_dict['png_name']) # used for compatibility with old inifile
logging.info(f' configuration {conf} loaded')
return partext, par_dict, res_dict, fits_dict, opt_dict
# ------------------------------------------------------------------------------
def write_configuration(conf, par_dict, res_dict, fits_dict, opt_dict):
"""
writes configuration to conf
:param conf: filename with ext .ini
:param par_dict: parameters for m_calib, partly used in m-spec
:param res_dict: results of m_calib, used for distortion
:param fits_dict: content of fits header
:param opt_dict: options, used in m_calib and m_spec
:return: None
"""
def configsetbool(section, option, boolean):
if boolean:
config.set(section, option, '1')
else:
config.set(section, option, '0')
# for compatibility with old versions
pngdir, opt_dict['png_name'] = path.split(opt_dict['png_name'])
config = configparser.ConfigParser()
cfgfile = open(conf, 'w')
if 'Lasercal' not in config: config.add_section('Lasercal')
if 'Calib' not in config: config.add_section('Calib')
if 'Fits' not in config: config.add_section('Fits')
if 'Options' not in config: config.add_section('Options')
for key in par_dict.keys():
k = key[0]
if k == 'b':
configsetbool('Lasercal', key, par_dict[key])
elif k == 'i':
config.set('Lasercal', key, str(par_dict[key]))
elif k == 'f':
config.set('Lasercal', key, str(par_dict[key]))
elif k == 's':
config.set('Lasercal', key, par_dict[key])
else:
print('unknown key in writeconf: ', key)
for key in res_dict.keys():
config.set('Calib', key, str(res_dict[key]))
for key in fits_dict.keys():
config.set('Fits', key.upper(), str(fits_dict[key]))
for key in opt_dict.keys():
if key in ('debug', 'fit-report', 'scale_win2ima', 'scale_ima2win', 'colorflag', 'bob', 'show_images'):
configsetbool('Options', key, opt_dict[key])
else:
config.set('Options', key, str(opt_dict[key]))
config.set('Options', 'pngdir', str(pngdir))
config.write(cfgfile)
logging.info(f' configuration saved as {conf}')
cfgfile.close()
# -------------------------------------------------------------------
def write_fits_image(image, filename, fits_dict, dist=True):
"""
writes image as 32-bit float array into fits-file
:param image: np.array with image data, scaled to +/- 1.0, b/w or color
:param filename: filename with extension .fit
:param fits_dict: content of fits header
:param dist: True: distorted image; False: undistorted image
:return: 1 if error, else 0
"""
if len(image.shape) == 3:
if image.shape[2] > 3: # cannot read plots with multiple image planes
sg.PopupError('cannot convert png image, try bmp or fit')
return 1
image = np.transpose(image, (2, 0, 1))
hdu = fits.PrimaryHDU(image.astype(np.float32))
hdul = fits.HDUList([hdu])
fits_dict['BSCALE'] = 32767
fits_dict['BZERO'] = 0
fits_dict['COMMENT'] = str(fits_dict['COMMENT']) # [:20]
for key in fits_dict.keys():
if dist:
hdu.header[key] = fits_dict[key]
else:
if key not in ('D_SCALXY', 'D_X00', 'D_Y00', 'D_ROT', 'D_DISP0', 'D_A3', 'D_A5'):
hdu.header[key] = fits_dict[key]
hdul.writeto(filename, overwrite=True)
hdul.close()
return 0
# -------------------------------------------------------------------
def get_png_image(filename, colorflag=False):
"""
reads png image and converts to np.array
:param filename: with extension 'png
:param colorflag: True: colour image, False: image converted to b/w
:return: image as 2 or 3-D array
"""
image = np.flipud(img_as_float(ios.imread(filename)))
if not colorflag and len(image.shape) == 3:
image = np.sum(image, axis=2) / 3
return image
# -------------------------------------------------------------------
def extract_video_images(avifile, pngname, bobdoubler, binning, bff, maxim):
"""
creates png images from AVI file
:param avifile: filename of avi file (full path, with extension)
:param pngname: filebase of png images, e.g. tmp/m for series m1.png, m2.png,...
:param bobdoubler: if True: interlaced frames are separated into fields of half height,
default: False, frames are read
:param binning: integer:
:param bff: if True: bottom field first read for interlaced video, else top field first
:param maxim: integer, limit for converting images
:return:
nim: number of converted images, starting with index 1
dattim: date and time of video, extracted from filename created in UFO Capture
sta: station name, extracted from filename created in UFO Capture
out: full path filebase of extracted images, e.g. data/out/mdist
"""
# extract dattim and station from filename (for files from UFO capture)
def tfits(p):
# f = Path(p).name
f, ext = path.splitext(path.basename(p))
t = Time(datetime(int(f[1:5]), int(f[5:7]), int(f[7:9]), int(f[10:12]), int(f[12:14]), int(f[14:16]))).fits
sta = f[17:22]
return t, sta
# -------------------------------------------------------------------
# subprocess is os specific
sys = platform.system()
if sys == 'Windows':
cshell = False
else:
cshell = True
logging.info(f'Platform: {sys}')
out = pngname
pngdir, tmp = path.split(pngname)
nim = 0
dattim = ''
sta = ''
if avifile:
avifile = '"' + avifile + '"' # double quotes needed for filenames containing white spaces
# path name for png images
if pngdir:
if not path.exists(pngdir):
os.mkdir(pngdir)
try:
if bobdoubler:
# read bottom and top fields
command = f"ffmpeg -i {avifile} -frames {maxim / 2} -vf field=top {pngdir}/top%d.png -loglevel quiet"
subprocess.call(command, shell=cshell)
command = f"ffmpeg -i {avifile} -frames {maxim / 2} -vf field=bottom {pngdir}/bot%d.png -loglevel quiet"
subprocess.call(command, shell=cshell)
nfr = 0
n = 0
end = False
# sort and rename fields
while not end:
try:
n += 1
nfr += 1
if bff:
os.rename(f'{pngdir}/bot' + str(nfr) + '.png', out + str(n) + '.png')
n += 1
os.rename(f'{pngdir}/top' + str(nfr) + '.png', out + str(n) + '.png')
else:
os.rename(f'{pngdir}/top' + str(nfr) + '.png', out + str(n) + '.png')
n += 1
os.rename(f'{pngdir}/bot' + str(nfr) + '.png', out + str(n) + '.png')
except:
end = True
nim = n - 1
elif binning > 1:
# binning bin*bin for reducing file size
command = f"ffmpeg -i {avifile} -frames {maxim} -vf scale=iw/{binning}:-1 {out}%d.png -loglevel quiet"
subprocess.call(command, shell=cshell)
nim = check_files(out, maxim)
else:
# regular processing of frames
command = f"ffmpeg -i {avifile} -frames {maxim} {out}%d.png -loglevel quiet"
subprocess.call(command, shell=cshell)
nim = check_files(out, maxim)
if debug:
print(f'last file written: {out}' + str(nim) + '.png')
# get dattim from filename
dattim, sta = tfits(avifile)
except:
sg.PopupError('problem with ffmpeg, no images converted', title='AVI conversion')
return nim, dattim, sta, out
# -------------------------------------------------------------------
def create_file_list(file, n, ext='.png', start=1):
"""
create a file series according to IRIS convention
:param file: filebase
:param n: number of files
:param ext: extension, default = .png
:param start: index of first file
:return: file_list
"""
result = []
for a in range(start, start + n):
filen = file + str(a) + ext
result.append(filen)
return result
# -------------------------------------------------------------------
def check_files(file, n, ext='.png'):
"""
check if files in file series file+index+ext exist, starting with index 1
:param file: filebase
:param n: last index to check
:param ext: file extension, default = .png
:return: number of files found, 0 if no file exists
"""
filelist = create_file_list(file, n, ext=ext)
index = 0
for i in range(len(filelist)):
if path.exists(file + str(i + 1) + ext):
index = i + 1
else:
index = i
return index
return index
# -------------------------------------------------------------------
def delete_old_files(file, n, ext='.png'):
"""
delete files in order to clean up directory before new calculation
:param file: filebase
:param n: last index to check
:param ext: file extension, default = .png
:return:
number of files found
number of deleted files
"""
oldfiles = check_files(file, n, ext)
deleted = 0
answer = ''
if oldfiles:
answer = sg.PopupOKCancel(f'delete {oldfiles} existing files {file}, \nARE YOU SURE?', title='Delete old Files')
if answer is 'OK':
for index in range(oldfiles):
os.remove(file + str(index + 1) + ext)
deleted = oldfiles
return oldfiles, deleted, answer
# -------------------------------------------------------------------
def create_background_image(im, nb, colorflag=False): # returns background image
"""
creates background image from first nb png images extracted from
video with VirtualDub
Parameters:
im: filebase of image without number and .png extension
e.g. m_ for series m_1.png, m_2.png,...
nb: number of images, starting with index 1,
for calculation of background image
n = 0: zero intensity background image
colorflag: True: color image, False: b/w image output
Return:
background image, average of input images, as image array
"""
if nb > 0:
# create list of filenames
image_list = create_file_list(im, nb)
# open a series of files and store the data in a list, image_concat
first = True
for image in image_list:
ima = get_png_image(image, colorflag)
if first:
image_sum = ima
first = False
else:
image_sum += ima
ave_image = image_sum / nb
else:
# zero background
ave_image = 0 * get_png_image(im + '1.png', colorflag)
return ave_image
# -------------------------------------------------------------------
def apply_dark_distortion(im, backfile, outpath, mdist, first, nm, window, fits_dict, graph_size, dist=False,
background=False, center=None, a3=0, a5=0, rotation=0, yscale=1, colorflag=False,
show_images=True, cval=0):
# subtracts background and transforms images in a single step
"""
subtracts background image from png images and stores the result
as fit-images
(peak image from series, sum image from series)
Perform a dist transformation
Parameters:
im: filebase of image without number and .bmp extension
e.g. m_ for series m_1.bmp, m_2.bmp,...
backfile: background fit-file created in previous step without extension
outpath: path to mdist (output files)
mdist: file base of output files, appended with number starting from 1
(IRIS convention) and .fit
first: index of first image converted
nm: number of images created (if exist)
dist: flag, if True the distortion is calculated,
with additional parameters
background: flag, if True the background image (backfile) is subtracted
center : (column, row) tuple or (2,) ndarray, optional
Center coordinate of transformation, corresponds to optical axis.
If None, the image center is assumed
a3 : float, optional
The cubic coefficient of radial transformation
a5 : float, optional
The quintic coefficient of radial transformation
(the linear coefficient is set equal 1 to preserve image scale
at center, even order coefficients are equal zero due to the
symmetry of the transformation
rotation : float, optional
Additional rotation applied to the image.
yscale : float, optional
scales image by a factor in y-direction to compensate for non-square
pixels. The center coordinate y0 is scaled as well
colorflag: True for colour images, False for b/w images
fits_dict: dictionary with fits-header info
cval : float, optional
Used in conjunction with mode 'constant', the value outside
the image boundaries.
Return:
actual number of images created
peak image from series, sum image from series
disttext: multiline info about success
The distortion was adapted from skimage.transform.swirl.
Instead of a swirl transformation a rotation symmetric radial transformation
for converting tangential projection to orthographic projection and/or to
correct lens distorsion described by
r =rp*(1+a3*rp^2 +a5*rp^4)
Other parameters, as used in swirl
----------------
# output_shape : tuple (rows, cols), optional
# Shape of the output image generated. By default the shape of the input
# image is preserved.
order : int, optional
The order of the spline interpolation, default is 1. The order has to
be in the range 0-5. See `skimage.transform.warp` for detail.
0: Nearest-neighbor
1: Bi-linear (default)
2: Bi-quadratic
# 3: Bi-cubic
# 4: Bi-quartic
# 5: Bi-quintic
mode : {'constant', 'edge', 'symmetric', 'reflect', 'wrap'}, optional
Points outside the boundaries of the input are filled according
to the given mode, with 'constant' used as the default.
clip : bool, optional
Whether to clip the output to the range of values of the input image.
This is enabled by default, since higher order interpolation may
produce values outside the given input range.
# preserve_range : bool, optional
# Whether to keep the original range of values. Otherwise, the input
# image is converted according to the conventions of `img_as_float`.
Also see
http://scikit-image.org/docs/dev/user_guide/data_types.html
"""
def _distortion_mapping(xy, center, rotation, a3, a5, yscale=1.0):
"""
the original images are converted to square pixels by scaling y
with factor yscale
if yscale is omitted, square pixels are assumed
Calculate shifted coordinates: xs,ys =x',y' – x0,y0
Calculate r', phi': r' =sqrt(xs^2+ys^2)
phi' =phi = arctan2(ys,xs)
Calculate r: r =r'*(1+a3*(r'/f)^2 +...)
Calculate x,y: x=x0+r*cos(phi)
y= y0 + r*sin(phi)
(Pixel value at x',y': I'(x',y') = I(x,y) in the original image)
"""
x, y = xy.T
x0, y0 = center
y0 = y0 * yscale # the center in the original image has to be scaled as well
# y has been scaled in a previous step with resize image
rp = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
phi = np.arctan2(y - y0, x - x0) + rotation
r = rp * (1 + rp ** 2 * (a3 + a5 * rp ** 2)) # 8sec, 2.9217, 2.906 for single image png
xy[..., 0] = x0 + r * np.cos(phi)
xy[..., 1] = y0 + r * np.sin(phi)
return xy
idg = None
dattim = ''
sta = ''
# scale image
back, header = get_fits_image(backfile)
# notice order of coordinates in rescale
if center is None:
center = np.array(back.shape)[:2][::-1] / 2
warp_args = {'center': center,
'a3': a3,
'a5': a5,
'rotation': rotation,
'yscale': yscale}
# warnings.filterwarnings('ignore') # ignore warnings for cleaner output
with warnings.catch_warnings():
warnings.simplefilter("ignore")
# create list of filenames
image_list = create_file_list(im, nm, start=first)
a = 0
if dist:
if len(back.shape) == 3:
multichannel = True
else:
multichannel = False
ima = tf.rescale(back, (yscale, 1), multichannel=multichannel) # scale sum and peak image start
if debug:
print('imy imx , x00 y00: ', ima.shape, center)
else:
ima = back
imsum = 0 * ima
impeak = imsum
t1 = time.time()
fullmdist = outpath + '/' + mdist
for image in image_list:
if path.exists(image):
a += 1 # create output filename suffix
fileout = fullmdist + str(a)
idist = get_png_image(image, colorflag)
if background:
idist = idist - back # subtract background
# calculate distortion
if dist:
if abs(yscale - 1.0) > 1.0e-3: # scale image if yscale <> 1.0
idist = tf.rescale(idist, (yscale, 1), multichannel=multichannel)
if len(idist.shape) == 3:
for c in [0, 1, 2]: # separate color planes for faster processing
idist2 = idist[:, :, c]
# use bi-quadratic interpolation (order = 2) for reduced fringing
idist2 = tf.warp(idist2, _distortion_mapping, map_args=warp_args, order=2,
mode='constant', cval=cval)
idist[:, :, c] = idist2
else:
idist = tf.warp(idist, _distortion_mapping, map_args=warp_args, order=2,
mode='constant', cval=cval)
write_fits_image(idist, fileout + '.fit', fits_dict, dist=dist)
if show_images:
image_data, im_scale = get_img_filename(fileout + '.fit', opt_dict)
# if idg: window['-D_IMAGE-'].delete_figure(idg)
idg = window['-D_IMAGE-'].draw_image(data=image_data, location=(0, graph_size))
# create sum and peak image
imsum = imsum + idist
file = path.basename(fileout + '.fit')
impeak = np.maximum(impeak, idist)
disttext = f'{file} of {nm} done\n'
window['-RESULT2-'].update(value=disttext, append=True)
window.refresh()
# write sum and peak fit-file
write_fits_image(imsum, fullmdist + '_sum.fit', fits_dict, dist=dist)
write_fits_image(impeak, fullmdist + '_peak.fit', fits_dict, dist=dist)
nmp = a
# print(nmp, ' images processed of ', nm)
logging.info(f'{nmp} images processed of {nm}')
tdist = (time.time() - t1) / nmp
disttext = f'{nmp} images processed of {nm}\n'
if dist:
info = f'process time for single distortion: {tdist:8.2f} sec'
logging.info(info)
disttext += info + '\n'
# print(f'process time background, dark and dist {t2:8.2f} sec')
if 'DATE-OBS' in fits_dict.keys():
dattim = fits_dict['DATE-OBS']
sta = fits_dict['M_STATIO']
else:
logging.info('no fits-header DATE-OBS, M-STATIO')
disttext += '\n!!!no fits-header DATE-OBS, M-STATIO!!!\n'
logging.info(f"'DATE-OBS' = {dattim}")
logging.info(f"'M-STATIO' = {sta}")
info = f'Bobdoubler, start image = {im}{first}'
if int(fits_dict['M_BOB']):
logging.info(f'with ' + info)
else:
logging.info('without ' + info)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
ios.imsave(outpath + '/' + mdist + '_peak.png', np.flipud(impeak * 255).astype(np.uint8))
return a, imsum, impeak, disttext
# -------------------------------------------------------------------
def register_images(start, nim, x0, y0, dx, dy, infile, outfil, window, fits_dict, contr=1, idg=0, show_reg=False):
"""
:param start: index of first image (reference) for registering
:param nim: number of images to register_images
:param x0: x-coordinate of reference pixel (int)
:param y0: y-coordinate of reference pixel (int)
:param dx: half width of selected rectangle
:param dy: half height of selected rectangle
:param infile: full filebase of images e.g. out/mdist
:param outfil: filebase of registered files, e.g. out/mdist
:param window: GUI window for displaying results of registered files
:param fits_dict: content of fits-header
if the procedure stops early, nim = index - start + 1
:return:
index: last processed image
sum_image: !average! of registered images
regtext: multiline text of results
dist: if True, distorted images, else False
outfile: filename of sum-image, e.g. for sum of 20 added images: out/r_add20
fits_dict: updated values of fits-header
"""
def _shift(xy):
return xy - np.array(dxy)[None, :]
index = start
sum_image = []
outfile = ''
dist = False
fits_dict.pop('M_NIM', None) # M_NIM only defined for added images
logging.info(f'start x y, dx dy, file: {x0} {y0},{2 * dx} {2 * dy}, {infile}')
regtext = f'start x y, dx dy, file: {x0} {y0},{2 * dx} {2 * dy}, {infile}' + '\n'
image_list = create_file_list(infile, nim, ext='', start=start)
regtext += f' file peak x y wx wy\n'
try:
for image_file in image_list:
im, header = get_fits_image(image_file)
if 'D_X00' in header.keys():
dist = True
if 'M_BOB' in header.keys():
fits_dict['M_BOB'] = header['M_BOB']
if len(im.shape) == 3:
imbw = np.sum(im, axis=2) # used for _fit_gaussian_2d(data)
data = imbw[y0 - dy:y0 + dy, x0 - dx:x0 + dx]
shifted = im
# selected area
else:
data = im[y0 - dy:y0 + dy, x0 - dx:x0 + dx]
params, success = _fit_gaussian_2d(data)
(height, x, y, width_x, width_y) = params # x and y reversed
width_x = 2 * np.sqrt(2 * np.log(2)) * np.abs(width_x) # FWHM
width_y = 2 * np.sqrt(2 * np.log(2)) * np.abs(width_y) # FWHM
# full image
x = x + y0 - dy # y and x reversed
y = y + x0 - dx
imagename = os.path.basename(image_file)
info = f'{imagename:12s} {height:7.3f} {y:6.1f} {x:6.1f} {width_y:5.2f} {width_x:5.2f}'
regtext += info + '\n'
window['-RESULT3-'].update(regtext)
window.refresh()
logging.info(info)
if index == start: # reference position for register_images
x00 = y
y00 = x
# register_images
dxy = [x00 - y, y00 - x]
if len(im.shape) == 3:
for c in [0, 1, 2]: # separate color planes for faster processing
im2 = im[:, :, c]
coords = tf.warp_coords(_shift, im2.shape)
sh2 = map_coordinates(im2, coords) # / 255
shifted[:, :, c] = sh2
else:
coords = tf.warp_coords(_shift, im.shape)
shifted = map_coordinates(im, coords) # / 255
if index == start: # reference position for register_images
sum_image = shifted
else:
sum_image += shifted
# write image as fit-file
write_fits_image(shifted, outfil + str(index - start + 1) + '.fit', fits_dict, dist=dist)
if show_reg:
image_data, idg, actual_file = draw_scaled_image(outfil + str(index - start + 1) + '.fit',
window['-R_IMAGE-'], opt_dict, idg, contr=contr,
resize=True, tmp_image=True)
window.set_title('Register: ' + str(actual_file))
window.refresh()
# set new start value
x0 = int(y)
y0 = int(x)
index += 1 # next image
index += -1
except:
# Exception, delete last image with error
if path.exists(outfil + str(index - start + 1) + '.fit'):
os.remove(outfil + str(index - start + 1) + '.fit')
index += -1
info = f'problem with register_images, last image: {image_file}, number of images: {index}'
logging.info(info)
regtext += info + '\n'
nim = index - start + 1
if nim > 1:
if index == nim + start - 1:
outfile = outfil + '_add' + str(nim)
sum_image = sum_image / nim # averaging
fits_dict['M_STARTI'] = str(start)
fits_dict['M_NIM'] = str(nim)
write_fits_image(sum_image, outfile + '.fit', fits_dict, dist=dist)
return index, sum_image, regtext, dist, outfile, fits_dict
# -------------------------------------------------------------------
def _gaussian(height, center_x, center_y, width_x, width_y):
"""Returns a _gaussian function with the given parameters"""
width_x = float(width_x)
width_y = float(width_y)
return lambda x, y: height * np.exp(
-(((center_x - x) / width_x) ** 2 + ((center_y - y) / width_y) ** 2) / 2)
# -------------------------------------------------------------------
def _moments(data):
"""Returns (height, x, y, width_x, width_y)
the _gaussian parameters of a 2D distribution by calculating its
_moments """
height = x = y = width_x = width_y = 0.0
total = data.sum()
if total > 0.0:
xx, yy = np.indices(data.shape)
x = (xx * data).sum() / total
y = (yy * data).sum() / total
col = data[:, int(y)]
width_x = np.sqrt(np.abs((np.arange(col.size) - y) ** 2 * col).sum() / col.sum())
row = data[int(x), :]
width_y = np.sqrt(np.abs((np.arange(row.size) - x) ** 2 * row).sum() / row.sum())
height = data.max()
# print('h: %5.1f'%height,'x: %5.1f'%x, 'y: %5.1f'%y, 'wx: %5.1f'%width_x, 'wy: %5.1f'%width_y)
return height, x, y, width_x, width_y
# -------------------------------------------------------------------
def _fit_gaussian_2d(data):
"""Returns (height, x, y, width_x, width_y)
the _gaussian parameters of a 2D distribution found by a fit
(ravel makes 1-dim array)"""
params = _moments(data)
success = 0
if params[0] > 0:
errorfunction = lambda p: np.ravel(_gaussian(*p)(*np.indices(data.shape)) - data)
p, success = optimize.leastsq(errorfunction, params)
return p, success
# -------------------------------------------------------------------
def get_fits_keys(header, fits_dict, res_dict, keyprint=False):
"""
gets fits-header from image and appends or overwrites the current fits-header
dictionary. It also converts the specially coded D_ keys to res_dict keys
and attributes floating type values
updates res_dict
:param header:
:param fits_dict: dictionary of fits-header values
:param res_dict: dictionary of distortion parameters (result of m-calib)
:param keyprint:
:return: updated fits_dict
"""
for key in fits_dict.keys():
if key in header.keys():
fits_dict[key] = header[key]
if keyprint:
print(key, header[key])
for key in res_dict.keys():
fkey = 'D_' + key.upper()
if fkey in header.keys():
res_dict[key] = np.float32(header[fkey])
fits_dict[fkey] = np.float32(header[fkey])
if keyprint:
print(key, fits_dict[fkey])
return fits_dict
# -------------------------------------------------------------------
def get_fits_image(fimage):
"""
reads fits image data and header
fimage: filename with or without extension
converts 32-bit floating values and 16-bit data to Python compatible values
reads also color images and transposes matrix to correct order
(normalizes images to +/- 1 range)
returns: image as np array, header
"""
fimage = change_extension(fimage, '.fit')
im, header = fits.getdata(fimage, header=True)
if int(header['BITPIX']) == -32:
im = np.array(im) / 32767
elif int(header['BITPIX']) == 16:
im = np.array(im)
else:
print(f'unsupported data format BITPIX: {header["BITPIX"]}')
exit()
if len(im.shape) == 3:
im = np.transpose(im, (1, 2, 0))
return im, header
# -------------------------------------------------------------------
def show_fits_image(file, imscale, image_element, contr=1.0, show=True):
"""
not needed at present, left in place for further use
loads fits-image, adjusts contrast and scale and displays in GUI as tmp.png
replaced by draw_scaled_image
:param file: fits-file with extension
:param imscale: scale for displayed image
:param image_element: where to display image in GUI
:param contr: image contrast
:param show: if True, image_element isupdated, otherwise only 'tmp.png' is created
:return:
"""
imbw, header = get_fits_image(file)
if len(imbw.shape) == 2:
im = tf.rescale(imbw, imscale, multichannel=False)
else:
im = tf.rescale(imbw, imscale, multichannel=True)
im = im / np.max(im) * 255 * contr
im = np.clip(im, 0.0, 255)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
ios.imsave('tmp.png', np.flipud(im).astype(np.uint8))
if show:
image_element.update(filename='tmp.png')
return
# -------------------------------------------------------------------
def select_rectangle(infile, start, res_dict, fits_dict, wloc, outfil, maxim):
"""
displays new window with image infile + start + 'fit
a rectangle around the selected line can be selected with dragging the mouse
:param infile: filebase of image
:param start: index of selected image
:param res_dict: dictionary
:param fits_dict: "
:param wloc: location of displayed window for selection
:param outfil:
:param maxim:
:return:
Ok if rectangle selected,
x0, y0: center coordinates of selected rectangle (int)
dx, dy: half width and height of selected rectangle (int)
"""
im, header = get_fits_image(infile + str(start))
im = im / np.max(im)
get_fits_keys(header, fits_dict, res_dict, keyprint=False)
# #===================================================================
# new rect_plt
# first get size of graph from tmp.png and size of image
# graph coordinates are in image pixels!
(imy, imx) = im.shape[:2]
image_file = 'tmp.png' # scaled image
imbw = np.flipud(ios.imread(image_file)) # get shape
(canvasy, canvasx) = imbw.shape[:2]
wlocw = (wloc[0] + 300, wloc[1] + 50)
# check for old files
delete_old_files(outfil, maxim, ext='.fit')
image_elem_sel = [sg.Graph(
canvas_size=(canvasx, canvasy),
graph_bottom_left=(0, 0), # starts at top, set y-scale here
graph_top_right=(imx, imy), # set x-scale here
key='-GRAPH-',
change_submits=True, # mouse click events
drag_submits=True)]
layout_select = [[sg.Text('Start File: ' + infile + str(start), size=(50, 1)), sg.Text(key='info', size=(40, 1)),
sg.Ok(), sg.Cancel()],
image_elem_sel]
# ---------------------------------------------------------------------------
winselect_active = True
winselect = sg.Window(f'select zero order or spectral line',
layout_select, finalize=True, location=wlocw,
keep_on_top=True, no_titlebar=False,
disable_close=False, disable_minimize=True)
# get the graph element for ease of use later
graph = winselect['-GRAPH-'] # type: sg.Graph
graph.draw_image(image_file, location=(0, imy)) if image_file else None
winselect.refresh()
dragging = False
start_point = end_point = prior_rect = None
x0 = y0 = dx = dy = 0
while winselect_active:
event, values = winselect.read()
idg = graph.draw_rectangle((0, 0), (imx, imy), line_color='blue')
if event == "-GRAPH-": # if there's a "Graph" event, then it's a mouse
x, y = (values["-GRAPH-"])
if not dragging:
start_point = (x, y)
dragging = True
else:
end_point = (x, y)
if prior_rect:
graph.delete_figure(prior_rect)
if None not in (start_point, end_point):
prior_rect = graph.draw_rectangle(start_point,
end_point, line_color='red')
elif event is not None and event.endswith('+UP'):
# The drawing has ended because mouse up
xy0 = [int(0.5 * (start_point[0] + end_point[0])),
int(0.5 * (start_point[1] + end_point[1]))]
size = (abs(start_point[0] - end_point[0]),
abs(start_point[1] - end_point[1]))
info = winselect["info"]
info.update(value=f"grabbed rectangle at {xy0} with size {size}")
start_point, end_point = None, None # enable grabbing a new rect
dragging = False
if min(size[0], size[1]) > 1: # rectangle
info.update(value=f"rectangle at {xy0} with size {size}")
x0 = xy0[0]
y0 = xy0[1]
dx = int((size[0] + 1) / 2)
dy = int((size[1] + 1) / 2)
elif event in ('Ok', 'Cancel'):
graph.delete_figure(idg)
winselect_active = False
winselect.close()
return event, x0, y0, dx, dy
# -------------------------------------------------------------------
def add_rows_apply_tilt_slant(outfile, par_dict, res_dict, fits_dict, opt_dict,
contr, wloc, restext, regtext, window):
"""