© The Author(s), published by EDP Sciences, 2026

1 Introduction

Structural colors have emerged as a promising inkless approach for generating durable and versatile colors with high spatial resolution [1, 2]. They originate from the interaction of light with subwavelength architectures that shape the spectral response through well-identified physical mechanisms. Depending on the considered material and structure, coloration can arise from resonant absorption and scattering in plasmonic or dielectric metasurfaces [3–5], diffraction from periodic patterns [6], Bragg reflection in photonic crystals [7], or interference effects in multilayer and cavity-based systems [8–10]. Fabrication strategies for structural colors patterning typically rely on top-down approaches such as electron beam lithography or nanoimprinting, which provide precise control over the geometry. However, these approaches remain intrinsically limited in terms of throughput, scalability, and compatibility with large-area or non-planar substrates. In this context, laser processing offers an effective route toward structural coloration compatible with industrial-scale applications. By inducing localized and controlled material transformations, it enables the direct formation of nanostructures in a single-step, non-contact process. This approach is compatible with a wide range of materials, large areas, complex geometries, and in-volume structuring, while maintaining high processing speeds and flexibility.

Laser-induced structural colors rely on several physical mechanisms depending on the initial material and structure. Exposing discontinuous metallic layers, or dielectric matrices containing metallic precursors, to a laser can induce the formation of metallic nanoparticle assemblies with controlled statistical properties. These random plasmonic metasurfaces exhibit optical responses governed by localized surface plasmon resonances and their collective interactions [11–18]. In parallel, laser-induced periodic surface structures can form through the interaction between the incident field and surface electromagnetic waves leading to diffractive effects on a wide range of materials [19, 20]. In thin film configurations, laser-induced changes in thickness or composition, for instance through oxidation or phase transformation, modify the optical path and generate interference-based coloration [21, 22]. More advanced approaches combine these mechanisms to further tailor the spectral response [23, 24]. These developments have enabled large-area color printing, high-resolution patterning, and the fabrication of functional optical surfaces for applications in decoration, optical encoding, and anti-counterfeiting. In particular, laser-induced multiplexed image printing, where different images are revealed under different observation conditions, provides a powerful approach for secure document personalization [25–27].

Despite these advances, laser processing remains limited in its ability to fully control the optical response. The physical parameters governing color generation cannot be tuned independently, and the resulting color gamuts are generally restricted compared to standard color spaces such as sRGB. This results in missing hues and reduced chromaticity. While gamut mapping strategies can be used to adapt images to the available color palette [28], they do not compensate for the intrinsic limitations of the process. Maximizing the accessible color gamut therefore remains a key objective. Achieving this requires identifying optimal laser processing parameters, such as power, scan speed, or repetition rate, within a high-dimensional parameter space. However, the relationship between these parameters and the resulting colors is highly nonlinear and sensitive to material variability. Advances in ultrafast laser processing have highlighted the critical role of precise parameter control in achieving reproducible and high-quality material structuring [29]. For industrial applications, the ability to rapidly adapt to new materials or to variations between samples is therefore essential. In this work, we address this challenge by developing an efficient method to identify, with a minimal number of experiments, the optimal laser parameters that maximize the achievable color gamut for a given material. While similar optimization problems have been addressed using genetic algorithms [30, 31], we show that Bayesian optimization provides a more efficient framework for this task.

While Genetic Algorithms (GA) and Bayesian Optimization (BO) are well-established global optimization techniques, their application to laser material processing has gained increasing attention in recent years. This is largely due to their ability to work on high-dimensional spaces and to efficiently explore parameter spaces without relying on gradient information, making them well suited to the nonlinear nature of laser-matter interactions. Recent studies have demonstrated the efficacy of surrogate-based optimization methods for complex laser processing tasks. For instance, hybrid frameworks combining physical models, such as the two-temperature model, and machine learning have been proposed to predict ultrafast laser ablation depths under data-scarce conditions [32]. Similarly, Gaussian Process-based surrogate models have been successfully employed for the inverse design of complex metasurfaces, significantly reducing the number of computationally expensive electromagnetic simulations [33, 34]. Beyond simulation-driven studies, surrogate-based models are increasingly applied to laser manufacturing processes, including the identification of process windows for laser-cutting [35], the exploration of parameter spaces for laser-induced graphene synthesis [36] or the optimization of photonic surfaces texturing [37]. However, maximizing the structural color gamut presents a unique challenge: the mapping between laser parameters and the resulting colors is highly nonlinear and strongly affected by microscopic variability, which prevents reliable analytical or predictive modeling. This limitation is particularly pronounced in random plasmonic metasurfaces, where the optical response from statistical distributions of nanostructures. The optimization process must therefore rely on an iterative experimental feedback, where each laser inscription provides information to guide the exploration of the parameter space. In this article, we introduce a grey-box Bayesian optimization framework tailored for rapid gamut expansion. In contrast to conventional black box formulations, where objective function is modeled directly, our approach uses a different structure. The goal here is not to optimize a single color, but to maximize the hypervolume of a set of printed colors in a perceptual space. To address this, we propose a grey-box, set-based Bayesian optimization framework. Instead of modeling the volume directly, a Gaussian Process is used as a surrogate model for intermediate physical quantities, namely the color coordinates. The acquisition function then analytically evaluates the expected contribution of a candidate point in terms of its ability to expand the volume of the current color gamut. This formulation simplifies the optimization problem and enables a more efficient exploration of the color space compared to genetic algorithms.

2 Material and methods

2.1 Material

