## Reference calculation of light propagation between parallel planes of different sizes and sampling rates |

Optics Express, Vol. 19, Issue 1, pp. 32-39 (2011)

http://dx.doi.org/10.1364/OE.19.000032

Acrobat PDF (1631 KB)

### Abstract

The article deals with a method of calculation of off-axis light propagation between parallel planes using discretization of the Rayleigh-Sommerfeld integral and its implementation by fast convolution. It analyses zero-padding in case of different plane sizes. In case of memory restrictions, it suggests splitting the calculation into tiles and shows that splitting leads to a faster calculation when plane sizes are a lot different. Next, it suggests how to calculate propagation in case of different sampling rates by splitting planes into interleaved tiles and shows this to be faster than zero-padding and direct calculation. Neither the speedup nor memory-saving method decreases accuracy; the aim of the proposed method is to provide reference data that can be compared to the results of faster and less precise methods.

© 2010 OSA

## 1. Introduction

2. E. Lalor, “Conditions for the validity of the angular spectrum of plane waves,” J. Opt. Soc. Am. **58**(9), 1235–1237 (1968). [CrossRef]

3. L. Onural, “Exact analysis of the effects of sampling of the scalar diffraction field,” J. Opt. Soc. Am. A **24**(2), 359–367 (2007). [CrossRef]

4. L. Onural, A. Gotchev, H. Ozaktas, and E. Stoykova, “A survey of signal processing problems and tools in holographic three-dimensional television,” IEEE Trans. Circ. Syst. Video Tech. **17**(11), 1631–1646 (2007). [CrossRef]

## 2. One-dimensional case

*Source*and

*Target*.

*Target*(e.g. a camera sensor) is affected by the shining of all points

*Source*(e.g. a SLM, see Fig. 1a ). As mentioned in the introduction, we will assume scalar approximation of coherent light, i.e. the light source can be completely described by amplitude

*A*and phase

*ϕ*, or complex amplitude

*p*. Light propagation is space invariant; therefore, it is just the mutual position of points on

*Source*and

*Target*that matters, not their absolute positions. It follows that the calculation will have a convolution form:In the equation there is no distance along

*z*axis between

*Source*and

*Target*because it is a constant and can be a part of the function

*p*.

*Source*and

*Target*into

*M*and

*N*basic elements

*s*is zero-padded to the size

*C*). Notice that the important values of

*C*is an arbitrary number bigger than or equal to

*C*, we can speed up the computation significantly, becausewhere

*t*,

*s*and

*p*are

*C*-dimensional vectors (arrays) of complex numbers,

12. M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proc. IEEE **93**(2), 216–231 (2005) (Special issue on “Program Generation, Optimization, and Platform Adaptation”). [CrossRef]

*C*is important because the actual calculation time is highly sensitive to its character.

*Target*is twice as fine as the sampling rate of

*Source*. Then the subset of even samples from

*Target*has the same sampling rate as

*Source*. Consequently, we can easily calculate the propagation of

*Source*to even samples of

*Target*. Obviously, we can do the same with odd samples. It means we can split the calculation of light propagation into two calculations; they are denoted in Fig. 1c by black and magenta arrows. It follows that the same principle can be applied when the sampling rate of

*Target*is

*τ*-times finer than the sampling rate of

*Source*, where

*τ*is a natural number; we have to split

*Target*into

*τ*“

*interleaved tiles*”, i.e. the calculation has to be split into

*τ*calculations.

*Source*has

*σ*-times finer sampling than

*Target*. The idea can be generalized: if the sampling rates of

*Source*and

*Target*are in a ratio

*σ*and

*τ*are coprime natural numbers (i.e.

*σ*elements of

*Source*have the same size as

*τ*elements of

*Target*), we can split the calculation into

*σ*-times calculate the propagation for the sampling rate ratio

*σ*and

*τ*that approximate the desired sampling fairly well.

*s*and

*t*have to be padded to size

*C*. This means that for

*M*and

*N*we can expect a lack of memory very soon, especially when using special hardware for

*Source*to an extended detector

*Target*could require hundreds of gigabytes.

