Neuromorphic Shack-Hartmann wave normal sensing
Chutian Wang Shuo Zhu Pei Zhang Jianqing Huang Kaiqiang Wang Edmund Y. Lam





Department of Electrical and Electronic Engineering, The University of Hong Kong, Pokfulam, Hong Kong SAR, China

*elam@eee.hku.hk 


The Shack-Hartmann wavefront sensor is widely employed in adaptive optics systems to measure optical aberrations. However, simultaneously achieving high sensitivity and large dynamic range is still challenging, limiting the performance of diagnosing fast-changing turbulence. To overcome this limitation, we propose neuromorphic Shack-Hartmann wave normal sensing (NeuroSH). NeuroSH is a unifying framework that harnesses the computational neuromorphic imaging paradigm to extract the high-dimensional wave normal from temporal diversity measurements. Both numerical analysis and experimental verification demonstrate the feasibility of NeuroSH. To the best of our knowledge, the proposed NeuroSH is the first scheme to surpass the optical dynamic range limitation under challenging dynamic scenarios, thereby advancing ultra-fast turbulence mitigation technology for cutting-edge imagers.





§ INTRODUCTION

Adaptive optics (AO) is a technology extensively applied in various fields, including astronomy <cit.> and microscopy <cit.>. Its primary objective is to enhance the performance of high-contrast imagers by integrating a wavefront sensor (WFS) for wavefront detection with a deformable mirror for wavefront compensation. Among the available WFSs, the Shack-Hartmann (SH) WFS stands out as the most widely adopted owing to its robustness and high optical efficiency. A SH WFS utilizes a microlens array (MLA) to sample the wavefront gradient as a combination of spots on the sensor plane, with their displacement linear to the corresponding local gradient. Hence, by comparing the registered reference position of SH spots and the aberrated one in the updated measurement, the map of the wavefront gradient can be recovered <cit.>.

Due to the fixed physical parameters in the utilized MLA and the pixelated sensor, the maximum and minimum magnitudes of the measured local tip/tilt term in the SH WFS are coupled <cit.>. The maximum value corresponds to the WFS dynamic range (DR), representing the range of detectable wavefront aberrations. While the minimum reflects the resolution of the sensor to detect subtle variations in the wavefront. This trade-off imposes significant constraints on the current design specifications, as achieving either high sensitivity or large DR becomes a compromise depending on the desired level and scale of aberrations of interest. However, current advanced AO systems pose a significant challenge if they necessitate the simultaneous achievement of large DR and high sensitivity. For instance, extreme AO systems aim to tackle fast-changing, large-scale, highly aberrated atmospheric turbulence under photon-starved and adverse conditions <cit.>. Consequently, extending the DR while preserving its sensitivity is a highly sought-after endeavor in the field.



Numerous DR enlargement methods have been proposed to enhance wavefront sensing, which can be broadly categorized into two groups: hardware-based and algorithm-based. The former requires considerable hardware modifications to the utilized MLA mask <cit.>, leading to configuration bulkiness, manufacturing complexity, or higher expenses. Alternatively, it involves additional lateral or axial scanning strategies at the cost of sacrificing a substantial amount of temporal resolution <cit.>. On the other hand, algorithm-based methods, such as sorting <cit.>, unwrapping <cit.>, Fourier modulation <cit.>, optical flow-tracking <cit.>, and deep learning <cit.>, offer advantages including cost-effectiveness and plug-and-play compatibility, without requiring additional modifications to its classical setup. However, their principle primarily revolves around avoiding the overlap of SH spots. This strategy imposes inherent restrictions on their DR enlargement capability, confining it within the limits of the optical dynamic range dictated by the wavefront curvature <cit.>.

To avoid the inherent ambiguity in the optical dynamic range without sacrificing temporal resolution, one prospective strategy is to leverage previous time-series shots. However, the conventional SH WFS employs an active pixel sensor (APS) that records the entire scene through a synchronous response mechanism, accompanied by large uninformative portions. The perceptual and statistical redundancies associated with the frame-based paradigm would require a significant waste of hardware storage and computational resources. In contrast, the emerging computational neuromorphic imaging (CNI) paradigm <cit.> utilizes a bio-mimic dynamic vision sensor (DVS) as the data acquisition terminal. The DVS asynchronously responds to logarithmic intensity changes rather than capturing the complete scene <cit.>. It offers microsecond temporal resolution and a high dynamic range exceeding 120dB <cit.>, which can potentially enhance the SH WFS performance <cit.>.