The first type of sample is a thin mesoporous film of titanium dioxide elaborated by sol-gel process containing silver nanoparticles. The fabrication process and the laser-induced optical properties have been described in several articles or our group so far [38–40]. This type of sample is semi-transparent and will be observed in two observation modes transmission and backside reflection. The second type of sample is a metal-dielectric-metal-dielectric thin film stack elaborated by magnetron sputtering. For more detailed descriptions, the reader can refer to the reference [41]. This sample does not transmit light and will be observed in the following two observation modes in front-side reflection: diffuse and specular. These two types of samples contain silver nanoparticles in thin films, whose shape and spatial arrangement can be modified by laser treatment. The observed colors arise from the coupled interaction between localized plasmonic resonances of the silver nanoparticles, strongly influenced by their morphology, spatial distribution and near- or far field coupling, and the interference and propagation effects in thin films, which vary under different viewing conditions [17, 42]. The material transformation at the nanometer scale and the color variations induced by changing the laser parameters are not predictable using current physical models, as these variations do not follow a simple or linear trend in any of the observation modes. Deep neural networks have successfully modeled the complex relationships between the laser parameters and the colors [26], however their reliance on large, sample-specific training datasets limits their rapid deployment across new materials.

2.2 Laser processing

Laser inscriptions were carried out with a nanosecond pulsed laser coupled to a scanner head. This setup, described in Appendix A, enabled the printing of micrometric square pixels under well-defined processing parameter sets. In all the experiments shown here, the only laser parameters used to explore the color space were only the laser power, the laser repetition rate, and the beam scanning speed. The laser fluence ranges between 40 mJ/cm² and 457 mJ/cm², the repetition rate ranges from 10 kHz to 600 kHz, the scanning speed from 50 mm/s to 2000 mm/s. Other parameters, such as the laser polarization, distance between hatched lines in a pixel, and defocus, could also generate different colors, but they were not varied during this study. The parameters are printed on a grid-like matrix where a given laser processing parameter set is selected for each square. The squares that are produced to record the laser-induced colors are made of 5 × 5 identical printed pixels, as exemplified in Figure 1.

Figure 1 Example of inscription matrix on a sol-gel sample.

2.3 Color measurement

Color images of the laser-processed samples were recorded using a standard high-resolution commercial RGB camera integrated into the experimental setup for in line acquisition. Different positions of the white light source and the camera relative to the sample enable the recording of images in different observation modes including transmission, front-side diffuse and specular reflection, as described in Appendix B. Since the optimization process operates in a closed loop where both the color measurements and the final image generation are performed using the same camera and illumination conditions, additional calibration of the camera or lighting was deemed unnecessary. Any biases would be consistent among all iterations and even for the final image printing. The color measurement is derived by calculating the mean value of the pixel color within an area centered on the square center that encompasses 60% of the square size. To ensure repeatability, a filtering algorithm is implemented to discard squares whose color appears inhomogeneous (inscription defects such as scratches or color gradients for the same parameter indicating instabilities). This is achieved by first calculating the standard deviation value of the RGB channels within the square. If the calculated value exceeds a predefined threshold, established through empirical means, the square is filtered accordingly.

2.4 Optimization methods

Given the unpredictable nature of the color in relation to the laser parameters, the only option is to measure them after inscription. This makes an iterative approach suitable for optimizing the color gamut on any type of sample. This approach is well-suited to the nature of the experiment, as the measurements can be taken in line on the setup by adding an RGB camera. Thus, the optimization step is integrated between the first inscription and parameter selection to print the image. After every new inscription on the sample, an algorithm suggests new parameters to inscribe, aiming at improving the gamut volume (Fig. 2). This continues until an exit criterion is met, such as a limit on the number of inscriptions (iterations), achieving a certain gamut value, or a limited increase in gamut.

Figure 2 Algorithmic flowchart of the iterative laser parameter optimization and image printing framework.

Both the GA and BO approaches share a common iterative experimental loop (as illustrated in Fig. 2): (1) Inscription of an initial random set of parameters, selected using Latin hypercube sampling (LHS), (2) measurement and filtering of the resulting colors, (3) algorithmic suggestion of new parameters based on previous evaluations, and (4) re-inscription. The process is repeated until the gamut volume improvement drops below 10% over three consecutive iterations. The fundamental difference between the two methods lies in how the new parameters are suggested, which is detailed in the following subsections.

Using a GA constitutes a well-established approach for this complex problem. It is a robust method for global exploration of the parameter space without assumptions regarding the parameter-to-color function. However, GA does not fully exploit previously measured color information. Therefore, a BO that is more deterministic with the use of a surrogate model might increase the convergence speed.

2.4.1 Genetic algorithm

This approach draws extensively from Cucerca’s work [30], adhering to the standard guidelines set forth in the GA framework. We included the hue diversity criterion as a metric for evaluation, which is a key takeaway for the algorithm to consider even the hues where the colors are the least saturated. The main difference was to scrap the DSD (Design Space Diversity) and PSD (Performance Space Diversity) metrics, as they were found to be ill-suited for the materials assessed in this study. Instead, an emphasis was placed on saturation and hue diversity. Another key difference is the replacement of the repeatability criterion by a more straightforward filtering process. Figure 3 schematizes the process:

Figure 3 GA steps.

A detailed description of the metrics, the multi-objective NDSA procedure, and parameterization choices are provided in Appendix D.

2.4.2 Bayesian optimization

The second approach employed Bayesian Optimization (BO) to guide the search for optimal parameters that would maximize the gamut. The proposed BO framework is:

