-
Notifications
You must be signed in to change notification settings - Fork 1
/
omnifocuslib.py
2263 lines (2040 loc) · 89.4 KB
/
omnifocuslib.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
# -*- coding: utf-8 -*-
#-------------------------------------------------------------------------------
# Name: omnifocuslib.py
# Purpose: helper functions for simulation of omnifocus image synthesis
# using Zemax and PyZDDE. In particular, the results presented in
# the paper: "Omnifocus image synthesis using lens swivel,"
# I. Sinharoy, P. Rangarajan, and M. Christensen, in Imaging and
# Applied Optics 2016, OSA.
#
# Author: Indranil Sinharoy, Southern Methodist university, Dallas, TX.
#
# Copyright: (c) Indranil Sinharoy 2015, 2016
# License: MIT License
#-------------------------------------------------------------------------------
'''utility functions for simulation of omnifocus image synthesis using Zemax and
PyZDDE. Several functions are included for geometric optics computations,
automating the simulation, storing and retrieving image stack data into and
from tagged HDF5 format, etc.
Assumed wavelength for tracing rays in Zemax is 0.55 mu
'''
from __future__ import division, print_function
import os
import sys
import numpy as np
import cv2
import collections as co
import pyzdde.zdde as pyz
import pyzdde.arraytrace as at
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colorbar import make_axes
from IPython.core import display
import ipywidgets as widgets
from ipywidgets import interactive, interact, fixed
import h5py as hdf
import time as time
from scipy.misc import imsave
# global variable
_initTime = 0.0
#%% File, directory management and other housekeeping utility functions
def get_directory_path(dirbranch=None):
'''returns absolute path of the leaf directory
If directory branch is not present, the function creates the directories under
current file directory
Parameters
----------
dirbranch : tuple or None
tuple of strings representing the directory branch. If `None`
the current directory is returned
Returns
-------
dirpath : string
absolute path to the directory
Example
-------
>>> get_directory_path(['data', 'imgstack'])
C:\Projects\project_edof\data\imgstack
'''
wdir = os.path.dirname(os.path.realpath(__file__))
if dirbranch:
dirpath = os.path.join(wdir, *dirbranch)
if not os.path.exists(dirpath):
os.makedirs(dirpath)
return dirpath
else:
return wdir
def set_hdf5_attribs(hdf5Obj, attribDict):
"""helper function to set attributes of h5py objects
Parameters
----------
hdf5Obj : HDF5 object
group or dataset object (including an hdf5 file)
attribDict : dict
attribute dict
"""
for key, value in attribDict.items():
hdf5Obj.attrs[key] = value
def save_to_IMAdir(arr, imagename):
"""helper function to save a numpy ndarray as image into the Zemax's Image
directory
Parameters
----------
arr : ndarray
ndarray
imagename : string
image filename, e.g. 'gridofdots.png'
Returns
-------
None
"""
usr = os.path.expandvars("%userprofile%")
IMAdir = os.path.join(usr, 'Documents\Zemax\IMAFiles')
filename = os.path.join(IMAdir, imagename)
imsave(name=filename, arr=arr)
def get_time(startClock=False):
"""returns elapsed time
The first time this function is called, `startClock` should be True.
The function doesn't print the time for the initiating call.
Returns
-------
time : string
elapsed time from the previous call
"""
global _initTime
if startClock:
_initTime = time.clock()
else:
_cputime = time.clock() - _initTime
m, s = divmod(_cputime, 60)
h, m = divmod(m, 60)
return "%d:%02d:%02d" % (h, m, s)
#%% Helper functions for plotting data
def show_pixel_grid(ax, pixX, pixY, gAlpha=None):
"""plots pixel grid on the given axes
Parameters
----------
ax : axes object
figure axes object
pixX : integer
number of pixels along x/columns
pixY : integer
number of pixels along y/rows
gAlpha : float (0 < gAlpha <= 1)
alpha for the pixel grid
Returns
-------
None
"""
gAlpha = gAlpha if gAlpha is not None else 1.0
xl = np.arange(1, pixX) - 0.5
yl = np.arange(1, pixY) - 0.5
ax.xaxis.set_minor_locator(plt.FixedLocator(xl))
ax.yaxis.set_minor_locator(plt.FixedLocator(yl))
ax.xaxis.set_tick_params(which='minor', length=0)
ax.yaxis.set_tick_params(which='minor', length=0)
ax.xaxis.set_tick_params(which='major', direction='out')
ax.yaxis.set_tick_params(which='major', direction='out')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
ax.grid(which='minor', color='r', linestyle='-', linewidth=0.75,
alpha=gAlpha)
def show_around(img_data, pixX=20, pixY=None, ctrX=None, ctrY=None,
pixSize=None, pixGrid=False, retData=False, ax=None):
"""display pixX x pixY pixels around (ctrX, ctrY).
Parameters
----------
img_data : list of list
data returned by ln.zGetImageSimulation()
pixX : integer
number of pixels along x
pixY : integer, optional
number of pxiels along y. If `None` or 0, `pixY` is assumed equal
to `pixX`
ctrX : integer, optional
center pixel along x
ctrY : integer, optional
center pixel along y
pixSize : real, optional
length of the side of a square pixel.
if passed, a second axes is shown to show the physical dimension
pixGrid : bool, optional
if True, show pixel grid
retData : bool, optional
if True, the cropped data is returned
and not plotted
ax : axes object, optional
if an axes object is passed the plot will be rendered
in the given axes and `plt.show()` will have to be invoked
by the caller. Use `0` if you want to use this function
with ipywidgets' interactive function.
Assumptions
------------
pixX, pixY, ctrX, ctrY are not 0
"""
img = np.array(img_data, dtype='uint8')
M, N = img.shape[:2]
ctrX, ctrY = ctrX or N/2, ctrY or M/2
pixX = pixX or 20
pixY = pixY if pixY else pixX
sRow, eRow = int(ctrY - pixY/2), int(ctrY + pixY/2)
sCol, eCol = int(ctrX - pixX/2), int(ctrX + pixX/2)
# ensure bounds are right
sRow = sRow if sRow >= 0 else 0
eRow = eRow if eRow <= M else M
sCol = sCol if sCol >= 0 else 0
eCol = eCol if eCol <= N else N
data = img[sRow:eRow, sCol:eCol, :]
if retData:
return data
else:
if not ax:
fig, ax = plt.subplots(1, 1, figsize=(5.5, 5.5))
iAlpha = 0.8 if pixGrid else 1.0
ax.imshow(data, interpolation='none', zorder=15, alpha=iAlpha)
ax.set_xlabel('pixels', fontsize=10)
ax.set_ylabel('pixels', fontsize=10)
if pixGrid:
show_pixel_grid(ax=ax, pixX=pixX, pixY=pixY)
if pixSize:
# create a secondary axes
ax2 = ax.twiny()
ax3 = ax2.twinx()
wideX = pixSize*(eCol-sCol)
wideY = pixSize*(eRow-sRow)
ax.xaxis.set_ticks_position('bottom') # required
ax.yaxis.set_ticks_position('left') # required
ax2.set_xlim(0, wideX)
ax3.set_ylim(wideY, 0)
ax2.set_xlabel('phy. dim.', fontsize=10)
ax3.set_ylabel('phy. dim.', fontsize=10)
ws = 'Width: {:2.3f}'.format(wideX)
hs = 'Height: {:2.3f}'.format(wideY)
ax.text(x=0.99, y=0.055, s=ws, transform = ax.transAxes,
color='w', zorder=16, alpha=0.65, horizontalalignment='right')
ax.text(x=0.99, y=0.02, s=hs, transform = ax.transAxes,
color='w', zorder=16, alpha=0.65, horizontalalignment='right')
if not ax:
ax.set_aspect('equal')
plt.show()
def imshow(image, fig=None, axes=None, subplot=None, interpol=None,
xlabel=None, ylabel=None, figsize=None, cmap=None):
"""Rudimentary image display routine, for quick display of images without
the spines and ticks
"""
if (subplot == None):
subplot = int(111)
if(fig==None):
if figsize and isinstance(figsize, tuple) and len(figsize)==2:
fig = plt.figure(figsize=figsize)
else:
fig = plt.figure()
axes = fig.add_subplot(subplot)
elif(axes==None):
axes = fig.add_subplot(subplot)
# plot the image
if len(image.shape) > 2:
imPtHandle = plt.imshow(image, interpolation=interpol)
else:
cmap = cmap if cmap is not None else cm.gray
imPtHandle = plt.imshow(image, cmap=cmap, interpolation=interpol)
# get the image height and width to set the axis limits
try:
pix_height, pix_width = image.shape
except:
pix_height, pix_width, _ = image.shape
# Set the xlim and ylim to constrain the plot
axes.set_xlim(0, pix_width-1)
axes.set_ylim(pix_height-1, 0)
# Set the xlabel and ylable if provided
if(xlabel != None):
axes.set_xlabel(xlabel)
if(ylabel != None):
axes.set_ylabel(ylabel)
# Make the ticks to empty list
axes.xaxis.set_ticks([])
axes.yaxis.set_ticks([])
return imPtHandle, fig, axes
def get_imlist(filePath, itype='jpeg'):
"""returns a list of filenames for all images of specified type in a
directory
Parameters
----------
filePath : string
full path name of the directory to be searched
itype : string, optional
type of images to be searched, for example -- 'jpeg', 'tiff', 'png',
'dng', 'bmp' (without the dot(.))
Returns
-------
imageFiles : list of strings
list of image filenames with full path.
"""
imlist = []
opJoin = os.path.join
dirList = os.listdir(filePath)
if itype in ['jpeg', 'jpg']:
extensions = ['.jpg', '.jpeg', '.jpe',]
elif itype in ['tiff', 'tif']:
extensions = ['.tiff', '.tif']
else:
extensions = [''.join(['.', itype.lower()]), ]
for ext in extensions:
imlist += [opJoin(filePath, f) for f in dirList if f.lower().endswith(ext)]
return imlist
#%% Basic geometric-optics functions
def gaussian_lens_formula(u=None, v=None, f=None, infinity=10e20):
"""return the third value of the Gaussian lens formula, given any two
Parameters
----------
u : float, optional
object distance from first principal plane.
v : float, optional
image distance from rear principal plane
f : float, optional
focal length
infinity : float
numerical value to represent infinity (default=10e20)
Returns
-------
glfParams : namedtuple
named tuple containing the Gaussian Lens Formula parameters
Notes
-----
1. Both object and image distances are considered positive.
2. Assumes frontoparallel configuration
Examples
--------
>>> gaussian_lens_formula(u=30, v=None, f=10)
glfParams(u=30, v=15.0, f=10)
>>> gaussian_lens_formula(u=30, v=15)
glfParams(u=30, v=15, f=10.0)
>>> gaussian_lens_formula(u=1e20, f=10)
glfParams(u=1e+20, v=10.0, f=10)
"""
glfParams = co.namedtuple('glfParams', ['u', 'v', 'f'])
def unknown_distance(knownDistance, f):
try:
unknownDistance = (knownDistance * f)/(knownDistance - f)
except ZeroDivisionError:
unknownDistance = infinity
return unknownDistance
def unknown_f(u, v):
return (u*v)/(u+v)
if sum(i is None for i in [u, v, f]) > 1:
raise ValueError('At most only one parameter can be None')
if f is None:
if not u or not v:
raise ValueError('f cannot be determined from input')
else:
f = unknown_f(u, v)
else:
if u is None:
u = unknown_distance(v, f)
else:
v = unknown_distance(u, f)
return glfParams(u, v, f)
def pupil_centric_lens_formula(u=None, v=None, f=None, mp=1, infinity=10e20):
"""return the third value of the pupil centric lens formula, given any two
Parameters
----------
u : float, optional
object distance from entrance pupil position
v : float, optional
image (geometrical focus) distance from exit pupil position
f : float, optional
focal length
mp : float, optional
pupil magnification
infinity : float
numerical value to represent infinity (default=10e20)
Returns
-------
pclfParams : namedtuple
named tuple containing the Pupil Centric Lens Formula parameters
Notes
-----
1. Both object and image distances are considered positive.
2. Assumes frontoparallel configuration
3. The pupil magnification (mp) is the ratio of the paraxial exit pupil size
to the paraxial entrance pupil size
Examples
--------
>>> pupil_centric_lens_formula(u=1016.0, v=None, f=24)
pclfParams(u=1016.0, v=24.580645161290324, f=24.0)
>>> pupil_centric_lens_formula(u=1016.0, v=24.58064516129, f=None)
pclfParams(u=1016.0, v=24.58064516129, f=23.999999999999694)
>>> pupil_centric_lens_formula(u=504.0, v=None, f=24.0, mp=2.0)
pclfParams(u=504.0, v=49.170731707317074, f=24.0)
>>> pupil_centric_lens_formula(u=535.6319, v=None, f=24.0, mp=0.55)
pclfParams(u=535.6319, v=14.370742328796812, f=24.0)
"""
pclfParams = co.namedtuple('pclfParams', ['u', 'v', 'f'])
if sum(i is None for i in [u, v, f]) > 1:
raise ValueError('At most only one parameter out of u, v & f can be None')
if f is None:
if not u or not v:
raise ValueError('f cannot be determined from input')
else:
f = (mp*u*v)/(mp**2 * u + v)
else:
if u is None:
try:
u = (1.0/mp)*(v*f)/(v - mp*f)
except ZeroDivisionError:
u = infinity
else:
try:
v = (mp**2 * f * u)/(mp * u - f)
except ZeroDivisionError:
v = infinity
return pclfParams(u, v, f)
#%% Helper functions to draw the cardinal and pupil planes
def set_surface_semidia(ln, surf, value=0):
"""set the surface `surf` semi-diameter to value.
This requires a special function because we need to first set
the solve on the semi-diameter to 'fixed', else Zemax will automatically
modify it.
Parameters
----------
ln : object
pyzdde object
surf : integer
surface number
Returns
-------
None
"""
ln.zSetSolve(surf, ln.SOLVE_SPAR_SEMIDIA, ln.SOLVE_SEMIDIA_FIXED)
ln.zSetSurfaceData(surfNum=surf, code=ln.SDAT_SEMIDIA, value=value)
def insert_dummy_surface(ln, surf, thickness=0, semidia=0, comment='dummy'):
"""helper function to insert dummy surface
Parameters
----------
ln : object
pyzdde object
surf : integer
surface number at which to insert the dummy surface
thickness : real, optional
thickness of the surface
semidia : real, optional
semi diameter of the surface
comment : string, optional
comment on the surface
"""
ln.zInsertSurface(surf)
ln.zSetSurfaceData(surf, ln.SDAT_COMMENT, comment)
ln.zSetSurfaceData(surf, ln.SDAT_THICK, thickness)
set_surface_semidia(ln, surf, semidia)
def draw_plane(ln, space='img', dist=0, surfName=None, semiDia=None):
"""function to draw planes at the points specified by "dist"
Parameters
----------
ln : pyzdde object
active link object
space : string (`img` or `obj`), optional
image space or object space in which the plane is specified. 'img' for
image space, 'obj' for object space. This info is required because
Zemax returns distances that are measured w.r.t. surface 1 (@LDE) in
object space, and w.r.t. IMG in image space. See the Assumptions.
dist : float, optional
distance along the optical axis of the plane from surface 2 (@LDE) if
`space` is `obj` else from the IMG surface. This assumes that surface 1
is a dummy surface
surfName : string, optional
name to identify the surf in the LDE, added to the comments column
semiDia : real, optional
semi-diameter of the surface to set
Returns
-------
None
Assumptions (important to read)
-------------------------------
The function assumes (for the purpose of this study) that surface 1 @ LDE is
a dummy surface at certain distance preceding the first actual lens surface.
This enables the rays entering the lens to be visible in the Zemax layout
plots even if the object is at infinity. So the function inserts the planes
(and their associated dummy surfaces) beginning at surface 2.
"""
numSurf = ln.zGetNumSurf()
inSurfPos = numSurf if space=='img' else 2 # assuming that the first surface will be a dummy surface
insert_dummy_surface(ln, inSurfPos, dist, 0, 'dummy')
ln.zInsertSurface(inSurfPos+1)
ln.zSetSurfaceData(inSurfPos+1, ln.SDAT_COMMENT, surfName)
if semiDia:
set_surface_semidia(ln, inSurfPos+1, semiDia)
thickSolve, pickupSolve = 1, 5
frmSurf, scale, offset, col = inSurfPos, -1, 0, 0
ln.zSetSolve(inSurfPos+1, thickSolve, pickupSolve, frmSurf, scale, offset, col)
def get_cardinal_points(ln):
"""Returns the distances of the cardinal points (along the optical axis)
Parameters
----------
ln : object
PyZDDE object
Returns
-------
fpObj : float
distance of object side focal point from surface 1 in the LDE,
irrespective of which surface is defined as the global reference
fpImg : float
distance of image side focal point from IMG surface
ppObj : float
distance of the object side principal plane from surface 1 in the
LDE, irrespective of which surface is defined as the global
reference surface
ppImg : float
distance of the image side principal plane from IMG
Notes
-----
1. The data is consistant with the cardinal data in the Prescription file
in which, the object side data is with respect to the first surface in the LDE.
2. If there are more than one wavelength, then the distances are averaged.
"""
zmxdir = os.path.split(ln.zGetFile())[0]
textFileName = os.path.join(zmxdir, "tmp.txt")
sysProp = ln.zGetSystem()
numSurf = sysProp.numSurf
ln.zGetTextFile(textFileName, 'Pre', "None", 0)
line_list = pyz._readLinesFromFile(pyz._openFile(textFileName))
ppObj, ppImg, fpObj, fpImg = 0.0, 0.0, 0.0, 0.0
count = 0
for line_num, line in enumerate(line_list):
# Extract the Focal plane distances
if "Focal Planes" in line:
fpObj += float(line.split()[3])
fpImg += float(line.split()[4])
# Extract the Principal plane distances.
if "Principal Planes" in line and "Anti" not in line:
ppObj += float(line.split()[3])
ppImg += float(line.split()[4])
count +=1 #Increment (wavelength) counter for averaging
# Calculate the average (for all wavelengths) of the principal plane distances
# This is only there for extracting a single point ... ideally the design
# should have just one wavelength define!
if count > 0:
fpObj = fpObj/count
fpImg = fpImg/count
ppObj = ppObj/count
ppImg = ppImg/count
# # Delete the temporary file
pyz._deleteFile(textFileName)
cardinals = co.namedtuple('cardinals', ['Fo', 'Fi', 'Ho', 'Hi'])
return cardinals(fpObj, fpImg, ppObj, ppImg)
def draw_pupil_cardinal_planes(ln, firstDummySurfOff=40, cardinalSemiDia=1.2, push=True, printInfo=True):
"""Insert paraxial pupil and cardinal planes surfaces in the LDE for rendering in
layout plots.
Parameters
----------
ln : object
pyzdde object
firstDummySurfOff : float, optional
the thickness of the first dummy surface. This first dummy surface is
inserted by this function. See Notes.
cardinalSemiDia : float, optional
semidiameter of the cardinal surfaces. (Default=1.2)
push : bool, optional
push lens in the DDE server to the LDE
printInfo : bool, optional
if True (default), the function also prints information about the
locations of the cardinal planes and pupils
Assumptions
-----------
The function assumes that the lens is already focused appropriately,
for either finite or infinite conjugate imaging.
Notes
-----
1. 'first dummy surface' is a dummy surface in LDE position 1 (between the
OBJ and the actual first lens surface) whose function is show the input
rays to the left of the first optical surface.
2. The cardinal and pupil planes are drawn using standard surfaces in the LDE.
To ensure that the ray-tracing engine does not treat these surfaces as real
surfaces, we need to instruct Zemax to "ignore" rays to these surfaces.
Unfortunately, we cannot do it programmatically. So, after the planes have
been drawn, we need to manually do the following:
1. 2D Layout settings
a. Set number of rays to 1 or as needed
2. For the pupil (ENPP and EXPP) and cardinal surfaces (H, H', F, F'),
and the dummy surfaces (except for the dummy surface named "dummy 2
c rays" go to "Surface Properties" >> Draw tab
a. Select "Skip rays to this surface"
3. Set field points to be symmetric about the optical axis
3. For clarity, the semi-diameters of the dummy sufaces are set to zero.
"""
ln.zSetWave(0, 1, 1)
ln.zSetWave(1, 0.55, 1)
# insert dummy surface at 1 for showing the input ray
ln.zRemoveVariables()
# before inserting surface check to see if the object is at finite
# distance. If the object is at finite distance, inserting a dummy
# surface with finite thickness will change the image plane distance.
# so first decrease the thickness of the object surface by the
# thickness of the dummy surface
objDist = ln.zGetSurfaceData(surfNum=0, code=ln.SDAT_THICK)
assert firstDummySurfOff < objDist, ("dummy surf. thick ({}) must be < "
"than obj dist ({})!".format(firstDummySurfOff, objDist))
if objDist < 1.0E+10:
ln.zSetSurfaceData(surfNum=0, code=ln.SDAT_THICK, value=objDist - firstDummySurfOff)
insert_dummy_surface(ln, surf=1, thickness=firstDummySurfOff, semidia=0, comment='dummy 2 c rays')
ln.zGetUpdate()
# Draw Exit and Entrance pupil planes
expp = ln.zGetPupil().EXPP
draw_plane(ln, 'img', expp, "EXPP")
enpp = ln.zGetPupil().ENPP
draw_plane(ln, 'obj', enpp - firstDummySurfOff, "ENPP")
# Get and draw the Principal planes
fpObj, fpImg, ppObj, ppImg = get_cardinal_points(ln)
draw_plane(ln,'img', fpImg, "F'", cardinalSemiDia)
draw_plane(ln,'obj', fpObj - firstDummySurfOff, "F", cardinalSemiDia)
draw_plane(ln,'img', ppImg, "H'", cardinalSemiDia)
draw_plane(ln,'obj', ppObj - firstDummySurfOff, "H", cardinalSemiDia)
# Check the validity of the distances
ppObjToEnpp = ppObj - enpp
ppImgToExpp = ppImg - expp
focal = ln.zGetFirst().EFL
v = gaussian_lens_formula(u=ppObjToEnpp, v=None, f=focal).v
ppObjTofpObj = ppObj - fpObj
ppImgTofpImg = ppImg - fpImg
if printInfo:
print("Textual information about the planes:\n")
print("Exit pupil distance from IMG:", expp)
print("Entrance pupil from Surf 1 @ LDE:", enpp)
print("Focal plane obj F from surf 1 @ LDE: ", fpObj, "\nFocal plane img F' from IMA: ", fpImg)
print("Principal plane obj H from surf 1 @ LDE: ", ppObj, "\nPrincipal plane img H' from IMA: ", ppImg)
print("Focal length: ", focal)
print("Principal plane H to ENPP: ", ppObjToEnpp)
print("Principal plane H' to EXPP: ", ppImgToExpp)
print("Principal plane H' to EXPP (abs.) calc. using lens equ.: ", abs(v))
print("Principal plane H' to rear focal plane: ", ppObjTofpObj)
print("Principal plane H to front focal plane: ", ppImgTofpImg)
print(("""\nCheck "Skip rays to this surface" under "Draw Tab" of the """
"""surface property for the dummy and cardinal plane surfaces. """
"""See Docstring Notes for details."""))
if push:
ln.zPushLens(1)
def insert_cbs_to_tilt_lens(ln, lastSurf, firstSurf=2, pivot='ENPP', offset=0, push=True):
"""function to insert appropriate coordinate break and dummy surfaces
in the LDE for tilting the lens about a pivot.
Parameters
----------
ln : object
pyzdde object
lastSurf : integer
the last surface that the tilt coordinate-break should include. This
surface is normally the image side principal plane surface H'
firstSurf : integer, optional, default=2
the first surface which the tilt coordinate-break should include.
Generally, this surface is 2 (the dummy surface preceding the
object side principal plane H)
pivot : string, optional
indicate the surface about which to rotate the lens group. Currently only
ENPP (with an offset from the ENPP) has been implemented
offset : real, optional
offset (in lens units) to offset the actual pivot point from the `pivot`
push : bool
push lens in the DDE server to the LDE
Returns
-------
cbs : 2-tuple of integers
the first and the second coordinate breaks. Also, if `push` is True,
then the LDE will be updated
Notes
-----
1. A dummy surface will be inserted at Surface 2 to move the CB to the
pivot point.
2. Check "skip rays to this surface" for the first and last dummy surfaces
inserted by this function.
3. The layout will display all the surfaces. This function is needed in spite
of the `ln.zTiltDecenterElements()` function because it is designed to tilt
all the the cardinal and associated dummy surfaces in between that may appear
before or after the actual lens surfaces in the LDE. Also, the pivot point is
generally not about the lens surface.
4. Use this function to rotate the lens group about a pivot point (PIVOT). To
rotate the image plane use a CB infront of the IMA surface.
Assumptions (weak)
-----------------
The following assumptions need not be strictly followed.
1. Surface 0 is the object surface which may or may not be at infinity
2. Surface 1 is a dummy surface for seeing ray visibility.
"""
ln.zRemoveVariables()
gRefSurf = ln.zGetSystem().globalRefSurf
assert ln.zGetSurfaceData(gRefSurf, ln.SDAT_THICK) < 10e10, 'Change global ref'
if pivot=='ENPP':
# estimate the distance from the firstSurf to ENPP
enpp = ln.zGetPupil().ENPP # distance of entrance pupil from surface 1
enppInGRef = ln.zOperandValue('GLCZ', 1) + enpp
firstSurfInGRef = ln.zOperandValue('GLCZ', firstSurf)
firstSurfToPivot = enppInGRef - firstSurfInGRef + offset
cbFirstSurf = firstSurf + 1
# insert dummy surface to move to the pivot position where the CB will be applied
cmtstr = 'Move to pivot (offset from ENPP)' if offset else 'Move to ENPP'
insert_dummy_surface(ln, surf=firstSurf, thickness=firstSurfToPivot,
semidia=0, comment=cmtstr)
lastSurf +=1
# insert dummy surface to show pivot ...
insert_dummy_surface(ln, surf=firstSurf+1, thickness=0, semidia=1.0, comment='PIVOT')
cbFirstSurf +=1
lastSurf +=1
# insert coordinate breaks
cb1, cb2, _ = ln.zTiltDecenterElements(firstSurf=cbFirstSurf, lastSurf=lastSurf,
cbComment1='Lens tilt CB',
cbComment2='Lens restore CB',
dummySemiDiaToZero=True)
# set solve on cb1 surface to move back to firstSurf
ln.zSetSolve(cb1, ln.SOLVE_SPAR_THICK, ln.SOLVE_THICK_PICKUP,
firstSurf, -1.0, 0.0, 0)
else:
raise NotImplementedError("Option not Implemented.")
if push:
ln.zPushLens(1)
return (cb1, cb2)
def show_grid_distortion(n=11):
"""Plots intersection points of an array of real chief rays and paraxial
chief rays from the object field with the image surface.
Rays are traced on the lens in the LDE (and not in the DDE server). The
paraxial chief ray intersection points forms the reference field. Comparing
the two fields gives a sense of the geometric distortion produced by the lens.
Parameters
----------
n : integer, optional
number of field points along each x- and y-axis. Total number
of rays traced equals n**2
Returns
-------
None
A matplotlib scatter plot displays the intersection points, against
the paraxial intersection points of the chief ray on that specified
surface. Therefore, the plot shows the distortion in the system.
Notes
-----
1. The reference field in Zemax's Grid Distortion plot is generated by
tracing a real chief ray very close to the optical axis, and then
scaling appropriately.
2. TODO:
1. Show rays that are vignetted.
"""
xPara, yPara, _, _, _ = get_chief_ray_intersects(n, False)
xReal, yReal, _, err, vig = get_chief_ray_intersects(n, True)
# plot
fig, ax = plt.subplots(1, 1, figsize=(6,6), dpi=120)
ax.scatter(xPara, yPara, marker='o', s=10, facecolors='none', edgecolors='b',
alpha=0.8, zorder=12)
ax.scatter(xReal, yReal, marker='x', s=40, c='r', alpha=0.9, zorder=15)
# these two lines is dependent on the way the hx, hy grid is created in the
# function get_chief_ray_intersects()
xGridPts = xPara[:n]
yGridPts = yPara[::n]
ax.vlines(xGridPts, ymin=min(yGridPts), ymax=max(yGridPts),
zorder=2, colors='#CFCFCF', lw=0.8)
ax.hlines(yGridPts, xmin=min(xGridPts), xmax=max(xGridPts),
zorder=2, colors='#CFCFCF', lw=0.8)
ax.set_aspect('equal')
ax.axis('tight')
plt.show()
def get_chief_ray_intersects(n=11, real=True, surf=-1):
"""Returns the points of intersection of the chief rays from the
object field (normalized) with the specified surface.
Parameters
----------
n : integer, optional
number of field points along each x- and y-axis. Total number
of rays traced equals n**2
real : bool, optional
If `True` (default) then real chief rays are traced,
If `False`, then paraxial chief rays are traced
surf : integer, optional
surface number. -1 (default) indicates image plane
Returns
-------
x : ndarray float32
x-intersects of the chief rays with `surf`
y : ndarray of float32
y-intersects of the chief rays with `surf`
z : ndarray of float32
z-incersects of the chief rays with `surf`
err : ndarray of int16
error code
vig : ndarray of uint8
vignetting code
Notes
-----
TODO:
1. How to handle errors that occurs during the ray tracing.
2. The function seems to work only if the field normalization is set
to "radial". Need to fix it for rectangular field normalization.
"""
# create uniform square grid in normalized field coordinates
nx = np.linspace(-1, 1, n)
hx, hy = np.meshgrid(nx, nx)
hx = hx.flatten().tolist()
hy = hy.flatten().tolist()
# trace the ray
mode = 0 if real else 1
rayData = at.zGetTraceArray(numRays=n**2, hx=hx, hy=hy, mode=mode, surf=surf)
# parse ray traced data
err = np.array(rayData[0], dtype=np.int16)
vig = np.array(rayData[1], dtype=np.uint8)
x = np.array(rayData[2], dtype=np.float32)
y = np.array(rayData[3], dtype=np.float32)
z = np.array(rayData[4], dtype=np.float32)
return x, y, z, err, vig
def plot_chiefray_intersects(ln, cb, tiltXY, pushNewLens=True):
"""plot chief-ray intersects for various rotations of the lens about a pivot point.
Parameters
----------
ln : object
pyzdde object
cb : integer
Element tilt coordinate break (the first CB) surface number
tiltXY : list of 2-tuples
each element of the list is a tuple containing the tilt-about-x
angle, tilt-about-y angle in degrees of the lens about a
pivot point.
e.g. tiltXY = [(0,0), (10, 0), (0, 10)]
Max length = 6
pushNewLens : boolean, optional
see notes
Returns
-------
None : matplotlib figure
Notes
-----
1. This function uses the array ray tracing module functions for generating
the chief-ray--image-plane intersection points. To do so, it needs to push
the lens in the DDE server into the LDE. By default, following the ray
tracing, a new lens is loaded into the LDE. Please make sure that there
are no unsaved/work-in-progress lens files in the LDE before invoking this
function.
2. This function always traces "real" rays. This is because paraxial ray
tracing doesn't seem to be affected by lens tilts, i.e. the ray-intersect
pattern is constant for all rotations of the lens.
"""
tiltAbtXParaNum = 3
tiltAbtYParaNum = 4
cols = ['b', 'r', 'g', 'm', 'c', 'y']
mkrs = ['o', 'x', '+', '1', '2', '3']
n = 11
fig, ax = plt.subplots(1, 1, figsize=(8, 8), dpi=120)
for i, angle in enumerate(tiltXY):
ln.zSetSurfaceParameter(surfNum=cb, param=tiltAbtXParaNum, value=angle[0])
ln.zSetSurfaceParameter(surfNum=cb, param=tiltAbtYParaNum, value=angle[1])
ln.zGetUpdate()
# push lens into the LDE as array tracing occurs in the LDE
ln.zPushLens(1)
x, y, _, err, vig = get_chief_ray_intersects(n=n, real=True)
legTxt = r'$\alpha_x, \alpha_y = {}^o,{}^o$'.format(angle[0], angle[1])
if i:
ax.scatter(x, y, marker=mkrs[i], s=40, c=cols[i], alpha=0.9, zorder=15,
label=legTxt)
else:
ax.scatter(x, y, marker=mkrs[i], s=10, facecolors='none',
edgecolors=cols[i], alpha=0.8, zorder=12, label=legTxt)
xGridPts = x[:n]
yGridPts = y[::n]
ax.vlines(xGridPts, ymin=min(yGridPts), ymax=max(yGridPts),
zorder=2, colors='#CFCFCF', lw=0.8)
ax.hlines(yGridPts, xmin=min(xGridPts), xmax=max(xGridPts),
zorder=2, colors='#CFCFCF', lw=0.8)
# finally load a new lens into the LDE
if pushNewLens:
ln.zNewLens()
ln.zPushLens(1)
ax.set_aspect('equal')
ax.axis('tight')
ax.set_ylabel(r'$\bf{\acute{y}}\,\it{(mm)}$', fontsize=15)
ax.set_xlabel(r'$\bf{\acute{x}}\,\it{(mm)}$', fontsize=15)
ax.legend(fontsize=14, scatterpoints=1, markerscale=1., scatteryoffsets=[0.5], mode='expand',
ncol=len(tiltXY), loc='upper left', bbox_to_anchor=(0.09, 0.965, 0.84, 0.005),
bbox_transform=fig.transFigure, handletextpad=0.5, handlelength=0.9)
plt.show()
#%% Functions for focal / angular sweep to EDOF simulation using Zemax
# Meta-data schema for storing frames in hdf5 file
# ------------------------------------------------
# / (root Group)
# Attribute: zemax file (string)
# Attribute: obj_sharp_focus_distance (real)
# Attribute: lens_tilt_about_x
# Attribute: lens_tilt_about_y
# Attribute: img_sharp_focus_distance (real)
# Attribute: img_tilt_about_x
# Attribute: img_tilt_about_y
# Attribute: img_sim_field_height
# Attribute: img_sim_oversampling
# Attribute: img_sim_wavelength
# Attribute: img_sim_pupil_sampling
# Attribute: img_sim_image_sampling
# Attribute: img_sim_psfx_points
# Attribute: img_sim_psfy_points
# Attribute: img_polarizatoin
# Attribute: img_sim_aberrations
# Attribute: img_sim_relative_illumination
# Attribute: img_sim_fixed_apertures
# Attribute: img_sim_reference
# Attribute: img_sim_pixel_size
# Attribute: img_sim_xpixels
# Attribute: img_sim_ypixels
# attributes for each frame
# 1. defocus_waves (only makes sense for fronto-parallel imaging)
# 2.
def get_image_plane_shifts(nearObj=800, midObj=1000, farObj=1200, fl=24, num=10):
"""Returns the image plane delta shifts in both forward and
backward directions from the base position of the image plane.
The "base position" is the position of the image plane for which
the middle object is in geometrical focus.
Parameters
----------
farObj : real
distance along the optical axis for furthest object
midObj : real
distance along the optical axis for the middle object
nearObj : real
distance along the optical axis for the nearest object
fl : real
effective focal length
num : integer