A Regularized Inverse Holographic Volume Reconstruction (RIHVR) GUI Program for Digital Inline Holography (DIH)

1 Introduction

Digital inline holography (DIH) as an experimental technique finds its origins in the works of the Hungarian physicist Dennis Gabor [1]. At its core, it records the information of a three-dimensional (3D) scene onto a two-dimensional image plane. This is achieved by exploiting the wave nature of light to record the two-dimensional interference pattern generated by the interaction of light with the propagation medium and the object field [1]. Due to its cost-effectiveness and simplicity, DIH has rightfully found its use in multiple disciplines such as fluid mechanics for particle tracking velocimetry (PTV) [2, 3], microbiology for visualization of microbes [4], particulate size measurements [5, 6], optical phase retrieval [7, 8], and so on. The reader is referred to Refs. [9, 10, 11] for a comprehensive treatment of the subject.

Due to its wide applicability, an accurate quantitative reconstruction of the 3D object field from its two-dimensional hologram is a fundamental problem in experimental physics. However, this reconstruction is an ill-posed inverse problem (IPIP), as the camera sensor only records the intensity of the interference pattern, leading to the loss of its phase information. During reconstruction, this inherent ill-posedness therefore manifests as the so-called twin-image problem, further details of which can be found in Stoykova et al. [12]. Typically, in a DIH setup such as the one shown in Figure 1, a light beam (typically assumed to be a plane wave), emitted from a coherent light source (usually a laser), interacts with the object field to produce the object wave O. The object wave O and the reference wave R interfere with each other. The recording media, such as a film or a camera, capture the intensity field I, which can be mathematically written as

Here, R^⋆ and O^⋆ represent the complex conjugate of the reference and the object wave, respectively. The first term in the sum |R|² is the amplitude of the reference wave and is typically subtracted from the recorded hologram. The second term |O|² is the amplitude of the recorded object wave and is typically small. The third term OR^⋆ contains the most crucial information as it is this cross term, which, when illuminated by the reference wave, that is $R(O{R}^{\star })={|R|}^{2}O$, results in a real image of the object field. At the same time, the fourth term O^⋆R is the main culprit behind the ill-posedness inherent to holographic inversion, as it leads to the conjugate image, also known as the twin-image of the object.

Figure 1

The presence of this twin image makes the object field reconstruction sensitive to noise, while also leading to artifacts in the reconstruction. To remedy this problem, the inversion of the hologram to compute the object field is typically performed iteratively while subjected to various physics-based constraints to suppress the artifacts and noise [13, 14]. This is done by reformulating the inversion as an optimization problem in the variational framework as described by Soulez [13], and Mallery and Hong [14]. The variational problem is given as follows.

Here, ${\lambda }_{l_1}$ and λ_TV are the case-dependent Lagrange multipliers for the optimization problem and are set by the user. Increasing ${\lambda }_{l_1}$ allows the user to promote sparsity in the solution, while increasing λ_TV reduces the total variation (TV), resulting in a smoother reconstruction. $\widehat{O}$ is the best estimate of the reconstructed object field. I is the recorded hologram and $\tilde{I}(O)$ is the hologram generated using the object field O. To generate a hologram $\tilde{I}(O)$ from an object field O, a forward model for hologram generation is required. A common approach to propagate the reference plane wave through the object field is via the Green’s response function for the homogeneous Helmholtz equation [15]. To ensure solution convergence and computational efficiency, this propagation operation is linearized, resulting in Eq 3 for computing $\tilde{I}(O)$.

In Eq. 3, N_z denotes the number of axial planes to be reconstructed. z denotes the axial depth. k₀ is the wavenumber, and h_z is the Rayleigh-Sommerfeld kernel. Note that the parity change of z to –z for h_z in Eq. 3 is due to the traditional sign convention used when writing h_z. As Eq. 3 contains a convolution operation, it is computationally efficient to perform it in the Fourier domain. The discrete Fourier representation of h_z is given in Eq. 4, and is known as the optical transfer function (OTF) for the Rayleigh-Sommerfeld kernel.

Here, m and n denote the discrete coordinates, s is the pixel size, and M and N is the size of the hologram in pixels. For a given holographic recording setup, the OTF can therefore be pre-computed for later use. Once the OTF is known, $\tilde{I}\left(O\right)$ is readily computed in the Fourier domain via the operation in Eq. 3. It should be noted that for backpropagation (reconstruction), and forward propagation (generation), only the sign of z changes in h_z. Hence, the same OTF can be used for both forward and back propagation.

Intuitively, the problem posed in Eq. 2 proposes that instead of reconstructing the object field directly from the hologram, the problem can be reformulated as a minimization problem that finds the object field $\widehat{O}$ that would produce the observed hologram I while subjected to the specified physical constraints. In this work, the physical constraints are sparsity and smoothness, also known as the fused lasso. These constraints are chosen as in a typical holographic reconstruction scenario, the object field is sparse. However, mere sparsity can lead to over-segmentation of the object field, hence it is coupled with a TV penalty that promotes smoothness in the object field [14]. The introduction of these additional physical constraints enforces a unique solution by regularizing the inversion. Further details on the solution to the problem can be found in Mallery and Hong [14].

Computationally, solving the minimization problem posed by Eq. 2 has its unique challenges. The first challenge is mainly due to the non-differentiability of the l₁ norm (the absolute value operation). Hence, most of the efficient optimization algorithms, such as the limited-memory Broyden–Fletcher–Goldfarb–Shanno (l-BFGS) algorithm or other quasi-Newton approaches, are rendered inapplicable as they require the gradient of the functional [16, 17]. One is therefore forced to rely on non-gradient-based methods that can be prohibitively slow for large problems [18]. The second major challenge, which is shared by all iterative solvers, is that solving Eq. 2 requires a single hologram to be reconstructed hundreds or thousands of times before converging to the optimal solution. This poses a major computational speed bottleneck. The third major challenge is the significant storage and working memory requirements for volumetric processing. For a typical one megapixel image (1024×1024 pixels), if the reconstruction is performed at 1024 planes, it would require approximately 16 GB to store it in memory (complex double-precision floating-point values, 16 bytes), while also requiring several additional variables of this size during the optimization process. One would therefore be forced to either limit the object field size to take advantage of the speed increase by GPU-acceleration or to rely on the much larger RAM available on most desktop computers while settling for a much slower CPU processing [14]. While it is true there has been a surge in the GPU memory recently, especially for data center use, most of the DIH applications are at a lab-scale, which are still limited in their DRAM. For instance, a consumer-grade RTX 5080 has a VRAM of 16 GB.