• Set-Based: the parameter optimization concerns the entire gamut range; each color is evaluated in comparison with the rest of the palette, i.e. the objective function leverages all preceding predictions to make each new evaluation.

• Grey-Box(8): the value of an intermediate quantity (the color) is leveraged here, contrary to typical Bayesian optimization functions that aim at predicting the objective function directly.

Thanks to its probabilistic nature, the BO weights both exploration and exploitation by targeting areas that have already been identified as performing (exploitation) or areas with high uncertainty (exploration). This approach prevents the algorithm from becoming trapped in a certain region, which is a common issue with GA. It is important to note the structural difference in the problem formulation between the two algorithms. The objective being to maximize $V\left(\mathcal{P}\right)$ the volume of (the convex hull of the colors produced by) a set of parameters, standard GA and BO need to be adapted. Indeed, if one considers the complete set $\mathcal{P}$ as the input space for GA or BO, then it means many laser inscriptions are to be done at each iteration.

For GA, it relies on individuals P_i to perform mutations and crossovers, and relies on individual fitness for selection. Because the volume is a property of a set, not an individual, the challenge is to design an individual fitness function that eventually correlates with maximizing the collective volume. For this purpose we followed Cucerca et al. [30] and converted the problem to a multi-objective one. The fitness of an individual is evaluated with respect to the current population based on its contribution to hue diversity (in the CIE LCh color coordinates) and saturation. This multi-objective framework is necessary to maintain a diverse Pareto front; a standard single-objective GA would quickly converge to a single dense region of the color space, failing to span a wide gamut.

For BO, as the approach is more flexible, we propose a new grey-box approach (i.e. where part of the objective function is known). The surrogate model (GP) models the intermediate physical property by predicting the color of individual points, but the objective function is evaluated at the set level. The acquisition function analytically estimates the expected gamut improvement by averaging over all possible outcomes, defined as the added volume a predicted color brings to the existing set (Eq. (5)). This objective function has the added benefit of direct applicability in scenarios involving multimodes colors. It has been shown already that predicting the resulting spectra obtained using any processing parameter set for those metasurface was doable [26] using a deep neural network. However our approach aims to reach the optimal gamut configuration as rapidly as possible while relying on a limited number of experimental data points, making Bayesian optimization more adapted than deep learning model that requires a large amount of data and will not adapt easily to a new initial material.

A gamut hypervolume can then be computed even for higher dimension (6D for 2 modes, 9D for 3 modes, etc.), which makes it easy to introduce the concept of image multiplexing.

The first two steps are the same as with GA. The last step looks similar as well, but it is here possible to only reinscribe one color in a new iteration, unlike the GA that required to inscribe a batch of colors each time.

The objective of the proposed framework is to maximize the diversity of structural colors generated on the material surface by intelligently exploring the available laser settings. The input space $\mathcal{X}\subset {ℝ}^3$ is strictly bound by the physical limits of the laser system: power (40–457 mJ/cm²), scanning speed (50–2000 mm/s), and repetition rate (10–600 kHz).

∀i,$p_i\in \mathcal{X}$ is a parameter vector with a corresponding output that is the 3-dimensional color coordinate vector $c_i=\mathit{col}\left(p_i\right)\in {ℝ}^3$ in the CIELAB color space.

Hypervolume V(P) is computed as the volume of the convex hull of the colors obtained from the finite vector parameter set P = {p₁, p₂,…} in the L*a*b* space. We choose this space for defining both the state space and the objective function due to its perceptual uniformity.

The hypervolume associated to a set of laser parameters can be expressed as:$$V\left(P\right)≔V\left(\text{ConvexHull}\left({\left\{\mathit{col}\left(p_i\right)\right\}}_{p_i\in P}\right)\right).$$(1)

Formally, for a fixed size n, the problem is to find the optimal set of laser parameters $P_n^{\ast }$ that maximizes the gamut hypervolume V, which is defined as the volume of the convex hull formed by the colors induced by P*, namely:$$P_n^{\ast }=\text{argma}{\mathrm{x}}_{P\in {\mathcal{X}}^n}V\left(P\right).$$(2)

However, because the true physical mapping col(p) is a highly complex unknown function, finding the optimal set P* is unrealizable. Instead, the optimization must be performed sequentially in order to acquire information about col, the function linking the laser parameters and the produced color on the current sample. In its sequential form, where the set of parameters is constructed one point after the other and the true color is obtained, the problem can be written, at each iteration n as:$$P_{n+1}^{\ast }=\left\{p_{n+1}^{\ast }\right\}\cup P_n^{\ast }$$(3)

However, as pointed out, the true physical color col(p) is a black-box outcome before the actual physical experiment. Consequently, the physical volume equation cannot be solved directly. The optimal sequential solution $P_n^{\ast }$ actually ignore the colors measured (as it supposed the col function is known), but any algorithm will use the measured colors ${\mathcal{C}}_n$ as a (very partial in the first iterations) information about the col function and thus will fail to achieve the optimal. The surrogate model consists in a Gaussian Process (GP), that is fitted using the (n − 1) data points and is fully defined by a vector-valued mean function and a multi-output covariance kernel. The GP implementation is further detailed in Appendix F.

Using this GP we substitute the physical col(p) by a predictive probability distribution $\mathbb{P}\left(c\right|p)$, the practical optimization step solved by the algorithm is therefore:$$p_n^{\ast }=\text{argma}{\mathrm{x}}_{p\in \mathcal{X}}\mathit{EV}\left(p\right).$$(4)

