Patent ID: 12191911

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

Referring toFIG.1Aan apparatus for transmission and reception of an image through a disordered medium is generally identified by reference numeral100. The image may represent a binary sequence of data (such as a digital bit stream or byte stream for example) and the apparatus may be adapted to transmit a plurality of different images over time whereby the apparatus is useful for parallel communication; that is to say two or more binary digits bits of data may be transmitted simultaneously in each image. The image transmission may take place through a disordered medium such as an optical fibre, and may be through a multi-mode optical fibre.

As used herein a disordered medium may be any quasi-static medium in which light is scattered multiple times as it passes through the medium. The disordered medium may be any medium that comprises a heterogenous spatial distribution of optical refractive index or speed of light in the medium. Examples of such disordered medium include, but are not limited to, multi-mode optical fibres (both graded and step-index), optical diffusers, and tissues of the human or animal body.

The apparatus comprises a transmitter101and a receiver102. The transmitter101may be adapted for inputting, or projecting, coherent light into one end a multi-mode optical fibre118. As such the transmitter may comprise a plurality of controllable electromagnetic radiation sources. The receiver103may be adapted receive light leaving the opposite end of the multi-mode optic fibre118. As such it may comprise a plurality of electromagnetic radiation detectors.

The transmitter101may comprise a laser103. The laser103may be an actively Q switched diode pumped solid state laser such as the SPOT-10-200-532 laser available from Elforlight. The SPOT-10-200-532 is adapted to produce a beam of coherent light at 532 nm with a minimum pulse width of <1.8 ns and a maximum pulse energy of 10 μJ (at ≤10 kHz). The spatial mode of the laser is TEM00with a beam diameter ˜1 mm.

A first tube lens104may be positioned to receive the beam of light from the laser103. The first lens104may be an achromatic doublet-type having a focal length of 30 mm, a diameter of 25.4 mm, SM1-threaded mount and anti-reflective coating for the range 400-700 nm. A suitable lens is available from Thorlabs. Inc. under code AC254-030-A-ML. A second tube lens106may be positioned after the first tube lens104. The second tube lens104may be an achromatic doublet-type having a focal length of 75 mm, a diameter of 25.4 mm SM1-threaded mount and anti-reflective coating for the range 400-700 nm. A suitable lens is available from Thorlabs. Inc. under code AC254-075-A-ML.

The transmitter101may further comprise a digital micromirror device (DMD)108that may be arranged to receive light from the lens106. The DMD108comprises an array of micromirrors110for spatially modulating the laser light from the laser103, and which may act as the aforementioned a plurality of sources of electromagnetic radiation. The DMD108may be a DLP7000 available from Texas Instruments. Inc. The DMD108comprises an array of 1024×768 micromirrors. Other array sizes are possible.

The transmitter101may further comprise a third tube lens112that may be positioned to receive spatially modulated light from the DMD108and pass the beam of light to a first objective lens114. The third tube lens112may be an achromatic doublet such as an AC254-050-A-ML available from Thorlabs which has a focal length of 50 mm. The purpose of the third tube lens112is to focus the spatially modulated light onto the first objective lens114. The first objective lens114may be an infinity-corrected plan achromat such as an RMS20X available from Thorlabs. The purpose of the first objective lens is to focus the beam of light onto the proximal end116of a multi-mode fibre (‘MMF’)118. As explained in greater detail below, three different MMFs were tested in the apparatus100.

A distal end120of the MMF118may be positioned to guide the laser light to the receiver102. The receiver102may comprise a second objective lens124which receives the laser light from the distal end120. The second objective lens122may be the same type of lens as the first objective lens114, such as the RMS20X. The second objective lens may be positioned to improve the collimation of the laser light leaving the MMF118. A second tube lens123may be positioned to receive laser light from the second objective lens122and to further collimate the laser light. The second tube lens may be an achromatic doublet such as an AC254-100-A-ML available from Thorlabs which has a focal length of 100 mm.

The receiver101may further comprise a CCD camera124that may be positioned to capture laser light from the second objective lens122. The CCD camera comprises a plurality of output pixels which may act as the aforementioned plurality of electromagnetic radiation detectors. The CCD camera may be a CMOS device such as model C11440-22CU01 available from Hamamatsu. The C1440-22CU01 provides 4.0 megapixels resolution at100frame/s with 37,000:1 dynamic range. Output data from the CCD camera124may be captured and stored by a computer processor126and a memory127.

In use an input image representing a binary sequence of data may be transmitted from the DMD108to the CCD camera124via the MMF118. In particular, the laser103may generate a beam of coherent light128having a wavelength of 532 nm. The lenses104and106may spread and collimate the beam128so that its diameter increases from 1.0 mm to 2.5 mm.

Each micromirror of the DMD108may controlled by a controller129(which may comprise a computer processor and memory) so that a portion of the beam128is either reflected toward or away from the third tube lens112. By reflecting a portion of the beam128toward or away from the third tube lens112, the DMD108may be used to send a binary sequence of data through the MMF118. For example, it may be that reflecting a portion of the beam128toward the third tube lens112indicates a ‘1’ and reflecting a portion of the light away from the third tube lens112indicates a ‘0’, or vice-versa. As each micromirror of the array110is independently controllable in this way, parallel communication of the binary sequence is possible.

However, the MMF118supports a number of propagation modes that cause the image to become spread in time as the beam128travels along the MMF118. This modal dispersion affects the beam128so that a speckle pattern is seen by the CCD camera124rather than the image transmitted by the DMD108. The speckle pattern is a seemingly random variation in the intensity of the beam128across its diameter, and it appears that the modulation applied by the DMD108to the beam128is lost by the time the beam128is received by the CCD camera128.

The electric field component Emof the coherent light field at the mthoutput pixel of the CCD camera124received from the distal end120of the MMF118can be expressed as:
Em=Σn=Ntmn·En(1)
where Enis the electric field component of the light field at the nthinput pixel of the DMD108, with a total of N input pixels. In other words, the electric field Emat the mthoutput pixel is the sum of the electric field Enfrom each of the N input pixels. However the electric field from each input pixel is modified by a complex-valued intensity transmission constant, tmn=Amneiθmn, that links Enwith Em. In this way there are N complex-valued transmission constants for each output pixel m.

Emand Encan be expressed as Em=Ameiθmand En=Aneiθnwith amplitude A and phase θ. As explained above, in binary modulation each micromirror of the DMD108can be switched between two states (‘ON’ or ‘OFF’) independently, with the ‘ON’ micromirrors deflecting a portion of the beam128onto the proximal end116of the MMF118. Since the beam128is coherent, the light fields at all input pixels of the DMD108are assumed to have the same phase and amplitude. Thus, the phase θnis 0 whilst the amplitude Anis either 1 (‘ON’) or 0 (‘OFF’). Hence, the light intensity at then nthinput pixel is also either 1 (‘ON’) or 0 (‘OFF’), and the light intensity Imat the mthoutput pixel can be expressed as:
Im=|Σn=Ntmn·In|2(2)
When all the micromirrors are switched ‘ON’ (In=1), each micromirror appears to the receiver (CCD camera124) to produce a specific output light field with the same phase and amplitude as those of tmn.FIG.2Ashows a schematic illustration of the CCD camera124having an array of output pixels127. It will be recalled that each micromirror of the DMD108is ‘ON’ and is therefore directing a portion of the light, beam128toward the MMF118. The electric field Emreceived by each output pixel127ais schematically shown as a phasor127b. It can be seen that the phase of Emis variable from output pixel to output pixel, even though each micromirror (input pixels) of the DMD108has the same phase. This difference is due to the modal dispersion in the MMF118, and is approximated by the complex-valued intensity transmission constants described above.

Referring toFIG.2Bthe contributions127bfrom all the micromirrors of the DMD108form a total output light field Rmwith phase Ømand amplitude ARmat the distal end120of the MMF118. The amplitude ARmcan be considered as the superposition of all of the transmission constants tmnthat are projected on Rm, and the intensity IRmcan also be expressed as:
IRm=ARm2=ARm|Σn=NAmncos(θmn−Ømn)|·In(3)

Now considering the case when a binary pattern is input to the DMD108(seeFIG.2C), only a portion of the micromirrors of the DMD108are ‘ON’. The total output field Emat the mthoutput pixel is generated with a phase θmand an amplitude Am(seeFIG.2D). The intensity. Im, of the total output field can be expressed as:
Im=AmΣm|Σn=NAmncos(θmn−θm)|·In(4)
We can define a ratio, αm, of the amplitude of the total light field when all micromirrors are ‘ON’ on the DMD108to the amplitude of the total light field when fewer than all micromirrors are ‘ON’ as:

αm=ARmAm(5)
A ratio, βm, of the phase difference of the total light field when all micromirrors are ‘ON’ on the DMD108to the phase difference of the total light field when fewer than all micromirrors are ‘ON’ can be expressed as:

βm=❘"\[LeftBracketingBar]"∑n=N⁢Amn⁢cos⁡(θmn-∅mn)❘"\[RightBracketingBar]"❘"\[LeftBracketingBar]"∑n=N⁢Amn⁢cos⁡(θmn-θm)❘"\[RightBracketingBar]"(6)
Finally a parameter, θm, which represents the ratio of the intensity of total light field when all micromirrors are ‘ON’ on the DMD108to the intensity of the total light field when fewer than all micromirrors are ‘ON’ can be expressed as can be defined as θm=αmβm. Following that, the intensity of the total output field when some fraction of the micromirrors is on can be re-written as:

Im=1ηm⁢ARm⁢❘"\[LeftBracketingBar]"∑n=N⁢Amn⁢cos⁡(θmn-∅mn)❘"\[RightBracketingBar]"·In(7)
Interestingly, it was found that when the number of input pixels that are switched ‘ON’ (J) is sufficiently large compared to the total number of input pixels (N), the value of the parameter ηmremained mostly consistent across all output pixels m, with a mean value of η and a small standard deviation (described in greater detail below).

In other words, when J is sufficiently large, there is a pseudo-linear (or approximately linear) relationship between the intensity of the input image and the intensity of the output image. This pseudo-linearity enables the system (i.e. input image, disordered medium and output image) to be approximated with a set of linear equations, and the constants of the intensity transmission matrix may be determined using a compressive sensing technique. A particular advantage of this is that the intensity transmission constants become real-valued, rather than complex-valued. Another particular advantage is that the real-valued intensity transmission constants may be determined quickly and with low computational overhead compared to the model-based and deep-learning methods mentioned elsewhere herein.

As such it is possible to approximate the output intensity at each pixel m from the input intensities due to the number of ‘ON’ input pixels. J, using a matrix containing real-values of the intensity transmission constants tmn. This matrix is herein called an intensity transmission matrix (ITM) and the output light intensity, Im, at the m°ioutput pixel may be approximated as:
Im≈|Σn=Nitmmn)|·In(8)
where

itmmn=(1η)⁢ARm⁢Amn⁢cos⁡(θmn-∅mn)(9)
are the real-valued intensity transmission constants (the elements in the ITM) linking Inwith Im. Hereinafter the real-valued intensity transmission constants will be called the ‘RVITCs’.

As the variation (e.g. expressed in terms of standard deviation) of η across all output pixels decreases as J increases, when the variation of η is sufficiently small compared to the mean value of η, an input intensity pattern (i.e. image130). Iin, may be reconstructed from the output intensity pattern134. Iout, by inverting the ITM as:
Iin≈(ITM)−1*Iout(10)

In some embodiments the ITM may be inverted using other equivalent techniques, such as solving for the inverse linear problem using any available method including through matrix factorisations or iterative (potentially matrix-free) solvers.

By ‘high quality’ it is meant that the reconstructed image has a correlation coefficient with the input image of 90% or higher.

As described above, we have realised that it is possible to recover the image generated by the DMD108from the intensity speckle pattern at the CCD camera124.

In this way a receiver may determine the image sent by a transmitter through a MMF118(or other disordered medium) whereby the binary sequence of data may be determined. In order to recover the binary sequence, a characterization process is performed on the MMF118to determine the RVITCs itmmn.

The characterization process may comprise the use of compressive sensing to determine the RVITCs. In particular the characterization process may comprise steps of generating a series of known input images, each input image comprising a pattern representing a plurality of binary values. Each input image may be independent of each other input image (for example, the input images may be such that no input image is a linear combination of any other image). The pattern may be stored as an input matrix (or other computer-processable equivalent). Each input image may be transmitted from the transmitter101into the disordered medium (e.g. MMF118). An output image may be received at the receiver102(e.g. CCD camera124) from the disordered medium in the form of an intensity-only output speckle pattern (i.e. with no phase information). The values of each output intensity pattern may be stored in the computer memory126as an output matrix (or other computer-processable equivalent). The values of each intensity pattern may for one column of the output matrix. Each value may correspond to a pixel of the receiver102or may be derived from a plurality of pixels of the receiver102. The input and output matrices may be used as a system of linear equations to determine an intensity transmission matrix comprising a plurality of intensity transmission constants. The intensity transmission matrix may be stored in a computer memory for later use. Alternatively an equivalent representation or any representation of its inverse may be stored in a computer memory for later use.

The system of linear equations to be solved may be expressed as:

[I11…I1p⋮⋱⋮Im1…Imp]=ITM·[H1,H2](11)
In terms of compressive sensing, the left-hand side of this equation represents the measurements. The ITM132is the sparse matrix and the measurement matrix is [H1, H2], which is generated from a Hadamard matrix. The generation of [H1, H2] is explained in greater detail below.

Although the matrix on the left-hand side of Eq. (11) may appear to have a large number of values, it is in fact a small number of the possible measurements of the system (i.e. transmitter, multi-mode fibre118, and receiver). If the ITM132were to be found using traditional linear algebra techniques, it would be necessary to input all possible binary input patterns and record all corresponding output images. This is not possible for binary patterns of any appreciable size. For example, an 8×8 binary pattern has 264possible images. In contrast, by using compressive sensing, 2N input images can be used and the RVITCs determined in seconds. Recalling that N is the number of input pixels of the DMD108, in this example 2N=2×(32×32)=2048 images) which is a much smaller number. As mentioned elsewhere, it is not essential 2N images are used, and this number could be smaller or larger.

Other measurement matrices that are used in compressive sensing can be used to determine the ITM132, as long as these matrices have the restricted isometry property. Examples of other measurement matrices include random matrices that are generated to follow a certain type of distributions such as Gaussian. Bernoulli. and random Fourier ensembles, and deterministic matrices such as second-order Reed-Muller codes. Chirp sensing matrices, binary Bose-Chaudhuri-Hocquenghem codes, and quasi-cyclic low-density parity-check code matrix. Particular reference is made to Arjoune Y, Kaabouch N. El Ghazi H, Tamtaoui A.A performance comparison of measurement matrices in compressive sensing, International Journal of Communication Systems. 2018 Jul. 10; 31(10):e3576, which is herein incorporated by reference.

An embodiment of the characterization process is illustrated inFIGS.1B and1D. At step S1-1a plurality of input images130may be generated as described above. In this embodiment there may be 2N input images130. At step S1-2the 2N input images130may be transmitted sequentially into the disordered medium, such as the MMF118. Each input image130is affected by modal dispersion within the MMF118, and as described above, this is approximated by multiplication between each input image and the ITM132. At step S1-3a corresponding plurality of output images134is received by the CCD camera124. Each output image134may comprise an intensity-only speckle pattern in which the corresponding input image130cannot be recognized by the human eye. Phase information is not captured by the CCD camera124.