The work of Mallery and Hong provides solutions for these challenges [14]. Firstly, they employ a GPU-accelerated version of the fast iterative shrinkage-thresholding algorithm (FISTA), which can efficiently handle a non-differentiable cost function via computing the so-called proximal gradient of the non-differentiability to guide the optimization process [14, 18]. This solves the challenges for computational efficiency and convergence. Secondly, by exploiting the inherent sparsity of the problem, the computations are performed in the COO format (coordinate sparse matrix), which only requires the indices and values of the non-zero elements to be stored. Therefore, they store only the non-sparse data to enable GPU processing on consumer-grade GPUs that have a typical VRAM of 12 GB. Thus, memory usage is reduced as long as the data sparsity (number of zero elements divided by total) exceeds 67% [14].

In this work, we present the numerical details of the CUDA-accelerated implementation of the regularized inverse holographic volume reconstruction (RIHVR) algorithm developed by Mallery and Hong [14] along with multiple background subtraction routines and a C++ implementation of the Crocker-Grier tracking algorithm [19]. All the routines are implemented in C++ and CUDA for maximum performance and cross-platform portability. For ease of use and distribution, all the routines are packaged into a coherent package with a graphical-user interface (GUI), making it a one-stop solution for DIH-based particle sizing and tracking.

2 Implementation and architecture

The method of Mallery and Hong, as described in Sec. 1, is implemented in CUDA and C++. A GUI is developed in C++ using QT6 [20]. The core workflow of the software is illustrated in Figure 2, and further details are provided in the subsequent discussion in Sec. 2.3. The green boxes represent the executable binaries. The GUI is used to handle the loading of the data, specification of the parameters, and management of the files. Once all the required settings are specified, the GUI calls the processing routines for each section as specified by the user. For the background computation, there are two processing executable binaries: ccbg and make-background. The first executable, ccbg, computes the background image using a temporal correlation-based strategy, while make-background is adapted to compute a mean, median, or a moving average background. Once the background image is computed, the sparse-inverse-recon executable is called to solve the optimization problem posed in Eq. 2, to compute the object field for each hologram. If selected, the object field is passed to the segmentation executable to extract the location of the centroids of the objects. These extracted centroids are then used by the particletracking routine to compute particle tracks via the Crocker-Grier algorithm. Hence, the whole package serves as a one-stop solution for determining particle tracks from a sequence of holograms. An in-depth, detailed discussion on all the parameters and the codebase is provided in subsequent sections.

Figure 2

2.1 End-user options

In this section, we describe the effects and behavior of multiple parameters and options that can be specified in the program. The section is divided into different sub-sections such that each sub-section corresponds to a different tab in the GUI of the program. The GUI consists of six different tabs. Namely: Background, Hologram Processing, Tracking, Visualization, Output Console, and Batch Automation. Each tab can be used independently or in a sequence with others using the batch automation process implemented in the Batch Automation tab.

2.1.1 Background

To increase the signal-to-noise (SNR) of the hologram and to get rid of the first term |R|² from Eq. 1, a reference background is typically subtracted from the recorded raw hologram [2, 3]. This background can either be recorded during the hologram capture process or can be generated from the recorded holograms, depending on the data. There are multiple strategies that can be used to generate a background. In the developed software, we have implemented four different methods: mean, median, moving average, and correlation-based. In the following sections, the differences and similarities among them are discussed, along with the effect of the end-user parameters on the final results.

2.1.1.1 Directories

The program uses a C-style naming convention for naming the files and the directories. To start, the user should browse to a sample image from an image sequence in the Image Path field. Once the sample image is selected, the name of the image is automatically turned into its C-style format and is displayed in the Image File Name field. Start Image and End Image fields are also populated by using the information from the selected folder. By default, the Start Image and End Image also populate the Background Images Range. An identical directory structure and file browse format is used for all the different tabs in the GUI.

2.1.1.2 Mean

The simplest choice for generating a background from a sequence of holograms is to compute a mean of the entire hologram sequence and use the resulting image as the reference background. To compute the mean background, the user can set the background type to ‘mean’ in the background tab of the GUI. The background images range, specified under the Background Computation Options group, decides the number of images that would be considered in the background computation. Once the computation is done, the computed background is stored as background.png in the specified folder. The mean background is an appropriate choice when a fixed pattern noise is present in the images, which can be removed by subtraction and normalization. This is the case for Gaussian noise, as a mean of the images containing Gaussian noise gives the lowest-variance estimate of the true background [21]. However, this requires that the holographic signal of the objects to be reconstructed moves enough over the sequence such that only the true background contributes significantly to the mean. This is easily achieved by holograms captured for the purpose of PTV, as the particle motion is significant enough that only the true background contributes to the mean.

Figure 3 shows a contrived example that illustrates the effect of mean background subtraction. Figure 3(a) shows an image from a sequence of synthetic holograms of random spherical particles containing synthetic fixed pattern additive Gaussian noise, (b) shows the mean of the sequence (colored for better visualization), and (c) shows (a) after mean background subtraction. It is clear from (b) that a high-fidelity estimate of the fixed pattern noise in a sequence of holograms can be generated by computing the mean of the image sequence. This is so because a fixed pattern noise stays the same, while the holographic signal moves over a sequence. By performing mean background subtraction, the SNR of the hologram has clearly increased, as evident by comparing the quality and the structure of the hologram fringes in (a) and (c). However, converging to the true background requires enough uncorrelated hologram samples in the recorded sequence such that the signal due to the particle motion converges to zero in the mean. In the illustration shown in Figure 3, 200 images were found to be sufficient to obtain a reasonable estimate of the fixed background pattern. It should be noted that the mean background requires the optics and illumination to be temporally constant, as any transient changes in the illumination result in a background that might be ill-suited for the entire sequence.