Let f_gam be a function that associates to a point c of the color space the corresponding contribution to the added volume in the gamut gam:$$f_{\mathit{gam}}\left(c\right)=\mathit{Vol}\left(\text{ConvexHull}\left(\mathit{gam}\cup \left\{c\right\}\right)\right)-\mathit{Vol}\left(\text{ConvexHull}\left(\mathit{gam}\right)\right).$$(5)

The EV is meant to operate by computing the volume that every parameter is expected to bring by summing the added volume of any given possible color the parameter set could produce times the density of probability of said color occurring:$$\mathit{EV}\left(p\right)={\int }_{{ℝ}^3}‍\mathbb{P}\left(c\right|p)V\left(\left\{c\right\}\cup {\mathcal{C}}_n\right)\mathit{dc}.$$(6)

While equation (6) defines the theoretical integral over the entirety of ${ℝ}^3$, computing this unbounded integral numerically is intractable. Because the probability density function $\mathbb{P}\left(c\right|p_n)$ of a Gaussian Process decays exponentially away from its predicted mean, the numerical integration is practically restricted to a localized bounding box in the color space where the probability mass is computationally significant.

Within this truncated domain, the integral is approximated numerically using a discrete step size of (ΔL^*, Δa*, Δb*) = (0.1, 0.1, 0.1). This expected improvement evaluation is systematically performed across a dense grid of candidate parameters with a resolution of (ΔPower, ΔSpeed, ΔFreq) = (1%, 10 mm/s, 5 kHz). The candidate parameter p_n that maximizes this numerical EV is ultimately selected for the next physical inscription.

The EV function is designed to consider the weight of any possible color within a given parameter range of colors available that lies outside of the already measured gamut.

Figure 4 shows how two different parameters will be evaluated by the algorithm: the parameter $P_1$ will most likely produce a color that is already printable as it is most likely inside the convex hull of the already measured colors gamut, $P_2$ and $P_3$ on the other hand have more chances to produce a color outside which would increase the volume.

Figure 4 Example of parameter selection based on predicted color distributions in the CIELAB color space. The a* and b* axes represent chromaticity coordinates, while the shaded blue polygon denotes the color gamut already obtained without adding new parameters. The colored ellipses (red for $P_1$, green for $P_2$, and purple for $P_3$) indicate the 95% confidence regions of the predicted color positions, showing where the resulting color is most likely to lie for each parameter. Dots within the ellipses represent individual simulated color coordinates, and the black-edged circles correspond to the mean predicted color for each parameter. The outer black line represents the current achievable gamut, while the green and purple outlines illustrate possible extensions of the gamut if the actual color corresponds to predicted points outside the existing region. This example suggests that $P_2$ is the most promising choice (most likely to expand the gamut), followed by $P_3$ (large uncertainty) and $P_1$ (unlikely to expand the gamut).

The maximum value for the expected improvement across the whole parameters set space is selected to be inscribed and the process carries on towards the next iteration.

2.5 Gamut mapping

When reproducing images containing colors that fall outside the achievable color gamut of a given printing process, it is necessary to adapt the image colors to the available palette. This process, known as gamut mapping, follows the approach introduced by Chosson and Hersch [28]. The procedure starts with the definition of a set of primary colors, including black, white, and chromatic primaries, which define the accessible color palette. To ensure consistent color appearance under varying viewing conditions, a chromatic adaptation transform is applied [43]. In practice, the tristimulus values of the palette are converted from their original white point to the D65 standard illuminant using the CAT02 chromatic adaptation model [44], as implemented in the colour-science library. This transformation maps the palette white (scaled by an adaptation factor) to D65 and adjusts all other colors accordingly. A lightness soft-compression is then performed to match the luminance range of the image to that of the available palette [45]. In parallel, the colors of the palette are transformed such that the gray axis defined by its black and white points aligns with the gray axis of the CIE 1976 L*a*b* color space. This alignment ensures a consistent treatment of achromatic components during the mapping process. Finally, a chroma soft-compression is applied to the image colors to constrain them within the available gamut while preserving, as much as possible, their perceived saturation. These successive transformations enable the selection of colors that provide the most faithful visual reproduction of the original image within the limitations of the accessible color space.

However, gamut mapping does not overcome the intrinsic limitations of the achievable color space. Improving the native color gamut therefore becomes a critical objective, which requires identifying laser processing conditions that maximize the accessible color space.

3 Results

In this section, a comparison on how much the GA, BO as well as random exploration strategies increase the gamut size and the printing quality. The comparison were done using “replay experiments” analysis, “live analysis” and by comparing the final images with both methods.

3.1 Replay experiments

A postmortem analysis was carried out to evaluate the speed with which each technique reached the total gamut. To do so, as many parameter combinations as possible were printed on the two types of samples and the colors were measured in different observation modes (described earlier). To ensure a fair comparison between the explorations, the GA, BO, and random explorations start with the same set of parameters, determined by a set seed. The use of LHS ensures a well-distributed initial sampling of the parameter space, which explains why the starting gamut volumes are relatively consistent and not drastically different across various seeds. The exploration was simulated according to each method. Random selection sampled new candidates uniformly at random. GA suggested new parameters and the closest non inscribed parameter (using Euclidian norm over power, speed, repetition rate space) in the database was selected. For the BO, as the proposition is done through an acquisition function, only the value for the parameters present in the dataset was evaluated. The performance was assessed by evaluating the rate of variation in the gamut volume (Fig. 5), and the number of points expected to reach a certain proportion of the total gamut. The left plot corresponds to a sample made by magnetron sputtering in scattering reflection mode, and the right one to the mesoporous film sample in backside reflection mode. Similar results are observed in other modes for the same samples as shown in Appendix C. BO method did outperform GA, achieving a faster convergence toward the maximum attainable gamut volume with fewer of the inscribed parameters were required to reach the maximum gamut.