*Target*in the middle into two parts (let us call them “

*common tiles*”) with

*Source*to

*Target*and sum the results. As in the “different sampling” case, even this idea can be generalized:

*Source*can be split into

*S*parts,

*Target*into

*T*parts and then we have to calculate

## 3. Two-dimensional case

*Source*to a part of

*Target*, it is necessary to carefully calculate the convolution kernel, i.e. 2D array

*p*. A practical aid is the fact that its element

*Source*to element

*Target*. The equations for sampling rates, samples counts and offsets are simple but technically demanding, so we will not show them here.

2. E. Lalor, “Conditions for the validity of the angular spectrum of plane waves,” J. Opt. Soc. Am. **58**(9), 1235–1237 (1968). [CrossRef]

*λ*is a wavelength), and

*r*is the distance between points

*x*,

*y*(

*z*is a constant). Alternatively, we can assume that a sample of

*Source*in fact expresses – using some pixel spread function – the shining of a particular non-zero-area element of

*Source*and alter the kernel accordingly. If we take non-zero area of

*Target*(sensor) elements into account, we can pre-filter the kernel. If

*Source*is a mathematical model of a real display, we can even measure the kernel. The proposed method therefore does not depend on particular features of the kernel; its only assumption is the description of light propagation as a convolution. In our implementation, we did not deal with advanced methods of kernel calculation, however, and we have chosen plain sampling.

## 4. Theoretical time of computation

*Source*and

*Target*have different sampling rates. Let us assume that

*Target*has

*τ*-times finer sampling than

*Source*, i.e.

*Source*must have the same sampling rate as

*Target*, i.e. it must have

*Target*into

*Source*to a part of

*Target*of size

*τ*for

*M*leads to a similar graph.

*Target*has more samples than

*Source*. That led to saving nearly one-third of all

*Source*and

*Target*of the same sampling rate into common tiles in a complicated ratio leads rather to slowdown if we do not mind memory restrictions.

## 5. Implementation notes

*Source*and

*Target*sampling rates, and to split them into interleaved tiles of common sampling. In case we have enough memory to calculate light propagations between these tiles, we can calculate them (but see later). In another case, we are facing a still unresolved problem.

*C*are available. The interleaved tiles created by splitting

*Source*and

*Target*due to the requirement of common sampling have to be split into

**DFT**of special array sizes up to 6 × faster than for comparable prime sizes. This is the first expected reason for scattered look of the measured data in Fig. 2. One can see in Fig. 2b that measured speedups are bigger than theoretical ones. This behavior starts to appear for

*Source*side size

*τ*-times, so

**DFT**has to work with a big array. Using interleaved tiling,

**DFT**works with smaller arrays that fit into cache memory more easily, and therefore bigger speedup can be expected.

*i*. From these times we have picked “friendly-size” ones that satisfy condition

*friendly-size*side sizes will be “friendly” too. Next, we have measured calculation times for those arrays. The propagation itself is then calculated in 2D “friendly-size” arrays. In case we need to split into common tiles due to memory restrictions, we choose tiles of maximum width (tiles are then in fact horizontal stripes); height of

*Source*and

*Target*tiles is chosen to be approximately equal and as big as possible. We have compared this heuristic to the optimal solution found using brute force; it suggests at most approximately

*Source*to the

*Target*of the same size (for several different sizes) using different tiling schemes. We have found that the relative difference of the calculated complex amplitudes is in the order of 10

^{−10}or smaller; the difference was calculated as the absolute value of complex amplitudes differences. This difference is so small that we have not tried to find where it comes from.

## 6. Conclusion

## Acknowledgments

## References and links

1. | J. W. Goodman, |

2. | E. Lalor, “Conditions for the validity of the angular spectrum of plane waves,” J. Opt. Soc. Am. |

3. | L. Onural, “Exact analysis of the effects of sampling of the scalar diffraction field,” J. Opt. Soc. Am. A |

4. | L. Onural, A. Gotchev, H. Ozaktas, and E. Stoykova, “A survey of signal processing problems and tools in holographic three-dimensional television,” IEEE Trans. Circ. Syst. Video Tech. |