Figure 3

2.1.1.3 Moving Average

The moving average background is a modified version of the mean background. Instead of computing the mean for all the images, a moving average of the window size specified by the user in Window Size field is used. Once the computation is done, the computed backgrounds are stored in the format of background_%04d.png in the specified folder. The moving average background is particularly useful in the presence of temporal background shifts, as the background is adapted for each image based on its surrounding images. However, sudden changes in the background are not remedied by a moving average background. It also typically requires a longer sequence of holograms to be effective.

2.1.1.4 Moving Median

To compute the median background, the user can set the background type to median in the background tab of the GUI. Currently, a temporal filter of length 5 is used around the image, leading to N–3 backgrounds for a sequence of N images. Once the computation is done, the computed backgrounds are stored in the format of background_%04d.png in the specified folder. The median background subtraction is an appropriate choice to increase the SNR of the moving objects in a hologram sequence, as it acts as a temporal edge preserver while also removing noise. The moving median is preferable to the related moving mean in situations where outliers, such as large dark objects, make the average unstable.

2.1.1.5 Correlation

Correlation-based background is an advanced version of the moving average background. Instead of using a fixed window size like a moving average or using images in a consecutive sequence, correlation-based background subtraction first computes a temporal correlation between all the images in the sequence. Once the temporal correlation matrix is computed, for every image, only the most temporally correlated images are used to generate a background. Further details on the working of the method are provided in Sec. 2.3. As computing the correlation matrix could be computationally expensive, the user can choose to down-sample the images for correlation computation by specifying a Resize Factor, which defaults to 0.25. By default, the correlation-based background uses an entropy-based estimate to determine the optimal number of temporal samples to be considered for per-image background generation. However, it is possible to provide an upper and a lower bound by specifying the Min frames and Max frames options in the Correlation-based Background Options group. The default minimum and maximum frames are set to 0 and 100 respectively. Once the computation is done, the computed backgrounds are stored in the format of background_%04d.png in the specified folder. If desired, the user can also save the background-subtracted images by selecting the Save Enhanced checkbox. Correlation-based background subtraction strategy generates the best possible per-image background by exploiting the temporal correlation between different images. This method is well-suited for hologram sequences with temporal fluctuations. However, a large sequence of samples is desired to compute a strong enough correlation that can be exploited for optimal background generation.

2.1.2 Hologram Processing

The object field from the holograms is computed using the parameters specified on the Hologram Processing page. Directories are automatically populated using the directories specified in the Background tab. The default background type is selected based on the options selected in the Background tab. The user also has an option to skip a certain number of images at a specified uniform interval by specifying the Increment, which is set to unity by default.

2.1.2.1 Reconstruction Parameters

The parameters needed to optimize the functional in Eq. 2 are specified in the Reconstruction Parameters group. Firstly, z0 is specified in µm, which sets the starting axial plane location for the reconstruction volume. dz specifies the axial step size in µm, that is, the distance between two consecutive axial planes. Nz is the total number of planes in the axial direction and is a positive integer. dx is the magnification, hence the size of a pixel in µm in the captured image. By specifying the four parameters: z0, dz, NZ, and dx, the full reconstruction volume is therefore declared. As the OTF given in Eq. 4 depends on the wavelength of light, the parameter lambda can be used to set the wavelength in µm. It should be noted that lambda specified should be divided by the index of refraction of the propagation medium. As the computations are performed in the Fourier domain and the holographic data is almost never periodic, computations can lead to ringing artifacts at the edges of the domain due to the Gibbs phenomenon [22]. To mitigate these artifacts and edge effects, the user can specify the size of zero-padding in the 0-pad field, which increases by domain size by padding zeros to all sides of the image domain. It is recommended to pad the image to get close to 2ⁿ type domain size, as the Fourier transform algorithms are much more efficient when the domain size can be written as 2ⁿ. Please note that all input images must use the unsigned 8-bit integer data type.

In the same group box, the two regularization parameters are also set by the user. Both of these parameters directly dictate the quality of the results obtained. The values specified for the Sparsity and TV field set ${\lambda }_{l_1}$ and λ_TV in Eq. 2 respectively. Figure 4 illustrates the effect of the choice of these parameters on the volume reconstruction by showing the maximum phase projections computed along the axial direction at different values of ${\lambda }_{l_1}$ and λ_TV of the background-subtracted hologram shown in (a). Maximum phase projections are generated by computing the voxel-wise maximum of the 3D reconstructed phase field along a chosen axis. This approach provides a clear, planar visualization of inherently volumetric data, making structural features easier to interpret at a glance. In these projections, color encodes the phase value: brighter regions indicate areas that induce greater phase delay, corresponding to optically denser media. Unless stated otherwise, for all the results shown in this section, the following parameters were used: z0 = 0 µm; dz = 8000 µm; dx = 10.4 µm; Nz = 10; lambda = 0.633 µm; 0-pad = 64; Sparsity = 0.001; TV = 0.02; Number of FISTA Iterations = 100; Number of TV Iterations = 5. The hologram in (a) is of two dandelion stems at different axial distances and is adapted from the work of Brady et al. [23]. When ${\lambda }_{l_1}$ and λ_TV are set to their optimal values of 0.001 and 0.1, respectively, the maximum phase projection in (b) shows that a clear stem and branch structure is recovered with very low noise. However, if the value of λ_TV is chosen to be less than optimal (0 for the case shown), the results show the appearance of noisy structures in the background, and the sharp boundary of the branches is clearly lost, resulting in under-smooth results. Furthermore, if the value of λ_TV is set to be greater than the optimal value (0.5 for the results shown), as shown in (d), the sharp boundaries of the branches are washed out. This leads to the appearance of a smoothened halo around the object. A similar behavior is observed when the value of ${\lambda }_{l^1}$ is chosen to be less than optimal, as shown in (e). When the solution is obtained with an underestimated sparsity value (0 for the results shown), film-like noisy structures are observed between the stems. This leads to an under-segmentation of the branches, as clear from the results shown in (e). On the contrary, if the sparsity value is overestimated (0.1 for the results shown), the solution field is oversegmented, leading to the breakdown of real coherent objects in the volume field as shown in (f). This is because the minimization problem over-prioritizes the sparsity constraint, leading to a large number of zero values in the volume reconstruction. The regularization parameters should therefore chosen with utmost care for optimal volume reconstruction.