In this paper, we propose neuromorphic Shack-Hartmann wave normal sensing (NeuroSH), which is a unifying CNI-based wavefront framework for directly perceiving the high-dimensional wave normal. NeuroSH is the first to surpass the current optical dynamic range limitation in a classical SH optical setup. We emphasize the advantage of CNI in terms of selective responses and low data redundancy, thanks to its asynchronous triggering mechanism. This enables the acquisition of a highly efficient data modality that represents the temporal diversity measurements within a complete AO servo-loop period. NeuroSH is therefore designed to perceive the spatio-temporal trajectory that approaches the temporal increment of wave normal.  The verification of NeuroSH is established on an adaptable SH setup for flexible MLA generation and dynamic phase ground truth overlaying. We further demonstrate its application in instantaneous flame front diagnosis, showcasing its substantial advantage over the DR limitations in ultra-fast scenarios.









§ PRINCIPLE




 §.§ SH WFS Geometry and the DR limitation


The wavefront is a surface of equal phase in Euclidean space ℝ^3. As depicted in Fig. <ref>(a), considering the Cartesian coordinates of the wavefront plane as ζ = (ζ,η), we can express the distribution of the wavefront on this plane as a 2D function denoted by W(ζ). The wave normal, which is a 3D vector perpendicular to the wavefront, indicates the propagation direction of the light wave. Given θ and ϕ as the polar and azimuth angles for an incidence point, we can represent the wave normal as (θ, -1)^𝖳, where θ = (θ, ϕ)^𝖳=∇ W(ζ) denotes a 2D angular incident vector.


A classical SH WFS comprises a MLA and a pixelated sensor placed at its focal plane, with its coordinates denoted by x = (x, y). The basic working principle of the SH WFS involves positioning the MLA at the desired wavefront plane to divide it into individual patches using the lenslet boundaries. In each patch, the regional wavefront is linearized to keep propagating along the direction of wave normal and is finally focalized to form a SH spot. The local principal ray vector is aligned with wave normal and can be therefore represented as N = (θ, -1)^𝖳f, with the location of the SH spot being the end point of N on the sensor at the focal plane. Consequently, the linear relationship between the spot displacement and the angular incident vector can be written as

    Δx=x - x_ref=θ f,

where f represents the focal length, and x_ref serves as the registered position at initial stage. If assigning a plane wavefront as a reference, it corresponds to the lenslet's vertex within each patch.






As depicted in Fig. <ref>(b), the fixed focal length and pitch size of the MLA used in SH WFS impose limitations on the detectability of the SH spots, referred to as the dynamic range (DR) of the SH WFS. Previously, two commonly adopted types of DR have been identified <cit.>, they are inherited in this work as Type 1 and Type 2 DR, respectively.  Herein, we further enlarge the DR to Type 3 from the perspective of sensing the wave normal. The following  explicates each of the three types of DR in detail:


     
  * Type 1. Spots can only be detected within their corresponding lenslet boundaries. As the detection of each spot is completely independent, it directly restricts the local maximum of the wavefront gradient, ,  |θ|.

     
  * Type 2. Spots that stray outside the lenslet boundary can to some extent be detected based on avoiding their overlapping. This detection scheme is highly dependent on the association between adjacent spots. It restricts the local maximum of the wavefront curvature, , the difference between θ_i +1 and θ_i.

   
  * Type 3.  Spots can be detected even when their displacement paths overlap and cross over with each other in this work. The proposed scheme is motivated by recovering the high dimensional wave normal rather than its projection, thereby solving the ambiguity. The specific limitations are contingent upon the specific strategies employed to recover N.







 §.§ Strategies of Perceiving Wave Normal

To perceive the high-dimensional 𝐍, the fundamental step is to increase the dimensionality of the measurements. One commonly employed strategy is the implementation of axial scanning, which involves acquiring images from multiple planes along the z-axis. Several methods have been systematically investigated that leverage axial diversity measurements, including defocused SH WFS <cit.> and the known phase diversity <cit.>. Nevertheless, the introduced axial scanning scheme will inevitably sacrifice temporal resolution, significantly restricting its applicability in dynamic scenarios that require ultra-fast wavefront diagnosis.

