version 2.3.2.0 (367 KB) by
Yury

Optometrika MATLAB library implements analytical and iterative ray tracing approximation to optical image formation using Snell’s and Fresne

OPTOMETRIKA

Developer: Yury Petrov

Optometrika library implements analytical and iterative ray tracing approximation to optical image formation using Snell’s and Fresnel’s laws of refraction and reflection.

Currently, the library implements refractive and reflective general surfaces, aspheric (conic) surfaces with astigmatism, Fresnel surfaces, cones and cylinders (elliptic too), planes, circular and ring-shaped apertures, rectangular flat screens, spheroidal screens, and a realistic model of the human eye with accommodating lens and spheroidal retina. See example*.m files for examples of ray tracing in general (user-defined shape) lenses, aspheric lenses, Fresnel lenses, prisms, mirrors, and human eye.

The library traces refracted rays, including intensity loss at the refractive surface. Reflected rays are currently traced for mirrors and also for a single total internal reflexion or double refraction, if it happens. Note that the Bench class object is not a real physical bench, it is only an ordered array of optical elements, and it is your responsibility to arrange optical objects in the right order. In particular, if you need to trace rays passing through the same object multiple times, you have to add the object multiple times to the bench array in the order the object is encountered by the rays. For example, double refraction/reflection for cylindrical and conical surfaces can be calculated by adding the surface twice to the bench.

The library is very compact and fast. It was written using Matlab classes and is fully vectorized. It takes about 2 seconds to trace 100,000 rays through an external lens and the human eye (8 optical surfaces) on a 3 GHz Intel Core i7 desktop. Fresnel lens tracing is somewhat slower due to looping through the Fresnel cones describing the lens surface. Tracing through user-defined (general) surfaces is significantly slower due to iterative search of ray intersections with the surface.

Thank you for downloading Optometrika, enjoy it!

————————————————————

List of examples:

example1.m: tests the basic functionality of the Optometrika library

example2.m: demonstrates the Optometrika's optical model of the human eye

example3.m: demonstrates accommodation of the human eye by minimizing the retinal image

example4.m: tests a ring lens with the cosine surface profile defined in coslens.m

example5.m: tests planar mirrors

example6.m: tests planar and parabolic mirrors (a Newtonian refractor telescope)

example7.m: tests a Fresnel lens

example8.m: tests a lens with polynomial aspheric terms

example9.m: tests cone mirrors

example10.m: tests cylinder and cone surfaces with double refraction

example11.m: demonstrates ray tracing for rays originating inside the human eye

example12.m: draws a lens and determines its front surface, back surface, and total height. Makes an animated gif of the lens and an engineering drawing of the lens.

example13.m: tests refraction through the lens edge and backward rays refraction (sub-aperture Maksutov-Cassegrain telescope)

example14.m: tests refraction through a lens with astigmatism (different vertical and horizontal radii of curvature)

example15.m: simulates a hexagonal array of spherical micro lenses

example16.m: demonstrates STL export of various lenses

Yury (2021). Optometrika (https://github.com/caiuspetronius/Optometrika), GitHub. Retrieved .

Created with
R2018b

Compatible with any release

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Michele MarroHello Yury Petrov,

Just a question: eye model is validated?

Alexander SHello Lucas, Noah and Uwe,

The original code did not implement consecutive rotation. On this branch you can find a patch which implements continuos rodrigues rotations.

https://github.com/alexschultze/Optometrika/tree/rotation_fix

Regards, Alex

Alaster McDonachLiked it and it did what I needed. Looking to use it more.

Lucas FaberHello Yuri,

great programm, it really helps me a lot.

I ran into a problem while applying serveral rotations around different axis to a plane and noticed that severall people had the same problem(Noah Birge,Uwe Apel).

The rotation is not visible in the .draw() function.

Is the program designed to do this? And do you have a workaround for this issue?

Myrtle42 proposed a solution on 22 Aug 2016, can I implement this without disturbing any other calculations?

Regards Lucas

Shuiqin Zhengvery very very beautiful code!! Not only useful in raytracing but also coding!

Alexander SGood job Yury. Thanks for starting this easy to understand ray tracer.

Joseph SmalleyHi Yury,

When you say it takes 2 seconds to trace 100,000 rays through 8 surfaces, are you referring to 100,000 source rays, or a total of 100,000 ray segments? Approximately how many times slower is ray tracing through user-defined surfaces?

Tanner FaderoHi Yury,

I love using Optometrika, thanks so much for putting so much time into writing it!

I work with light sheet fluorescence microscopy, where I often need to ray trace through plano-convex cylindrical lenses. Your CylinderLens.m classdef is very useful, but I was wondering if there was a way to introduce an argument that would "crop" (i.e. introduce a bounding box or maximum outer diameter) the CylinderLens.m object in the bench container, similar to how you can specify an "outer diameter" argument for the Lens.m object. This would help me with figure generation, since the cylindrical lens wouldn't be "intersecting" the rays at different points in the optical path (even though they don't refract through the lens at that point).