5. | V. Katkovnik, J. Astola, and K. Egiazarian, “Discrete diffraction transform for propagation, reconstruction, and design of wavefield distributions,” Appl. Opt. |

6. | N. Delen and B. Hooker, “Free-space beam propagation between arbitrarily oriented planes based on full diffraction theory: a fast Fourier transform approach,” J. Opt. Soc. Am. A |

7. | J.-L. Kaiser, E. Quertemont, and R. Chevallier, “Light propagation in the pseudo-paraxial fresnel approximation,” Opt. Commun. |

8. | E. Sziklas and A. Siegman, “Diffraction calculations using fast Fourier transform methods,” Proc. IEEE |

9. | R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, “Shifted Fresnel diffraction for computational holography,” Opt. Express |

10. | F. Zhang, G. Pedrini, and W. Osten, “Reconstruction algorithm for high-numerical-aperture holograms with diffraction-limited resolution,” Opt. Lett. |

11. | K. Matsushima and S. Nakahara, “Extremely high-definition full-parallax computer-generated hologram created by the polygon-based method,” Appl. Opt. |

12. | M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proc. IEEE |

**OCIS Codes**

(070.2025) Fourier optics and signal processing : Discrete optical signal processing

(070.7345) Fourier optics and signal processing : Wave propagation

**ToC Category:**

Fourier Optics and Signal Processing

**History**

Original Manuscript: November 2, 2010

Revised Manuscript: December 10, 2010

Manuscript Accepted: December 12, 2010

Published: December 20, 2010

**Citation**

Petr Lobaz, "Reference calculation of light propagation between parallel planes of different sizes and sampling rates," Opt. Express **19**, 32-39 (2011)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-1-32

Sort: Year | Journal | Reset

### References

- J. W. Goodman, Introduction to Fourier Optics (Roberts & Company Publishers, 2004), 3rd ed.
- E. Lalor, “Conditions for the validity of the angular spectrum of plane waves,” J. Opt. Soc. Am. 58(9), 1235–1237 (1968). [CrossRef]
- L. Onural, “Exact analysis of the effects of sampling of the scalar diffraction field,” J. Opt. Soc. Am. A 24(2), 359–367 (2007). [CrossRef]
- L. Onural, A. Gotchev, H. Ozaktas, and E. Stoykova, “A survey of signal processing problems and tools in holographic three-dimensional television,” IEEE Trans. Circ. Syst. Video Tech. 17(11), 1631–1646 (2007). [CrossRef]
- V. Katkovnik, J. Astola, and K. Egiazarian, “Discrete diffraction transform for propagation, reconstruction, and design of wavefield distributions,” Appl. Opt. 47(19), 3481–3493 (2008). [CrossRef] [PubMed]
- N. Delen and B. Hooker, “Free-space beam propagation between arbitrarily oriented planes based on full diffraction theory: a fast Fourier transform approach,” J. Opt. Soc. Am. A 15(4), 857–867 (1998). [CrossRef]
- J.-L. Kaiser, E. Quertemont, and R. Chevallier, “Light propagation in the pseudo-paraxial fresnel approximation,” Opt. Commun. 233(4-6), 261–269 (2004). [CrossRef]
- E. Sziklas and A. Siegman, “Diffraction calculations using fast Fourier transform methods,” Proc. IEEE 62(3), 410–412 (1974). [CrossRef]
- R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, “Shifted Fresnel diffraction for computational holography,” Opt. Express 15(9), 5631–5640 (2007). [CrossRef] [PubMed]
- F. Zhang, G. Pedrini, and W. Osten, “Reconstruction algorithm for high-numerical-aperture holograms with diffraction-limited resolution,” Opt. Lett. 31(11), 1633–1635 (2006). [CrossRef] [PubMed]
- K. Matsushima and S. Nakahara, “Extremely high-definition full-parallax computer-generated hologram created by the polygon-based method,” Appl. Opt. 48(34), H54–H63 (2009). [CrossRef] [PubMed]
- M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proc. IEEE 93(2), 216–231 (2005) (Special issue on “Program Generation, Optimization, and Platform Adaptation”). [CrossRef]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.