Figure 4

From experience, the authors note that for most of the experimental cases, λ_TV = 0.1 and ${\lambda }_{l^1}=0.1$ is typically a good starting point. In the context of particle tracking applications, it is typical to decrease ${\lambda }_{l^1}$ with an increasing particle density to allow for a larger number of non-zero voxel intensities in the reconstruction volume. This may produce overly irregular solutions, which can be mitigated by a slight increase in λ_TV. In most of the cases, it is expected that both λ_TV and ${\lambda }_{l^1}$ remain within the 0.01–0.5 interval.

2.1.2.2 Optimization Parameters

Two important parameters, namely Number of FISTA iterations and Number of TV iterations, control the convergence of the solution. These can be specified by the user in the Optimization Parameters group. Figure 5 shows the dandelion hologram from Figure 4(a) reconstructed at different numbers of iterations. (a) shows the maximum phase projection results when Number of FISTA iterations is set to be less than the optimal; 10 iterations instead of 100 in the results shown. In such a scenario, the solution has not converged to its minimum value, and hence, the presence of noisy structures is observed in the reconstruction. If the number of FISTA iterations is raised above the optimal of 100 to 200, as shown in (b), the results appear to be identical to the ones shown in Figure 4(b) at 100 iterations. Therefore, adding more iterations than necessary only increases the computational time without improving the results once the minimum is achieved. Hence, the user is recommended to track the evolution of the functional in the Output Console to set an optimal value for Number of FISTA iterations. In the same spirit, when Number of TV iterations is set to be less than optimal (set to unity in the results shown), as done for (c), the results appear to be lacking detailed, sharp boundaries around the structures. However, once the value is raised above the optimal from 5 to 20, only minor gains in SNR of the reconstruction and structural integrity are observed, as evident from comparing the results in Figure 5(d) with Figure 4(b), but at a higher computational cost. The user is therefore recommended to choose these values such that acceptable results are obtained without an unnecessary increase in the computational cost.

Figure 5

2.1.2.3 Cropping Parameters

The user can also crop the input hologram images to get rid of the unnecessary regions. This can be achieved by using two different ways of image cropping: Center and Rectangle, in the Cropping/Resizing group. The parameters for each of the cropping options can be set in their respective group boxes, namely, Center Cropping Options and Rectangle Cropping Options. Both of these options require the user to select the reference pixel location and the size of the crop. Figure 6 illustrates the coordinate system for both the cropping settings. For center crop in Figure 6(a), the user is required to set the value of the center pixel (marked as red cross) of the final region (boxed as purple) in x0 and y0 fields, along with the side length of the box in L. For the rectangle crop in Figure 6(b), the user is required to set the top-left pixel (marked as red cross) of the final region (boxed as purple) in x0 and y0 fields, along with the size of the rectangle in Lx and Ly.

Figure 6

2.1.2.4 Saving Options

The selections made by the user in the Save Options group box decide the output of the hologram processing. All the data is saved in the specified Output Directory. If Phase Projections is checked, a folder named Projections is created inside the specified Output Directory to store the maximum phase projections along all three axes. These files are named in the C-Style format as cmb_xy_%04d.tif, cmb_xz_%04d.tif, cmb_yz_%04d.tif, for XY, XZ, and YZ projections, respectively. If the user wishes to save the voxel data in the sparse matrix format (which contains the indices and values of the non-zero elements), the Voxel Data checkbox can be enabled. Voxel data is stored inside the VoxelData directory, which is created inside the Output Directory. The naming convention for the voxel data is inverse_voxels_%04d.csv, with columns of the .csv file being $x,y,z,O$ values. This data can later be used for reconstructing the full 3D volumes. The user can also choose to save the background-subtracted images by checking the Enhanced Images checkbox, which ensures that the background-subtracted images are saved inside the Enhanced directory in the specified Output Directory. The enhanced images are named as enhanced_%04d.tif. It is also possible to save the volume reconstruction at every single plane as a TIFF image by enabling the Phase Reconstruction option. However, this could be very storage-heavy. By default, the images are stored inside the Output Directory and are overwritten for every image. If the user has selected the Create Subfolders option, the images are not overwritten and are stored inside a sub-folder corresponding to the image number of the image. To visualize the effect of different numbers of FISTA iterations on the results, the user can enable the Each Iteration option that outputs the results at each iteration inside the Output Directory. However, these are also overwritten by every image processed to save storage.

2.1.2.5 Segmentation Options

To compute particle tracks, it is necessary to extract the centroids of the particles from the reconstructed volume field. This can be done by enabling the Extract Centroids checkbox. Once enabled, the Segmentation Options group box is activated. This group box contains three core parameters for segmenting the volume field. The first parameter, Minimum Voxel Size, is used to neglect the oversegmented small regions. It is recommended to be 1 by default. The second parameter, Minimum Voxel Intensity, can be used to threshold the volume such that only the intensities above the specified threshold count as a detected particle. Note that the voxel intensities are between 0–255. The third parameter, Segment Close Size, sets the kernel size for a morphological closing operation performed prior to connected component labeling to identify each segmented object. The extracted centroids are saved as .csv files in a folder named Centroids inside the Output Directory. There are different files that are generated for centroid data. The first file, centroids_%04d.csv, contains the centroids data in the $x,y,z,\text{Size}$ format (Size is the number of voxels in each object). The second file, centroids_weighted_%04d.csv also contains the centroids data in the $x,y,z,\text{Size}$ format, but the stored coordinates are intensity weighted. The third file, centroids_verbose_%04d.csv, contains the coordinates with additional information such as the minimum and maximum intensity, along with weights for each axis.