Thanks again!

Tanner

YuryWilfred Mason, Optometrika is not my main line of work, so I didn't get to writing a tutorial. But I feel that if you understand how Matlab works in general (e.g., understand its data structures an know how to debug a program) you should find it fairly easy to understand Optometrika by following the provided examples. Any time you are in a doubt about what is going on just put a breakpoint and dive into the Optometrika function that you want to understand. Otherwise use "help function_name" to read its description. I tried to put enough information there.

Wilfred MasonCan be very useful but the lack of documentation makes it difficult to use. Please consider incorporating a documentation package.

Wilfred MasonAwesome package! Does anyone know how to modify the input rays so that they have a Gaussian distribution. I see the functions, but I don't know how to incorporate them.

Johannes StoerkleTristan DavisHello,

I am having trouble getting this to work with freemat? Any chance someone knows which parts of these would hold up the code? It has to be formatting of some kind. In a general sense, is there specific items that people have ran into with this?

Nick StuartAlaster McDonachEnjoyed using it.

Jin YangNice job! I was wondering is that easy to modify your codes to trace light rays pass through non-spherical (parametric function is given) transparent lens ?

Noah BirgeHi Yury,

I'm having some trouble making my aperture rotate.

a=Aperture(positionVector, rectangularDimensions);

a.rotate([0 1 0],pi/2);

bench.append(a)

^ when I use "a.display()" in the command line, I get the correct information, but bench.draw() shows the aperture in the un-rotated position. The rays seem to act as though the aperture is in the correct location, but the visualization of the aperture location seems to be different. Am I doing something incorrectly?

Thanks!

Bradley FrankStefi AnguDear Yury, I have a question on how to get number of rays in a desired coordinates/screen to represent the light intensity.Thank you.

Aarón Cruz RamírezYuryStefi, reflections are implemented, see examples 5, 6, 9, 13. You can build a retroreflector, e.g., a corner reflector or a refractive sphere using the existing objects like Planes arranged into a 3D corner.

Stefi AnguDear Yury, I'm trying to simulate a confocal system which requires Retro-reflection function in Optometrika. Is there anyway i can capture the reversed reflection because i cant seem to make it work.

Yuryf d, export_stl() works with Plane surfaces in Optometrika 2.2. Also, the internal diameter (hole) was added to Plane. Thanks for pointing this out!

YuryKonstantina, beamsplitters are not implemented. A possible workaround is to run the same beam of rays through a bench twice: the first time with a mirror, the second time - without. Then you can combine the two sets of Rays using Rays.append() function and send the resulting bundle to a screen or elsewhere. Just note that no phase information is currently taken into account in Optometrika, it is a purely geometric optics approximation.

YuryUwe Apel, you can use the Retina object as a spherical screen.

YuryWick, thanks a lot for catching this bad bug! This problem was fixed in Optometrika 2.2.

YuryPeter Cheimets, the bug affecting backward rays and hyperbolic surfaces was fixed in Optometrika 2.2. Thanks for finding it!

Peter CheimetsI'm trying to simulate a Cassegrain telescope. I started with example 13 (Maksutov-Cassegrain), removed the lenses and tweaked the element positions so they focus. Then I put in a -1 conic constant for both mirrors. This worked fine. However, in a Cassegrain, the secondary has a <-1 conic constant. When I changed the conic constant to less than -1 (even as small a change as going to -1.00000001), the beams coming off the secondary no longer converged. Am I missing something, or is this incapable of dealing with conic constants less than -1? (I downloaded my version in April.