Therefore, we develop a novel strategy to fully leverage the long-term information of the changing wavefront. It replaces the conventional active axial scanning scheme with a passive temporal-sequential measurement approach. To be specific, considering a temporally varying wavefront W(ζ, t),  as shown in Fig. <ref>(c). The changing rate of the wavefront will lead to the corresponding velocity of the spot for each lenslet, thus yielding a map of optical flow on the sensor plane

    ∇(∂/∂ tW)   =∂/∂ t(∇ W ) 
       =1/f∂/∂ t(θ f) =1/f∂/∂ t(x -x_ref) 
       =1/f∂x/∂ t - 0  = v(t)/f,

where the interchangeability of partial derivative order in the first line is established due to Schwarz’s theorem, and the equality in the second line holds according to <Ref>. 

Consequently, given two consecutive timestamps separated by an interval of Δ t (, t_1 and t_2), the total increment of the angular incident vector θ can be precisely reconstructed by directly integrating the instantaneous optical flow within Δ t


    1/f∫_t_1^t_2𝐯(t)  t =  
    ∫_t_1^t_2∂θ/∂ t t=. θ(t) |_t_1^t_2.


If assigning t_1 as the reference timestamp, the measuring obtained by <Ref> is simplified as an instantaneous estimation of its increment within Δ t, , Δθ(t).

For a sufficiently small Δ t or uniform motion, <Ref> will trace a straight 3D motion trajectory T within a spatio-temporal neighbourhood, represented by 

    T = Δ t ( , 1)^𝖳 =((∂θ/∂ tΔ t)f,Δ t)^𝖳≈(Δθf, Δ t)^𝖳,

where the approaching on the right side is due to the finite difference approximation.  Note that <Ref> yields an expression similar to that of the scaled wave normal N, , (Δθf, -f)^𝖳, with the third dimension now representing the temporal axis t rather than the axial axis z. 

However, as shown in the top row of Fig. <ref>(c), the classical WFS employs an active pixel sensor (APS) for data acquisition. The acquired temporal-varying Hartmanngrams form a stack of frames denoted as I(x, t), with the majority of values being non-informative zeros. Consequently, due to limited output channel capacity and computational memory, the significant redundancy in the frame modality inevitably restricts the maximum achievable sampling rate. Illustrated in the same row of Fig. <ref>(d), the limited sampling rate diminishes the effectiveness of the temporal diversity scheme, reducing it back to a conventional 2D displacement tracking between adjacent frames through iterative looping. Consequently, the accurate reconstruction of the temporal derivative becomes challenging within the frame-based modality, ultimately leading to tracking failure in the case of Type 2 and 3 DR.

Hence, we propose leveraging the emerging CNI paradigm to facilitate a cheap temporal diversity strategy, , NeuroSH. As depicted in the bottom row of Fig. <ref>(c) and (d), NeuroSH employs a DVS to transform the redundant frame stack as a highly efficient neuromorphic data modality, termed the neuromorphic Hartmanngram (NHG). Thanks to the asynchronous triggering mechanism in CNI, the NHG benefits from inherent data sparsity and low redundancy. By approaching the temporal-varying scaled wave normal N as the trajectories T within a spatio-temporal neighborhood, NeuroSH is capable of perceiving the high-dimensional wavefront normal directly.








 §.§ Forward Model
 
To derive a model for extracting the wave normal in utilizing the CNI paradigm, we begin by examining the fundamental event-triggering mechanism. As the DVS output, an event e_k = (x_k, p_k, t_k) is triggered independently at the pixelated location x_k = (x_k, y_k) if the logarithmic intensity change exceeds a predetermined threshold C at the timestamp t_k <cit.>. Herein, k represents the counting index of each event point, and the symbol p_k∈{-1,+1} indicates whether there is a positive (p_k=+1) or negative (p_k=-1) change. As such, the asynchronous event-triggering mechanism can be formulated as


    Δ L(x_k, t_k) ≜ L(x_k,t_k+1)- L(x_k, t_k)  = p_kC,