2.1.3 Tracking

Once the centroid data has been extracted from the reconstructed volumes, the user may choose to perform particle tracking via the Crocker-Grier algorithm. The Centroid Data is selected according to the settings on the Hologram Processing page. By default, the tracks are saved to Tracks directory created inside the Output Directory specified in the Hologram Processing page. As the extracted centroid .csv files contain a column specifying the size of the centroids, it is possible to apply a further size filter by specifying the Min Particle Size and Max Particle Size in the Coordinate Filtering group box. The user can specify the maximum displacement considered for tracking in the Max Displacement field. It is possible to remove trajectories shorter than a certain length by specifying Minimum Trajectory Length. This is typically used to denoise the tracks, as noise is temporally incoherent. In principle, it therefore cannot lead to long tracks. It is also possible to set the particle memory of the algorithm by setting the Memory field. This decides the number of frames for which the ID of a particle is remembered. If the particle disappears and reappears, it would be treated as the same particle as long as it reappears within the specified memory window. As the tracking algorithm is implemented for any general data, it can be applied to two-dimensional data by setting Dimensions to 2. The final tracks are saved as a file named tracked_particles.csv which contains $x,y,z,\text{Frame},\text{ID}$ for each particle.

2.1.4 Visualization

In the Visualization tab, the user may choose to visualize the computed results. There are three visualization panes. By default, the directories are set according to the settings specified in the Hologram Processing tab. The first pane is for hologram visualization. Using the horizontal slider in the hologram pane, the user can swiftly scroll through all the holograms in the specified folder. If the backgrounds have already been computed for the data, it is possible to check the Subtract BG option to display the background-subtracted image. In the Contrast group box, the user can specify the highest and lowest values for the image coloring. By default, the range is set to 0–1 as the images displayed are normalized to be in the unit interval. In the same spirit, if the maximum phase projections have been computed, it is possible to visualize them in the maximum phase projection display pane. By enabling the Sync checkbox, the sliders for the hologram and the maximum phase projections can be synchronized. The projection plane being displayed can also be changed by selecting the desired choice in the Projection Axis group box. If the tracks are computed, they can be displayed in the tracks pane. The end time and the number of tracks displayed can be controlled by their respective sliders. To update the tracks displayed, the Plot button can be used. The number of tracks displayed and other information about the tracks is displayed in real-time in the Output Console page.

2.1.5 Console

For all the processes, the text output is displayed in the console present in the Output Console tab. If the user wishes to stop the processes, the console also contains a button to stop all processes. It is also possible to save the console output log by clicking on the save log button. By default, the log file is saved in the input image directory. The console can also be cleared using the clear console button. Whenever a process is initiated, the program automatically switches to the Output Console tab to display the console output to the user.

2.1.6 Batch Automation

Typically, one wants to perform all three steps: background subtraction, hologram processing, and tracking on the data. It is possible to perform these operations in a sequence automatically by selecting all the desired options on their respective page, and then selecting the steps to be performed in the Processes to Execute group box present in the Batch Automation tab. It is also possible to automatically process various datasets. This can be done by creating a .rihvr file for each dataset by selecting the appropriate settings in the GUI and saving the state. Once all the .rihvr files are created, the user can simply provide the path to the folder containing all the files in the directory of .rihvr files field. All the cases specified by the .rihvr files in the specified folder will be automatically processed in sequence. This can be of use when the user wants to test the effect of different parameters or process different datasets. Further information regarding the .rihvr files can be found in Sec. 2.2.1.

2.2 Miscellaneous Routines

The program uses multiple different routines to save and load data, as well as interface handling. These routines are not a part of the core computation algorithm. Therefore, these routines are described in this section.

2.2.1 Save/load .rihvr files

For convenience, the program can save and load state from a .rihvr file. These files are JavaScript Object Notation (JSON) type files that store information about all the fields specified in the GUI. Therefore, one can share the same JSON file with different users to exactly reproduce all the results. These files can be saved and loaded using the Save State and Load State buttons in the Menu. There is a default .rihvr file that is loaded at the startup of the program. By default, it is saved in the installation directory and is named as defaultsettings.rihvr. If desired, one can set a different .rihvr file for the startup. This can be done by changing the <USER>\Documents\RIHVR_FFIL\startuppath.txt or inside the settings dialog box that can be accessed through the Menu.

The routines that handle the saving and loading of the .rihvr files are defined inside the rihvr.cpp file. Code 1 is a pseudo-code for the save state function defined as saveStateToFile. Shortly, the settings from the settings dialog box are queried. These are stored in the registry of the operating system. Then the JSON object is built and populated using all the fields and parameters from all the different pages of the GUI.

Pseudo-code 1

Similarly, to load the state, the loadStateFromFile function is used. The pseudo-code for loadStateFromFile is given in Code 2. The function is called automatically to load the default file at startup, and can be called manually by clicking on the load state button in the Menu.

Pseudo-code 2

2.2.2 Track visualization

The visualization of the tracks in the Visualization page is produced using OpenGL. Firstly, the tracks are loaded from a .csv file with columns $x,y,z,\text{Frame},\text{ID}$ using the routine named loadTracksFromTextField, which calls the loadTracksFromCSV function to sort the tracks into an array. The visualization of the tracks is also filtered based on the selected time and concentration value via sliders using the updateTrackDisplay function. These three functions are defined inside rihvr.cpp. Once data is sorted and filtered, it is sent to the tracksPlotter object and the setTracks operation is called, which plots the tracks. Further details on the plotting logic can be found in the tracksplot.cpp file, which contains the OpenGL implementation of the tracks plotting operation. The tracks are colored according to the velocity magnitude. A bounding box is drawn based on the maximum extent of the tracks being plotted. A corner of the bounding box displays its size for an easier interpretation.