Thanks Peter

Kashmira NakhodaHi Yuri,

Thank you for developing this toolbox, it is extremely powerful.

I am trying to use the Plane element in your package to transform it as a coded aperture of 0 and 1 and add a source that is an extended object.

Have you ever tried introducing such elements before?

KonstantinaHi Yury! Are beamsplitters implemented in the optical elements and if not, can they be implemented? Thank you for your answer! Konstantina

Stefan SchuberthIs there a solution to Wick's problem?

WickI have found a bug in using an aspheric lens with rays.intersection. Consider this snippet of code:

clearvars

R = 50; k = -1.3; a4 = 5e-7;

bench = Bench;

bench.append(AsphericLens([0 0 0],50,R,k,[0 a4],{ 'air' 'n-bk7'}));

bench.append(Screen( [20 0 0 ], 50, 50, 512, 512 ));

rays_in = Rays( 2, 'collimated', [ -5 0 0 ], [ 1 0 0 ], 2*8.789, 'linear', 'air', 532e-9); % just two rays in the XY plane at +/- 8.789 above the centerline

rays_out = bench.trace( rays_in );

pass_angle_relative_to_centerline = atand((rays_out(3).r(:,2) - rays_out(2).r(:,2)) ./ (rays_out(3).r(:,1) - rays_out(2).r(:,1))) % arctan (deltaY / deltaX)

The angle relative to the centerline after passing into the lens surface at +/- 8.789 is -/+ 6.4 degrees. However, at that height, a horizontal ray has an angle of incidence of 10 degrees. For n = 1.519 glass, the transmission angle relative to the internal surface normal is 6.56 degrees. So the transmitted angle relative to the centerline should be -/+ 3.44 degrees.

The points where the rays intersect the aspheric surface appear to be correct so the problem must be with the calculation of the angle.

I see you had a check for the Optimization toolbox in that section of the code. Disabling the Optimization toolbox did not change the result. Thank you in advance for taking another look at your aspheric surface calculations.

WickThanks for the submission. I'm just now getting started. Is there any built-in feature to determine wavefront error (optical path difference, OPD) shy of me evaluating the optical path length for each ray? Thanks.

f dDear Yury, I encounter a problem, the front lens is aspheric, and the back len is planar, when I use export_stl() to export, cause error below,

error: export_stl (line 79)

p2 = asphlens( Y2 * cs, Y2 * sn, pars2, 0 );

no variable pars2.

I found is the problem of planar back len whose D is scalar not vector, how to solve this question?

Peter CheimetsI downloaded the program last week and have been playing with it. It looks very powerful. Has anyone written up some notes on how to use it? I have been digging into the software itself, and looking at the examples, but I was hoping for something more manual-like.

Thanks

Peter

Zhaozhong ChenUwe ApelHi Yury - After switching from Matlab2013b to Matlab2018 now almost all works fine. However, I would like to use a spherical screen, but I wasn't able to find any example how the default setting R=inf could be changed successfully. Could you kindly advise? Many thanks.

Ann OmanHi Yury - What are some of the anticipated use cases for the eye model? We are curious to test focusing on near and far objects and what is seen by the eye as well as potential near and farsighted cases.

Is this possible?

Pratul NijhawanThanks a lot for such a wonderful library and so much respect for the amount of help you did to the society. Is the 3D tracing of rays through an array of some points, say 3*3 matrix of points in a 3D space and obtaining the final outputs as zones on image screen possible with this library, say a set of some rays are passing point by point to lens's 3D surface, then the obtained image on screen would be in form of zones, so I needed some suggestions regarding that. In the designed library, which examples I should refer for the same?

Kind regards!

YuryDear 123, run "help coslens" to see a description of how to define a freeform surface to be used with Optometrika. It doesn't have to be described by a formula, like in example4, you just need a function that provides the lens surface's height and normal at any given (x,y) point in the lens plane. It could be obtained by spline interpolation of the surface data from a look-up table, for example.

xuYu XiaoZhijie Ma123Hello Yury,

This program is very usefull. Can this program use for freeform lens with the coordinates (X, Y, Z) that is produced by the user? Not just use for the 'GeneralLens' 'AsphericLens' 'FresnelLens', etc. I have sent an email to you in this web from the email Contact. I hope you will see it.

