Simulating Speckle with Mathematica ® Joseph W. Goodman Simulating Speckle with Mathematica ® Joseph W. Goodman SPIE PRESS Bellingham, Washington USA Library of Congress Cataloging-in-Publication Data Names: Goodman, Joseph W., author. Title: Simulating speckle with Mathematica / Joseph W. Goodman, Stanford University. Description: Bellingham, Washington : SPIE, [2022] | Includes bibliographical references. Identifiers: LCCN 2022030726 | ISBN 9781510656543 (paperback) | ISBN 9781510656550 (pdf) Subjects: LCSH: Speckle--Computer simulation. | Speckle--Mathematical models. | Mathematica (Computer file) Classification: LCC QC427.8.S64 G659 2022 | DDC 535/.470113--dc23/eng20220919 LC record available at https://lccn.loc.gov/2022030726 Published by SPIE P.O. Box 10 Bellingham, Washington 98227-0010 USA Phone: +1 360.676.3290 Fax: +1 360.647.1445 Email: [email protected] Web: www.spie.org Copyright © 2022 Society of Photo-Optical Instrumentation Engineers (SPIE) All rights reserved. No part of this publication may be reproduced or distributed in any form or by any means without written permission of the publisher. The content of this book reflects the work and thought of the author. Every effort has been made to publish reliable and accurate information herein, but the publisher is not responsible for the validity of the information or for any outcomes resulting from reliance thereon. Wolfram Mathematica and Mathematica are registered trademarks. Wolfram Language, MathLM, Computable Document Format, Wolfram Desktop, Wolfram|One, and Wolfram Cloud are trademarks of Wolfram Research, Inc. Wolfram|Alpha is a registered trademark of Wolfram Alpha LLC. The cover image was simulated from the painting “The Audience” by Gene Walch. Printed in the United States of America. First printing 2022. For updates to this book, visit http://spie.org and type “PM355” in the search field. Simulating Speckle with Mathematica® Joseph W . Goodman Stanford University Table of Contents 1. Introduction 1.1. Mathematica Background 1.2. Speckle Background 1.3. Methods for Simulating Speckle 2. First-Order Statistics of Speckle Amplitude 2.1. Speckle as a Sum of Independent Random Complex Phasors 2.2.Amplitude Statistics of the Sum of Many Random Phasors with Unit Lengths and Random Phases 2.3.Amplitude Statistics of the Sum of a Large Number of Unit-Length Random Phasors and One Large Phasor 2.4.Amplitude Statistics of the Sum of a Small Number of Unit-Length Random Phasors 3. First-Order Statistics of Speckle Intensity 3.1. Intensity Statistics of the Sum of Many Random Phasors with Unit Lengths and Random Phases 3.2. Intensity Statistics of the Sum of a Large Number of Unit-Length Random Phasors and One Large Phasor 3.3. Intensity Statistics of the Sum of a Small Number of Unit-Length Random Phasors 3.4. Intensity Statistics for Sums of Independent Speckle Patterns 3.5. Intensity Statistics of Partially Developed Speckle 4. Simulation of Speckle in Optical Imaging 4.1. Generation of a Discrete DiffuserArray with Correlated Phases 4.2. Speckle in an Imaging Geometry 4.3. TheAutocorrelation Function of Speckle Intensity vi Table of Contents 4.4. The Power Spectrum of Speckle Intensity 4.5. Effect of Speckle on Resolution 4.6. Speckle in Color Images 5. Simulation of Speckle in Free-Space Propagation 5.1. The Diffuser 5.2. The Fresnel Transfer FunctionApproach 5.3. The Fresnel Transform Approach 6. Speckle at Low Light Levels 6.1. Photocount Image of a Uniform Intensity 6.2. The Negative-Binomial Distribution 6.3. Photocount Image for a Speckle Pattern with Uniform Statistics 6.4. Photocount Image for a Diffuse Structured Object 7. Speckle Phase Vortices 7.1. Generating the Speckle Pattern 7.2. Finding the Zeros of Intensity 7.3. Phase Behavior in the Vicinity of the Zeros of Intensity 8. Polarization Speckle 8.1. The Polarization Ellipse and the Degree of Polarization 8.2. Generating the Two Polarization Components 8.3. Generating a Filtered Speckle Pattern 8.4. Generating the Polarization Ellipses 8.5. Visualizing Polarization Speckle 9. Speckle Simulation for Metrology 9.1. Measurement of In-Plane Displacement 9.2. Electronic Speckle Pattern Interferometry 9.3. Phase-Shifting Speckle Interferometry Appendix A – Some Subtleties in Speckle Simulation With the 4f System A.1. Effects on the Speckle Contrast A.2. Simulation With a Smoothed Phase A.3. Simulation With an Unsmoothed Phase Appendix B – Some Subtleties in Dealing With Mathematica Images B.1. Dimensions of DataArrays and Images B.2. Images When the Data Range Exceeds (0, 1) Table of Contents vii B.3. Effect of Using ImageAdjust[] on an Image B.4.Arrays With Bipolar Values B.5.A Problem Encountered When Starting With an Image Acknowledgements References 1. Introduction Thespeckle phenome nonisubi quitousi nman yfields ofsciencea ndt echnology. Specklep heno menacanbeseen inma nydifferen t im aging modalities,includingacousticalimaging(e.g.,medicalultrasound)andmicrowaveimaging(e.g.,synthetic-apertureradarimaging).Thisbook focusesonsimulatingopticalspecklewithMathematica®,butthesamemethodsusedcaninmanycasesbeappliedtootherimagingmodalities. The r eade r may wonder w hy Mathema tica has be en chosen as the soft ware packa ge f or this b ook.There a re sev eral reasons f or this c hoice. First,andmos timp ortant, Mat hematicaall ows the intersper sin gof bothc ontinuous anddis cretec alculationsun der one umbrel la.Secon d,using Mathematic a, text, c ode, and illus trations ca n be included i n the same docum ent. Th ird, using Mat hem ati ca we c an crea te d yna mic fi gur es, param etersof whichtheus ercan cha ngea twill.H owe ver, such man ipu lationc an notbe perform edinth epr intedv ersion ofthebook, so we haveavoid edm anipulablefigures i nw hatfo llows and replac ed themby arr ays ofstaticfig ure s. L astly, this autho rl ovesM ath ematica for its flexibil ity a nd c ompre hensivene ss. It seems the re a re a lm ost an infini te set of c apabi lities of t he pr ogram, many of w hich lie hidde n for t he novi ceus erb utw hichgrad uallyarereve al edas the use ofthepro graminc reases. Thisbook ha sbeenwri ttene ntir elyinMa them atica .Itcanbe read with the fu llprogram M ath ematic ao rw iththefree progr amWolframPlayeravailablefordownloadfromtheWolfram site.T h eMa thema t- icafile s fo rallchapte rs can befo undatth efollowing U RL:http ://spie. org/S amples/Press boo k_Suppl emental/P M3 55_su p.zip This book is meantas aco mp aniontoth ebo ok Speck le Phe nomenainOptics:TheoryandApplications,2ndEdition,publishedbySPIEPress(Ref.[1]).An extensivelistofreferencescanbefoundinthatbook. 1.1 Mathematica Background Forsomeback grou ndt ha twillbehe lpfultotheread erne wto Math em ati ca, seeR ef.[2].T hereare man ybooksth atdescri be thecapa bilitiesof Math ematica.Ref .[3] ise spe ciallycom preh ensive .O urg oalhereist on otonlypresen tm ethod sforsim ula tingspec klein variou ssituatio ns and applicati ons, but also to introduc e t he reader to the capabilities of Mathematica. In what follows we present some salient features of Mathematicathatwillhelpthenovicegetstarted. 1. B uilt-in funct ionso fM athematica alwaysbeginwithupper-caseletters.Ifthecommandconsistsofaconcatenationoftwoseparatewords, bo t hw ordsm us tbegi ninupper- casel etters. 2. In order to avoid confusion with built-in commands, user-defined functions should usually begin with a lower-case letter. However, it is permissibletobegincommandswithanupper-caseletterinafontthatisdifferentthanSource Code Pro,theusualfontforMathematica input.Asanexample,isnotthesameasN. 3. Math ema ticaredu ces aninputtha tisarat ionaln um bertot h eequiva lentsim plest ra tionalnu mber p ossi ble. Th eratio 32 /6yields1 6 /3, butdoesnotproduceanapproximatedecimalresult.To obtainadecimalresult,placeadecimalpointattheendofeitherthenumeratororthe denominator.Thus,32./6yieldstheresult5.33333. 4. The arguments of functions must be enclosed in square brackets [ ]. Curly brackets { } are reserved for containing lists, which can represent vectors and matrices, and are also used for variable ranges in, for example, plot commands. Parentheses ( ) are reserved for groupingmathematics. 5. When a fun ction is first defined, e ac hinde pen dentvari ab leont he lefto fthedefinition mustbefollowedbyanunderbar_.Underbarsare notusedontherightofthedefinition,orwhenthefunctioniscalledbylatercode. 6. Thesymbol:=isadelayedequality,forwhichtheoperationofequalityismadeonlywhenthesymbolontheleftiscalledinlaterprogram- ming. 7. The symbol %, when used as the argument of a function, represents the last previous output, which is now used as the input to this new function.%%representsthesecond-to-lastoutput,etc. 8. Oftencertainsymbolsmaybeusedmorethanonceinaprogram,withtheirmeaninghavingchangedbetweenuses.Ifhisthesymbol,then thecommandClear[h]erasesthedefinition ofthesymbolhandallowsanewdefinition tobemade.To clearmorethanonesymbol,for example, h and g, use Clear[h,g]. To clear all symbols and functions that have been previously defined, use the command Clear["Global`*"]. 9. The sy mbol ; at the e nd o f an ex press ion to be evaluated indicates that the evaluation should take place, but the result of the evaluation shouldbesuppresseduntilitisneededlaterintheprogram. 10. The symbol I (or)isu sed inMathemati cafor the imagin ar yco nstant - 1.Th esy mbolD[] in dicatesa deriv ati veo perat ion . 11 . T o executea com mand in Math ema tica ,pl acethe cur sorinthesamecell(cellsareindicatedbyverticallinesontherightofthenotebook) asthecommand,andpresstheshiftkeyandthereturnkeysimultaneously. 12. Thesymbols&&meanalogicaland;thesymbols||meanalogicalor. 2 Chapter 1 13. Multiplicationcanbeindicatedbyeither*orbyaspacebetweenthetwoquantitiestobemultiplied. 14. When dealingwithacomplexnumberc,thecommandAbs[c]yieldsthemagnitudeofc,andArg[c]givesthephaseangleofcin radians. 15. Ma thematic ain put co m mand sar ealw ays writt en inSou rc e Co de Pro sem i-boldf ont. Intheev en tth ata Ma the maticacom mand extends beyo ndtheen dof aline, you mayfindthat a hyp henhasb eeninse rtedinthe co mm and.Th ishy phe nis no tpartoft hec omm and.W hen viewing this mat erialw ith either theMath ematic aortheWolfram Playerprograms,ifyouchangethesizeofthewindow, thelinebreakswill change,withtheresultthatthesehyphensshoulddisappear. At the start we d efine a cons iderable numbe r of init ializ ation c omma nds that dete rmine cer tai n parame ters, su ch as th e s ize of images, and certainop eratorst ha twewill d escr ibe when we first use them.T he se i nitial izatio nce llsarei na hidden cel latthe start ofe achM athem at ica notebo okcha pter.InC hapte r1yo uca nopenthisc ellb ycli ckingo nthegray linet hat ap pears jus tabov eth echa ptertit le,an dthe ng oing to the menu above, clic ki ng Cell, the n C ell Proper tie s, and then O pen . Howev er, the re is no need to ope n that c ell unle ss y ou wish to s ee it.The initi aliza tio nce llisrepe at ed att hebegin nin go fea chchaptertha tfoll ow ss othatyo ucana cce ss anych apterdirectlyandevaluatetheinitializa- tioncellinthatchapter.Itisnotpossibletoseetheinitializationcellsintheprintedversionofthebook. The astute rea der may notic e tha t som etimes we bun dle severa l M ath emati ca co mm ands in on e in put cell, wh ile in oth er c ases we list the comm andsin separ atece lls .Th e sam ecomm an ds areexecu tedinthesam eor der inbot h cases, but some tim esbre aking on ecellinto se veral helpsreducewhitespaceinthePDFversionofthechapters.MathematicawillnotbreakasinglecellanditsoutputwhenitpaginatesthePDF version. Fourse tsofuser-d efined M a them aticaco m mandsareus edma ny timesin th isbo ok.T hefirst arec omman dsforforwardandinversediscrete Fouriertransforms(DFTs)inanynumberofdimensions,whichareappliedtolists,i.e.,discretedataarrays: dft[arg_]:=Fourierarg, FourierParameters{1,-1} idft[arg_]:=InverseFourierarg, FourierParameters {1,-1}. There are cor responding use r-defin ed commands for Fourier transforms and inverse Fourier transforms of two-dimensional (2D) continuous functions(bycontinuous,wemeannotdiscrete): ft2[arg_]:= FullSimplifyFourierTransformarg,{x,y}, {u,v},FourierParameters 0,-2*Pi ift2[arg_] := FullSimplify InverseFourierTransformarg,{u,v}, {x,y},FourierParameters 0,-2*Pi. Anadditionalpairofcommandsconvertsa2Drealnumericalarraytoimageformat:Image[array],andacommandtoconverta2Dimage (inimageformat)backtonumericaldata:ImageData[image].Hopefully, thereaderwillrememberandunderstandthesecommandswhen theyareencounteredinthechaptersthatfollow.SeeAppendixBforadiscussionofsomesubtletiesoftheselattertwocommands. Lastlyaresomedefinitionsoffunctionsthatareusedanumberoftimesinthechaptersthatfollow.Theseinclude: rect[x] = UnitBox[x] whichisunityforx< 1 and0otherwise, 2 circ[r]= UnitBox[r - 1/2] whichisacircularlysymmetricfunctionthatisunityforradiuslessthan 1 andzerootherwise,and 2 gaus[x]= Exp[-π*x^2]. ForthereaderwhoisusingeithertheMathematicaorWolframPlayertoreadthenotebooks,theauthorrecommendsthatyousetthemagnifica- tion(intheupperright-handcornerofthenotebook)to150%foreasierreading.IntheMathematicanotebooks,theTable ofContentsprovides linkstoeachofthechapters. Theauthor wis hestoe mphasizeth ath eisno tan “exper t” M athemat icaus er,b uthasm aster edthe prog ram wellenoughto per form almost any calcula tion hew ishes .Hope fully, t here ad erw ill arr iveat as imilarlevelafterreadingthisbook.ButwithMathematica,onecanalwaysdive deeperandlearnnewtricks.That’spartofthefunofusingtheprogram. 1.2 Speckle Background Speckle ari ses fu ndamen tally wheneve r a sum of randomly phased complex qua ntities (i .e., p ha so rs, which a re qua ntiti es h aving bot h a magnitude a nd a p hase) are summ ed tog ether. In co her ent imag ing ap plicati ons, spec kl e can be a si gnifi cant no ise. A s we shall se e, the fluctuat ion s ofthe ima gein tensitycanbecomparableinstrengthtothemeanintensityateachpointintheimage.Foramoredetailedback- groundonspeckle,seeRef.[1]. Itisimport antt odistingui sh betw een theampli tudea nd the intensity ofas um ofphasors be cau seth ese twoqua ntitiesha vediff eren tstatistic al properties.Theamplitudeofthesumgenerallyreferstotheabsolutevalueormagnitudeofthesumofcomplexphasors,whiletheintensityof