2.2.3 Settings Dialog

The user can access some important settings for the program via the Settings button in the Menu. For proper functioning of the program, it is crucial to select the correct paths for all the settings. As the processing executable and the GUI make use of OpenCV, a path to the OpenCV dynamic libraries is required. By default, the OpenCV path is set to be inside the installation directory, as the opencv_world4130.dll is shipped with the Windows redistributable of the application. The path to CUDA libraries is set to the CUDA folder inside the application folder, as the required CUDA libraries are also shipped with the redistributable for Windows. By default, all the executable files are stored inside the lib folder inside the program directory. There are five executable files in total. Namely, ccbg, particletracking, make-background, segmentation, and sparse-inverse-recon. If the user decides to compile the executable files from the source, these five files are generated. The path to these executables is set in the settings dialog. By default, all the settings are saved as a key inside the registry of the operating system. Code 3 shows the implementation of the registry writing inside the settings.cpp file. By default, the installer writes the correct registry values according to the installation settings. A user should therefore never need to change anything in the settings. It should be noted that as the .rihvr file stores the values from the settings dialog as well; these settings will change if the user loads a .rihvr file generated on a different machine with a different installation directory. If a user desires to load .rihvr files from a different machine, it is recommended to open the .rihvr file in a text editor and delete the JSON entries corresponding to the path settings. Then the settings on the user’s machine would remain unchanged when that modified .rihvr is loaded.

Pseudo-code 3

2.3 Algorithm

In this section, the full code workflow is described. The core steps are illustrated in Figure 2, and the details are discussed below. The computations for each step are conducted by the executable files listed under each block in green boxes in Figure 2.

2.3.1 Background

Once all the parameters are set inside the GUI, and the Compute Background button is clicked, one of the three functions computeBackground, computeCCBackground, or computeMovingBack ground is called based on the selected background type. These functions are defined inside the rihvr.cpp file and act as monitors and callers for two executables: make-background and ccbg. When any background computation needs to be performed, a .params file is generated inside the params folder in the specified Background Path. These .params files vary according to the background type selected and the image number. The executable binaries take these .params files as inputs to perform computations. The mean, median, and moving average background computations are performed using the same executable make-background by making appropriate changes to the .params file. Code 4 shows the logic defined inside the main_create_background.cu file that is at the heart of the make-background executable. The code executes one of the two branches: median or mean (standard), based on whether the user has selected the median background as an option. If the median background is selected, a temporal median of 5 frames is computed and stored as a background. Therefore, for N frames, only N–3 backgrounds are computed. During hologram processing, for the last three images, the (N–3)th background is used. If the median filter is not selected, an average of all the images in the sequence is computed. For the moving average background, the size of the image sequence input to the function is varied by changing the input .params file in a loop by the computeMovingBackground routine.

Pseudo-code 4

As described in Sec. 2.1.1.5, the correlation-based background subtraction involves inherently different operations; this background generation is performed by the ccbg executable. When Correlation is selected as the background type, clicking the Compute Background calls the computeCCBackground routine inside the rihvr.cpp. This routine creates a .params file as an input for the ccbg executable. A pseudo-code for the correlation-based background computation of the logic followed by the ccbg executable is given in Code 5. The entire implementation can be found in ccbg.cpp. hlThe images are downsampled according to the user-specified Resize Factor and a temporal correlation matrix is computed between the downsampled images. For each image, the most temporally correlated frames (not downsampled) are added until the entropy of the computed background monotonically increases to decide the maximum number of images to consider for generating the background for a particular image. It is possible to overwrite the entropy-based estimate by providing a bound for the number of frames considered for background averaging by specifying Min Frames and Max Frames in the GUI as described in Sec. 2.1.1.5.

Pseudo-code 5

2.3.2 Processing

Once the background is generated, the holograms are processed according to the parameters set inside the Hologram Processing tab. Clicking the Process Holograms button calls the compute HolographicReconstruction routine defined inside the rihvr.cpp file. This routine generates the .params file for the sparse-inverse-recon binary. The main function for this binary is defined inside the main_sparse_inverse.cu file. This file essentially handles the overall process of loading the holograms, calling the reconstruction functions, and saving the results. The core logic of the reconstruction algorithm is defined inside the sparse_compressive_holo.cu. It relies on other core files such as reconstruction.cu for generating the OTF and its adjoint as shown in Eq. 4 for the forward and back-propagation. It also uses the sparse_gradient.cu for gradient computation, and other files for handling the volumes, object field, etc. The reader interested in further details is encouraged to read through the codebase itself. However, the core logic of reconstruction processing adapted from sparse_compressive_holo.cu and reconstruction.cu is shown in Code 6. The OTF and the adjoint functions are defined inside reconstruction.cu, and are used to compute the forward and gradient model for a hologram. The optimization is performed under the FISTA reconstruction loop, where the step-size is estimated using the non-monotone backtracking (golden-section search) method. The estimated hologram is changed according to the gradient. The number of iterations is specified by the user in the Number of FISTA iterations field in the GUI.

Pseudo-code 6

2.3.3 Segmentation

Once the sparse volume fields are reconstructed, the segmentation binary is executed (if selected) to extract centroids from the voxel field data. The binary takes in the inverse_voxels.csv file from the Output Directory and segments the volume field. The core logic is implemented inside main_segmentation.cu, which makes use of the functions defined inside sparse_segmentation.cu and sparse_volume.cu. The basic algorithm is shown in Code 7. The voxel data CSV file is loaded and binarized according to the parameters set in the GUI. After binarization, a morphological closing operation is performed to connect nearby components and to fill gaps. Once the connected components are identified and labeled, centroids are extracted. Before saving the extracted centroids, a blob-size filter is used to remove objects smaller than the specified size.

Pseudo-code 7

2.3.4 Tracking