To construct the input images130, a generating matrix such as a Hadamard matrix H ∈(−1, +1) was constructed with dimension N×N using Sylvester's method (for further details attention is directed to J. J. Sylvester.Thoughts on inverse orthogonal matrices, simultaneous sign successions, and tessellated pavements in two or more colours, with applications to Newton's rule, ornamental tile-work, and the theory of numbers. Philosophical Magazine, 34:461-475, 1867, and which is incorporated herein by reference). Using this method, a first binary matrix H1was generated by replacing ‘−1’ with ‘0’, and then H was used to generate a second binary matrix H2by changing ‘−1’ to ‘1’ and ‘+1’ to ‘0’. A new matrix was generated from these two matrices as [H1, H2]. Since each matrix H1and H2has a size N×N, the matrix [H1, H2] has a size N×2N. i.e. N rows and 2N columns. The input images130were generated using the columns of the matrix [H1, H2]. In particular, the first input image130was generated using the first column of the matrix, the second input image using the second column, and so on to generate 2N input images. Each column of the matrix [H1, H2] was converted into a square matrix of sire √{square root over (N)}×√{square root over (N)} (recalling that each column of [H1, H2] has N elements. Each column of a binary matrix [H1, H2] was displayed as a square pattern on the DMD.

An advantage of generating the input images in this way is that the input images are independent of one another, so that no input image is a linear combination of any other image. This helps to ensure that the maximum information is obtained about the transmission matrix of the MMF118, and that there are no repeat measurements with the same input image.

As mentioned above, the output intensities may be expressed as:

[I11…I1p⋮⋱⋮Im1…Imp]=ITM·[H1,H2](12)
where Imprepresents the intensity value at the mthoutput pixel in the pthoutput image134, where p=1, 2, . . . , 2N. In other words, in the matrix on the left-hand side of this equation, the intensity values for each output image are placed in a respective column of the matrix. The number p of input images130may be lower or higher than 2N (recalling that in an embodiment N is the number of input pixels of the DMD108). At step S1-4, the output intensity values are stored in a multi-column output matrix (or other computer-processable equivalent), where each column of the matrix contains intensity values from one output speckle pattern.

As all micromirrors were switched ‘ON’ for the first input pattern, a standard Hadamard matrix [H, −H] was constructed as:

[2⁢I11-I11…2⁢I1p-I11⋮⋱⋮2⁢Im1-Im1…2⁢Imp-Im1]=ITM·[2⁢H1-1,2⁢H2-1]=ITM·[H,-H](13)
In this equation, the measurement matrix [H1, H2] has been expressed in terms of the original matrix H and the remaining terms adjusted accordingly. According to the properties of Hadamard matrices, the RVITCs itmmna of the ITM132can be obtained by multiplying both sides of this equation by [H, −H]T(the transpose of the matrix) to yield:

itmmn=(12⁢N)[2⁢Im1-Im1,2⁢Im2-Im1,…,2⁢Imp-Im1,2⁢Imp-Im1,2⁢Im2⁢N-Im1]*[hn1,hn2,…,hnp,…,hn2⁢N]T(14)itmmn=(12⁢N)⁢∑p=2⁢N⁢(2⁢Imp-Im1)⁢hnp⁢itmmn=(1N)⁢∑p=2⁢N⁢Imp⁢hnp
where hnp∈(−1, +1) represents the values at the nthinput pixel of the pthinput image130in the Hadamard matrix [H, −H].

Thus each RVITC of the ITM132may be found as follows:

for a first pair, mn, of output pixel m (n=1, 2 . . . m) and input pixel n (n=1, 2, . . . N):(a) take the first input and output image pair (p=1) and determine the product of (i) the measured output intensity or amplitude at output pixel m, (Imp), and (ii) a binary value, (hnp), indicating whether the corresponding input pixel n of the pair nm was on or off for that input and output image pair p (p=1);(b) repeat step (a) for each input and output image pair p (p=2, 3, . . . , P);(c) sum the products obtained in steps (a) and (b);(d) divide said sum by the number of input pixels N and store the result as the mnthRVITC in the ITM132;

repeat steps (a) to (d) for each other pair of output pixel n and input pixel n to generate m×n RVITCs and store as the ITM132.

This process may be performed at step S1-5, and the RVITCs stored in step S1-6. It is possible to further process the RVITCs into another equivalent form or version (e.g. an inverse of the ITM132) and to store such equivalent version instead of the ITM132itself.

An advantage of this is that the ITM132can be calculated comparatively quickly by a processor (since it involves only multiplication and addition of real-valued numbers).

It is noted that it is not essential to use binary patterns based on Hadamard matrices to generate the input images. Although Hadamard matrices provide some computational advantages (for example, the RVITCs of the ITM132may be obtained using the transpose of the original Hadamard matrix which is equivalent to the inverse), it may be that random binary pattern % are used as the input images for the characterization process. In that case, the Hadamard matrices. H, above would be replaced with a random binary matrix, B. and the inverse of B would be used to determine the RVITCs of the ITM132.

Once the ITM132has been generated and stored (either directly or in some other equivalent version or form), it is possible to use it to generate a reconstructed image136from an output speckle pattern134.FIGS.1C and1Eillustrate the reconstruction process. An input image130is transmitted by the DMD108into the MMF118. At step S2-4an output image134in the form of an intensity-only speckle pattern is received by the CCD camera124. At step S2-2the intensity values of the output image are stored in the memory126as an output matrix. At step S2-3the computer processor125processes the intensity-only values of the output image134by multiplying it with an inverse of the ITM132to generate the reconstructed image136. It is noted that the step of matrix multiplication may be performed using any equivalent computer-implemented operation (e.g. with a combination of linear operators) to achieve the same result.

To study the relationship between the intensities of input images130and output speckle patterns134, numerical simulations were performed with a custom MATLAB program. In order to generate the output speckle patterns134, a complex-valued transmission matrix TM with 8192 output pixels and 1024 input pixels was generated. The phases and amplitudes of the TM were randomly generated to obey uniform and Gaussian distributions between 0 and 2×, and 0 and 1, respectively. The characterization process described above was used to obtain the ITM132. In this case the ITM132has m×n=8192×1024=8,388,608 elements.

To investigate the effect of the number of switched ‘ON’ input pixels, a series of binary images with varying J (from 32, 64, 96 . . . to 1024) were generated as input images130. With each J, 64 different input patterns were generated by setting the values at J random pixel positions as ‘1’ and the rest of pixels with ‘0’. For each J 64 reconstructed images were determined from the output intensity speckle patterns. The standard deviations of αm, βmand ηmacross all output pixels were calculated and compared when J was varied. Correlation coefficients between reconstructed images and their corresponding input images (also called ‘ground truths’) were calculated for the evaluation of the image reconstruction performance. The correlation coefficients were calculated as a percentage of the reconstructed pixels correctly determined. It is worth noting that the correlation coefficients for the input and output images with all the mirrors switched ‘ON’ were calculated by changing the value of the first pixel of the ground truth from 1 to 0.999999 so that it is not undefined.

Several physical experiments were also performed to study the impact of a variety of fibre parameters on the performance of the image reconstruction retrieval algorithm. Firstly, to study the impact of input pixel counts (N) of the input images130, the number of input pixels on the DMD108used to generate each binary Hadamard pattern (used in each set 2N) was varied from 8×8, 16×16, 32×32 to 64×64, thereby producing ITMs132based on varying input pixel count. In other words, whilst the number of pixels in the binary pattern was kept constant (8×8), the number of illuminated micromirrors on the DMD108used the generate that binary pattern was varied.

After the characterisation process and the RVITCs were stored in the ITM132, a set of random binary patterns of 8×8 size were projected onto the DMD108as the input images130, or ground truths. As such, although the same set of binary patterns were used as ground truths, the input pixel count (N) of the DMD108varied and the reconstructed images136were based on different values of N. Correlation coefficients between the reconstructed images and their ground truths were calculated to evaluate the image reconstruction.

Secondly, to study the impact of the number of supported transverse modes of the MMF118, three fibres with different core diameters and numerical apertures (NA) were tested (see table300inFIG.3). Each fibre was 1 m in length. The mode count shown in the table ofFIG.3is calculated by:

Mfibre=12⁢(π·D·NAλ)2

where D is the diameter of fibre, NA is the numerical aperture and λ is the wavelength of light beam128from laser103.

Thirdly, to study the impact of variability of the input patterns, binary images of different types, including handwritten figures, schematic plants, animals. Chinese characters and random patterns, were used as ground truths for image reconstruction with a step-index multimode fiber (diameter 200 μm. NA=0.22, length=1 m). In addition, the reconstructed images were binarized with the Otsu threshold method, which is available in Matlab. The accuracies of reconstructed binary images, which represented the percentages of pixels with correct values, were calculated.

Finally, to study the impact of the fibre length, three step-index fibres with the same diameter (Ø200 μm) and NA (NA=0.22), but different lengths 0.1 m, 1 m and 10 m) were used for the retrieval of the same input image. After the fibre characterisation process, the input image130(ground truth) was displayed on the DMD108while output speckle patterns134were captured at different times. In order to evaluate the output decorrelation over time caused by fibre drift, correlation coefficients were calculated between (i) each output speckle pattern134and the first output speckle pattern, and (ii) between each reconstructed image136and the input image130(ground truth).

Results

Pseudo-Linearity

FIG.4Ais a graph400of the number of ‘ON’ input pixels J (x-axis) versus the standard deviation (y-axis) of the parameters αm, βmand ηmas determined in the numerical simulation described above. As the number of ‘ON’ input pixels J increases from 32 to 1024, the standard deviations of both αmand ηmdeclined rapidly, while the change in the standard deviation of βmwas much smaller. It is recalled that αmis the ratio of the amplitude of the total light field when all micromirrors are ‘ON’ on the DMD108to the amplitude of the total light field when fewer than all micromirrors are ‘ON’, whilst ηmis the same ratio in terms of intensity of the light field. The decline in the standard deviation of αmand ηmindicates that the amplitude and intensity of the light field at each output pixel127atends to a more constant value as J increases (i.e. the number of input pixels switch ‘ON’ increases). The standard deviation of the parameter βm(the ratio of the phase of the output light field to the input light field) is a much smaller value, ranging from about 0.3 when J=32 to about 0.25 when J=640. This indicates that there is a very small effect on the ratio of the output and input light fields J increases.

FIG.4Bis a graph402of the number of ‘ON’ input pixels J (x-axis) versus the correlation coefficient (y-axis) determined by comparing each reconstructed image136to its corresponding input image130. As such the graph402indicates the performance of the image reconstruction algorithm. Two trends are shown: a first trend404indicating performance of the image reconstruction algorithm when performed in the numerical simulation environment, and a second trend406indicated the performance of the image reconstruction algorithm in the physical experiments. A correlation coefficient of 1.0 indicates that all output pixels match the input pixel values (0 or 1), whereas a correlation coefficient of 0.0 indicates that none of the output pixel values match the input pixel values.

With reference to the first trend404(numerical simulation) the correlation coefficients between reconstructed images136and their corresponding input images130rapidly increased from 0.5 to 0.9 as J increased from 32 to 320 and then gradually increased from 0.9 to 1 as J increased from 320 to 1024. With reference to second trend406(physical experiments), the correlation coefficient also increased rapidly from 0.051 to 0.89 as J increases from 32 to 384, and remained largely consistent as J increases from 384 to 896. As J increased from about 986 to about 1024, the correlation coefficient decreased rapidly from 0.85 to 0.1. This discrepancy between the simulations and measurement results may be attributed to the loss of low-spatial-frequency information in the input images130due to diffraction of light from the micromirror array110of the DMD.

Referring to bothFIGS.4A and4Bit can be seen that the ratio between the light intensity of the input image and the output speckle pattern (indicated by αmand ηm) becomes roughly constant, or pseudo-linear, in region408when Nis greater than about 30% of the total number of input pixels. An alternative way of expressing this is that the correlation coefficient becomes roughly constant when N is greater than about 30% of the total number of input pixels. The starting point of this pseudo-linearity is not precisely defined, and could be said to start somewhere in the range 25%-40%. The starting point may also be dependent on the acceptable level of correlation coefficient. For example, some methods may need only a lower correlation coefficient compared to some other methods. e.g. in biomedical endoscopy, lower correlation coefficients can be tolerated if the user can see the necessary detail of the target in the reconstructed image.

Input Pixel Count and Number of Transverse Modes

Turning now to the physical experiments, three different multi-mode fibres were tested.FIG.3shows properties of those multi-mode fibres. Each multi-mode fibre took the position of the MMF118inFIG.1, and the fibre characterisation and image reconstruction methods described above were used to assess the effect of the different fibres. Prior to transmission of any input image130, each fibre was characterised and a corresponding ITM132stored in memory. After characterisation of each fibre, a series of input images were transmitted into the fibre and received as output speckle patterns134at the CCD camera124.FIG.5Aillustrates the binary pattern500that was passed through each a fibre. However, it is noted that the number of pixels (on the micromirror array 11 of the DMD108) making up the binary pattern500was 8×8, 16×16, 32×32 and 64×64. Thus a series of four input images was transmitted, with each input image constructed with a different number of input pixels on the DMD108.

Once each output speckle pattern was received by the CCD126, a corresponding reconstructed image was generated from each output speckle pattern. An example of the reconstructed images is shown inFIG.5A: reconstructed images136a-136dwere generated from the output speckle patterns134received at the CCD camera124from the fibre called “Fiber-200-0.22” inFIG.3.

FIG.5Bis a graph502of the number of input pixels used to form the binary input pattern500(x-axis) versus the correlation coefficient between the input image130and the reconstructed image136. It is noted that the number of input pixels inFIG.5Bis the total number of input pixels (e.g. 8×8, 16×16, etc.) used on the DMD108to transmit the binary input pattern500, and is not the same as the number of ‘ON’ pixels, J, mentioned above. First trend504indicates the results for the “Fiber-200-0.50”, second trend506indicates the results for the “Fiber-200-0.22”, and third trend508indicates the results for the “Fiber-105-0.22”. It can be seen that as the number of input pixels increases, the quality of the reconstructed image (defined by the correlation coefficients) decreased. For example, second trend506shows that the correlation coefficient declined from 99.44% to 76.81% as input pixel count increased from 8×8 to 64×64.

The core diameter and NA of each MMF had significant impact on the quality of the reconstructed image. This can be explained by the varying number of supported transverse modes in the fibres (seeFIG.3). The fibres with larger numbers of supported transverse modes (Fiber-200-0.22. Fiber-200-0.50) were able to transmit images with higher input pixel counts.

The computation time for the characterisation process (in order to estimate the RVITCs of the ITM132) increased with the increase of both the input and output pixel counts (N and M). For example, with a desktop PC (Intel i7 8700, 3.2 GHz, 16 GB RAM), when N and M increased from 32×32 and 360×360 to 64×64 and 500×500, respectively, the computation time for ITM estimation increased from ˜8 s and ˜240 s, respectively.

Variability of Binary Input Patterns

The apparatus ofFIG.1was also used to investigate the impact of changes to the binary input pattern used to generate the input image130. In this part of the experiment the MMF118was the optical fibre “Fiber-200-0.22” (seeFIG.3). The characterisation process was used to generate and store an ITM132for the optical fibre.

