1. Introduction

Tomographic Volumetric Additive Manufacturing (TVAM) [1,2] is a novel technique allowing for rapid 3D printing. As shown in Fig. 1(h), it is based on a computed tomography principle (CT) which consists of projecting 2D light patterns from many angles into a photosensitive resin. The light patterns build up a cumulative 3D energy dose distribution within the photosensitive resin. Above a specific light dose received at a point in space, enough conversion of the monomer into polymer occurs, thereby solidifying the resin. To compute this set of 2D light patterns, most recent works have built upon the Radon transform [3]. A common approach is to voxelize the 3D model, calculate the tomographic projections using an implemented Radon transform function, and then filtering these projections. The Radon transform assumes un-attenuated, un-distorted, straight ray propagation, which differs from the physical implementation. The naive approach can be corrected to compensate for lack of telecentricity and limited etendue [4]. To include absorption of the ray upon propagation, the attenuated Radon transform [5] is the correct model. Fig. 1 shows an example of the computed tomographic patterns to construct the 3DBenchy boat (3DBenchy.com) for a slice. The object shown in Fig. 1(e)), after thresholding, consists of solidified pixels.

 

Fig. 1. The general principle behind TVAM. a) a set of 2D projection patterns is propagated into space. b) shows how a slice of the pattern propagates through the volume and c) how the incoherent sum results in a total energy dose. d) the object polymerizes if it reaches an energy threshold. e) polymerization threshold results in a printed slice. f) is the intensity histogram of b). g) is the 3D view of the Benchy boat. h) is the general setup.

Download Full Size | PDF

Finding the optimal set of patterns for the projector has been a problem that is being tackled since the introduction of TVAM. Gradient descent-based methods have been proposed in several different works [1,6,7], with various loss functions used to achieve different properties for the prints. The object-space optimization algorithm proposed by Rackson et al. [8] is a conceptually fast algorithm and can be implemented without explicit gradient descent formulation. The forward model in this work is the attenuated Radon transform. Orth et al. presented a deconvolution approach for TVAM [9] that accounts for diffraction and for the effect of chemical diffusion of polymerization inhibitors as a convolutional kernel in the optimization scheme. Although this method introduces a diffraction model, it is still based on ray optics by treating a single ray, as not spatially sharp but instead describing it as a Gaussian intensity function over a spatial region.

Further, a general loss function for TVAM has been proposed to incorporate material specific effects or grayscale dose into the model [10]. The same work also mentions the possibility of using a wave optical operator but does not provide a methodology. In addition to this work, a multi-beam optimization wave optical framework was recently introduced [11] which has also not been experimentally verified. The methodology and the optimization approach is similar to our work, but it is restricted to phase optimization and does not reveal details about the implementation. Their simulated volume is discretized with $N_\text {pixels}=256$ and $N_\text {angles}=64$. Hence, the simulated volume is around 10 times smaller and the total amount of optimized pixel values is 40 times smaller than in this work. Moreover, larger volumes and more beams is not possible in their implementation because of memory issues.

Here we introduce a scalar wave optical model which includes wave optical effects to TVAM that bears many similarities with diffraction tomography [12,13] which introduced the wave nature of light to CT problems at micrometer resolution. Instead of tracing straight rays through the volume, we propagate the 2D projections step-wise through the volume with an angular spectrum method of plane waves (AS). This step is repeated for each angular projection pattern, and the intensities are then summed incoherently. This is similar to multi-slice diffraction tomography [14] except that we do not take into account the variation of the refractive index during printing. Also we do not consider material properties such as the chemical diffusion of inhibitors; which other works have modeled computationally [15]. Based on efficient algorithms with custom gradient rules written in the programming language Julia [16], we are able to simulate volumes as large as $550 \times 550 \times 550$ voxels with $600$ projection angles on a single NVIDIA A100 GPU. The comparison between wave optical and ray optical methods reveals that at 20 µm the wave optical nature impacts the printing results noticeably.

2. Theory

2.1 Angular spectrum of plane waves method

Free-space propagation allows for propagating vectorial components of the electrical field $\psi$ in a uniform medium with refractive index $n$ by solving the homogeneous Helmholtz equation [17]

2.2 Pattern optimization

Given the wave optical forward model, we can search for a set of patterns which minimizes a certain loss function for optimal prints. The loss function we use is inspired by [8]. Additionally, we introduce an extra term to avoid overexposure of certain voxels. Mathematically, we seek a optimal set of patterns to minimize