Figure 5 Gamut size evolution for BO, GA and random, (a) for magnetron sputtering sample in scattering reflection mode, (b) for mesoporous film in backside reflection mode.

3.2 Real time optimization

A “live” experiment was also performed to complement the postmortem analysis. Contrary to the post-mortem method, this allows GA and BO to actually improve the gamut of inscribed colors by suggesting never-inscribed-before parameters that will be marked experimentally on the sample. The volume improvement on the gamut was tracked for both methods: the same initial random selection of 40 laser parameters combinations was used for both methods. This experiment was done on several variations of the reflective sample elaborated by magnetron sputtering, all elaborated with different thicknesses of the intermediate dielectric layer. Measurements were done here in diffuse reflection observation mode.

The optimization results are compiled in Table 1. They show that, irrespective of the sample, the BO optimization leads to the largest gamut volume after the inscriptions of 400 colors. It is also the most effective optimization method after reaching the exit criterion, with the exception of one sample (120 nm) where it is marginally surpassed by GA. The study also demonstrates that the optimal thickness to achieve the widest color gamut is 140 nm. However, if we were to rely solely on the GA optimization, the conclusion would have identified the 120 nm thickness as the most effective one to achieve a broader color palette. For the sake of illustration of the gain brought by the optimization process, Figure 6 compares images simulated by using either the color gamut obtained by considering randomly chosen laser parameter sets (pre-optimization), or the gamut optimized with the GA or BO. The corresponding gamuts can be found in Appendix E. This figure clearly illustrates the gain brought by the optimization processes. The first row corresponds to the gamut in transmission of a semi-transparent sample elaborated by sol-gel process. The second row corresponds to a gamut in specular reflection on a reflective sample elaborated by magnetron sputtering with a thickness of 140 nm. The third row corresponds to the same sample as in the second row, but for a gamut observed in diffuse reflection. In order to assess the precision of the reproductions, a metric designed for highly distorted and restricted gamut was employed, as defined in [46]. The score of this metric is indicated on the bottom right side of each simulated image. The better the reproduction, the lower the metric score. The proposed image-quality score is not calibrated as an identity metric. It is a linear model fitted to predict relative rankings from metrics such as δL* range or a*b* hull area, after ANOVA/AIC selection. Since each basic computable metric is normalized using fixed bounds and several terms depend jointly on the image content and the gamut mapping (chromatic adaptation, gray-axis alignment, lightness/chroma soft-compression), the model is not constrained to return exactly 1 for the original image and small deviations around unity are expected. These images indicate that both optimization methods are useful in enhancing the perceptual fidelity to the original image by simply improving the gamut. BO optimization appears even better as expected from the larger increase in the gamut. The fourth row of Figure 6 shows the images recorded in diffuse reflection of the printed samples before and after GA and BO optimizations of the gamut. The image quality metric value is provided but should be taken carefully as the metric was designed with simulated images only. Except some defects on the samples, the experimental results confirm the predictions.

Figure 6 Comparison of optimization outcome with 100 points for reproducing a portrait, the picture of a dog and a palm tree doing simulation with the measured colors. The framed palm tree images were laser inscribed. The number on the bottom right image is the image quality metric (IQM) score. A lower IQM value corresponds to a better reproduction. The IQM evaluation is performed all images but is more relevant for the simulated images

3.3 Summary of optimization performance

Based on the outcomes of both the post-mortem analysis and the live experiments, the BO approach proposed in the paper exhibits several advantages over the GA approach:

• Simplified problem formulation: The BO approach operates with a single-objective formulation, whereas the GA approach, which requires reformulating the problem as a multi-objective optimization. The objective is defined as the improvement of the gamut volume alone, without complexifying it with criteria such as hue diversity, because BO does not persist in trying parameters in the same zone repeatedly. BO reduces the complexity of the optimization and the associated computational overhead.

• Faster expansion of the color gamut: By focusing on enhancing the gamut volume, the method rapidly increases the number of colors that can be achieved with each inscription. This efficiency stems from the Bayesian optimization framework’s ability to intelligently select the subsequent candidate points based on predicted improvement, as opposed to attempting to balance multiple competing objectives, which often slows convergence.

• Better gamut coverage at fixed cost: Given a “printing” budget (a fixed number of inscriptions), the method can produce a larger color gamut than GA, which struggles to efficiently explore the parameters and often results in redundant sampling in dense regions. The single-objective Bayesian optimization model naturally avoids this problem by targeting regions of the color space that maximize volume gain.

• Effective without explicit diversity criteria: In GA, a hue diversity criterion is required to prevent the algorithm from repeatedly sampling in the same color regions. The Bayesian optimization approach balances exploration and exploitation to ensure coverage of new regions without the need for explicit diversity enforcement. This makes the method both faster and more effective, achieving broader gamut expansion with a reduced number of evaluations.

Ultimately, the proposed Bayesian optimization framework proves to be highly generalizable and can be applied to any material system for which color can be modified by laser processing. Its versatility has been demonstrated here on two distinct types of samples, including plasmonic random metasurfaces. In all cases, the framework successfully identified the laser processing parameters required to maximize the achievable color gamut in a given observation mode.

4 Conclusion