where L(x_k, t_k) ≜log(I(x_k, t_k)) is defined as the logarithmic scale of the pixel-wise intensity I(x_k, t_k), and Δ L(x_k, t_k) represents the recorded temporal increment from t_k to t_k+1.




Without losing generality, considering one moving SH spot on the sensor, it will trigger a cluster of events on the pixels it traverses within Δ t.  Such a motion-induced cluster is defined as ℰ≜{e_k}_k=1^N_e, containing N_e events.  Under static ambient illumination, the resulting temporal increment Δ L is completely caused by the moving of the spot pattern or the point spread function (PSF) of the sub-aperture. More details of theoretical derivation in Fourier Optics are included in Supplement 1. For a short Δ t, the resulting event cluster ℰ will enjoy a spatio-temporal correspondence due to the basic optical flow constraint <cit.>. From that, we successfully link the measured increment to the physical prior of the reference target

    Δ L(x) = -∇ L(x) ·Δ t,

where Δ L(x) represents temporal increment in a short Δ t, {·} indicates the dot product of two vectors, ∇ L(x)=(∂ L/∂ x,∂ L/∂ y)^𝖳 represents the spatial gradient of the target distribution near t_ref. 



In some naive approaches, the left-hand side of <Ref> is approximated by accumulating the individual increment contributed by each event to form an “event frame”, ,

    Δ L(x; Δ t) = ∫_t_ref^t_ref+Δ t e_k(x_k, p_k, t_k)  t,

where e_k(x_k, p_k, t_k)=Cp_kδ(x-x_k)δ(t-t_k) serves an alternative representative of an individual event that describes the amount of the carried increment. 






However, the event frame is just a 2D image and will suffer from the loss of dimensionality. The direct accumulation along the t-axis will cancel the spatio-temporal correspondence and cause motion blur.  Therefore, we propose to manipulate the event clusters in NHG from a neural manifold perspective, as is discussed in <cit.>.  Each event cluster is considered to describe a simple manifold embedding of the motion within the spatio-temporal neighborhood.  As shown in Fig. <ref>(a), if consider a long period, the triggered event cluster ℰ will contain replicas of the spatial gradient along the motion trajectory T.   That is, the triggered events by the i^th SH spots will “grow” along the motion trajectory T that starts from its own lenslet boundary Ω_i to form a distinct “tubular” morphology <cit.>. As depicted in Fig. <ref>(a), the tubular event clusters provide a distinction from each other in the whole NHG in the cases of Type 1, 2, and 3 DR. 


To extract the simple manifold embedding within NHG, a modified state-of-the-art contrast maximization framework is utilized to tackle the circumstance of cluster occlusion <cit.>. As is illustrated in Fig. <ref>(b), the conventional contrast maximization applies an affine transformation 𝒜 to the event clusters ℰ in the NHG to produce a geometrically warped event cluster ℰ≜{e_k}_k=1^N_e, with each point representing the warped event e_k'. The warping transformation can be defined as

    x_k () = 𝒜(x_k; t_k, ) = x_k - (t_k -t_ref) ,

where the flow  of the SH spot directly serves as the warping parameter.



Aggregating all the events in ℰ yields an accumulation image, namely an image of warped events (IWE), ,

    H(x; ) = ∫_t_ref^t_ref+Δ te_k (x_k, p_k, t_k) t,

where e_k (x_k, p_k, t_k)=Cp_kδ(x-x_k (v))δ(t-t_k) represents the discrete warped event by applying the affine transformation 𝒜, with δ (·) being the Kronecker delta function. 

By substituting <Ref> and  <Ref> into <Ref>,  the IWE can be interpreted as

    H(x; ) =-∇L̂(x) ·Δ t,

where ∇L̂(x) represents the estimated spatial gradient.





The contrast maximization assumes that if all the events triggered by the moving edges are properly aligned, the replica of these edges will be warped back to the reference state. This process is referred to as motion compensation, as all the moving edge replicas after warping now share the same pixel coordinates. Hence, as depicted in Fig. <ref>(c), if =0, then no motion is compensated, the producing IWE, , H(x;  =0) is equivalent with the definition of event frame, according to <Ref>. This is a histogram that merely estimates a blurred gradient. On the other hand, if the candidate flow is the optimal one ^*, the yielded IWE will produce the clearest estimation of the PSF spatial gradient, , ∇L̂^*(x) = ∇ L (x), given a specific t_ref. The expression in <Ref> can be therefore rewritten as

    H(x; =^*) = - ∇ L (x)·^*Δ t,