After the characterisation process was completed, a series of input images were transmitted through the MMF118.FIG.6Ashows different input images600a-600h(also called ground truths) that were generated using 32×32 input pixels of the DMD108. Each input image600a-600hwas transmitted into the multi-mode fibre and corresponding output speckle patterns602a-602hwere received at the CCD camera124. The image reconstruction process was used to generated reconstructed images604a-604hfrom the output speckle patterns. Correlation coefficients606a-606hwere calculated for each pair of input image600a-600hand reconstructed image604a-604h. Finally, the reconstructed images604a-604hwere binarized to generate binary output images608a-608h.

The quality of the reconstructed images604a-604hwas weakly dependent on the binary input patterns in the input images600a-600h. In particular the correlation coefficients between the reconstructed images604a-604hand the input images600a-600hvaried from 91.64% for a handwritten digit (600a/604a) to 97.56% for a random binary pattern (600h/604h). A further experiment was performed with the same set up, except that a higher number of input pixels (64×64) was used on the DMD108for each input image600a-600h. In that case the correlation coefficient varied from 76.22% for the handwritten digit (600a/604a) to 90.43% for the random binary pattern (600h/604h), respectively.

The accuracy610a-610hof the binary output images608a-608hdemonstrates that there is a weak dependency on quality of the reconstructed images604a-604h(as defined by the correlation coefficient606a-606h). In particular, the accuracy610a-610hof each binary output image608a-608hwas almost 100% irrespective of the correlation coefficient. This indicates that the apparatus and methods described herein would be especially useful for transmission and reception of binary data across a disordered medium, such as a multi-mode optical fibre.

FIG.6Bshows the results of a further numerical simulation (similar to that reported with reference toFIG.4B) to investigate whether the methods described herein can be used to reconstruct grayscale input images, rather than binary input images. The input images had 32×32 input pixels and each pixel was randomly assigned a value between 0 and 255. As can be seen, grayscale input images may be reconstructed with an accuracy of greater than 98%. This demonstrates the plausibility of the methods described herein for transmitting grayscale input images.

Fibre Length

The apparatus ofFIG.1was used to investigate the effect of optical fibres of different lengths. Three MMFs118of length 0.1 m, 1.0 m and 10 m were investigated. Each MMF118had a diameter 200 μm and numerical aperture 0.22. Prior to transmission of an input image through each MMF118, the characterisation process was used to generate and store an ITM132for that optical fibre.

Referring toFIG.7Athe input image700was a binary pattern of the letters ‘KCL’. 32×32 input pixels of the DMD108were used to transmit the input image into the MMF118. The CCD camera124received output speckle patterns702a,702band702ccorresponding to the 0.1 m, 1.0 m and 10 m MMFs118respectively. Reconstructed images704a,704b,704cwere determined for each input image700and correlation coefficients706a,706b,706cwere determined between each reconstructed image and the input image. As shown inFIG.7A, all MMFs118produced high correlation coefficients of 97%, 96% and 94%, respectively. These results suggest that the characterisation and image reconstructions processes are insensitive to fibre length.

However, it was found that the 10 m fibre suffered from fibre drift (causing decorrelation of the output speckle patterns). Referring toFIG.7Ba graph710of time (x-axis) versus correlation coefficient (y-axis). In this particular part of the experiment, the correlation coefficient of both the output speckle pattern and the reconstructed image (each compared to the input image) was monitored over time. First trends712aand712bshow the change in correlation coefficient of the reconstructed image and output speckle pattern using the 0.1 m fibre. Second trends714aand714bshow the change in correlation coefficient of the reconstructed image and output speckle pattern using the 1.0 m fibre. Third trends716aand716bshow the change in correlation coefficient of the reconstructed image and output speckle pattern using the 10 m fibre.

The first trends712band second trends714bshow that the correlation coefficients of the output speckle patterns using the 0.1, m and 1.0 m length fibres remained relatively stable (˜99%) over a 5-minute period. Accordingly the first trends712aand714aof the correlation coefficient of the reconstructed images also remained at a stable level (˜97% and ˜96% respectively). However, for the 10 m fibre, the output speckle pattern captured 5 minutes after fibre characterization process had degraded from 100% to ˜92%, whilst the correlation coefficient of image retrieval degraded from ˜94% to ˜75%. The faster degradation of the output speckle pattern from the 10 m fibre was mainly caused by two factors: first, the longer length suffered more serious fibre drift; and second, both 0.1 and 1.0 m fibres were in cable suits and fixed on an optical table, whilst the 10 m fibre was twined on a mount and hence suffered more vibration.

Referring toFIG.8an apparatus for transmission and reception of an image through a disordered medium is generally identified by reference numeral800. The apparatus may comprise a transmitter801and a receiver.

The image may represent a binary sequence of data (such as a digital bit stream or byte stream for example) and the apparatus may be adapted to transmit a plurality of different images over time whereby the apparatus is useful for parallel communication; that is to say two or more binary digits bits of data may be transmitted simultaneously in each image. The image transmission may take place through a disordered medium such as an optical fibre, and may be through a multi-mode optical fibre.

The transmitter801may comprise a plurality of controllable electromagnetic radiation sources. In an embodiment the plurality of controllable electromagnetic radiation sources may comprise a light transmitter array unit803optically coupled with a light modulator array unit804. The light transmitter array unit803may comprise a plurality of coherent light sources that are optically coupled via respective a bundle of optical fibres (not shown) to the light modulator array unit803. These coherent light sources may be laser diodes, or solid-state lasers. The light transmitter array unit803may be a 2D array of vertical-cavity surface-emitting lasers (VCSEL).

The light modulator array unit804comprises a plurality of input pixels806. In an embodiment the light modulator array unit804may be a spatial light modulator. In an embodiment the light modulator array unit804may be an array of electrical circuits that modulate the current or voltage supply for driving the light transmitter array unit803. In another embodiment the spatial light modulator may be a deformable mirror. Each input pixel806is controllable to either transmit or not transmit a portion of the light received from the light transmitter array unit803. In this way the light modulator array unit804may indicate a binary pattern on its output side. For example a binary ‘1’ may be indicated by light being allowed to pass through an input pixel806, and a binary ‘0’ may be indicated by light not be allowed to pass through an input pixel806, or vice-versa. The control of each input pixel is performed by a first computer processor808in conjunction with a first memory810(such as RAM and/or non-volatile memory). The first computer processor808and first memory810may be in the form of an ASIC, system-on-a-chip or photonic-integrated circuit.

The computer processor808and memory810are adapted to cause a binary input pattern812to be displayed on the output side of the light modulator array unit806. The binary input pattern812shown inFIG.8is merely exemplary. The binary input pattern812may also be a pattern indicating a sequence of binary data.

Light from the light modulator array unit804passes to a first lens814that focuses the light to the input of a first objective lens816. The first objective lens816focuses the light onto the proximal end817of an optically disordered medium818. The optically disordered medium818may be a multi-mode optical fibre, such as a step-index multi-mode optical fibre. The light passes through the optically disordered medium818and out of a distal end819. The light is scattered as it passes through the optically disordered medium818.