In this study, we introduced a Bayesian Optimization (BO) model for laser-induced color printing and compared it with the Genetic Algorithm (GA) method, which had not yet been applied to this type of sample. Both approaches resulted in a substantial expansion of the accessible color gamut. BO demonstrated faster convergence and a more systematic exploitation of the available data, leading to a larger gamut achieved. The perceptual comparison of generated images further confirmed the benefit of the optimization.

The proposed BO model, relying on a single and easily defined acquisition function, facilitates scalability to higher-dimensional problems involving multimodal colors. This framework can be directly extended to the optimization of a hyper-gamut spanning multiple viewing conditions, where the same principle applies. The algorithm not only enhances image quality but also favors a better control over the inscription process, enabling the simultaneous optimization of multiple constrained images on a single sample by adjusting their respective hypercolors.

As a next step, this approach could be extended toward image multiplexing applications, thereby facilitating simultaneous optimization across multiple output modes. Although this was not explored in the present work due to the lack of a setup allowing simultaneous color measurement in both modes, it remains a promising direction for future research. Another possible next step would be to return to a multi-objective Bayesian Optimization (MOBO) framework, introducing reproducibility as an additional optimization criterion to complement the current single-objective formulation.

Abbreviations

GA: Genetic Algorithm

BO: Bayesian Optimization

DSD: Design Space Diversity

PSD: Performance Space Diversity

NDSA: Non-Dominated Sorting Algorithm

IQM: Image Quality Metric

Funding

This work is funded by the ANR project SLICID (ANR-23-CE39-0006).

Conflicts of interest

The authors declare no conflicts of interest.

Data availability statement

The replay gamut exploration data are available at https://github.com/robinmermillod3-fileshare/Bayesian-optimization-of-laser-process-for-color-gamut-data.gitgithub in the form of image files containing the evolution of the gamut volume with different methods.

Author contribution statement

Conceptualization: Rémi Emonet and Nathalie Destouches;

Methodology: Robin Mermillod-Blondin and Rémi Emonet;

Software and Model Development: Robin Mermillod-Blondin;

Experimental Setup and Data Acquisition: Robin Mermillod-Blondin, Nicolas Dalloz and Nathalie Destouches;

Validation and Analysis: Rémi Emonet;

Visualization: Robin Mermillod-Blondin;

Writing – Original Draft Preparation: Robin Mermillod-Blondin;

Writing – Review and Editing: Jessica Pellegrino, Nicolas Dalloz, Nathalie Destouches, Rémi Emonet;

Supervision: Jessica Pellegrino, Rémi Emonet and Nathalie Destouches;

Acknowledgments

The authors thank Dr. Francis Vocanson (Laboratoire Hubert Curien) and Nicolas Jacquot (HEF) for kindly providing samples used in this study.

References

Appendix A: Laser processing setup

The primary source is a nanosecond IPG laser operating at a wavelength of 532 nm, with a pulse duration of 1.3 ns, a Gaussian intensity profile and a circular beam cross-section. It is focused on the sample surface by a 10 cm focal length f-theta lens that is mounted on a scanner head enabling the precise control of the spot location in the focal plane (plane of the sample surface) over a maximum field of 6 × 6 cm². The intensity profile is transformed into a top-hat profile and the circular shape is converted into a square shape by the CANUNDA beam-shaping module from Cailabs, which is located before the scanner head. The focused square spot has a nearly flat intensity profile over a width of 13.5 μm (full width at 90% of the maximum) Figure A1 and a Rayleigh length of 150 μm. The beam is linearly polarized. The laser power and repetition rate, the beam path and the scanning speed are computer controlled. Each printed pixel is a square of 50 μm width and is obtained by drawing hatched lines oriented at 45° from the x and y directions on the sample plane and separated by a distance of 10 μm.

Figure A1 Square top hat beam shape.

Appendix B: In line measurement setup.

The laser-induced colors were measured in line using a Nikon Z7 II camera mounted above the scanner head. The camera was equipped with a macro lens providing a working distance of approximately 20 cm and a field of view covering the entire sample area (≈ 5 cm). This configuration allows the acquisition of high-enough resolution images where each inscribed square (200 μm × 200 μm) is captured with a sufficient number of pixels to be able to perform an averaging of the color. A schematic of the color measurement setup is shown in Figure C1.

Figure B1 Measurement setup using a basic camera equipped with macro objective and two light sources, one for specular reflection (dark blue ray), one for the scattered reflection (red ray), the last one in the other side is for transmission in the case of the semi-transparent mesoporous film

Figure B2 Focus bracketing performed on 4 images taken with the focus being on different locations of the sample, the red rectangle indicates where the image is.

To ensure uniform sharpness across the entire sample surface, focus bracketing was employed during image acquisition as illustrated in Figure B2. It means taking several images focused at different depths and combining them to get one image that’s sharp everywhere. The use of the same camera and illumination setup for both measurement and final image production ensures consistency across all iterations of the optimization process.

Appendix C: Replay experiment with other sample and mode

Gamut evolution with different optimization method for the modes-sample combinations not shown in Figure 5.

Figure C1 Gamut size evolution for BO, GA and random, (a) for magnetron sputtering sample in specular reflection mode, (b) for mesoporous film in frontside reflection mode, (c) for the same sample in transmission mode.

Appendix D: Genetic Algorithm full description