where ∇ L (x) represents the motion-resolved spatial gradient.


This problem is thereby transformed into an optimization problem by selecting an appropriate cost function to evaluate the contrast. According to the performance comparison list offered in <cit.>, the variance of IWE is the lowest time-consuming (only 16.9μ s ) and the most robust type of loss. The calculation of variance is also independent of the specific pattern, which indicates that it is a unifying metric. For each flow v starting from its own Ω, the selected cost function can be defined individually 

    g(v)=-σ^2(H (x;)) = - 1/|Ω|∫_Ω(H - μ(H (x;)))^2  x,

where |Ω| is the area within the lenslet boundary, μ(·) represents the statistical mean.



<Ref> indicates that there is no need to pre-separate the NHG into different clusters ℰ. To be specific, for a flow allocated to the i^th lenslet, denoted as v_i, the sub-cost function g(v_i)  can be computed by first applying a global warping transform to the entire NHG and then examining the cropped accumulated IWE in its own Ω_i. This derives an overall evaluation by enumerating all N_MLA lenslets and summing their sub-costs. By incorporating it with <Ref>, we link the flow to the angular incident vector, thus formulating the forward model, , 


    minimize_θ_i,v_i     ∑_i = 1^N_MLA g(θ_i)
    subject to      
        f∂θ_i /∂ t=v_i.



Consequently, since the model that links the physical optics to the neuromorphic modality is differentiable, the optimal solution can be easily extracted using well-established optimization algorithms in PyTorch <cit.>. To further accelerate the iterative optimization, other advanced convex optimization variants <cit.> and learning-based algorithms can be integrated <cit.>.





§ EXPERIMENT AND RESULTS



 §.§ Experimental Verification


To verify the feasibility of NeuroSH, we conducted an experimental validation of NeuroSH using a self-built adaptable SH setup, as depicted in Fig. <ref>. The adaptable SH setup offers two significant advantages. The prominent advantage is the ability to generate a series of dynamic known patterns on the same plane as the MLA, thereby providing a reliable phase ground truth (GT) for the temporal diversity scheme. Additionally, the setup allows for flexible overlaying of the Fresnel diffractive lens array pattern on the spatial light modulator (SLM), enabling a high degree of freedom in arranging the simulated MLA configuration <cit.>.





 




As depicted in Fig. <ref>(a), the adaptable SH setup incorporates a reflective phase-only SLM (LETO-3, 1920 × 1080 pixels, 6.4μ m) as the core component. In the setup, a He-Ne laser (Thorlabs, HNL100L, 632.8nm) is utilized for coherent illumination. The emitted laser beam passes through a 10 × beam expander (Thorlabs, BE10M-A) and a linear polarizer to ensure collimated and linearly polarized illumination. To account for the spatial limitations imposed by the beam splitter used in the reflection arm, a relay lens L1 with a focal length of 15cm is positioned to relay the Hartmanngram to an intermediate reflector (R2). Subsequently, another camera lens L2 is employed to further relay the light onto the sensor plane of an event camera (iniVation, DAVIS 346). The event camera has a resolution of 346 × 260 pixels and supports both APS and DVS output modes. 

In practice, the diffractive efficiency of the reflective SLM is typically less than 67%, resulting in the formation of a significant amount of zeroth-order light from the dead areas of the SLM. To address this problem, a steer modulation technique is employed <cit.>. This technique involves superposing a 2D linear phase pattern, also known as a blazed grating, onto the entire SLM screen. The added steer modulation aims to induce a lateral frequency shift on the Fourier plane, separating the diffractive reflective beam from the zeroth-order reflective beam (non-diffractive). To select the diffracting part of the reflective beam and exclude the higher-order terms caused by the periodic SLM pixels, a blocker is introduced on the Fourier plane. Additionally, due to the special logarithmic response curve characteristic of the event-triggering mechanism in CNI, the output from the DVS channel exhibits high sensitivity to ambient noise as well as the adjustment of the liquid crystals in the utilized SLM.  Therefore, an additional step is introduced to further relay the Hartamanngram to an intermediate diffusive reflector R2. This allows for flexible adjustment of background illumination, thus leveraging the optimal segment of the logarithmic curve.