2.3 Implementation details

Pseudocode for the full wave optical optimization algorithm is sketched in Algorithm 1. The basis is a gradient descent optimization scheme. Starting with a set of 2D amplitude patterns for each angle $\varphi$, the patterns are propagated into space. After propagation of all angular patterns, the loss function given by Eq. (5) is calculated. Then, the gradient of $\mathcal {L}$ with respect to the amplitude patterns $A$ is calculated. With the gradient, we can perform one optimization step. This optimization is restricted to amplitude optimization currently. Optimizing phase and amplitude at the same time is left for future work.

Algorithm 1. Gradient descent based wave optical optimization algorithm.

View Table | View all tables in this article

Equation (4) is the key equation for the wave optical forward model. For the AS method, in fact we use the bandlimited AS [19] implementation to prevent any wrap-around artifacts. The computational complexity to propagate a single 2D field $\psi _\varphi$ of size $(N \times N)$ into the 3D volume ($N$ different $x$ values) is given by $\mathcal {O}(N^3 \cdot \log N)$. This is because we propagate the field (size $N\times N$) $N$ steps into the volume and each of the FFTs evaluates at $\mathcal {O}(N^2 \cdot \log N)$.

Algorithm 2. Implementation of the wave optical forward model $\mathcal{P}_W$.

View Table | View all tables in this article

The amount of angles $N_\text {angles}$ covering $2\pi$ should be chosen approximately as $N_\text {angles} \approx N \cdot \pi$ (to avoid angular subsampling) which results in a total complexity of $\mathcal {P}_W$ equal to $\mathcal {O}(N^4 \cdot \log N)$. This is promising because the Radon based optimization method has a computational complexity of $\mathcal {O}(N^4)$. It is obvious that a Radon based implementation will be faster in practice since the AS method in case of $\lambda \rightarrow 0$ evaluates to an identity matrix and is computationally trivial. Apart from $\log N$, the computational complexity is the same between wave optical and ray optical forward model.

In Algorithm 2 we show the pseudocode for the full wave optical propagator. The most time consuming operation is line 5. First, a 2D FFT of the field is calculated and then in Fourier space the 3D kernel is multiplied and finally an inverse batched 2D FFT is calculated. Each resulting volume is rotated with the rotation operator $\mathcal {R}_\varphi$ and added to the total intensity in the volume.

Algorithm 1 requires to adjoint the gradient of $\mathcal {P}_W$ with respect to the amplitude patterns $A$. However, most reverse mode automatic differentiation (AD) engines such as the common ones in Julia (Zygote.jl [20]) or PyTorch [21] fail to calculate the gradient of the loss function in memory efficiently. In line 6 of Algorithm 2, the pullback (backpropagator for the reverse mode AD) of $|\psi _\varphi (x,y,z)|^2$ is $2 \cdot \psi _\varphi (x,y,z) \cdot \Delta I_\varphi ^{R}$. In consequence, the AD needs not only to store the sum of all angular backprojections (memory consumption $\mathcal {O}(N^3)$) but also the result of each individual backprojection. This would require a memory of $\mathcal {O}(N^4)$. This memory scaling is beyond the memory limit of modern GPUs. For example, in our naïve PyTorch implementation volumes larger than $N\gtrsim 100$ were not possible. Instead of keeping all $\psi _\varphi (x,y,z)$ in memory, we recalculate those fields in the gradient pass. This method is called memory checkpointing. For the pseudocode of the whole adjoint of $\mathcal {P}_W$ see Algorithm 3.

Algorithm 3. Adjoint of the wave optical forward model $\mathcal{P}_W$.

View Table | View all tables in this article

All algorithms were implemented in Julia and made use of ChainRules.jl to register a custom adjoint for the critical $\mathcal {P}_W$ operation. Further, all routines feature CUDA acceleration in Julia [22] and the bilinear differentiable image rotation routine is implemented in KernelAbstractions.jl [23] which results in fast CUDA executable code too. As optimizer, we used the L-BFGS algorithm [24] implemented in Optim.jl [25]. L-BFGS is a quasi-Newton method which builds up an estimate of the inverse Hessian which helps minimizing the objective in less steps. We encourage to use Fabio Crameris color-blind friendly color maps [26].

3. Simulation results