Thank you very much.

123123Hello Yury,

There is bug in example16.m when use matlab 2018b. It does not work.

When the example is used in freeform lens which only has coordinate information, how to modify it? In other words, how to use them to simulate freeform lens and reflectors?

Thank you very much.

Uwe ApelHello Yury,

Just another question: All example solutions that use function "quadric_intersection" do not work, failure is mismatch of matrix dimensions at both last formulas:

t1 = 0.5 * ( -b - sqrt( D ) ) ./ a;

t2 = 0.5 * ( -b + sqrt( D ) ) ./ a;

r1 = r0 + n .* t1;

r2 = r0 + n .* t2;

t1, t2 are only one-dim vectors while r0 and n are 127x3 matrices. I'm working with Matlab13b,

might this be a version issue? Any hint? Thanks.

Uwe ApelHelly Yury,

Very useful approach for visualization of optical test benches. Some more basic information would be appreciated.

I observed some not-understood behaviour. So I modified example5 by rotating the 2nd mirror additionally a small amount around the y-axis. The course of the rays behaves in plausible way, but the orientation of the displayed mirror plane is completely different:

mirror2 = Plane( [ 60 50 0 ], 40, 40, { 'mirror' 'air' } ); % pay attention to the glass order here!

mirror2.rotate( [ 0 0 1 ], -pi / 4 );

mirror2.rotate( [ 0 1 0 ], -pi / 18 );

bench.append( mirror2 );

Any idea?

Jongmin KimThis is very good example!!

Joshua BoneWhat is the purpose of the attenuation property of the Ray object?

Stephan SuckowThis was very helpful for me.

MPOlivier Ripollyou may want to replace the for loop code adding the polynomials in the aspheric.m by a call to polyval, or rewrite it using Horner's algorithm to make it faster and more reliable (polyval uses the Horner algorithm)

Hugo FerreiraHello Yury. First of all, great work with this toolbox! I read the comments here and it seems there was a grating simulation tool in a previous version. Why was it removed, and can you provide me with that tool? Thanks!

Edgar BruckeEdvin AtteboR SHey Yury - thanks for the quick reply. I'm on R2016a so that might not be the issue, but I'll keep an eye out for your next update. Thanks!

YuryIf you have problems similar to the one described by Justin, R S, and Pascal, please update your Matlab to 2016 or later. I'll fix this in the next update, so it will work for earlier Matlab versions.

R SI have the same problem as Pascal & Justin, any help would be appreciated. Thanks!

Pascal BertoHello Yury,

I have the same problem than Justin Little...

Please, do you know where it could come from?

Justin LittleHello Yury,

I get the following error when attempting to trace rays in example1 (and other examples). I am using Matlab 2015B. Any thoughts?

Error using ./

Matrix dimensions must agree.

Error in Rays/intersection (line 712)

nrms = nrms ./ sqrt( sum( nrms.^2, 2 ) );

Error in Rays/interaction (line 732)

[ rays_out, nrms ] = self.intersection( surf );

Error in Bench/trace (line 327)

rays( i + 1 ) = rays( i ).interaction( self.elem{ i }, out_fl );

Error in example1 (line 42)

rays_through = bench.trace( rays_in ); % trace rays

Evgeniy AndreevEvgeniy AndreevHello, Yury

First I want to say thank you for the great toolbox! It really helped me in studying and scientific work.

What do you think about adding the project to GitHub? I'd like to contribute to it, like many others I believe.

Thanks!

Tobias SimmetdsdeocsHello Yury,

I am trying to add an additional rotation to mirror1 in example5:

(ZIP file downloaded on 9/21/2017)

% planar mirror

mirror1 = Plane( [ 60 0 0 ], 40, 40, { 'air' 'mirror' } ); % pay attention to the glass order here, it defines the mirror orientation!

mirror1.rotate( [ 0 0 1 ], -pi / 4 );

mirror1.rotate( [ 0 1 0 ], 0.00000001 ); % add this line to example5.m

bench.append( mirror1 );

The rendered image draws mirror1 at an unexpected angle.

What is the correct way to add additional surface rotations?

Thanks