The light passes through a second objective lens820that collimates the light leaving the optically disordered medium818. A second lens822may further collimate the light before it arrives at a plurality of electromagnetic radiation detectors. In this embodiment a plurality of electromagnetic radiation detectors may be a light detector array unit824. The light detector array unit824may be a focal plane array. The light detector unit may be an array of optical fibres arranged such that light is collected at the ends of the fibres nearest the second lens822and is delivered to a photodetector, such as a photodiode. The light detector array unit824may comprise a plurality of output pixels825. Each output pixel may be a photodetector, such as a photodiode. The light detector array unit824may be an avalanche photodiode (APD) array, such as an 8×4 Si APD array available from Hamamatsu (product code S8550-02). An analogue-to-digital converter (ADC array unit826for readout of output signals from the output pixels825. The ADC array unit may also digitize the output signals and make the digitized signals available to a second computer processor828. The computer processor828may store the digitized signals in a second memory830. The digitized output signals represent an output speckle pattern832received at the light detector array unit824resulting from the scattering of the light as it passes through the optically disordered medium818. The second computer processor828and second memory830may be in the form of an ASIC, system-on-a-chip or photonic-integrated circuit.

In order to implement the characterisation process to obtain image intensity constants of the ITM for the optically disordered medium818, the first memory818may stow computer executable instructions for implementing the input image transmission steps of the characterisation process steps described above in conjunction withFIG.1D. In particular steps S1-1and S1-2may be stored as computer executable instructions in the controller129. Furthermore, the second computer processor828and second memory830may store computer executable instructions for implementing the receiver steps of the characterisation process described with reference toFIG.1D. In particular, steps S1-3, S1-4, S1-S and S1-6be stored in the memory126and executed by the processor125.

In order that the second computer processor828may determine the image intensity constants, the second memory830may have stored the set of input images transmitted under control of the first computer processor808. The set of input images may be transferred from the first memory810to the second memory830by a separate transmission method (not shown), such as use of the Internet, at some point during the characterisation process or beforehand.

Once the image intensity constants have been determined for the optically disordered medium818, these may be stored by the second processor828in the second memory830, and the apparatus800may be used for parallel communication of binary data. In order to implement the image reconstruction process described above in conjunction withFIG.1E, the second memory830may store computer-executable instructions for implementing the steps of the reconstruction process.

The first computer processor808and first computer memory810may be adapted to transmit a sequence of input images representing binary data. For example, each input image may comprise a 2×2, 4×4, 8×8, etc. image of binary data using the light modulator array unit814. The second computer processor828and second memory830may receive each output speckle pattern corresponding to the input image and use the image reconstruction process to recover a reconstructed image. The reconstructed image may be binarized to recover the binary data of the input image.

In some embodiments, input images may be displayed by the light modulator array unit804at a rate of 22,000 frames per second. However, higher or lower speeds are envisaged. The light detector array unit824may be able to capture the output images at a rate of 250 frames per second. However, higher or lower speeds are envisaged. The display and detection rates may be matched.

It may be desirable to repeat the characterisation process from time-to-time. In this way a new set of RVITCs (or any computer-processable equivalent) is generated and stored by the second memory830for use in the next image reconstruction process. The repetition of the characterisation process may take place periodically (for example every 30 mins, although other time periods are envisaged which may be dependent on the length of the fibre—a longer fibre may require more frequent characterisation), or as often as necessary desired. It may be that known input images are transmitted at certain time intervals or every nth input image. In this way the second computer processor828and second memory830may check reconstructed images remain accurate (for example with a correlation coefficient above a certain threshold, e.g. 99%). If the accuracy has degraded, the characterisation process may be triggered.

Referring toFIG.9an apparatus for obtaining an image inside a body is generally identified by reference numeral900. The apparatus900comprises an endoscope902and a control console904. The endoscope902comprises a tube905. The tube905may be substantially rigid or substantially flexible. The tube905contains a single-mode optical fibre906and multi-mode optical fibre908, both of which terminate at a distal end910of the tube904. The endoscope902may comprise a grip (not shown) via which an operator may support the weight of the endoscope and assist its entry into a body912. The fibres906and908extend from a proximal end914of the tube905and into the control console904. Mounted inside the control console904is a transmitter914, a receiver916and an illuminator918. The transmitter914is similar to the transmitter part of the apparatus800described above. That is the transmitter914comprises a laser, spatial light modulator, and a computer processor and memory (not shown inFIG.9) for controlling the spatial light modulator to display known input images. However, in this embodiment the transmitter914is used only for the characterisation process. In particular, the control console904comprises a port920for receiving the distal end910of the tube905. When the distal end905is connected to the port920, the end of the fibre908is brought to a position where it may receive input light from the transmitter914as will be described in greater detail below.

The receiver916is similar to the receiver part of the apparatus800described above. That is the receiver comprises a light detector array unit922, and optics924for receiving a computer processor and memory for processing received light intensity data and reconstructing images using the image reconstruction process. The control console904may be connected for a display926so that reconstructed images may be displayed to a used as the endoscope902is in use.

The illuminator comprises a laser928arranged to transmit light into the fibre906. In use, light travels along the fibre906and leaves the distal end910of the tube905to illuminate inside the body921.

Before the endoscope902can be used to view an internal part of the body912, the multi-mode optic fibre908may be characterised in order to determine and store an ITM containing the RVITCs. To do that, the distal end910of the tube905is inserted into the port920. That brings the end of the fibre908into alignment with the transmitter914. Once in place, the characterisation process can be performed on the fibre908and the RVITCs stored by the receiver916.

It is noted that, since the transmitter914and receiver916are in the same location (in the control console904), they may share computing resources. For example, the transmitter914and receive916may share one or more computer processors. In another embodiment, the transmitter914and916may have dedicated computing resources. For example, the transmitter914and receiver916may have one or more dedicated computer processor (e.g. in the form of an ASIC or ASIP).

After the characterisation process is completed, the endoscope902is ready for use. If the endoscope is being used to look inside a patient (e.g. inside the human or animal body), the distal end910of the endoscope902is inserted through an opening or cavity in the body. Light from the laser912illuminates an interior portion of the body, for example an imaging target930. Light is scattered and reflected by the imaging target930and a portion of the light is received by the end of the multi-mode fibre908at the distal end910. This light is an input image into the multi-mode fibre908, in a similar way to input images generated by the digital micromirror device described in embodiments above.

The input image is scattered inside the fibre909as it travels toward the receiver916. At the receiver916the input image has become an output speckle pattern, as described above. The receiver916may take samples of output speckle pattern (e.g. at a certain number of frames per second), and may use the image reconstruction process described above to generated a reconstructed image for each sample of the output speckle pattern. The number of frames per second may be high enough so that a video may be displayed on the display926.

Focusing Light at a Target Via a Disordered Medium

Referring toFIG.10an apparatus for focusing electromagnetic radiation (e.g. laser light) at a target via a disordered medium is generally identified by reference numeral1000. As used herein the phrase “focusing light at a target via a disordered medium” is intended to include a number of different scenarios including, but not limited to: (a) light transmitted from a first side to a second side of the disordered medium with focusing taking place beyond the second side; (b) light transmitted into the first side, scattered inside the disordered medium and then focused inside the disordered medium; and (c) light transmitted toward the first side, diffusely reflected from the first side, and then focused after that diffuse reflection. In all cases, including (a), (b) and (c), the focusing is achieved by controlling transmission of the light before it enters or is reflected from the disordered medium.

The electromagnetic radiation may be collimated, substantially spatially and temporally coherent light (e.g. laser light). The electromagnetic radiation may be at optical wavelengths, for example between about 100 nm and about 1000 nm. The electromagnetic radiation used by the apparatus1000may be pulsed or may be intensity-modulated continuous wave.

The disordered medium may be any disordered medium as described anywhere in this document. Examples of such disordered media include, but are not limited to, in vitro and in vivo biological tissue belonging to a human or other animal, an optical waveguide such as a multi-mode optical fibre, a surface from which optical light is reflected such as a wall, floor or ceiling in a building, and any aerosol such as fog, mist or dust.

The apparatus1000comprises a transmitter1002and a receiver1004. The transmitter may be adapted for inputting, or projecting, collimated and spatially and temporally coherent electromagnetic radiation at optical wavelengths into an optical system and may be adapted to cause the light to be focused onto a target1006even though it passes via a disordered medium1008before reaching the target. As such the transmitter1002may comprise a plurality of electromagnetic radiation sources. The receiver1004may be adapted for receiving a direct or indirect indication of the intensity of electromagnetic radiation at the target. In one embodiment the receiver1004receives an indirect indication in the form of ultrasound waves generated by the photoacoustic effect, generates signals representing those ultrasound waves and then processes those signals, as described in greater detail below.

The transmitter1002may comprise a laser1010. The laser1010may be an actively Q switched diode pumped solid state laser such as the SPOT-10-200-532 laser available from Elforlight. The SPOT-10-200-532 is adapted to produce a collimated beam of spatially coherent light at 532 nm with a minimum pulse width of <1.8 ns and a maximum pulse energy of 10 μJ (at ≤10 kHz). The spatial mode of the laser is TEM00, with a beam diameter ˜1 mm.

The laser1010is positioned so that, in use, the beam may be directed onto an achromatic doublet comprising a first convex lens1012(f=30 mm) and a second convex lens1014(f=50 mm). In an embodiment the achromatic doublet may be model AC254-030-A-ML available from Thorlabs, Inc. Similarly to the apparatus100, the apparatus1000comprises a digital micromirror device (DMD)1016. The DMD1016is identical to the DMD108and all details of the DMD108described above are incorporated into this embodiment. The purpose of the achromatic doublet is to expand the beam diameter to cover the DMD1016. A first computer1017may be provided for generating and storing DMD patterns, and for sending those DMD patterns and instructions to the DMD1016.

The transmitter1002may comprise a synchronization device1018for synchronizing intensity variation (e.g. as individual pulses, or as intensity modulation of a continuous wave) of the laser beam with the operation of the DMD1016as described in greater detail below. In this embodiment the synchronization device1018is an arbitrary waveform generator (‘AWG’) such as a model33600A available from Key sight Technologies, Inc.

The transmitter1002may comprise a third convex lens1019(f=30 mm, AC254-030-A-ML available from Thorlabs) positioned to received spatially modulated light from the DMD1016and to direct that light onto a diffuser1020. The diffuser1020may be a ground glass diffuser such as model N-BK7 with a 220-grit polish available from Thorlabs. The purpose of the diffuser1020is to act as a disordered medium and provide an effect on the spatially modulated light similar to scattering by biological tissue for example, or diffuse reflection from a surface such as a wall, floor or ceiling.

The target1006may be placed about 5 mm behind the diffuser1020to receive light that has passed through it. The target1006acts as an optical absorber and may be a piece of black insulation tape, although any other material that is capable of absorbing laser light and generating ultrasound can be used. These materials include absorbing tissue chromophors including, but not limited to, haemoglobin, myoglobin, melanin, bile, collagen, deoxyribonucleic acid, ribonucleic acid, lipid and water. The purpose of the target is to absorb at least a portion of the incident light from the diffuser1020so that it may be detected by a detector1022. The detector1022may be an ultrasonic transducer such as a flat single-element piezoelectric ultrasonic transducer. In the present embodiment, the ultrasonic transducer is model V358 from Olympus, which has a central frequency of 50 MHz and a diameter of 6.4 mm (0.25 inches). A focusing element1024may be provided in the ultrasonic detector to provide a receiver focus1025. In one embodiment the focusing element1024may be silica piano-concave lens attached on the active surface of the ultrasonic transducer, such as model LC4210 from Thorlabs which has a focal length of −25 mm. Both the target1006and the detector1022were immersed in a water bath1026to provide acoustic coupling between them. The distance between the detector1022and the target1006was adjusted to maximise received ultrasound signals so that the target was situated at the receiver focus1025of the detector1022. In use, the detector1022receives ultrasound signals from the target1006and provides an output electrical signal representing the detected ultrasound.

An amplifier1028may be provided to receive and amplify the electrical signal from the ultrasonic detector. The amplifier may be a model SPA.1411 available from Spectrum Instrumentation. A second computer1030(e.g. a personal computer having a processor such as an Intel® i7, 3.2 GHz) may be connected to the amplifier1028to receive an amplified electrical signal from the amplifier1028. The computer1030may have a Data Acquisition (DAQ) card1032, such as model Mi.4420 available from Spectrum Instrumentation, to digitise the amplified electrical signal. The computer1030may also comprise a computer memory1034for storing digitised ultrasonic data generated by the DAQ1032. It is noted that in other embodiments the first computer1017and the second computer1030may take different forms, such as ASIC, ASIP and system on a chip, etc.

The computer memory1034also stores a set of computer-executable instructions that, when executed, perform steps to determine a set of real-valued intensity transmission constant (‘RVITCs’) for the diffuser1020. The steps are described in greater detail below. However, the purpose of the steps is to characterise the diffuser1020and determine and store a focusing spatial modulation pattern1036that will achieve an improved focus of the light through the diffuser1020onto the target1006. Once the focusing spatial modulation pattern1036has been determined, the second computer1030may communicate it to the first computer1017for storage by the first computer1017and subsequent use. The communication between the second computer1030and the first computer1017may take place using a data communication system1040(shown schematically by an arrow inFIG.10). The data communication system is not a part of the apparatus1000, although the transmitter1002and receiver1004may have one or more interface (not shown) respectively for sending and receiving data over the data communication system1040. The data communication system may take a number of forms, including by not limited to, wireless systems (e.g. WiFi. Bluetooth, cellular) and wired systems (e.g. Ethernet, power-line), or any combination of these.

In certain aspects, to provide a visual evaluation of the focusing, the target1006and detector1022were replaced with a CCD camera provided with a convex lens (e.g. AC254-050-A-ML from Thorlabs) for capturing the output speckle patterns at the focus of the ultrasound transducer.

Operation

In use, light is generated by the laser1010, spatially modulated by the DMD1016and transmitted to the diffuser1020. The light is scattered by the diffuser1020and portion is transmitted toward the target1006. At least some of the energy of the light which reaches the target1006is converted via the photoacoustic effect into ultrasound that is detected by the detector1022.

A difference between the operation of the apparatus1000and operation of the apparatus100is that during use of the apparatus100there are P output images, with each output image p comprising in output pixels and each pixel representing an intensity value. For each input DMD pattern k there are thus m output intensity values. However, during operation of the apparatus1000each output p is a single ultrasonic intensity signal taken by the detector1022. Since the detector1022is arranged to measure intensity of ultrasound at the receiver focus1025on the target1006, each ultrasonic intensity signal represents the degree to which light reaching the target1006is focused at the focal point: a higher valued measurement indicates better focus, whereas a lower valued measurement indicates that the light leaving the diffuser1020is dispersed over a wider area of the target1006. In this way it is possible to use the receiver focus1025of the detector1022as a guide to focus light onto the target1006.

It is possible for the second computer1030to determine the ultrasonic signal intensity by any number of techniques. For example within a given period of time: as the difference between maximum and minimum measured values; as the maximum measured value: as the absolute value of the minimum measured value: integration of the absolute value of the measured signal over time; as the average of the absolute values of the measured values: as the standard deviation of the absolute values of the measured values, to name but a few.

The method of focusing electromagnetic radiation at the target1006comprises two steps: a characterisation process step and a focusing process.

Characterisation Process

As explained above with reference toFIG.1D, the DMD1016may be controlled to provide a sequence of spatial modulations, in this embodiment 2N input patterns1038to characterise the diffuser1020. Since there is a single output ultrasonic intensity signal for each input pattern k some adaptation of the characterising process is required, as explained in the following.

As described above a transmission matrix comprising a plurality of real-valued transmission constants (‘RVITM’) can be used to approximately connect the input and output light intensities of a disordered medium (e.g. multimode optical fibre, diffuser1020), so that input binary and grayscale images can be retrieved from measured intensities of the output speckles. We have found that the RVITM can also be used for focusing light scattered by a disordered medium since the RVITM encodes both the phase and amplitude information of the light field changes from the DMD1016to the receiver focus1025. The present method enables a non-iterative and higher speed characterisation of the disordered medium than has been reported previously, and directly determines a preferred light input (e.g. certain DMD pattern) to achieve improved focusing of light scattered by the disordered medium.

Similarly to the characterisation process described above, the first step is to use a Hadamard matrix H ∈(−1, +1) with dimensions of N×N to construct two binary matrices H1=(H+1)/2 and H2=(−H+1)/2. Each column of the binary matrix, [H1, H2] was then converted to a square matrix that was used to spatially modulate the incident laser beam onto the diffuser1020using the DMD1016, whilst the corresponding light-generated ultrasound waves were recorded by the detector1022.

According to the principles of photoacoustic signal generation, the amplitude of the received ultrasound signal Q with the kth input pattern displayed at the DMD1016can be expressed as: Qk=aSTμaFk, where a is a constant account for the attenuation loss during ultrasound propagation, S is the sensitivity of the detector1022, Γ is the Grüneisen coefficient (a dimensionless constant defining the conversion efficiency of heat energy to pressure). μais the optical absorption coefficient, and Fkis the local optical fluence at the target1006. With Fk=∫Ik′dt/A, where ∫Ik′dt is the light intensity, Ik′, at the target1006when the kth pattern is displayed and A is the illumination area at the target1006. Since the pulse duration T is constant, Fkcan be further expressed as

Fk=Ik⁢TA
where Ikis the average light intensity over T, and we have

Qk=aS⁢Γμa⁢TA⁢Ik⁢where⁢aS⁢Γμa⁢TA
is a constant under the conditions of the experiment, which we define as

α=aS⁢Γμa⁢TA.
As Ikis linearly proportional to Qkincreasing the amplitude of Qkis equivalent to improving the focusing of light on the target1006. For example if the amplitude of Qkcan be maximised the light would be focused on the target1006. By ‘focused’ it may be meant that the intensity of the light at the target1006is increased over the light intensity at the target1006when a random DMD pattern is displayed on the DMD1016. In some embodiments. ‘focused’ may mean that the light intensity at the target1006is at a maximum or close to a maximum over the light intensity at the target1006when a random DMD pattern is displayed on the DMD1016

As explained above the approximate relationship between the intensities of the input and output light field through a disordered medium can be connected by the 24) RVITCs of the RVITM. Accordingly it is also possible to approximately connect the light intensity Ikat the target1006and the input patterns1038(represented by the matrix [H1, H2]) by a row of the RVITCs. RVITMr, in the intensity transmission matrix RTITM. The reason there is a single row of RVITCs is that each input pattern A is represented as a single column in the matrix [H1, H2]. The row of elements RVITMrcorresponds to the light transport from all the DMD input positions to the target1006, and may be expressed as:

[I1,I2,…,I2⁢N]=1α[Q1,Q2,…,Q2⁢N]=RVITMr·[H1,H2](15)

Equation (15) can be re-arranged as follows:

2×1α[Q1,Q2,…,Q2⁢N]=2×RVITMr·[H1,H2]1α[2⁢Q1,2⁢Q2,…,2⁢Q2⁢N]=2×RVITMr·[H1,H2]1α[2⁢Q1,2⁢Q2,…,2⁢Q2⁢N]-RVITMr·[D1,D1⁢…⁢D1]=2×RVITMr·[H1,H2]-RVITMr·[D1,D1⁢…⁢D1]
where D1is a column matrix of dimension N×1 with all elements having a value of 1. With further re-arrangement:

1α[2⁢Q1,2⁢Q2,…,2⁢Q2⁢N]-RVITMr·[D1,D1⁢…⁢D1]=RVITMr·{2[H1,H2]-[D1,D1⁢…⁢D1]}
The part of this expression {2[H1, H2]−[D1, D1. . . D1]} is equivalent to [H, −H], so that:

1α[2⁢Q1,2⁢Q2,…,2⁢Q2⁢N]-RVITMr·[D1,D1⁢…⁢D1]=RVITMr·[H,-H]
Recalling that Qk=αIkwe can also relate the real-valued intensity transmission matrix constants as RVITMPA=αRVITMrwhere RVITMPAis a row of matrix elements for the output photoacoustic signal (PA). Further re-arrangement of the equation yields:

1α[2⁢Q1,2⁢Q2,…,2⁢Q2⁢N]-1α⁢RVITMPA·[D1,D1⁢…⁢D1]=1α⁢RVITMPA·[H,-H](16)[2⁢Q1,2⁢Q2,…,2⁢Q2⁢N]-RVITMPA·[D1,D1⁢…⁢D1]=RVITMPA·[H,-H]

Since the matrix [D1, D1. . . D1] represents the input pattern with all micromirrors switched ON, and as that pattern was the first DMD pattern displayed in the experiment, the expression RVITMPA·[D1, D1. . . D1] corresponds to the measured output ultrasonic signal is Q1. Equation (16) can be re-written:

[2⁢Q1,2⁢Q2,…,2⁢Q2⁢N]-[Q1,Q1⁢…⁢Q1]=RVITMPA·[H,-H](17)[2⁢Q1-Q1,2⁢Q2-Q1,…,2⁢Q2⁢N-Q1]=RVITMPA·[H,-H]
Recalling that [H, −H]·[H, −H]T=[H, −H]−1, we can multiply equation (17) by [H, −H]Tto yield an expression for RVITMPA:
RVITMPA=[2Q1−Q1,2Q2−Q1, . . . ,2Q2N−Q1]·[H,−H]T(18)
The dimension of the first matrix on the right hand side of equation (18) is 1×2N and the dimension of the matrix [H, −H]Tis 2N×1. Accordingly the dimension of the matrix RVITMPAis 1×N. i.e. it comprises one row with N elements (recalling that there are N micromirror used by the DMD1016to form each input pattern). Thus each micromirror (or group of micromirrors acting as one micromirror) has a corresponding element in the matrix RVITMPA.

A positive element in the RVITMPAmeans that the corresponding micromirror of the DMD1016contributes positively (i.e. increases) the output ultrasonic signal, and therefore increases the light intensity at the target1006. Accordingly the aim of the characterisation step performed by the second computer1030is to determine a DMD pattern (which may or may not be not one of the DMD patterns1038displayed in the characterising process) that increases the number of elements in RVITMPAthat are positive. In some embodiments, the characterisation step may maximise the number of elements in RVITMPAthat are positive.

To further understand the physics of how the matrix elements of RVITMPAcan be made positive, we investigate below the effect of the nth micromirror of the DMD1016on the output ultrasonic signal by considering the contribution of the nth micromirror on the corresponding single element of the matrix RVITMPA, denoted RVITMPAn:
RVITMPAn=2Σk=12N(Qk−Q1)hnk=2 Σk=12NQkhnk(19)
where hnkis the element from [H, −H]∈(−1, +1) corresponding to the nth micromirror position (ON or OFF) in the kth DMD input pattern1038. In equation (19) the term −Q1can is cancelled out because [H, −H] has the same number of values +1 and −1 in each column. Substituting for Qk=αIkequation (19) can be written as:
RVITMPAn=2αΣk=12NIkhnk(20)

Based on conventional transmission theory the light intensity Ikat the target1006corresponding to the kth DMD input pattern1038can be expressed at Ik=|Σn=1NtnEkn|2where N is the total number of input micromirrors and tnrepresents the complex-valued transmission constants. Ekn∈(0,1) represents the light field at the nth micromirror (either ON or OFF) and is can be expressed in terms of hnk∈(−1, +1) as hnk=Enk−1.