The generation of the Fresnel diffractive lens array pattern onto the SLM follows the principle described in <cit.>. In our setup, we generated a 6 × 6 Fresnel lens phase pattern with a focal length of 10cm and a pitch size of 180 SLM pixels. This design allows us to make full use of the limited spatial resolution of the SLM. The system magnification is approximately 0.5, which resizes the SLM aperture with a size of 6.9mm onto our sensor with a size of less than 4.8mm. As a result, each MLA pitch corresponds to approximately 31 sensor pixels. While using a singlet to relay the Hartmanngram onto R2 introduces additional phase aberration, it is important to note that SH WFS measures the wavefront increment rather than its absolute value. Therefore, even if the position of r_ref slightly deviates from the center of each sub-aperture due to the introduced aberration, the linear relationship remains valid.



In terms of overlaying the dynamic phase GT, the 12^th order Zernike polynomial, , the spherical aberration Z_12, is chosen to simultaneously evaluate all three types of DR. As shown in Fig. <ref>(b), the normalized phase distribution of Z_12 is characterized in polar coordinates within a Zernike radius of R. We then generate a series of patterns corresponding to Z_12, with their amplitude increasing linearly in the range of [0, 5] π rad. To achieve optimal continuity, we ensure a maximum SLM refreshing rate of 60Hz and adjust the step to 0.2 π rad. This results in a loop period Δ t, which is nearly 400ms. In the top right of Fig. <ref>(b), we present a series of images corresponding to the Z_12 GT and its resultant simulated Hartmanngram within the same range but with a larger step of π rad for clearer demonstration. With an increase in mode amplitude, the edge region of Z_12 represents a sharper change, particularly around the "valley" region. This leads to contrary increasing directions of the local spatial gradient, gradually being dominated by Type 2 and 3 DR. The system calibration of NeuroSH is illustrated in the top right of Fig. <ref>(c). The specific procedures involve identifying a calibration matrix, denoted as D, by fusing the long exposure frame from the APS channel under static conditions and applying it to the dynamically triggered event streams. More details are included in Supplement 1.




 §.§ Algorithmic Implementation




Fig.<ref>(a) illustrates the performance of the event camera from different channels in practical dynamic scenarios. The APS channel preserves an exposure time of 300ms to adapt to the ambient illumination, however, it suffers from significant motion blur. Therefore, after completing the calibration procedure and registering the reference position x_ref, only the DVS output is utilized in the algorithmic implementation step. Additionally, to ensure a fair performance comparison with<cit.>, only positive events, , p_k = +1, are utilized.




The coarse alignment is realized by adopting conventional centroid methods through a temporal segmentation scheme, similar to that in <cit.>. Specifically, for any timestamp t_i of interest, the events within the spatio-temporal neighborhood are divided into several temporal segments using a smaller time window Δ t'. The events within adjacent segments are then accumulated according to <Ref> to produce the segmented event frames S_i and S_i+1. The displacement between the centroids of S_i and S_i+1 within each Ω serves as the starting point. As depicted in Fig. <ref>(b), for the reference timestamp, since the reference position has already been determined in the calibration step, only the centroid of S_1 needs to be calculated. A robust starting point, in addition to helping the optimizer avoid being trapped in a local minimum due to flow discontinuity, also aids in determining whether a reliable motion is occurring within each Ω. For instance, in this case, we have identified 7 sub-apertures to be excluded in the following optimization, as they contain events that fall below a certain threshold in the first segment. Consequently, these points are filtered out from further analysis.
 

The fine alignment procedure involves applying the proposed NeuroSH forward model to perform iterative optimization.  As shown in the bottom of Fig. <ref>(c), the combined figure contains the optimal candidate flow as well as the corresponding IWE. By examining the contrast of IWE, the effectiveness of NeuroSH has been demonstrated.  To showcase that NeuroSH has successfully tackled all three types of DR, we preserve an independent view of each extracted flow starting from each lenslet boundary. In the case of Type 2 and 3 DR, the fixed region of each Ω, colored in dark blue, is scaled to adapt to the displacement outside the lenslet boundary. The independent views of Type 1, 2, 3 DR are marked with green, orange, and red, respectively. Fig. <ref>(d) illustrates the optimization outcome for all the reliable starting points. Herein, the extracted optimal trajectory T^* is straight within the long period Δ t due to the linear increasing of the Z_12 mode. It is worth noticing that this verification procedure is still fundamental as well as effective in non-linear scenarios through basic temporal segmentation.