dsdeocsMarie LuHi Yury,

great software! What is the default unit for lengths ?

Thanks a lot!

Vinoin Devpaul VincelyHello everyone,

I am having a bit of trouble trying to grasp how to use the toolbox. Is there a simple user manual that I can use to better understand the toolbox?

Thanks,

Philip AndersonHello Yury,

The code Just what I was looking for to get solid angle sensitivity of a photosensor buried in clear plastic.

One question: I wish to trace rays that go through AND past a hemisphere lens, so I have been using

bench.trace( rays_in, 0 );

to allow missed rays t be traced. Works fine if I use a "conelens" (but this is not what I need) , but not with "lens" (which I use to build the 1/2 sphere.

Looking in Rays, the line

rays_out.I( in, : ) = self.I( in );

or its equivalent is not included for "lens", (as it is for "conelens" and others) so the intensity is not passed on.

Is there a fix for this please (not trivial because there is no obvious indices 'in', Otherwise the code is very close to ideal for me.

Thanks

Daniel AharoniExcellent optics package! I think there is an error in the focal point calculation in Rays.m. I think f = osr + +d * nav; should actually be f = osr + rav +d * nav;

Noam CohenGreat.

Could be even greater if will be expanded to NIR/SWIR/MWIR/LWIR.

Sotirios NousiasSotirios NousiasIbrahimGreat toolbox! Thank you, Yury!

Question for anyone; does anyone know how to add a custom lens like a Powell lens? I was able to generate one with a different script, but the rays don't interact with it when I incorporate it into an optical system. Thanks in advance!

Myrtle42In case anyone runs into the same problem, I corrected it by editing the rotate method in Surface.m:

Comment out lines:

41-43

Add in:

% Calculate new normal vector to the plane

self.n = rodrigues_rot(self.n, rot_axis, rot_angle);

% Calculate rotation between original and final vectors

rotvec = vrrotvec([1 0 0],self.n);

% Assign correct rotation axis and angle

self.rotax = rotvec(1:3);

self.rotang = rotvec(4);

I didn't check lines 44-49 are still correct since I don't rotate more than 180 deg

Clarification welcome if there is a better way to apply multiple rotations in 3d

Myrtle42How to create a prism? Can't find it in the examples. Thanks!

Max LowndesI'm not really clear on how to use this! Would love some help as it seems to do exactly what I need. Trying to run a ray tracing simulation of a fiber optic!

xizhi guHi,Yury, Thanks for the toolbox!

I have a question that why the "grating" is removed from the newest version? Because I can't find it anymore..

Robert ShahWhats the default length unit, mm?

Evgeniy AndreevValeriyYury, your work is impressive. Is it possible to treat influence on the image of diffraction grating placed between object and imaging lens?

juzimmVery impressive and useful, thanks for sharing!

Marshall BremerThank you. This is exceptional.

DarinYury, Thanks for sharing this AWESOME tool.

I was sitting down to write something for doing the same job... and would not have done it nearly this well. A nice education in OOP.

Juan PerezHi, this is a very useful and versatile tool, congratulations and thanks for sharing it! I am trying to get the coordinates of the intersection of the rays with the retina (in example2.m for example), can you give me some idea of how to do this?

Thanks!

MatthewAwesome ideas and code. Documentation could be much less sparse.

YuryIf you want to use a material which is not in the refrind.m file the easiest is to edit the file. Just add another 'case', i.e., case 'quartz' nref = [ 1.4585 1.4631 1.4564 ];

Just make sure that the order of the refractive indices corresponds to wavelengths of 587.6, 486.1, and 656.3 nm.

Gerald RothenhoeferHi, my first impression is very good. Great work!

But, it takes some round about to figure out which glass materials are available and which properties are assumed.

The links given under doc(refrindx) seem to be dead.

So, the only way to get the information is by edit(refrindx).

Would be nice to make the refractive index an input argument to the lens function (as an alternative to 'glass').

GR

AdamRichard YoungerBartExcellent program. I am wondering if there is a solution to create two source rays in one model. Tried to fix this by using two aperatures, however they block each other. Keep up the good work.

NachoExcellent. It would be nice to implement cylindrical lenses

IanExcellent and makes wonderful use of OOP, using a virtual optical bench object.