Once the centroids are computed, the user may choose to compute the particle tracks using the built-in Crocker-Grier algorithm [19]. When the Compute Tracks button is pressed, the computeTracks function defined inside rihvr.cpp is called. This function call generates the .params file to run the particletracking executable. The particletracking executable is defined inside the particletracking.cpp file, which handles the loading, size-based filtering, and saving of the particle data and tracks. The core of the tracking algorithm is implemented inside tracking.cpp routine. A pseudo-code is presented in Code 8. The algorithm assigns particle IDs based on the nearest neighbor logic. Firstly, it is checked if the current frame being processed is trivial or not. If it is found to be trivial, the trivial bond computation branch is activated; otherwise, the computation is done on the non-trivial branch, which requires a large number of permutations.

Pseudo-code 8

2.3.5 Summary

In summary, the GUI calls the make-background or ccbg executable to compute the backgrounds via computeBackground, computeCCBackground, or computeMovingBackground routines. Once the background computation is done, sparse-inverse- recon binary is executed by the compute HolographicReconstruction function, which also calls the segmentation binary if required. Once the centroids are extracted, the user can run the Crocker-Grier 3D tracking algorithm to obtain particle tracks. To do so, the GUI calls the computeTracks function, which generates the .params file to run the particletracking executable.

3 Quality control

The program has been extensively tested on various synthetic as well as experimental hologram data [14]. The method has been applied in the context of various scenarios relating to microbiology and fluid mechanics by Mallery in his thesis work [24]. For testing, multiple datasets, both synthetic and experimental, have been made available by the authors with the code that contains hologram images along with the ground truth to evaluate the accuracy of the algorithm. All processing in this section was performed on a Windows 11 machine with an NVIDIA RTX 2080 Ti (10 GB) GPU, an Intel i9-9820X CPU, and 64 GB of system memory.

We first look at the results obtained from processing a hologram of two dandelion stems placed at different axial distances [23]. For further details, please see Brady et al. [23]. Figure 7(a) shows the raw hologram and the maximum phase projection in the (b) XY plane, (c) YZ plane, and (d) XZ plane. For reconstruction, the following parameters were used: z0 = 0 µm; dz = 800 µm; dx = 10.4 µm; Nz = 100; lambda = 0.633 µm; 0-padding = 64; Sparsity = 0.001; TV = 0.02; Number of FISTA Iterations = 100; Number of TV Iterations = 5. The results shown in Figure 7 demonstrate that the developed software is able to reconstruct and resolve the sharp branches of the stems without any significant gain in noise. It is also observed that the axial position of the dandelion stems is well resolved, as evident from the sharp localization of the dandelion planes in Figure 7(c) and (d). The processing time was found to be approximately 12 s. for this single frame.

Figure 7

To quantify the accuracy of the reconstruction and tracking algorithm, the developed software is tested on synthetic holograms generated using the 3D homogeneous isotropic turbulence (HIT) flow from the Johns Hopkins Turbulence Database (JHTDB) [25]. The reader is referred to Mallery and Hong [14] for further details regarding the synthetic dataset. Figure 8 presents (a) a sample background-subtracted hologram from the synthetic 3D-HIT sequence, (b) a temporal maximum of the maximum phase projection in the XY plane with the tracks overlayed, and (c) 3D tracks computed for the first 50 frames from the dataset. Firstly, from (b), it is clear that the volumetric reconstruction algorithm is able to extract the position of the point particles consistently over time. It is observed that the identified tracks, colored lines in (b), overlay exactly on the phase projection shown in (b). The tracking algorithm is therefore able to connect all the tracks even in such a turbulent flow where the trajectories cross each other. Furthermore, it is evident from Figure 8(c) that the turbulent features of the flow are clearly captured. Therefore, the developed software can be used to investigate the turbulent flow field using holographic PTV data. Our machine took about 75 s to reconstruct a single volume field from these holograms.

Figure 8

To quantify the error in localization and tracking, we present (a) the percent error in the magnitude of the estimated root-mean-square (RMS) velocity from the HIT dataset and (b) the error in position localization via reconstruction, deconvolution, and RIHVR (presented software) in Figure 9. Please note that here reconstruction refers to finding the location of the waist of the particle fringes in the back-propagated hologram stack— as the fringes for a particle converge at its focal plane [2]. Deconvolution refers to deconvolving the hologram using the point spread function (PSF) of the system [2]. It is observed in Figure 9(a) that RIHVR provides the best axial velocity (w) accuracy in terms of the estimated RMS. For the in-plane motion (u and v), all the methods compared perform almost identically, as the in-plane location of the particles can be easily identified from a hologram, as evident from Figure 8(a). Furthermore, the accuracy of the localization performed is also evaluated and compared to the known ground truth for all three methods. Figure 9(b) presents the error in the estimated particle positions in all three directions. The solid bar represents the standard deviation around the mean, while the vertical error bar represents the maximum and minimum error values. For the in-plane directions (x and y), all three algorithms perform identically. However, for the axial direction (z), it is clear that not only RIHVR leads to more accurate results (i.e. a lower value of the mean error), but also more precise results as evident from the smaller standard deviation and the extrema of the error values.

Figure 9

In the context of fluid mechanics, the performance of the processing methods can significantly vary depending on the synthetic data and the experimental data [26, 27, 28, 29, 30, 31]. Therefore, we test the developed software on experimentally captured PTV holograms. The holograms were recorded in a T-junction flow, and for further details regarding the dataset, please see the work of Mallery and Hong [14]. Figure 10 presents (a) a background-subtracted hologram from the T-Junction sequence, (b) a temporal maximum across the first 100 computed maximum phase projections in the XY plane to visualize the motion of the particles with tracks overlayed, and (c) 3D volume reconstruction. Correlation-based background enhancement was used to compute Figure 10(a). For hologram processing, the following parameters were used: z0 = 2500 µm; dz = 2 µm; dx = 1.1 µm; Nz = 1000; lambda = 0.632 µm; 0-padding = 0; Sparsity = 0.15; TV = 0.07; Number of FISTA Iterations = 100; Number of TV Iterations = 5. From (b), it is clear that the software is able to extract the exact shape and structure of the seeded particles coherently across a sequence of holograms without any major noise. Consistent and coherent trajectories are clearly evident from Figure 10(b), which shows that the extracted trajectories exactly follow the temporal maximum phase projections. By examining the 3D volume field in Figure 10(c), the tumbling motion of the rods can also be observed. This demonstrates the potency of the algorithm in not only extracting the location, but also the 3D orientation of the seeded rods in the volume. Therefore, the developed software is able to extract the exact rod-like shape of the particles along with their orientation in 3D space. These results also demonstrate the potential use of the developed software in the context of particle diagnostics, where the exact shape, orientation, and size of particles need to be measured. On our machine, the holographic reconstruction for each frame was completed in about 40 s.