After the laser processing and measurement, each candidate solution is then evaluated compared to others. Since optimization requires yielding more saturated colors, selecting the most chromatic colors is important; however, since the gamuts are asymmetrical it would favor some gamut regions far too much. Cucerca [11] addressed this issue by defining another metric called hue diversity to ensure that even the gamut region with low saturation colors get some representatives. This metric was adapted to the present study where the gamut asymmetry is more dire. Other metrics introduced like Design Space Diversity (DSD), Performance Space Diversity were discarded as the high sensitivity of the material towards the laser parameters rendered them uninformative for the exploration. The DSD metric, which assesses the distance of each input to its closest neighbors, does not provide meaningful insight for Pareto front identification when the initial input set lies on a predefined grid. Likewise, the PSD metric offers limited relevance in our case, as the objective is to highlight configurations where the resulting color deviates significantly from the others, rather than simply maximizing dispersion in performance space. Same for the repeatability metrics which was overtaken by the parameter filtering during the measurement. Naturally, the criteria are not the same when tackling with best black (resp. white) research within the sample, what is more favored here is the closeness to the neutral axis and to minimize (resp. maximize) the lightness.

Figure D1 Non Dominated Sorting Algorithm example. (a) Determining the first front. (b) Determining the second front. (c) All points sorted in their front.

To compute chromaticity, the RGB channels value of the squares were converted to the CIE Lab color space, known for being perceptually uniform. Then it was converted to cylindrical coordinates LCh (Lightness, Chroma, Hue) where the chroma value is directly accessible.

The second criterion for color exploration is to determine every color’s contribution to hue diversity. It is intended to encourage the correct classification of colors that extend the gamut the most in hue ranges where the gamut is restricted, so that these colors are more likely to be selected subsequently and the gamut can be explored around them. In the case of the sol-gel and SLICID sample, the initial gamuts are often eccentric, which tends to penalize certain colors, which could nevertheless be interesting, on this criterion of hue diversity. It was therefore decided not to apply it to the colors as measured, but to the colors obtained after the gamut grey axis alignment (3), which enables them to be refocused around the grey axis. To compute this metric, the a*b* plane of the CIE Lab space is segmented into a random number of segments ranging from 4 to 72, with a random angular offset on the starting sector. Then in every sector the colors are ranked based on their chroma. For each different hue wheel, every color can have a different ranking. Every candidate is then characterized by its ranking frequency, appearing n times in rank 1, m times in rank 2, … The global ranking of the solutions is done using lexicographic ordering meaning number of rank 1 is given strict priority, then rank 2, …

For the search of black and white, the two criteria that were discussed for color exploration are discarded. For achromatic exploration, the chroma must be minimized as the colors closer to the grey axis should be prioritized. For the black colors, the lightness criteria in CIE Lab has to be minimized, for the white colors it should be maximized.

Following the framework to tackle multi-objective GA, the candidate solutions were sorted using the Non-Dominated Sorting Algorithm (NDSA), an algorithm that sorts entries according to the number of entries that are better than them in both criteria at the same time. The NDSA works by first identifying the pareto front, i.e. the set of solutions such that no other solutions constitute an improvement in both criteria, the pareto front is labeled as front 1 Figure D1a. Then every point that were not part of the pareto front go through the NDSA again to find the front 2 Figure D1b. The process continues until every point have been assigned to a front number Figure D1c. Once all points have been sorted, they are randomly selected using a Roulette-wheel selection, with selection probabilities depending on the front number, points belonging to lower order front having higher chance of being selected. The number of colors selected is fixed at about 20% of the initial population size.

After the selection, offspring parameters were generated using uniform crossover (crossover rate of 0.2) and mutation (mutation rate of 0.2 and parameter modification rate of 0.1), applied independently. Crossover promotes recombination of promising solutions and mutation maintains diversity. In the following paragraph, we explore the impact of those parameters. The stop criteria to interrupt the iteration was when the gamut has note increased significantly in the past iterations.

The adjustment of the mutation and crossover parameters was performed by a postmortem experiment. Using convergence rate for different parameters reported in Table D1, defined by the proportion of the number of points required to reach 95% of the gamut volume, done on a 14000 points exploration made on SLICID.

Appendix E: Gamuts obtained

The figures here show gamuts obtained and used to obtain the images on the first row of Figure 6 in CIE a*b* and L*C* planes:

For the first row (sol-gel in transmission mode) Figure E1:

For the second row (slicid in specular reflection mode) Figure E2:

For the third row (slicid in scattering reflection mode) Figure E3:

Figure E1 Gamuts obtained on sol-gel in transmission mode.

Figure E2 Gamuts obtained on slicid in specular reflection mode.

Figure E3 Gamuts obtained on slicid in scattering reflection mode.

Appendix F: Implementation details of the Gaussian process Surrogate model

Rather than directly modeling the macroscopic set volume $V\left(P\right)$, the surrogate model approximates the underlying physical black-box function $f:\mathcal{X}\to {ℝ}^3$. Because the output coordinates ($L^{\ast }$, $a^{\ast }$, and $b^{\ast }$) are inherently correlated in the CIELAB color space, we employ a Multi-Output Gaussian Process based on the Linear Model of Coregionalization (LMC). Unlike independent, single-task GPs, the LMC captures cross-channel dependencies.

Figure F1 Illustration of 4.5 δE76 color difference.

The multi-output surrogate model is defined by a vector-valued mean function $\mu \left(p\right)$ and a multi-output covariance kernel $K(p,p^{\prime })\in {ℝ}^{3\times 3}$:$$f\left(p\right)\sim \mathit{GP}(\mu \left(p\right),K(p,p^{\prime }))$$

where the coregionalized kernel is expressed as the Kronecker product of a coregionalization matrix B (which models the cross-channel correlations) and a spatial kernel $k(p,p^{\prime })$ (which models the parameter dependencies):$$K(p,p^{\prime })=B\otimes k(p,p^{\prime })$$