In this section we want to highlight three different simulations. Each of them was computed at $\lambda = 405 \text{nm}/n$ where $n=1.5$ represents the refractive index of the resin. The discretization is $N=550$ with $N_\text {angles}=600$. The absorption was set to 0, which implies that a set of angles in the range $[0, \pi )$ is sufficient. To compare the wave optical optimization with the classical ray tracing approach, we optimize total volumes of sizes $L_1=100\;\mathrm{\mu}\textrm{m}$, $L_2 = 1\;\textrm{mm}$ and $L_3 = 10\;\textrm{mm}$. Each optimization was performed with Julia 1.10 on a NVIDIA A100 with 80GB GPU RAM. The runtime was roughly 4 hours for the wave optical optimization. The ray tracing optimization takes 5 minutes. The L-BFGS optimizer requires usually around 40 iterations to achieve a well separated histogram comparable to Fig. 1(f). This work sets a computationally efficient and differentiable basis upon the rigorous pattern optimization can be based on. It is worth to note that, without the memory checkpointing the algorithm would not run at all because the hardware would run out of memory.

The results of the simulations are shown in Fig. 2. Each row represents a different resolution of the boat. The wave optical forward model results in the prints in column d). The Intersection over Union (IoU), a metric of fidelity between the target model and the resulting polymerized region, was always 1, so the print could be considered as perfect; in practice, absorbance from the photoinitiatior, chemical diffusion, and the time dependence of the refractive index upon polymerization would hinder resolution [9,27,28]. In Fig. 2(e)) we use the ray optical patterns as input for the wave optical forward model and apply the threshold of 0.75. For $L_3=10\;\textrm{mm}$ the wave optical and ray tracing optimization yield almost identical results. This is expected since the wave optical nature is not dominating at those resolutions. If we divide the physical dimension by 10, we see no substantial difference between the ray optical patterns and the wave optical patterns. Only the small feature in the top (purple arrow) is not printed by the ray optical patterns. However, at $L_1=100\;\mathrm{\mu}\textrm{m}$ the ray optical patterns start to misprint the object significantly. The wave optical patterns result in a perfect print. In Fig. 3 we plot the general trend of the IoU for different resolution scales of the boat. The wave optical simulation was carried out for each datapoint. Especially for $L<200\;\mathrm{\mu}\textrm{m}$ the IoU drops significantly. Quantifying the resolution limit of ray optical TVAM is not straightforward since it depends on the printed sample. But from this target we infer that below 20 µm wave optical effects impact the print.

 

Fig. 2. Overview of the simulation results for different printing sizes. a) is one 2D pattern. b) is the slice of this pattern propagated into the volume. c) indicates the incoherent sum of the patterns. d) after applying a relative energy threshold of 0.75. e) are the ray optical patterns propagated with the wave optical formalism and with the same threshold applied.

Download Full Size | PDF

 

Fig. 3. Intersection over Union for the different print results. The ray optical patterns are propagated with the wave optical model.

Download Full Size | PDF

4. Discussion

This presented scalar wave optical modeling shows that in principle an amplitude modulated TVAM is able to print micrometer features inside a single photon absorptive resin. Phase optimization only offers higher light efficiency but does not increase the resolution. Our simulations also reveal that the recently printed 20 µm features [29] are most likely limited by the ray optical assumptions underlying the existing optimization schemes. Of course, the next step is to experimentally demonstrate that wave optical patterns indeed allow for micrometer features. This poses several challenges. At the moment, any refraction of the glass vial (or an index matching bath) is ignored. It remains unclear, how severe the modulation of the glass vial on the electrical field is. Describing strong refraction of glass interfaces with a wave optical model at millimeter scale is a computationally very demanding problem. Further, the rotating vial in TVAM has to be digitally or mechanically stabilized to account for any motion blur. Also, as recently demonstrated [9], chemical diffusion can have a major impact on resolution too. During printing, the refractive index inside the volume changes according to the accumulated light dose. The scattering effect caused by the refractive index changes is not included in the wave model yet and is left for future research. However, it can be foreseen that the gradient rule requires special handling because otherwise it will run into similar memory issues as discussed in this manuscript. Gravity and buoyancy acting on the regions with different density during printing become more important the longer the printing time is. Also, a digital micromirror device (DMD) can not be used for experimental implementation. A DMD is not a true grayscale amplitude modulator but instead binary patterns are time multiplexed to achieve grayscale intensity patterns. Our model assumes continuous and coherent patterns which can be achieved with a true grayscale liquid crystal display (LCD) for example. Most likely, micrometer TVAM requires to take all those effects into account for reproducible printing.