Additionally, NeuroSH exhibits further superiority in reconstructing the motion-resolved Hartmanngram with low-light enhancement compared with the conventional methods based on the event frame. As shown in Fig. <ref>, given a selected timestamp of interest, say t_ref, the production of an event frame requires a manual selection of the segmentation window Δ t'. This is prone to be problematic in general cases where both fast and slow flows occur simultaneously within one single segment. On the one hand, to catch the fast flow,  Δ t should be chosen as short as possible. However, this leads to a low signal-to-noise (SNR) in such a thin segment. For example, setting Δ t' as 10ms makes the spots in the center of the event frame nearly invisible, since each event naturally conveys less information and is subject to noise. On the other hand, choosing a larger window, , setting Δ t to 200ms, will again result in the motion blur effect to the spots near the edge. The NeuroSH tackles this challenge by combining the compensated IWE of each Ω together to reconstruct a HDR image that is free from motion blur. The comparative outcome between the yielded IWE and the APS reference frame showcases the robustness of NeuroSH, owing to its weak reliance on the selection of Δ t.





 §.§ Results Comparison and Error Analysis

To analyze the performance of NeuroSH, we reproduce the conventional centroid tracking method to form the comparative group. Specifically, by accumulating events within a segmentation time window Δ t' of 10ms to form the event frame, the centroid in each Ω can be calculated. As shown in Fig. <ref>, the 3D trajectories of the SH spots can be extracted through this temporal tracking scheme. To highlight the centroid tracking outcomes in terms of different types of DR, we colored the trajectory of Type 1, 2, and 3 DR in green, orange, and red, respectively. By comparing the outcomes of the NeuroSH with that of GT, it is evident that the proposed NeuroSH model unequivocally provides rapid and dependable tracking for all types of DR. Conversely, the results obtained from classical centroid tracking cannot deal with  2, and 3 DR and exhibit a notable fluctuation attributable to the low SNR within each thin segment and cross-talk from adjacent sub-apertures.










In error analysis, we first examine the end-point measurement outcome, that is, to extract the displacement map directly from the last segment in Fig. <ref>(a). The centroid outcome calculated with a Δ t' of 10ms is extremely noisy. To guarantee a fair comparison,  the last segment window Δ t' is tuned to 80ms to exact the tracking of the spots as much as possible. However, as discussed before, the challenge of precisely determining the fast and slow flow simultaneously within the produced event frame is impossible. Furthermore, it is clear that nearly no displacement can be detected by the centroid method for spots allocated to Type 2 and 3 DR. Fig. <ref>(b) illustrates the mean squared errors on the average of the sub-apertures allocated to the 3 types of DR for all segments. It is obvious that once the points move outside their designated aperture Ω, the centroid tracking fails and exhibits significant tracking errors.  Consequently, the advantages of NeuroSH in tracking spots demonstrate its capability of achieving accurate and robust wavefront sensing that surpasses the previous DR limitations.







 §.§ Application: Instantaneous Flame Front Visual Odometry



Here we exhibit an application of the proposed NeuroSH for instantaneous flame front visual odometry. Flame front visualization provides valuable insights into the physics and chemistry of flames, such as flame structure, flame stability, flame propagation, and reaction kinetics. A flame front is defined as the boundary or interface separating the unburned fuel from the combustion products within a propagating flame. It represents an active area where the combustion reaction takes place, exhibiting a rapid emission of energy, heat, and light. That is, its specific sub-region encompasses a significant temperature gradient, resulting in a large-gradient wavefront. As a result, the observation and analysis of the flame front's intricate structure during the dynamic propagation process necessitate the wavefront sensing technology with a large DR. 