Since the matrix element RVITMPAnis the sum of the output power over the total number (2N) input DMD patterns1038, equation (20) can be re-written as:
RVITMPAn=2α|Σn=1NtnEkn|2hnk
RVITMPAn=2α[|tn|2Σk=12N(Enk)2hnk+Σi=1,i≠nN|ti|2Σk=12N|Eik|2hnk+Σi=1,i≠nN(tnti*+tn*ti)Σk=12N(EnkEik)hnk+Ei=1,i≠nNΣj=2,j≠n,j>iN(titj*+ti*tj)Σk=12N(EikEjk)hnk]  (21)
where * denotes the complex conjugate operator.

By noting the structure of the Hadamard matrices and the resulting combinations of mirrors that are ON and/or OFF in the 2N input DMD patterns1038it is possible to simplify equation (21). In particular, during the characterisation step (i.e. sequential display of 2N input patterns) the nth micromirror is ON in half of the input patterns, i.e. N times. For the N patterns where the nth micromirror is ON (hnk=1, Enk=1) the jth micromirror (i≠n) is ON (Eik=1) for half of the N patterns, i.e. N/2 times, and the ith and jth micromirrors (i≠j≠n) are simultaneously ON (Eik=1, Ejk=1) half as many times again, i.e., N/4 times. Similarly when the nth micromirror is OFF (hnk=−1, Enk=0), the ith micromirror (i≠n) is switch ON for N/2 times, and the ith and jth micromirrors (i≠j≠n) are simultaneously ON (Eik=1, Ejk=1) N/4 times.

With this in mind, some of the terms in equation (21) reduce as follows:

∑k=12⁢N(Enk)2⁢hnk=N∑k=12⁢N❘"\[LeftBracketingBar]"Eik❘"\[RightBracketingBar]"2⁢hnk=0∑k=12⁢N(Enk⁢Eik)⁢hnk=N/2∑k=12⁢N(Eik⁢Ejk)⁢hnk=0
so that equation (21) becomes:

RVITMPAn=2⁢α[N⁢❘"\[LeftBracketingBar]"tn❘"\[RightBracketingBar]"2+N2⁢∑i=1,i≠nN(tn⁢ti*+tn*⁢ti)](22)RVITMPAn=α⁢N[tn⁢tn*+tn*⁢tn+∑i=1,i≠nN(tn⁢ti*+tn*⁢ti)]RVITMPAn=α⁢N⁡(tn⁢R*+tn*⁢R)RVITMPAn=α⁢NAR⁢An⁢cos⁡(θn-∅R)
where R is the output light field at the target1006when all of the micromirrors are switched ON. ØRand ARare the phase and amplitude of R: θnand Anare the phase and amplitude of tn, respectively. According to equation (22), when the phase θnof the transmission constant tnis within the range [ØR−π/2, ØR+π/2] the value of RVITMPAnis positive and increases the intensity of ultrasonic signal at the target1006, which corresponds to an increase in constructive interference of the light field for focusing. If that is so, the nth micromirror corresponding to RVITMPAnshould be switched ON. On the other hand if the phase θnof the transmission constant tnis outside the range [ØR−π/2, ØR+π/2] the value of RVITMPAnis negative and decreases the intensity of the ultrasonic signal at the target1006, which corresponds to a decrease in constructive interference of the light field. In this case, the nth micromirror corresponding to RVITMPAnshould be switched OFF.

Continuing with the description of the characterisation process, reference is made toFIG.11. At step S11-1a sequence of spatial modulations is generated. In this embodiment the sequence of spatial modulation comprises the input patterns1038, as described above. There may be 2N input patterns1038. At step S11-2the 2N input patterns1038may be displayed sequentially on the DMD1016, thereby modulating the intensity of each pulse of laser light. There be one pulse of laser light per input pattern1038. From the DMD1016the pulse of laser light is transmitted into the disordered medium, such as the diffuser1020. The laser light is scattered within the diffuser1020, and then reaches the target1006. By means of the photoacoustic effect a portion of the energy in the pulse of laser light is converted to ultrasound waves. At step11-3an output ultrasonic signal may be received by the detector1022and a time sequence of digital values may be generated by the DAQ1032. The output ultrasonic signal represents the amplitude of the ultrasound at the target1006and does not include any phase information. At step1-4the digital values may be processed by the computer and a single value generated that represents the amplitude of the ultrasound signal at the target1006. The computer1030may store the single value for subsequent use. Following display of the 2N patterns, the computer1030will store 2N values representing the ultrasonic signal amplitude measured for each pattern. The 2N values may be stored as an output matrix of dimension 1×2N.

At step S11-5the computer begins determination of the focusing spatial modulation pattern using the matrix RVITMPAwhere each element n of the matrix is determined using equation (19):

RVITMPAn=2⁢∑k=12⁢NQk⁢hnk
recalling that Qkis the ultrasonic signal amplitude stored by the computer1030for the kth pattern input by the DMD1016, and hnkis either +1 or −1, i.e. whether the nth mirror is ON or OFF in the kth pattern. In this way, each element of RVITMPAis determined by the following steps:

1. for the mirror n=1:a. for input pattern k=1 retrieve the value hnk, i.e. whether the mirror n was on or off for input pattern k;b. retrieve the ultrasonic signal value Qkand determine and store the product Qkhnk;c. repeat steps (i) and (ii) for input patterns k=2 . . . 2N;d. sum the products Qkhnk, and store as element no in the matrix RVITMPAthe magnitude of the sum and an indication whether the sum is positive or negative; and

2. set n=2 and repeat steps (i) to (iv).

The outcome of these steps is that the computer stores1030the matrix RVITMPAcontaining N real numbers, each being either positive or negative.

At this point the method ofFIG.1Imay proceed in two different ways. In some embodiments, the method proceeds from step S11-5, to step S11-5aand then to step S11-6. In other embodiments (described later in this document), the method proceeds from step S11-5, to step S11-5a, then to step S11-5band then to step S11-6.

At step11-5athe computer1030identifies whether each element of the matrix RVITMPAholds a positive real number or a negative real number. At step S11-6the focusing spatial modulation pattern is generated in which the computer1030sets state indicators of the focusing spatial modulation pattern according to the determination made in step S11-5. For example, when a positive real number of matrix RVITMPAis identified, a state indicator in the focusing spatial modulation pattern is set so that the corresponding micromirror is to be switched to an ‘ON’ state. When a negative real number of matrix RVITMPAis identified, a state indicator in the focusing spatial modulation pattern is set so that the corresponding micromirror is to be switched to an ‘OFF’ state.

The state indicators may be stored in any computer-processable data structure including, but not limited to, a list and an array. The data structure may comprise any suitable data type acting as the state indicator, such as the Boolean data type. For example, the focusing spatial modulation pattern could comprise state indicators in the form of binary values e.g. +1 or −1, or +1 or 0 to indicate whether each mirror of the DMD1016should be set to an ON state or an OFF state in the focusing spatial modulation pattern

It is recalled that, owing to the properties of Hadamard matrices and the way the input patterns have been constructed from them, each mirror n is ON for half the number of input patterns1038and OFF for half the number of input patterns. The function of equation (19) is to compare the sum of the ultrasonic signal amplitude values when the mirror n is ON (directing light to the diffuser1020and target1006) to the sum of the ultrasonic signal amplitude values when the mirror n is OFF (not directing light to the diffuser1020and target1006). If the comparison is positive, this indicates that when the mirror n is ON there is constructive interference of light at the target1006, indicating that light is focused in the plane of the target. On the other hand if the comparison is negative, this indicates that when the mirror n is ON there is negative interference of light at the target1006, indicating that light is not focused.

Thus the elements of the matrix RVITMPAenable the focusing spatial modulation pattern1036to be determined as described above which will tend to increase focus of light at the target1006. The focusing spatial modulation pattern1036comprises an indication whether each mirror (or a group of mirrors acting in unison) should be set to ON or OFF in order to increase focusing at the target1006. At step S11-6the focusing spatial modulation pattern1036can be sent to the transmitter1002for subsequent use in focusing light onto the target1006through the diffuser1020. A number of different communication mechanisms could be used to send the focusing spatial modulation pattern1036to the transmitter1002, such as any wired or wireless communication system e.g. WiFi, Bluetooth, cellular, electrical cables, etc.

Focusing Process

Once the focusing spatial modulation pattern1036(or any computer-processable equivalent) has been received by the transmitter1002it may be stored (in a computer memory remotely or locally accessible to the transmitter1002) for subsequent use.FIG.12shows steps of a focusing method performed by the transmitter1002. At step S12-1the transmitter may configure the DMD1016with the focusing spatial modulation pattern1036. In other words, the DMD1016is configured so that its mirrors match the focusing spatial modulation pattern1036. At step S12-2the transmitter1002may cause the laser1010to transmit laser light toward the DMD1016. Mirrors that are ON direct light toward the diffuser1020and mirrors that are OFF direct light away from the diffuser1020. At step S12-3there is an increased intensity of light at the target1006compared to a random pattern on the DMD1016for example.

It is noted that the laser light need not be the same as the laser light used to perform the characterising process. For example the laser light used at this step may be a series of individual pulses, an intensity-modulated continuous wave or as a continuous wave. Furthermore the wavelength of light may be slightly different between the characterising process and the focusing process.

In some embodiments, the characterising process may be repeated for different wavelengths, and one or more focusing spatial modulation pattern determined and stored for each wavelength. In use, the different focusing spatial modulation patterns may be used with the respective wavelength to provide a multi-spectral imaging function.

Results

FIG.13shows various results obtained when using random patterns on the DMD1016compared to a focusing spatial modulation pattern1036obtained using the characterisation process described above.FIG.13Ashows a graph of time (μs) versus normalised amplitude for the ultrasonic signal measured by the detector1022. Trace1050shows the average of 64 different ultrasonic signals measured when 64 different random DMD patterns were sequentially displayed on the DMD1016. Trace1052shows the average of 64 different ultrasonic signals measured when the DMD1016was configured with a focusing spatial modulation pattern1036. As can be seen, the trace1052(focusing spatial modulation pattern) provides a greater signal amplitude than the trace1050(random DMD patterns). This greater signal amplitude corresponds to a greater light intensity at the target1006, and in this experiment the signal amplitude when using the focusing spatial modulation pattern1038was about 6.7 times greater than the signal amplitude when using the random DMD patterns.

FIG.13Bshows the optical speckle pattern1054obtained using one of the 64 different random DMD patterns when the detector1022was replaced with a CCD camera. A grayscale key1056indicates the scale of the normalised light intensity of the speckle pattern1054. InFIG.13Bthe optical speckle pattern was spread across an area measuring approximately 320 μm by 320 μm at the target1006.FIG.13Cshows the optical speckle pattern1058obtained using the focusing spatial modulation pattern1036. As can be seen light was concentrated on a smaller area than the optical speckle pattern1054and contained a region1060of saturated pixels indicating high light intensity in the region1060.FIG.13Dshows the same optical speckle pattern1058but with the use of a neutral density filter to reduce the light intensity reaching the target. The improved light focusing effect when using a focusing spatial modulation pattern1036compared to a random DMD pattern (seeFIG.13B) can be clearly seen inFIG.13D.

FIG.13Eis a graph of distance (in μm) versus normalised ultrasonic signal amplitude (arbitrary units) taken along the ‘lateral profile’ dashed line shown inFIGS.13B and13D, andFIG.13Fis a similar graph but taken along the ‘vertical profile’ dashed line inFIGS.13B and13D. Trace1064shows that the ultrasonic signal intensity (the average of the 64 different signals) is generally low across the full width of the CCD camera, with a low peak toward the centre. Traces1066are Gaussian fits of the intensity data points (small open circles in the graphs) along the lateral and vertical profiles. The traces1066clearly show a pronounced peak in intensity at and near the centre of each profile, indicating the improvement given when using the focusing spatial modulation pattern1036. The dimensions of the optical focus (indicated by the full width at half maximum of the Gaussian fits (traces1066) were approximately 56 μm by 40 μm (which is consistent with the ultrasonic transducer focus size of approximately 49 μm diameter). The ratio of the average values of the ultrasonic signal intensity in the focal region1062before (FIG.1313) and after (FIG.13D) focusing was found to be 6.89, indicating a considerable increase in light intensity in the region1062.

FIG.14is an illustration of the time taken to perform an embodiment of the method comprising the characterising and focusing process described above. In this embodiment 8.192 input patterns ([H1, H2])1038were uploaded into the memory of the DMD1016. This time was not recorded as a cost as, in practice, this step can be completed in advance of the need to characterise a disordered medium and focus light at a target within or beyond the disordered medium. During a first phase1070of the method the 8,192 input patterns1038were displayed sequentially on the DMD1016at a rate of 47 kHz, taking a time of 175 ms. During the same time period the laser1010emitted pulses of light at the same rate (synchronised by the synchronisation device1018), at least a part of each pulse reflected from the DMD1016, through the diffuser1020and onto the target1006. The detector1022measured an ultrasonic signal corresponding to the light intensity received at the target1006, output an analogue signal which was received and processed by the DAQ1032into digital data and stored. During a second phase1072the digital data representing all of the ultrasonic signals was sent from the DAQ1032to the memory1034of the computer1030, taking approximately 66 ms. During a third phase1074the digital data received by the computer was processed as described above to determine the elements of the matrix RVITMPA, i.e. the focusing spatial modulation pattern1036to achieve improved focusing at the target1036. The third phase1074took approximately 58 ms. During a fourth phase1076the focusing spatial modulation pattern1036was uploaded to the t) memory of the DMD1016and displayed, taking approximately 7 ms. Thus the total time taken by the method to characterise and focus was 306 ms. It is noted that the most significant time cost is the first phase1070. Time savings could be made here with a spatial light modulator that is able to display the input patterns1036faster for example. Another way to speed up the characterising and focusing process would be to reduce the number of mirrors of (or used by) the DMD1016, although it is expected that lower light intensities would be achieved at the focal spot on the target1006. The time cost of the second phase1072and the third phase1074could be reduced by using a field programmable gate array (FPGA) for example.

Example Application: Focusing Light Through Multi-Mode Optic Fibres

Referring toFIGS.15A and15Ban apparatus generally identified by reference numeral1100is like the apparatus1000shown inFIG.10, with like reference numerals indicating like parts. For the sake of brevity, the description of these parts will not be repeated, but these features are incorporated into this embodiment by reference to FIG. and the description above.

The transmitter1102is like the transmitter1002except for the following. The first lens1112may be a tube lens with a focal length of 30 mm, the second lens1114may be a tube lens with a focal length of 75 mm, and the third lens1118may be a tube lens with a focal length of 50 mm. Furthermore an objective lens1121may be positioned after the third convex lens1118. The purpose of the objective lens1121is to direct light from the laser1110into a first end1123aof a multi-mode optic fibre1123. The multi-mode optic fibre is a disordered medium, and may be like any multi-mode optic fibre described herein. In some embodiments the multi-mode optic fibre has a diameter in the range several tens of micrometres to several millimetres. In some embodiments the length of the multi-mode optic fibre is such that the speckle decorrelation time is longer than the time taken to characterise and generate the focusing spatial modulation pattern. It is expected that such length may be anywhere in the range several millimetres to several tens of kilometres. In use the multi-mode fibre1123scatters light from the laser1110, generating an output speckle pattern (i.e. unfocused light) at a second end1123bof the multi-mode optic fibre1123.