Figure 10

To show the versatility and the wide-scale applicability of the developed software, we also test it on holograms recorded in the context of microbiology. One thousand holograms of wild-type Vibrio cholerae C6706 cells swimming in tryptone broth in an open well at room temperature were captured at a frame rate of 120 Hz with 0.108 µm per pixel resolution. Figure 11 presents (a) a background-subtracted hologram from the sequence, (b) a temporal maximum of the maximum phase projection in the XY plane with the tracks overlayed, and (c) 3D tracks computed for the first 1000 frames from the dataset. For hologram processing, the following parameters were used: z0 = 0 µm; dz = 0.216 µm; dx = 0.108 µm; Nz = 150; lambda = 0.304 µm; 0-padding = 0; Sparsity = 0.1; TV = 0.1; Number of FISTA Iterations = 50; Number of TV Iterations = 5. For tracking, the following parameters were used: Min Particle Size = 110, Max Particle Size = 10000, Max Displacement = 10, Minimum Trajectory Length = 200, Memory = 20. From the maximum phase projections presented in Figure 11(b), it is evident that the method is able to accurately extract the shape and size of the bacteria coherently across time. Furthermore, it is also observed that the tracked trajectories exactly overlay the temporal maximum of the maximum phase projection. Therefore, the tracking algorithm is able to accurately extract the trajectory of the bacteria. The smoothness and the coherence of the trajectories are clear from the 3D trajectories presented in Figure 11(c). Hence, the developed software is able to extract the exact shape and structure of the microbes coherently across a sequence of holograms without any major noise, while also computing highly accurate trajectories of the swimming bacteria.

Figure 11

From the results presented in this section, it is clear that the developed software is capable of extracting both large and small-scale volumetric structures. It is also evident from the results presented for the synthetic 3D HIT dataset that the method outperforms traditional reconstruction and deconvolution in terms of both accuracy and precision. The dandelion results show the algorithm’s capability to reconstruct large, smooth, and sharp structures like dandelion stems, while the T-Junction and 3D-HIT results show that the algorithm can accurately reconstruct particle fields across time. Hence, it can be used for 3D PTV even in turbulent flows. The results presented for swimming Vibrio cholerae demonstrate that the software can also be used to track microbes and to investigate the flow fields around them. Our machine took about a minute on average to process each of the frames.

Although the developed software has been shown to be highly accurate, it has its limitations. For instance, as discussed by Mallery and Hong [14], at extremely high particle concentrations, the method can struggle with extracting all the particles. Even though the method makes use of sparse matrices, it still requires a large amount of GPU VRAM. For a $1024\times 1024\times 1000$ pixel reconstruction volume, the memory requirement can easily get over 10 GB. The tracking code implemented is based on a nearest-neighbor assignment, which can struggle with accurately linking particle trajectories when the particle concentration is very high. As evident from the discussion in Sec. 2.1.2.1, the solution to the minimization problem depends heavily on the user-specified Lagrange parameters λ_TV and ${\lambda }_{l^1}$. These parameters need to be carefully adjusted by the user, which can pose a barrier to entry. The developed program also assumes that the object field to be reconstructed is sparse and smooth, which may not be true in all cases, for example, if the medium has fluctuations in the index of refraction. Finally, the FISTA optimization routine assumes that the regularized minimization problem is convex. This assumption typically breaks down for large objects when a large number of reconstruction planes are used. Although these restrictions generally do not apply to PTV applications, the user must be aware of these limitations during experimental design and data processing. The future work would focus on quantifying the uncertainty in the obtained results along with the sources of error, so that the experiment design process can be standardized.

Operating system

Windows, UNIX/Linux

Programming language

Cuda 12.9, C++ standard 17, QT6, upward compatible

Additional system requirements

CUDA 12.9 compatible GPU with a compute capability of 7.0 or higher. 1 GB disk space, 4 GB RAM.

Dependencies for compilation

Nvidia CUDA toolkit 12.9 or higher; OpenCV 4.0 or higher; C++ compiler with standard 17 support; Microsoft Visual Studio 2017 or higher; glib v2.35 or higher; QT6.

Dependencies for executing the shipped binaries

Nvidia GPU drivers supporting CUDA 12.9 or higher; glib v2.35 or higher; Microsoft Visual Studio 2017 or higher runtime libraries.

Software location

Archive An archive version of the software is made available as an OSF repository.

• Name: OSF

• Persistent identifier: https://doi.org/10.17605/OSF.IO/QWRDP

• Licence: GNU-GPLv3 license

• Publisher: Gauresh Raj Jassal

• Date published: 27th January, 2026

Archive An archive version of the software is made available as a Zenodo repository.

• Name: Zenodo

• Persistent identifier: https://doi.org/10.5281/zenodo.18991876

• Licence: GNU-GPLv3 license

• Publisher: Gauresh Raj Jassal

• Date published: 12th March, 2026

Code repository The code is hosted in a GitHub repository.

• Name: Github

• Persistent identifier: https://github.com/GRJ-Satyakama/RIHVR—Regularized-Inverse-Holographic-Volume-Reconstruction

• Licence: GNU-GPLv3 license

• Date published: 27th January, 2026

Language

CUDA, C++, and QT