This architecture was implemented using the LMCVariationalStrategy provided by the GPyTorch library (https://docs.gpytorch.ai/en/v1.12/variational.html). Under the LMC framework, each output task $f_i\left(p\right)$ is expressed as a linear combination of Q independent latent Gaussian processes:$$f_i\left(p\right)=\sum _{q=1}^Q‍a_i^{\left(q\right)}{g}^{\left(q\right)}\left(p\right)$$

where $i\in \left\{123\right\}$ represents the specific output task (corresponding to the $L^{\ast }$, $a^{\ast }$, and $b^{\ast }$ color channels, respectively), $g^{\left(q\right)}\left(p\right)$ is the q-th independent latent GP, and $a_i^{\left(q\right)}$ are the scalar weighting coefficients. In our formulation, the model was configured using $Q=3$ independent latent variables.

These coefficients $a_i^{\left(q\right)}$ constitute the elements of a task projection matrix A. The coregionalization matrix B is $A{A}^T$.

Finally, to assess the model performance under different structural assumptions, several combinations of mean and covariance functions were tested. Each model was trained on the same subset of data, and the kernel hyperparameters (e.g., length-scale, output scale, noise level) were optimized by minimizing the mean CIE $\mathrm{\Delta }{E}_{76}$, defined as $\mathrm{\Delta }{E}_{76}=\sqrt{{\left(\mathrm{\Delta }{L}^{\ast }\right)}^2+{\left(\mathrm{\Delta }{a}^{\ast }\right)}^2+{\left(\mathrm{\Delta }{b}^{\ast }\right)}^2}$, between the prediction and the ground truth. As shown in Table F1, the best model achieved a mean $\mathrm{\Delta }{E}_{76}$ of 4.5, which is exemplified in Figure F1. While this value exceeds the typical Just Noticeable Difference (JND ≈ 2.3) making the error visually perceptible, this level of precision is entirely acceptable for the purpose of the surrogate model. The GP is not used to precisely predict final colors for high-fidelity image printing, but rather to identify high-probability regions of gamut expansion to guide the physical exploration phase. The results shown in Table F1 correspond to predictions obtained from models trained on 20 measured colors (similar to the initial points) and used to predict, one at a time, other 10000 measurement points. The discrepancies between predictions and true values were averaged to quantify the model error. This process was repeated 10 times with different fixed random seeds to ensure the robustness and the relevance of the comparison.

Based on minimizing the mean color prediction error (CIE ΔE₇₆) on the initial measurement points, a constant mean function and a Radial Basis Function (RBF) spatial kernel k(p, p′) were selected as they yielded the highest predictive accuracy for this physical system.

All Tables

All Figures

	Figure 1 Example of inscription matrix on a sol-gel sample.
In the text

	Figure 2 Algorithmic flowchart of the iterative laser parameter optimization and image printing framework.
In the text

	Figure 3 GA steps.
In the text

Figure 4 Example of parameter selection based on predicted color distributions in the CIELAB color space. The a* and b* axes represent chromaticity coordinates, while the shaded blue polygon denotes the color gamut already obtained without adding new parameters. The colored ellipses (red for $P_1$, green for $P_2$, and purple for $P_3$) indicate the 95% confidence regions of the predicted color positions, showing where the resulting color is most likely to lie for each parameter. Dots within the ellipses represent individual simulated color coordinates, and the black-edged circles correspond to the mean predicted color for each parameter. The outer black line represents the current achievable gamut, while the green and purple outlines illustrate possible extensions of the gamut if the actual color corresponds to predicted points outside the existing region. This example suggests that $P_2$ is the most promising choice (most likely to expand the gamut), followed by $P_3$ (large uncertainty) and $P_1$ (unlikely to expand the gamut).

In the text

	Figure 5 Gamut size evolution for BO, GA and random, (a) for magnetron sputtering sample in scattering reflection mode, (b) for mesoporous film in backside reflection mode.
In the text

Figure 6 Comparison of optimization outcome with 100 points for reproducing a portrait, the picture of a dog and a palm tree doing simulation with the measured colors. The framed palm tree images were laser inscribed. The number on the bottom right image is the image quality metric (IQM) score. A lower IQM value corresponds to a better reproduction. The IQM evaluation is performed all images but is more relevant for the simulated images

In the text

	Figure A1 Square top hat beam shape.
In the text

	Figure B1 Measurement setup using a basic camera equipped with macro objective and two light sources, one for specular reflection (dark blue ray), one for the scattered reflection (red ray), the last one in the other side is for transmission in the case of the semi-transparent mesoporous film
In the text

	Figure B2 Focus bracketing performed on 4 images taken with the focus being on different locations of the sample, the red rectangle indicates where the image is.
In the text

	Figure C1 Gamut size evolution for BO, GA and random, (a) for magnetron sputtering sample in specular reflection mode, (b) for mesoporous film in frontside reflection mode, (c) for the same sample in transmission mode.
In the text

	Figure D1 Non Dominated Sorting Algorithm example. (a) Determining the first front. (b) Determining the second front. (c) All points sorted in their front.
In the text

	Figure E1 Gamuts obtained on sol-gel in transmission mode.
In the text

	Figure E2 Gamuts obtained on slicid in specular reflection mode.
In the text

	Figure E3 Gamuts obtained on slicid in scattering reflection mode.
In the text

	Figure F1 Illustration of 4.5 δE76 color difference.
In the text