In order to characterise the multi-mode optic fibre1123for focusing at or near the second end1123b, the receiver1104of the apparatus1110may comprise a CCD camera1125for receiving light from the laser1110through a second objective lens1127and a fourth tube lens1129(which may have a focal length of 100 mm). A purpose of the objective lens1127is to provide the receiver focus1025in front of the CCD camera1125at the second end1123bof the multi-mode optic fibre1123where light is emitted. Accordingly the objective lens1127is positioned so that its focal point at the second end1123bof the fibre1123

In use, the characterisation process is repeated as described above in conjunction with the apparatus1000. However instead of a single ultrasonic signal (which provides an indirect measurement of light intensity at the target1006) the receiver1104provides an image indicating a plurality of light intensity values at the second end1123bof the fibre over the field of view of the CCD camera1125. In particular the CCD camera1125comprises a plurality of output pixels, each of which provides a signal indicating the intensity of light received by that output pixel and this data is provided to the computer1130in the form of an image, or as a dataset comprising a time series of intensity data values for each output pixel. For each input pattern1038displayed on the DMD1116the computer stores the image (or the time series of intensity data values) taken by the CCD detector1125.

So instead of a single ultrasonic signal intensity value (Qk) determined for each input pattern1038(e.g. in the embodiment ofFIG.10), there is now a plurality of light intensity values recorded for each input pattern1038, with a single intensity value determined for each output pixel. In particular, once all of the input patterns1038have been displayed, the second computer1130has a set of k images (or k datasets), and determines k light intensity values for each output pixel. For each output pixel the computer1130uses the k intensity values in equation (19) to determine a focusing spatial modulation pattern1036for that output pixel. After processing the data, the second computer1130will have generated and stored a plurality of focusing spatial modulation patterns1036(in this embodiment there would be N focusing spatial modulation patterns). These focusing spatial modulation patterns1036are sent to the transmitter1102using the data communication network1140. The number of focusing spatial modulation patterns1036may or may not be the same as the number of output pixels of the CCD camera1125. For example, if the field of view of the second objective lens1127is circular, the number of output pixels may corresponding to the circular are of light reaching the CCD camera1125from the second objective lens1127. In another example, output pixels could be treated in groups by the second computer1130with a single light intensity value for each group being determined based on the light intensity values recorded for each output pixel in the group, per image k.

Once the transmitter1102has received the focusing spatial modulation pattern1036it can be loaded into the memory of the DMD1116for display. When in use (for example as a photoacoustic endoscope) the second end1123band a portion of the multi-mode optic fibre1123may be inserted into a body cavity for viewing of the interior of the human or animal body. In use, the receiver1102ofFIG.15Ais replaced with a receiver1102as shown inFIG.15B. The receiver1102ofFIG.15comprises an ultrasonic detector1122for receiving pulses of ultrasonic sound generated by the photoacoustic effect as the pulses of laser light leave the second end1123bof the fibre1123and are scattered and absorbed by tissue.

In order to generate an image for viewing by a person, the DMD1116is configured to display the particular focusing spatial modulation pattern1036for each output pixel in sequence. When displayed, a pulse of laser light (or series of pulses, or an intensity-modulated continuous wave) is transmitted toward the DMD1116, a portion is reflected into the fibre1123(by ON mirrors in the pattern) and a portion directed away from the fibre1123(by OFF mirrors in the pattern). The pulse is scattered inside the fibre1123. An ultrasonic signal is recorded by the detector1122, sent to the DAQ1132and on to the computer1130. The second computer1130may store a value (or a time series of values) representing the ultrasonic signal in such a way that it associated with the corresponding output pixel of the CCD camera1125used in the characterising process. This process is repeated using each of the focusing spatial modulation patterns1136in sequence so that the computer ends up storing a set of ultrasound signal values, each value (or a time-series of values) associated with a corresponding output pixel of the CCD camera1125. Since this data represents a 2D array of ultrasound transducers (corresponding to the 2D array of pixels of the CCD camera1125) it can be processed by the computer1130to produce a variety of different image reconstructions, including (but not limited to) a volumetric (3D) image generated using maximum intensity projection with depth-resolved colour map, although other image reconstruction algorithms are possible (such as a 3D rendering algorithm). Image reconstruction may be performed by image reconstruction software on the second computer1130or the ultrasonic data stored and transmitted to a remote computer (e.g. accessible over a data communication network) for image reconstruction.

In order to illustrate some of the steps in the process described above, and to show some reconstructed photoacoustic images obtained using this process, reference is made toFIG.16.FIG.16Ais a speckle pattern1150take by CCD camera1125before any characterisation of the fibre1123has been conducted. As can be seen, the light intensity reaching the CCD camera1125from the second end1123bof the fibre has an approximately random distribution.FIG.16Bis a DMD mirror configuration1152of a focusing spatial modulation pattern determined by the characterising process for focusing light onto one pixel of the CCD camera1125as described above. White regions1154of the display are ‘ON’ micromirrors and black regions1156of the display are ‘OFF’ micromirrors.FIG.16Cshows an image1158taken by the CCD camera1125when the focusing spatial modulation pattern ofFIG.16Bwas displayed by the DMD1116and illuminated with continuous wave (CW) laser light.FIG.16Dis a mesh plot1159of the pixel array of the CCD camera1125in the X and Y axis (horizontal plane) with light intensity along the Z axis (vertical). As can be seen light was focused onto a small region of pixels1160and not onto any other pixels. A strong peak1161of the light intensity inFIG.16Dshows the degree of focusing achieved. This demonstrates that, by using different focusing spatial modulation patterns, it is possible to focus laser light onto different groups of pixels of the CCD camera1125(or even onto to different individual pixels). For example,FIG.16Eshows that by using another focusing spatial modulation pattern light is focused onto a different group of pixels compared toFIG.16C.

The size of focus is determined by a number of factors including the wavelength of the light, the distance between the focus and the fire distal end, and the numerical aperture (NA) of the fibre. There are different cases: (1) if light is focused at distal fibre tip to a small distance in front of the fibre tip facet so that the NA of the output focusing light is larger than fibre NA, the size of the focus is determined by fibre NA. (2) If the light is focused to a large distance from the fibre tip, so that the actual NA of the output focusing light is smaller than the fibre NA, the size of focus is determined by the actual output NA.

FIG.16Fis a 3D photoacoustic image1162obtained using the apparatus ofFIG.15B. It shows a carbon fibre sample placed at the second end1123bof the fibre1123using water or an ultrasonic coupling gel for transmission of ultrasound. The ultrasonic data obtained was processed using a maximum intensity projection reconstruction algorithm to generate a 3D image. It is recalled that a focusing spatial modulation pattern1036has been determined for each output pixel of the CCD camera1125. In this example, the CCD camera1125has 400×400 pixels. When the apparatus ofFIG.15Bis used to generate an image, it is noted that it is not necessary to use all of the focusing spatial modulation patterns1036corresponding to all 400×400 pixels of the CCD camera1125. In this way it is possible to obtain different fields of view and images with different spatial resolutions and imaging frame rates. In particular it is possible to select a combination of patterns1036from all of the patterns1036to make up an image. For an image with a lower spatial resolution and a higher imaging frame rate, a combination of non-adjacent patterns1036can be selected, for example patterns1036corresponding to every 2nd, 3rd, 4thpixel, etc. are used. For a narrower field of view a combination of patterns1036can be selected in which patterns1036can be selected corresponding to a block of adjacent pixels of dimension x and y (where both x and y are smaller than then horizontal and vertical numbers of pixels in the CCD respectively). For exampleFIGS.16F and16Gshow 2D photoacoustic images1162and1164respectively that were generated using non-adjacent focusing spatial modulation patterns1036, in this case corresponding to every other pixel of the CCD1125thereby generating a 200×200 pixel image. InFIG.16Ha 3D photoacoustic image1166was generated using a combination of focusing spatial modulation patterns1036corresponding to a block of adjacent pixels of the CCD camera1125, in this case corresponding to a 200×200 block of pixels in the centre of the CCD. In this way the photoacoustic image1166has a narrower field of view than the photoacoustic images1162and1164inFIGS.16F and16G.

FIGS.16I-16Kspeckle patterns1268,1269and1270respectively showing focusing of light at a plurality of areas or regions1271of the CCD camera simultaneously (note that only one area or region is referenced in each ofFIGS.16I-16Kfor illustration).FIG.16Kis a longer exposure version ofFIG.16J. To focus light in multiple areas or regions simultaneously, each area or region is characterised as described above to determine a focusing spatial modulation pattern that will improve the focus of light in that area or region. Once a focusing spatial modulation pattern has been determined for each one, the patterns are processed further by computer in order to determine a multi-focusing spatial modulation pattern that will focus light simultaneously in the plurality of areas or regions. In particular, the corresponding state indicators in each pattern are summed (as binary values 0 and 1) and then averaged (mean) to determine a value between 0 and 1. The average value is then subjected to a threshold to determine whether the corresponding state indicator of the multi-focusing binary modulation pattern will set the corresponding micromirror to ON or OFF. For example the threshold could be 0.5, so that if an average value is greater than 0.5 the state indicator of the multi-focusing spatial modulation pattern is set to ON and if an average value is less than 0.5 the state indicator of the multi-focusing spatial modulation pattern is set to OFF. By processing the patterns in this way, a balance is struck between them that permits light to be focused in multiple areas or regions.

FIGS.17A-17Dare schematic drawings showing some example tip embodiments of the second end1123bof the multi-mode fibre1123that can be used inFIGS.15A and15B. InFIGS.17A-17Da region of biological tissue is generally indicated by reference numeral1172. The biological tissue1172may be a tissue sample taken from a human or animal (for example for histology) or may represent some part of the body of the human or animal in vivo. Generally, when the apparatus is used in vivo, the second end1123bwould be inserted into a cavity or passage of the human/animal body for internal examination. The second end1123bwould be moved into contact with a target such as tissue. Laser light suitable for causing a photoacoustic effect (e.g. pulsed or intensity-modulated continuous wave) is sent along the fibre1123to reach the second end1123b(utilising the focusing process described above). Light is scattered and absorbed by the target, with a portion of the light generating ultrasonic waves1174by the photoacoustic effect (shown schematically at1176). The ultrasonic waves propagate1174in all directions from the target, and so the ultrasonic detector1122can be positioned anywhere in or on the biological tissue1172providing the ultrasound waves can be detected. InFIG.17a first tip embodiment1170comprises a separate fibre1123(which, as in any embodiment described herein, may be a bundle of multi-mode fibres or a single multi-mode fibre) and ultrasonic detector1122. As shown the fibre1123is within the biological tissue1172and the ultrasonic detector1122is external of the biological tissue1172, but sonically coupled with it to facilitate detection ultrasonic waves.

InFIG.17Ba second tip embodiment1180comprises a fibre1123similar toFIG.17A. However, in this embodiment, the ultrasonic detector1122comprises a piezoelectric transducer of diameter 12.7 mm that is mounted at the distal end of an electrical cable1182, and adjacent the second end1123b. The electrical cable1182is moveable in unison with the fibre1123(for example they may be provided within the same housing (not shown), or otherwise attached to one another) so that both the second end1123band the ultrasonic detector1122may be brought into proximity of the part of the biological tissue1172to be inspected. In use, the generation and processing of ultrasonic signals is the same as embodiments described above.

InFIG.17Ca third tip embodiment1190comprises a fibre1123similar toFIGS.17A and17B. However, in this embodiment, the ultrasonic detector1122comprises an optically transparent piezoelectric transducer having the same diameter as the optic fibre mounted at the second end1123bof the fibre1123. In some embodiments, this optically transparent piezoelectric transducer may comprise an active layer of polyvinylidene difluoride or lithium niobate with a transparent electrical conductive material coated on both sides served as the electrodes. An example of a suitable optically transparent piezoelectric transducer is described in Chen. H et al, “Optical-Resolution Photoacoustic Microscopy Using Transparent Ultrasound Transducer” 20 Sensors (2019) 19, no. 24: 5470, the contents of which is incorporated herein by reference. Electrical connections for the detector1122may be routed around the outside of the fibre1123and contained withing a sleeve (not shown). In use, the generation and processing of ultrasonic signals is the same as embodiments described above.

InFIG.17Da fourth tip embodiment1200is similar to the second tip embodiment1180except that the electrical cable1182and detector112are replaced with a fibre optic ultrasound sensor1202. An example of a suitable fibre optic ultrasound sensor is described in Guggenheim, James A., et al. “Ultrasensitive piano-concave optical microresonators for ultrasound sensing.” Nature Photonics 11.11 (2017): 714-719, the contents of which is incorporated herein by reference. In use, the generation and processing of ultrasonic signals is the same as embodiments described above

FIG.18is an embodiment of an apparatus1210which is generally similar to the embodiment ofFIG.15Awhen using a second end1123hof the fibre1123as shown inFIG.17. However, this embodiment provides a dual imaging modality: in addition to the photoacoustic imaging modality, there is a light detector1212(e.g. photodetector, photodiode, CCD camera) and a dichromatic beam splitter1214(or a beam splitter and wavelength filter) for provision of fluorescence microscopy. After characterising the fibre1123as described above, laser light can be focused ‘pixel-by-pixel’ over a particular area the biological tissue1172into which fluorophores (e.g. fluorescent dye) have been introduced (e.g. by injection). Longer wavelength photons produced by the fluorescence are directed by the dichromatic beam splitter1214to the photodetector1212. By recording the fluorescence signal indicated by the photodetector1212at each pixel, an array of data can be compiled and this data processed to produce a fluorescence image of the target. At the same time, the ultrasonic detector1122can measure the ultrasound waves produced by the photoacoustic effect at the target, and be processed to produce an ultrasound image as described above. Different laser light may be used for photoacoustic and fluorescent signal excitation respectively.

Example Application: Optical Wireless Communication

Another application of the characterising and focusing process described above is optical wireless communication in free space which there is no line of sight between a transmitter and receiver, such as indoor optical wireless communication. In this application, data is transmitted using laser light at optical wavelengths in free space (e.g. through the air). In many situations there is no line of sight between the transmitter and receiver. If the light is directed at a surface (e.g. wall, ceiling) for reflection toward the receiver, reflection is diffuse and only a small portion of the transmitted laser light reaches the receiver. Furthermore, light that reaches the receiver has a time-varying speckle intensity pattern making high speed data transmission difficult or impossible owing to a poor signal to noise ratio.

Laser light may be reflected from one or more existing surface within buildings (e.g. walls, floors, ceilings, furniture) but, as mentioned above, these surfaces are types of disordered media and cause a high degree of scattering and diffuse reflection. By treating the communication link between the transmitter and receiver as the disordered media (similar to the diffuser1020inFIG.10and the multi-mode optical fibre1123inFIGS.15Aand B) the characterising process and focusing process may be performed to determine a focusing spatial modulation pattern to be used for that communication link. In this way, there is improved focus of light at the receiver with an increase in signal to noise ratio.

FIG.19is a schematic diagram of an embodiment of a non-line-of-sight optical wireless communication system (‘OWS’)1220installed inside a building1222. The building1222may be any structure in which a data communication system is useful. e.g. home, office, school, factory, shopping mall, airport, stadium, station, etc. The OWS1220may comprise a transmitter1224having a laser1226, such as an external cavity laser, and a spatial light modulator (‘SLM’)1228such as a DMD, like the DMD described in embodiments above. Control of the laser1226and the SLM1228may be performed by a first computer1229(e.g. in the form of an ASIC). The OWS1220may further comprise one or more receiver1230which may be in the form of a CCD camera. CMOS or array of photodiodes. The receiver may comprise (or may be adapted to communicate with) a second computer1231for performing functions like the functions of the computer1130described above. Objects1232(e.g. furniture, building structures) within the building1222block line of sight communication between the transmitter1224and one or more of the receivers1230. To overcome that, the SLM1228may be arranged to direct laser light toward a reflecting surface1234within the building. Generally the reflecting surface1234is a part of the building structure (such as a wall, floor or ceiling) is not a dedicated reflection device (e.g. mirror). As described above most such reflecting surfaces tend to produce a diffuse reflection of laser light. By reflecting laser light from an existing surface within the building1222fewer component parts are required for the OWS1220and it is easier to install and maintain.