To achieve flame front characterization, the 6×6 sampling grid used in the GT verification experiment is inadequate for finer sampling. Due to the limited resolution of the phase-only SLM, generating a higher-dimensional MLA is a great challenge. In addition, The flame front visualization should necessitate a ten times larger field of view, which cannot achieved for the current SLM screen without scaling. Referring to <cit.>,  the background-oriented schlieren (BOS) technology shares a similar gradient-based forward model with the classical SH wavefront sensing, but it offers a larger FOV. We implemented the BOS experiment by generating a chessboard mask on the screen to simulate a densely hexagonally arranged MLA with a 36 × 48 sampling rate.  More implementation details about the BOS setup are contained in Supplement 1.


As illustrated in Fig. <ref>(a), the APS channel's fixed exposure time of approximately 74ms restricts its ability to capture anything beyond blurred aberrated screens in dynamic scenarios. Hence, the information provided by the APS channel regarding the flame front is limited, and even simple subtraction of the screens yields insignificant differences. On the other hand, Fig. <ref>(b) displays the comprehensive NHG obtained from the DVS channel, which records the detailed propagation behavior of the flame front. The accumulated event frame within a window of interest of 14 ms provides a rough outline of the flame front. By further reducing the window size to 2.8ms, the scatter plot reveals the rough structure of the flame front. A comparison between the outputs of the APS and DVS channels can be found in Video 1.

Subsequently, we compare the results in terms of the absolute velocity map and the directional pseudo-color velocity map using centroid tracking and NeuroSH methods. As depicted in Fig. <ref>(c), the centroid method extracts only a few flame front structures, which are also blurred, indicating that the tracking method is only effective in the region of Type 1 DR.  On the other hand, NeuroSH successfully captures the morphology of the large-gradient hemispherical flame front, where Type 2 and Type 3 DR dominate. Consequently, the results obtained from NeuroSH align with the flame front formation theory under laminar burning conditions. Extensive outcomes showcasing the dynamic flame front propagation behavior are included in the last section of Supplement 1 and Video 2.
 







§ CONCLUSION AND DISCUSSION

In conclusion, we develop NeuroSH, which is a physics-informed wavefront sensing framework based on the bio-inspired CNI paradigm. NeuroSH enables ultra-fast and HDR wave normal sensing, which successfully handles all three types of DR. To achieve this, NeuroSH first utilizes the proposed low-cost temporal diversity strategy based on the CNI paradigm to efficiently gain data dimensionality. Then it leverages a modified global warping scheme to extract the optimal flow that produces the clearest IWE of each sub-aperture. By approaching the 3D trajectory within the spatio-temporal neighborhood as the high-dimensional wave normal, NeuroSH can resolve the inherent ambiguity in displacement distinction. The conduction of an adaptable SH setup on a phase-only SLM provides reliable experimental verification of NeuroSH. Furthermore, NeuroSH is utilized for instantaneous flame front odometry in a BOS setup, showcasing its remarkable advantage of monitoring rapidly changing phase flows with a large DR.



In particular, we emphasize that the proposed differentiable forward model is an important innovation, which establishes a bridge between the neuromorphic modality and the SH optical configuration.  Due to the property of the utilized contrast maximization scheme, the implementation of  NeuroSH is applicable for all angle-based wavefront sensing models without specifying the mask pattern. Furthermore, the obtained motion-corrected image can facilitate the decoupling between the moving targets and the static background under photon-starved scenarios. As such, NeuroSH is widely applicable in enhancing  3D tracking and quantitative phase imaging applications, including holographic particle tracking <cit.>, speckle analysis <cit.>, and non-line-of-sight imaging <cit.>.









	
Funding
The work is supported by the Research Grants Council of Hong Kong (GRF 17201620), ACCESS — AI Chip Center for Emerging Smart Systems, sponsored by InnoHK funding, Hong Kong SAR.

Acknowledgments
We thank Haosen Liu, Zhou Ge, Dr. Hai Gong,  Prof. Evan Peng, and Prof. Wolfgang Heidrich for their technical support and useful discussions. J. H appreciates the partial financial support from Shanghai Jiao Tong University and the Hong Kong Scholars Program (XJ2022032).

Disclosures
The authors declare no conflicts of interest.

Data Availability Statement
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.  Codes will be released once published.

Supplemental document See Supplements 1, Video 1, and 2 for supporting content.
 


