Before the OWS1220is used to transmit data, the characterisation step (like that described in conjunction withFIG.10) is performed to determine a focusing spatial modulation pattern for the SLM1228in order to improve the focus of light at the receiver1230. In this case, the laser1226and SLM1228are controlled by the first computer1129to display the sequence of input patterns1038, thereby characterising the communication link. The second computer1231may be adapted to process signals representing the light intensity at the receiver1230as each input image is displayed on the SLM1228, and then utilise equation (19) above to determined the focusing spatial modulation pattern1036for the communication link. Once determined, the focusing spatial modulation pattern may be sent by the receiver1230to the transmitter using a different communication mechanism, such wireless (e.g WiFi), powerline, wired (e.g. Ethernet), and cellular. Once the transmitter1024has received the focusing spatial modulation pattern1036it may be stored and then displayed on the SLM1228and then data transmitted using the communication link. Once the focusing spatial modulation pattern is in use, there will be an improved focus of light at the receiver giving an increased light intensity and signal to noise ratio, and thereby improved data transmission. Since the characterisation step takes a comparatively short period of time (a few hundred milliseconds) may be that the first computer1129and the second computer1130are adapted to repeat the characterisation process periodically, non-periodically, under instruction of the user, and/or if the signal to noise ratio drops below a predetermined threshold.

FIGS.20A to20Dshow results obtained in an experimental set up of an OWS like the OWS1220.FIG.20Ais a reflected speckle pattern1240at the receiver1230when a random spatial modulation pattern was displayed by the spatial light modulator1228.FIG.20Bis the light field1242captured by the receiver1230in free space at a distance of approximate 0.6 m from the reflecting surface, before a focusing spatial modulation pattern is determined (i.e. a random selection of ON/OFF mirrors on the DMD).FIG.20Cis a reflected speckle pattern1244showing a spot of light1246focused at the receiver1230.FIG.20Dis the light field captured by the receiver1230in free space at a distance of approximate 0.6 m from the reflecting surface, after determination of a focusing spatial modulation pattern. As shown by these results, characterising the communication link enables the transmitter1024to improve the focus of light at the receiver1230. Some advantages of improved focus include increasing the signal to noise ratio at the receiver1230and an increase in the length over which the communication link is functional.

Improving Peak to Background Ratio (PBR)

As described above, micromirrors of the DMD1116are switched ‘ON’ when the phase θnof the transmission constant tnis within the range [ØR−π/2, ØR+π/2], where ØRis the phase of the light field at the focusing pixel. This increases the intensity of light at the focusing pixel. Although increasing (or even maximising) the intensity at the focusing pixel is desirable, it does not necessarily indicate best performance. An alternative measure of the effectiveness of the focusing of light via a disordered medium in imaging applications is the so-called Peak to Background Ratio (PBR). PBR may be defined as the ratio of the maximum intensity of light at a focusing pixel to the average light intensity of the background (i.e. light received by pixels other than the focusing pixel). This measure takes account of both the light intensity at the focusing pixel and the background light intensity, so that maximum PBR does not necessarily coincide with maximum light intensity at the focusing pixel.

In particular, as each switched ‘ON’ micromirror also contributes to the background light intensity, some of the ‘ON’ micromirrors may contribute even more to the light intensity of the background than to that of the focusing pixel. Thus, the highest PBR, which largely determines the achievable signal-noise-ratio (SNR) in imaging applications, is not necessarily achieved by this modulation approach. As such, switching ‘OFF’ those micromirrors that contribute significantly to the background light intensity can lead to a higher PBR

In the paper D. Wang. E. H. Zhou, J. Brake. H. Ruan, M. Jang. and C. Yang. “Focusing through dynamic tissue with millisecond digital optical phase conjugation.” Optica 2, 728-735 (2015), it was shown that the ensemble average of the peak output intensity with binary modulations using a DMD can be expressed as:

Ip=2⁢NA2⁢σ2(sin(2⁢φ2)+φ2⁢π)+π2⁢N⁡(N-1)⁢A2⁢σ2(sin⁢φπ)2(23)
in which the transmission constants tmnare assumed to obey a Rayleigh distribution and |tmn|2follows an exponential distribution |tmn|2˜e−1/2σ2with 2σ2as the ensemble average intensity of each element, φ is the upper bound of the absolute phase difference: 0≤Δ527≤φ, and N is the total number of input micromirrors, and A is the amplitude of the incident light field at each DMD mirror. The average background light intensity is expressed as:

Ip=2⁢NA2⁢σ2⁢φπ(24)
It can be shown that PBR can be determined as follows:

PBR=IpIb≈N⁢sin2⁢φ4⁢φ(25)
Equation (23) indicates that the maximum light intensity at the focusing pixel occurs when φ=π/2, whilst equation (25) indicates that PBR is a function of the absolute phase difference φ. In fact it was shown in the Wang paper that PBR can be maximised when φ=0.371π, although the Wang paper only considered phase information for producing a DMD pattern that maximised the PBR.

We have discovered that in some embodiments the PBR of the focusing spatial modulation pattern may be further improved by identifying those transmission constants tnfrom the matrix RVITMPAwith the highest real number values, and then switching ON the corresponding mirrors of the DMD (and leaving all other mirrors OFF). Furthermore, we have discovered that in some embodiments transmission constants tnmeeting this criterion can be identified as a top percentage or a top fraction of the all the transmission constants in the matrix RVITMPA. For example, it may be that the transmission constants tnin the matrix RVITMPAwith real number values lying within the top 30% are identified and the corresponding mirrors switched ON. Not only is there an improvement in the PBR, but the improvement may be achieved quickly compared to other methods, as will be explained below.

Referring again toFIG.11, in order to implement this aspect of the invention optional step S11-5bmay be performed. Once the matrix RVITMPAhas been generated at step11-5, the computer1030proceeds to step S11-5bin which the matrix RVITMPAis processed to identify n % of transmission constants with the highest real number values in the matrix (the percentage being of all the transmission constants in the matrix including positive and negative values, and where positive values are considered higher or greater than negative values). In one embodiment n % may be about 30% (although other percentages are envisaged as described below) Once those transmission constants have been identified the computer1030proceeds to step S11-6in which the focusing spatial modulation pattern is generated. To do that, the computer1030takes the top n % of transmission constants and may generate a new matrix which represents all of the ON/OFF state of the mirrors of the DMD1016. The matrix values are set so that mirrors are ON corresponding to the top n % of transmission constants, and OFF otherwise. This new matrix could contain only binary values e.g. +1 or −1, or +1 or 0 for that purpose. The new matrix is then used as the focusing spatial modulation pattern to be sent to the transmitter in step S11-7. It should be clear that the two methods that can be performed by the computer1030as illustrated inFIG.11will result in different focusing spatial modulation patterns sent to the transmitter. In particular, performing step S11-5then S11-6will generate a different focusing spatial modulation pattern than performing step S11-5, step S11-5band then step S11-6.

In order to demonstrate the improvement in PBR by generating the focusing spatial modulation pattern in this way, it was compared to a number of other DMD-based non-holographic algorithms that function to improve PBR. These other techniques, and an explanation of their implementation (both numerical simulation and experiment), are as follows:

(A) Real-Valued Intensity Transmission Matric (RVITM) Algorithm

This algorithm is as described herein, for example with reference toFIG.11. Results are presented below covering the two options inFIG.11, i.e. without step S11-5band with step S11-5b.

(B) Estimated TM-Based Algorithm (ETA) This algorithm is described in A. Drémeau. A. Liutkus, D. Martina, O. Katz, C. Schülke, F. Krzakala, S. Gigan, and L. Daudet. “Reference-less measurement of the transmission matrix of a highly scattering material using a dmd and phase retrieval techniques,” Opt. express 23, 11898-11911 (2015), and the algorithm. The algorithm was used for transmission matrix (TM) estimation and light focusing through a MMF. A total number of 6N random binary patterns with 50% micromirrors ‘ON’ were displayed on a DMD whilst the speckle intensities behind the MMF were captured by a camera. A Bayesian phase retrieval algorithm (described in A. Drémeau and F. Krzakala. “Phase recovery from a bayesian point of view: the variational approach,” in 2015IEEE International Conference on Acoustics, Speech and Signal Processing(ICASSP). (IEEE, 2015). pp. 3661-3665)] was then used to calculate the complex-valued TM from intensity-only input and output pairs via iterative optimizations. This algorithm was chosen as it benefits from a moderate computational cost. An open-source script of the phase retrieval algorithm described in the Drémeau paper was used in this work. As this algorithm provides phase values of a row of the TM corresponding to the mthoutput mode, optimal DMD patterns with both Re(tmn>0) and |arg(tmn)|<0.371π were used for focusing and their performance were compared.
(C) Conditional Probability-Based Algorithm (CPA)

The CPA is described in detail in T. Zhao. L. Deng, W. Wang, D. S. Elson, and L. Su, “Bayes' theorem-based binary algorithm for fast reference-less calibration of a multimode fiber,” Opt. express 26.20368-20378 (2018). Like the ETA, a total number of 6N random binary patterns were used as inputs, whilst the intensities of speckles at output of a MMF were captured. There were three steps involved in the generation of an input DMD pattern for focusing. First, an intensity threshold was used to divide the output intensities into two groups: a ‘focusing’ and a ‘non-focusing’ group. Second, Bayes' theorem was used to calculate the conditional probability that switching ‘ON’ each micromirror leads to light focusing at the target output position (‘focusing’ group). Finally, a threshold was used to produce the optimal DMD pattern for light focusing through the MMF by switching ‘ON’ micromirrors a conditional probability higher than a threshold. As reported in the Zhao paper to maximise the PBR a first threshold was set as the 80 percentile of all intensities at the target position, while a second threshold was set as the median value of all probability values.

(D) Genetic Algorithm (GA)

The method for implementing the GA for light focusing through a diffuser was described in D. B. Conkey, A. N. Brown, A. M. Caravaca-Aguirre, and R. Piestun, “Genetic algorithm optimization for focusing through turbid media in noisy environments.” Opt. express 20, 48404849 (2012) and in X. Zhang and P. Kner, “Binary wavefront optimization using a genetic algorithm,” J. Opt. 16, 125704 (2014). We used the same process but employ the PBR of the output light field as the feedback to be maximised. First, a total number of 20 random binary patterns with approximate 50% micromirrors ‘ON’ were used as the 1BC generation population. Each DMD pattern is considered as the chromosome of an individual, and the status of a micromirror was considered to be a chromosome code (‘1’ for ‘ON’ and ‘0’ for ‘OFF’). Output speckles intensities were recorded when displaying these binary patterns on a DMD and their PBR values were compared. Individuals in the 1BC generation population were ranked according to their corresponding PBRs in the outputs. Then individuals with larger PBRs were assigned larger probabilities to be selected as parents to produce the next generation population by crossing the parent chromosomes with a constant cross rate. Mutation was also introduced by randomly switching a small number of chromosome codes with a mutate rate to avoid locally optimal solutions. In the next step, the new generation was ranked according to the resulting PBRs and produced the next generation patterns through the aforementioned progress. After a large number of iterations, the chromosome codes leading to a high PBR were saved in a DMD pattern. In our numerical simulations, 30,000 generations were implemented with the cross rate setting to be 0.6 and the mutate rate 0.02. In our experiments, 4.000 generations were implemented, the mutate rate was set to be 0.1*e−G/600+0.02 to speed up the optimisation, where G is the index of the generation.

Numerical simulations were implemented in MATLAB to investigate the performance of different algorithms (A)-(D). A complex-valued simulated TM was generated with random phases and amplitudes following a uniform and a Rayleigh probability density function between 0 and 2π, and 0 and 1, respectively. The number of input micromirrors (N) was set to be 32×32 while the number of output pixels M was set to be 64×64. Output light intensities were calculated based on the simulated TM via Em=|Σn=1NtmnEn|2, which were fed to those algorithms for comparison. The resulting PBR was calculated as the ratio of the intensity at the focusing pixel over the average intensity in the background, for the evaluation and comparison of the focusing performance with different algorithms.

Numerical Simulation

Numerical simulations were performed with a MATLAB program. In order to generate the output speckle patterns134, a complex-valued transmission matrix (TM) with 8192 output pixels and 1024 input pixels was generated. The phases and amplitudes of the TM were randomly generated to obey uniform and Gaussian distributions between 0 and 2π, and 0 and 1, respectively. The TM can be used as a ground truth (or reference TM) for comparison with results achieved with DMD-based algorithms. As the TM elements follow a circular Gaussian distribution with phases obeying a uniform distribution in [−π, π], modulating phases of output light fields coming from all input modes to an ideal phase Ø=0 or to Ø=ØRleads to approximately the same constructive interference at the output position. For the case of comparison with the RVITM algorithm which employs ØRas the ideal phase, the focusing condition was chosen as θmn−ØR<φ rather than |θmn|<φ. Both φ=π/2 and φ=0.371π were used as the upper boundary for producing the DMD pattern to focus light in the complex-valued TM-based approaches.

Experimental Setup

The experimental setup was very similar to the apparatus ofFIG.15A(except a data communication system1140was not required). The light source was a collimated diode-pumped solid-state laser module (532 nm, 4.5 mW, CPS532. Thorlabs. New Jersey. USA). After beam expansion through two tube lenses (AC254-030-A-ML; AC254-075-A-ML, Thorlabs, New Jersey. USA), the light was spatially modulated using a DMD (DLP7000, 768×1080 pixels, Texas Instruments. Texas. USA) and then projected onto the proximal tip of a MMF (105 μm, 0.22 NA, 1 m. M43L01, Thorlabs, New Jersey. USA) via a tube lens (AC254-050-A-ML. Thorlabs. New Jersey. USA) and an objective (20×, 0.4 NA. RMS20X, Thorlabs, New Jersey. USA). The light illuminated area on the DMD included 32×32 independent micromirrors. An objective (20×, 0.4 NA. RMS20X. Thorlabs. New Jersey, USA) and a tube lens (AC254-0100-A-ML. Thorlabs. New Jersey, USA) were used to magnify the output light beam before it was captured by a complementary metal-oxide-semiconductor (CMOS) camera (C11440-22CU01. Hamamatsu Photonics. Shizuoka. Japan) with a frame rate of 200 frames per second (fps) for MMF characterisation.

Results

In the RVITM-based algorithm (A), each micromirror has a respective phase difference (θmn−ØR) and a transmission constant (rvitmn).FIG.21Ashows a graph1250of phase difference versus the normalised value of the transmission constants rvitmn1252obtained via equation (22) for each micromirror in the numerical simulations. The distribution of the transmission constants rvitmnfor different phases differences has an envelope of a cosine function1254. As the phases of the transmission constants rvitmnobey a uniform distribution in [−π, π], half of the transmission constants rvitmn1256have positive values and their corresponding phase differences (θmn−ØR) are in the range of [−π/2, π/2]. So, switching ‘ON’ the corresponding micromirrors leads to a linear sum of the transmission constants rvitmn1256with the constructive interference at the target output position and hence produces the light focusing. To further improve the PBR (FIG.11, step S11-5b), all values of the transmission constants rvitmn1252were ranked in descending order according to real number value, different proportions (P) of total number of the transmission constants rvitmnwith larger values were selected, the corresponding micromirrors were switched ‘ON’, and then the resulting PBR determined. By selecting a group of micromirrors with the highest real number values of the transmission constants rvitmn1252, not only is the phase difference of each micromirror taken into account, but also the amplitude contributed to the PBR by each micromirror at the target. In some embodiments this may result in a different selection of micromirrors than using phase difference only, the different selection producing an improved PBR compared to the selection by phase only.

FIG.21Bshows a graph1260of the proportion of ON micromirrors P versus normalised PBR. A first line1262shows the results obtained by simulation and a second line1264shows the results obtained by experiment. For example, when P=50% the micromirrors with transmission constants rvitmnabove a first threshold1258(seeFIG.21A) of rvitmn=0 are switched ‘ON’. When P=30% the micromirrors with transmission constants rvitmnabove a second threshold1259(seeFIG.21A) of rvitmn=0 are switched ‘ON’.

As seen inFIG.21Bwith both experiments and simulations the PBR was greater than one over the approximate range 20%≤P≤50%. In other words, switching ‘ON’ micromirrors with transmission constants rvitmndistributed in the dome region above the second threshold1259and below the cosine function1254inFIG.21Aproduced the maximum PBR with the RVITM-based algorithm. To investigate the maximum range of P for which the PBR is greater than one, reference is made toFIG.21Cwhich shows the results of numerical simulations. In particular a graph1265shows a line1266indicating variation of PBR with P. Lines1267and1268show the values of P for which PBR is greater than one, namely about 15%≤P≤50% (all figures ±1%). In practice it may be that the precise percentage range is slightly different to this range, but the experimental data ofFIG.21Bindicates this range is broadly justified. In some embodiments it is envisaged that the number of ON micromirrors having the greatest real-valued transmission constants rvitmnmay be selected to fall within any one of the following ranges: 15%≤P≤50%, 20%≤P≤45%, 22%≤P≤40%, 25%≤P≤35%, 25%≤P≤33%, and 29%≤P≤31%. In some embodiments maximum PBR can be obtained by switching ON the 30% of micromirrors having the highest real-numbered transmission constants. In some embodiments, it may be desirable to increase both intensity of light at the target and improve the PBR. That can be achieved by selecting a percentage of ON micromirrors between the peak PBR value (indicated by line1269) and highest percentage of micromirrors where PBR is equal to one (indicated by line1268). In an embodiment that range is 30%≤P≤50%. In other embodiments it is envisaged that the number of ON micromirrors having the greatest real-valued transmission constants rvitmnmay be selected to fall within any one of the following ranges: 30%≤P≤45%, 30%≤P≤40%, 30%≤P≤35%, 30%≤P≤33%, and 30%≤P≤31%.

One advantage of at least some methods according to the invention is that selecting micromirrors to switch ON according to the real number value of a real-valued transmission constant means that micromirrors are selected based on both amplitude and phase difference (of the corresponding complex-valued transmission matrix), rather than just phase difference only. This produces a different percentage of ON micromirrors than the same percentage of ON micromirror selected by phase difference only for example.

FIG.22is a table showing a performance comparison of the algorithms (A)-(D). The notes to the table are as follows:a, b with reference TM and estimated TM, the PBR refers to two different conditions for determining the patterns for focusing: |arg (tmn)|<π/2 (top) and |arg (tmn)|<0.371π(bottom).c in RVITM (D), the PBR refers to P=50% (top) and P=30% (bottom).S simulationsE experimentsd in experiments, the time cost includes the time for DMD pattern display during the fibre characterisation and the computation time for producing the optimal DMD patterns, while in simulations, only the latter was included in the time cost.

For the ease of comparison, the rvitmnand θmn−ØRvalues were calculated for all micromirrors using the ground truth TM (described above in conjunction with the numerical simulation), whilst those corresponding to the switched ‘ON’ micromirrors determined by the different algorithms am shown inFIGS.23A-23F.FIGS.23A-23Eshows a graphs1275,1280,1290,1300and1310respectively of phase difference versus the values of the transmission constants rvitmn.FIG.23Fis a graph1320of index of generation versus PBR for the genetic algorithm (D) as will described in the following. In the case of the ground truth (or reference) TM graph1275, all of the data points represent micromirrors that were switched ‘ON’ corresponding to values of rvitmnabove the first threshold1258of rvitmn=0. i.e. a phase difference φ in the range [−π/2, +π/2]. Data points shown with open circle symbols in regions1272represent additional micromirrors that were switched ‘OFF’ when the second threshold1259was used. i.e. the phase difference φ in the range [−0.371π, +0.371π]. Referring to the table shown inFIG.22a PBR value of 208.4 was achieved in the latter case, which is 13.9% higher than a PBR value of 183.0 achieved in the former case.

Still referring to the table inFIG.22the PBR values achieved with the ETA (B) were slightly smaller than that achieved with the reference TM (170.4 for phase difference φ in the range [−π/2, +π/2] and 190.3 for phase difference φ in the range [−0.371π, +0.371π], respectively). This slight reduction of PBR values can be attributed to the errors of the TM calculation, which is indicated by the different micromirrors chosen as shown by the rvitmnvalues inFIG.23B.

The CPA (C) produced lower PBR of 156.7 compared to the estimated TM method (see alsoFIG.23C).

With the RVITM algorithm (A), the first scenario with rvitmnvalues above the first threshold1258(i.e. P=50% and a phase difference φ in the range [−π/2, +π/2]) resulted in the same DMD pattern (as represented by all of the data points inFIG.23D) for light focusing as that obtained from the reference TM and hence the same PBR value of 183.0. In the second scenario when the second threshold1259was applied (i.e. P=30% and phase difference φ in the range [−0.371π, +0.371π]), those micromirrors that are represented by lighter shaded data points in the region1302were switched ‘OFF’, and the PBR value increased to 228.2, which is even higher than that achieved with the reference TM with phase difference φ in the range [−0.371π, +0.371π] (PBR=208.4).

The PBR value achieved with the GA (D) reached the highest value of 239.7 among all the algorithms but after the evolution of 30,000 iterations (seeFIGS.23E and23F).

To compare the focusing speed of different methods, the average times taken for different methods to compute an optimal DMD pattern for focusing over 100 output locations (e.g. different pixels of a CCD) were obtained on a PC with 2.3 GHz Dual-Core Intel Core15(see the table inFIG.22). Although providing the highest PBR, the GA (D) had the longest computation time of 400 s. In comparison, with 200 iterations the computation time for the ETA (B) and the CPA (C) was 15 s and 4 s, respectively. The RVITM-based algorithm (A) calculated the rvitmnvalues for all output positions at the same time (i.e. in parallel), while for each output position the computation time for focusing was 7.5 ms.

Performance Comparison in Experiments

FIGS.24A-24DFIG.5are plots1330,1340,1350and1360respectively of the output light intensity patterns of an optical focus generated at the distal end of the MMF using algorithms (A)-(D). The scale bars1332represent a dimension of 10 μm. Referring again toFIG.22the PBR of the ETA (B) was estimated to be 58.5 when switching ‘ON’ micromirrors with Re(tmn) >0, and increased to 64.2 with |arg(tmn)|>0.371π. The PBR achieved with the CPA (C) was the lowest, the GA (D) produced the largest PBR of 91 among all the algorithms after 4000 iterations (seeFIGS.24D and24E, the latter being a graph1370showing index of generation versus PBR). The relationship between the achieved PBR and the proportion of micromirrors P with the RVITM-based algorithm (A) is shown and described above with reference toFIG.21B. The highest PBR value of 89.4 was reached at P=30% (with the output light intensity shown inFIG.24C) compared to a PBR of 79.9 at P=50%. The profiles of foci are shown inFIG.24Fwhich is a graph of position (or distance) versus intensity (in arbitrary units). The focus achieved with different methods has a uniform diameter of 1.7 μm. Unlike the results obtained in the numerical simulation, the time cost in experiments includes the time for DMD pattern display. As a result, the GA (D) had the longest time of 2.2 h, while the time cost was 46 s, 36 s and 10.3 s for the ETA (B). CPA (C) and the RVITM (A) algorithms, including both input pattern display and data processing run-times.

In summary embodiments of the invention can improve the PBR obtained with the ETA (B) by switching ‘ON’ a group of micromirrors for which the corresponding real-valued transmission constants have the highest real number values and leaving ‘OFF’ the remaining micromirrors (i.e. those with lower values). The PBR in such embodiments is improved compared to PBR achieved with the commonly used criteria |arg(tmn)|>π/2 for DMD-based waveform shaping. With the RVITM-based algorithm (A), the maximum PBR was achieved at P=30%, although it is noted that various ranges of P around this value are expected to produce similar results. Although the precise phase values are not achievable with RVITM (A), we have found that the transmission constants rvitmnencode (or represent) phase and amplitude information. This enables a higher PBR than that of the ETA (B). Similarly, the higher PBR achieved with the RVITM algorithm (A) also suggests that the amplitude information should be considered for achieving the maximum PBR with a DMD. The GA (D), which does not require the phase and amplitude information of the TM as prior knowledge, achieved the highest PBR via a large number of iterations as expected but with a high time cost. Interestingly, a small number of micromirrors with negative rvitmnwere switched ‘ON’ with the GA for focusing (seeFIG.23E). This indicates that although these micromirrors had negative contributions to light intensity at the focusing position, they have substantially reduced the intensity of the background, leading to a higher PBR compared to those methods considering only micromirrors with positive rvitmn(seeFIGS.23A,23B,23D). To further increase the PBR, more independent micromirrors could be used in the future, however, at the expense of the focusing time.

Embodiments of the invention using employing improved PBR may be used in (but not limited to) any of the embodiments and applications described with reference toFIGS.10,15A,15B,17A-D,18,19and25(described below).

Referring toFIG.25Aan apparatus for discovering an out of sight object is generally identified by reference numeral1390. The apparatus1390comprises a transmitter1392and a receiver1394. The transmitter1392may comprise a first computer1396, laser1398and a spatial modulator1400, similar for example to the transmitter1102described with reference toFIG.15A. The receiver1394may comprise a camera1402(e.g. CCD camera, CMOS or array of photodiodes) and a second computer1402for receiving and processing data from the camera1402. The transmitter1392and receiver1394may have access to a data communication network1406so that they may exchange data with one another. The data communication network1406is not a part of the apparatus1000, although the transmitter1002and receiver1004may have one or more interface (not shown) respectively for sending and receiving data over the data communication network1406. The data communication network may take a number of forms, including by not limited to, wireless systems (e.g. WiFi, Bluetooth, cellular) and wired systems (e.g. Ethernet, power-line), or any combination of these.

It is noted that the apparatus1390is shown in use inFIG.25Awith a obstruction1408and a reflecting surface1410, all shown in plan view. The obstruction1408may be any temporary or permanent object or structure that prevents or obscures a direct line of sight between the transmitter1392and receiver1394. The reflecting surface may by a disordered medium that produces diffuse scattering and reflection of light, such as a surface of a building. In one embodiment the transmitter1392may be part of a vehicle (such as a car, van or truck), the obstruction may be the corner of building, and the reflecting surface may be another building on the opposite side of the street.

In a first step (shown inFIG.25A), a transmission path1412between the transmitter1392and receiver1394may be characterised in the same way as described with reference toFIG.15Afor example, to determine a focusing spatial modulation pattern for the spatial light modulator1400for each pixel or group of pixels of the camera1402. For the sake of brevity the description of these steps is not repeated here, but the steps from the description ofFIG.15Aare incorporated into the present embodiment. It may be necessary to temporarily place the camera1402behind the blocking object1408for the purposes of characterisation. Once the characterisation process is complete, the receiver1392may transmit the plurality of spatial modulation patterns (or any computer processable equivalent of them) to the transmitter1392for storage (either locally or remotely. e.g. a remote server accessible via the internet) and subsequent use. The camera1402may be removed at this point.

Referring toFIG.25Bthe transmitter (mounted on or part of the vehicle) may return to the same location at a future point in time. Once at the location, the vehicle may establish whether them is an object1414out of direct line of sight, obscured by the obstruction1408. The object1414could be another vehicle (stationary or moving), a cyclist or pedestrian, i.e. any object that could temporarily occupy the same location or a location close to the location of the camera1402during the characterising process. The transmitter may use the first computer1396to retrieve the plurality of focusing spatial modulation patterns for the transmission path1412. These focusing spatial modulation patterns may be used in a similar way as described with reference toFIG.15Bto generate an image, i.e. each focusing spatial modulation pattern is applied to laser light from the transmitter and spatially modulated light sent toward the object1414by reflection from the reflecting surface1410and focusing in an imaging plane previously occupied by the camera1402during the characterisation process. If the object1414is present, light is scattered and reflected by it. Some of the light will be reflected back towards the reflecting surface1410. The reflecting surface1410will scatter that reflected light and a small portion may reach a photodetector1416. All that is required is for the photodetector1416to detect some light reflected by the object1414and the reflecting surface1410, and for the first computer1410to store the reflected light intensity In some embodiments, the wavelength of the laser light from the transmitter may have a wavelength greater than about 800 nm and may be within, but not limited to, the near-infrared portion of the electromagnetic spectrum (e.g. about 900 nm to about 2,000 nm). A wavelength filter may be used to reduce the light intensity within a certain range (e.g. 400-800 nm) to distinguish the ambient light from the reflected light from the object.

This process is repeated for each focusing spatially modulating pattern (each one corresponding to a pixel of the camera1402). If an image is to be displayed to human (or otherwise required by an image analysis algorithm), the recorded intensity values from the photodetector can be used to generate an image in a way similar to that described in conjunction withFIG.15B. Alternatively the first computer1410may process the recorded intensity values to determine if the object1414is present in the location obscured by the obstruction1408. Once the first computer1410has made the determination, this may be used to cause a subsequent action or a further decision. For example, the first computer1410could cause a warning to be communicated to a user of the vehicle that an object is present behind the obstruction1408. Alternatively, the first computer1410could cause the vehicle to wait until the object1414has gone, or to turn around the corner if no object is detected. This may be useful in embodiments where the vehicle is an autonomous vehicle for example: the outcome of the determination may be used by the vehicle to make further control decisions (e.g. go, wait, etc.).

Whilst the embodiments above have been described with reference to silicon-based computer processors, it is envisaged that one or more of these may be replaced with an optical computing processor. For example, it may be that the process of matrix inversion may be performed using an optical network. One such example is described in Wu, K., et al. “Computing matrix inversion with optical networks” (see Wu K, Soci C, Shum P P. Zheludev N I. Computing matrix inversion with optical networks. Optics express, 2014 Jan. 13; 22(1):295-304, which is incorporated herein by reference).

Some embodiments have been described with reference to electromagnetic radiation at visible wavelengths. The invention is not limited to one or more wavelength in this band. Other embodiments may use electromagnetic radiation in other portions of the spectrum, for example, but not limited to ultraviolet and near-infrared.