Band selection for mid-wave infrared thermography through extended radiometric calibration and heuristic-based optimization

Non-contact radiation-based multispectral thermometry methods are inherently underdetermined due to the inability to decouple the effect of emissivity and temperature on the spectral exitance of the target. This work introduces a calibration method that incorporates the assumption of local emissivity linearity, leading to an extended radiometric system model with better regularization of the inverse problem. A heuristic-based multiobjective optimization scheme for designing optimal filter configurations per target scenario is developed, which considers the spectral shape of the band-pass filters. Through a per-scenario suggestion of filter combinations, the suggested procedure can use prior information about the spectral emissivity and the temperature range of interest. The optimization procedure is demonstrated for a four-channel thermographic system and numerically evaluated over a broad range of common aerospace materials. An optimally designed system, operated according to the suggested guidelines, is predicted to result in temperature recovery errors within 10 K over a range of 423–873 K.

Spectrally-resolved readings are a well-known data source for determining the temperature and emissivity of various types of targets [12][13][14]. Multispectral radiation thermometry (MRT) is a family of data acquisition and processing techniques for obtaining temperature based on radiometric measurements performed in several spectral bands. In these techniques, independent band-intensity datasets are used alongside a spectral emissivity model to reconstruct the temperature and apparent emissivity of the target. Thermography based on MRT may allow the acquisition of 2D surface temperature distributions using very weak assumptions on the target's emissivity [15][16][17].
Multispectral thermographic cameras are usually constrained by their small number of filter slots, typically up to 8. The limited spectral information makes it harder to differentiate similar targets, diminishing the resolution and the range of the resulting temperature measurements. These effects, in turn, limit the variety of targets or measurement conditions for which the device can produce acceptable outputs. For this reason, designing a universally applicable multispectral thermometer is extremely difficult, if at all possible. However, universality can be traded for better performance for a narrower range of targets [13,[18][19][20]; an extreme example of specialization is that of pyrometers built for specific materials in known temperature ranges [21][22][23][24]. Thus, a balanced MRT system would retain some target-independence, in the sense that it will not require ad-hoc calibration per target while providing acceptable temperature accuracy for all targets for which it was designed. A known difficulty with pyrometric methods is associated with their inherent under-determinedness, where there is always at least one more unknown than the number of available equations (one temperature and one emissivity at each measurement wavelength). This difficulty is addressed either through supplemental measurements (i.e., increasing the number of equations), or via assumptions about the targets (i.e., decreasing the number of unknowns). An example of the former is a class of methods known as "active", where the target is excited in order to obtain another quantity (like reflectivity) from which the emissivity is deduced [25,26]. An example of the latter are assumptions about the spectral dependence of emissivity, including constant, stair-wise, polynomial, exponential and other models [27,28].
Despite receiving much attention in prior literature, there appears to be no consensus about a universal modeling technique that applies to a broad range of targets [29]. There is an important tradeoff when choosing emissivity models, and specifically the number of degrees of freedom (DOF) in the model. A model with many parameters constitutes a "weak" assumption on the spectral emissivity, which makes the device applicable to a broader range of targets, yet this increases the risk of overfitting and reduces the redundancy of the system, making it more error-prone. A model with few parameters, like the single-parameter a gray-band model, constitutes a "strong" assumption, in that fewer materials can be truly considered gray, leading to a potentially more substantial modeling error. Significant errors in temperature might occur when the target deviates from the chosen emissivity model [28]. Typically, emissivity models are chosen such that the resulting equation system is overdetermined. However, such over-determinedness might be all but superficial if individual measurements are correlated to a large degree, resulting in an ill-conditioned system, making parameter recovery very error-prone [28,30]. This fact sparked investigations into optimal spectral configurations, attempting to achieve sufficient measurement heterogeneity without overextending the emissivity model, especially in scenarios of limited channels [31][32][33][34][35].
The outcomes of previous optimizations were sets of wavelengths (or spectral bands) most conducive to the considered measurement scenarios. In practice, these recommendations are realized in the form of narrow band-pass filters, owing to their affordability, diversity, and ease of replacement. However, commercial filters are not perfectly monochromatic, and neglecting this aspect results in a degradation of optimality and increased error when attempting to build such systems ( [36] reported a 200-fold increase in temperature uncertainty and a 30-fold increase in the standard deviation of the temperature). Thus, balancing the involved tradeoffs warrants a broader optimization scope, capable of suggesting not only the measurement wavelengths but also the spectral shapes of the filters.
The present work suggests a procedure that significantly expands the applicability of a linear emissivity model. This is achieved through an extended characterization of the measurement system and better utilization of prior information about the expected range of targets. The proposed calibration procedure empirically captures not only the relation between temperature and detected signal but also between temperature and the drift of the characteristic wavelength. Furthermore, the problem of filter design space exploration is formulated as a multi-objective optimization. The outcomes of this work can be used as a low-cost and effort way to enhance the temperature recovery capabilities of existing thermographic systems.

Outline
The presented methodology consists of three avenues that ultimately lead to the suggestion of optimal filter combinations, as shown in Fig. 1. Section 2.1 presents a model of the measurement system and the procedure for its calibration; Section 2.2 explains how simulated measurements are generated; and Section 2.3 shows how temperature is estimated. Section 2.5.1 introduces the heuristics which are used in a multi-objective optimization scheme to create a Pareto set of filter combinations. Finally, Section 2.5.2 discusses how optimal filter combinations are selected from the Pareto set by applying the temperature recovery scheme to simulated measurements using a database of real emissivities.

Optical system calibration
Relating measured signals to physical quantities requires radiometric calibration using a well-defined radiation source. This procedure is commonly done using a blackbody radiator whose temperature relates to spectral radiance via the Planck equation: where c 1 , c 2 , T and λ are the 1 st and 2 nd radiation constants, temperature and wavelength, respectively. However, equation (1) describes radiation emitted by the targetwhich is larger than the detected radiation. The comparably smaller signal put out by the measurement system, R BB , is given by: where R(λ) is the detector's spectral responsivity, τ lens (λ) is the spectral transmittance of the lens, τ filter (λ) is the spectral transmittance of the filter, τ atm (λ) is the atmospheric transmittance along the optical path, and R eff (λ) is the effective spectral responsivity of the combination of a detector, a lens and a filter. The main downside of equation (2) is the difficulty of using it with real optical setups, where the various spectral functions are only known approximately (including inaccuracies in the spectral transmittance of all elements, the optical density of filters, and the responsivity of the detector). Thus, the evaluation of the integral is inherently error-prone.
A common approach to describe the change of the measured photoresponse, R, with respect to the known temperature, T, is using approximate interpolation equations [37][38][39]. One of the interpolation formulae that mimic equation (1) well is a three-parameter model, where a i are parameters fitted to every distinct optical configuration, which take into account the detector response, the spectral transmittance along the optical path, and the form factor in the experiment [38,40]. The benefits of using an interpolation equation are twofold: experimentally, it sums up the calibration of an optical system to 3 coefficients; and computationally, it allows to evaluate (2) on a coarse temperature grid, which enables quick and accurate interpolation at all intermediate values. This approach was successfully used in the past to study thermographic systems [41]. Unlike blackbody radiators, which have a spectrally constant emissivity close to unity, real targets almost always exhibit spectral, temporal, and temperature-related variations in emissivity. Moreover, since the emissivity can have an arbitrary spectral shape, every spectrally distinct measurement comes with its unknown emissivity, making the problem underdetermined. Furthermore, the total incoming radiance, Γ Σ , also includes a potentially significant reflected radiance term (Γ ρ ) and an atmospheric emission term (Γ a ), where Γ is the emitted radiance term, ε is the target's emissivity, B ref is the weighted sum of the reflected radiance sources (where the weights are given by the configuration factors for each surrounding source), and T amb is the ambient temperature. Then, the total photo-response with filter k, R Σ,k , can be expressed by where R eff,k is the effective responsivity with filter k, R k , R ρ,k and R a,k are the photo-response components related to the target's exitance, the target's reflectance, and the atmospheric emission with filter k, respectively. The mentioned difficulties are addressed by making several as-sumptions. While it is not assumed that emissivity is independent of temperature, it is considered that image acquisition is instantaneous, such that radiations observed across all filters originate from a single set of temperature and surface conditions (such as oxidation) -thus, explicit reference to T is omitted henceforth [42]. Next, it is assumed that the considered wavelength band Δλ is sufficiently narrow that the spectral variation of the emissivity can be approximated by a linear model, A low-order polynomials model is reasonable from a mathematical standpoint as it constitutes the Taylor series expansion of the (continuous) emissivity function, known to provide a reasonable approximation within sufficiently small domains. This model allows the regularization of the otherwise ill-posed inverse problem of temperature recovery, a reason why such models are commonly used in MRT systems [29,43].
Another simplification comes from confining measurements to atmospheric windows, short measurement distances (up to several meters), and optically isolated conditions. Then, the atmospheric absorption and emission become negligible and therefore τ atm may be treated as unity and R a may be omitted [44].
Lastly, reflections are assumed to have been suppressed by the user, either through experimental means (physical blocking) or via postprocessing [40,[45][46][47]. Thus, reflections are treated as originating from a blackbody at room temperature, in which case they are negated during the dark frame subtraction stage of the measurement methodology employed in this work (detailed in [48]). For this reason, reflections are not explicitly treated in the scope of the present investigation.
As a result, there remain three unknowns in the system: one temperature and two parameters of the emissivity model (slope mand intercept n). Then, the photo-response for a linear-emissivity target, R k,lin , can be expressed as follows: The meaning of G k (T) can be elucidated by factoring the R BB,k term in Eq. (7), We note that the ratio of integrals appearing in (8) represents a continuous weighted average of λ k over the integration domain, Then, combining (8) and (9) yields the measurement equation for filter k: where λ k (T) is a temperature-dependent characteristic wavelength of the optical system employing filter k. This result follows an approach similar to [38,[49][50][51], in which the characteristic wavelength (typically the central wavelength, or CWL, of the filter) of an optical configuration is not considered constant, but insteadis a function of the target's temperature. Thus, by allowing the characteristic wavelength to change, this approach provides additional flexibility to multi-wavelength pyrometry models that assume measurement bands to be narrow enough such that incoming radiation is absorbed at discrete wavelengths [29].
During the present investigation, no single model for λ k (T) was found suitable for all the different filters. Instead, it was determined empirically that λ k (T) is well-captured by one of three models (the model with the lowest sum of squared errors (SSE) is chosen for each filter): where (11)-(13) are library models available in the MATLAB Curve Fitting Toolbox' fit function, under the designations 'power2', 'exp2', and 'rat11', respectively. Thus, equation (10) can be used to predict the optical system's response, R k , when observing a target of temperature T and emissivity ε lin . Interestingly, if targets are assumed to always exhibit a linear emissivity, it becomes possible to develop a system model that includes the temperature-related drift of the characteristic wavelength. This extended calibration process results in two continuous functions with a total of 6-7 parameters, providing a more detailed characterization of the measurement system. To the best knowledge of the authors, this is the first attempt that integrates the emissivity assumption into the system calibration procedure.

Simulated experiments
In an experimental scenario, R k,lin (T i ) and R BB,k (T i ) appearing in equation (10) would be measured directly as part of a calibration procedure. Calibrating the interpolation equation of the photo-response vs. temperature, R BB,k (T i ), is a well-established procedure (with specifics related to the present methodology discussed in [41]). Calibrating R k,lin (T i ) is similar-though with one important difference-instead of using a blackbody radiator, one must observe a target that exhibits a known spectrally-linear variation in emissivity (i.e. {m, n}), within the passband of the used filter. A suitable calibration target would be either a material with a known linear emissivity (such as certain metals) or a blackbody radiator observed through an optical element with a linear transmissivity function. If neither option is available, G k (T) (and in following, λ k (T) and R k,lin (T i )) may be computed numerically based on the available transmissivity and responsivity documentation 1 . This allows the recovery of λ k (T i ), and in turn, the fitting of R BB,k (T) and λ k (T) according to models (3) and (11)-(13), respectively.
Alternatively, if simulated data is sought, specific shapes for individual spectral transmittance functions appearing in equation (2) are assumed 2 , and the integrals appearing in equation (7) are evaluated to obtain G k (T) and R BB,k (T). In following, R BB,k is fitted with the interpolation equation (3), and the ratio G k /R BB,k is fitted with the models in (11). Finally, equation (10) can be evaluated for any desired combination of {m, n, T}. Thus, simulating signals allows the rapid evaluation of many different optical configurations. Fig. 2 exemplifies the R(λ)⋅τ lens (λ) term in equation (2) by showing the spectral responsivity of the InSb detector 3 [52] and the transmissivity of the L0118 lens [53] representative of the FLIR SC7600 camera. The figure also presents spectral bands where atmospheric absorption is significant (shown in light red) and negligible (shown in white), corresponding to τ atm (λ). The last spectral function, τ filter (λ), corresponding to the employed band-pass filter, is discussed in the following section.

Temperature recovery
Measurements conducted using N filters lead to a system of equations, (14), which contains N equations and 3 unknownsleading to an overdetermined system if N ≥ 4. Therefore, the procedure for recovering m, n, T involves solving an optimization problem. In line with previous works [27][28][29], the present study does so by minimizing the least-squares functional.
Since box constraints for the parameter ranges are easily identified and specified, a constrained minimization algorithm such as trust-regionreflective 4 can be used on top of Levenberg-Marquardt, as demonstrated in [31]. The residual of the system can be obtained by rearranging equation (14), dividing by the calibration functions R BB (T), and then equating to zero: The optimization objective chosen is minimal SSE of the residual, The subsequent optimization problem is defined as Suitable initial guesses for the temperature computation may be determined based on error statistics such as maximum absolute error, mean absolute error, standard deviation, interquartile range and number of diverged cases; then one can choose a guess that is most conducive to convergence. The present work uses the guesses: 1 While reasonable in a simulated study, this approach compromises accuracy when applied to temperature measurement in practice, as it defeats the purpose of a radiometric calibrationwhich implicitly assumes that a-priori information about the optical components is unavailable or unreliable. 2 It should be noted that assumed knowledge of the spectral functions does not constitute "an unfair advantage" in favor of a simulated analyses, because calibration is introduced at the same stage as with experimental datathus ensuring maximum compatibility. 3 Although the FLIR SC7600 camera contains an SCD Pelican-D detector, personal correspondence with SCD engineers pointed to another detector, Blue Fairy, whose responsivity was said to represent similar detectors by the company sufficiently well. 4 The exact implementation used in the present work is a vectorized adaptation of MATLAB's lsqnonlin solver (from the Optimization Toolbox) and SciPy's least_squares solver (from the Optimization and Root Finding module, optimize).
for nonpositive-sloped (hereinafter referred to as "negative" for easier differentiation), positive-sloped and uncategorized (i.e., containing a mix of profiles of the previous types) emissivities, respectively. It is pointed out that the guesses in (19) are suitable for the particular range of targets, filters, and solver configurations used, meaning they should be re-evaluated when the aforementioned aspects are modified.

Filter definition
The configuration of optical filters is perhaps the main DOF available to designers and users of multispectral thermography systems. Within the scope of the present work, symmetric trapezoids as shown in Fig. 3 are said to represent commercial filters sufficiently well, and are deemed to provide sufficient variability in the resulting filter combinations. In this work, the transmittance in the blocking regions is fixed to 1% (optical density of 2) and so four parameters are required to define the shape. Lastly, filters are always chosen such that their main lobe falls entirely within one of the atmospheric windows (of CO 2 and H 2 O) shown in Fig. 2 in white. The 2.2 − 2.4 μm band is also ignored due to low responsivity of the considered thermographic system there.
After choosing the parameters for a trapezoidal filter and validating its shape, the subsequent procedure consists of these steps: (1) Mimicking experimental data acquisition, the two integrals shown in equation (7) are evaluated at several temperatures T i to yield G k (T i ) and R BB,k (T i ).
(2) Mimicking experimental calibration, interpolation equations for R BB,k (T i ) and G k (T i )/R BB,k (T i ) are fitted according to models (3) and (11) When these steps are performed for the N filters of a combination, over a grid of K temperatures and P combinations of {m, n} (representing different linear emissivities), it results in the measurement matrix R of size (K⋅P) × N that represents the entire variety of measurements the said combination could produce: where the (⋅) lin subscript has been omitted from all R i in (20) for brevity. Each combination analyzed in this stage of the investigation was evaluated for a K = 50, P = 10 4 and N = 4.

Filter optimization framework
It appears that the only way to quantify the goodness of a filter combination is by simulating measurements (R shown in equation (20)), estimating temperature, and studying the resulting errors. Unfortunately, there is no straightforward way to associate the shapes of individual filters to the performance of the combination as a whole, nor can the smoothness of the aggregate error be guaranteed. Thus, it appears that an efficient gradient-based search cannot be performed, and one must resort to the more computationally intensive derivative-free algorithms. However, if one or more relevant heuristics can be identified, they could serve as a more efficient goodness metric in an automated optimization scheme. Ultimately, the choice of heuristics is arbitrary, yet an educated guess based on mathematical considerations is possible. Also, due to the arbitrary nature of the heuristic selection process, their relative importance cannot be determined in advance. Thus, an a-posteriori optimization approach is taken: an initial design space exploration and a formation of a Pareto set are performed using potentially fast-to-evaluate heuristics, followed by a sifting (or "decision") stage that employs the more computationally-heavy recovered temperature error criterion.  Fig. 3. Parameters required to define a trapezoidal filter. In a trapezoidal shape the blocking regions and the peak are horizontal and connected by linear segments. A symmetric (with respect to the central wavelength) trapezoidal filter with a fixed optical density has four degrees-of-freedom.

Heuristics
We note that although the measurement equation (5) is non-linear, it could still be useful as a first approximation, to use linear analysis tools, such as the singular value decomposition (SVD) algorithm. SVD allows one to compute the condition number(κ) of a matrix, where SV 1 is the first (and largest) singular value of R , SV N is the last (and smallest) singular value of R , and N corresponds to the number of measurement channels (i.e., filters). A small condition number indicates that information is more evenly distributed between different measurement channels and less sensitive to perturbations in any one channel-a useful metric even when explicit pseudo-inversion is unnecessary. However, this criterion possesses a couple of downsides: a) it gives no information about intermediate singular values (i.e., SV 2 , ⋯, SV N− 1 ), and thus the degree of ill-posedness [54]; b) its minimization produces narrow filters with low peak transmittance 5 . For these reasons, a mechanism to compensate for the effects of condition number minimization is required. By considering equation (2), one can deduce that filters with a larger integral area tend to produce correspondingly larger values of R. Thus, maximizing the elements of the measurement matrix could result in larger filters, in terms of both peak transmission and passband. The Frobenius norm aggregates the values of individual matrix elements and is thus useful as a maximization objective for the present problem.
One other effect that must be considered is the validity of the locally linear emissivity assumption. As the considered band grows wider, it becomes increasingly inaccurate to assume that the target's emissivity is indeed linear within the entire band. If this effect is neglected, filter combinations tend to span the entire allowable wavelength range in an attempt to optimize (21) and (22). This shortcoming of the previous heuristics is addressed by introducing another minimization objective. This objective indicates the "spectral bandwidth" of the whole combination, defined as the distance between the farthest edges of the passbands of the farthest-apart filters, is mathematically given by: A significant difference between previous heuristics and Δλ C is the fact it is computed based on a combination's parameter vector and not the measurement matrix. This reason makes it trivial to optimize (minimize) Δλ C manually, which indicates that Δλ C provides little information on its own. However, Δλ C is useful as a supporting objective, because its presence means that optimal κ and ‖M‖ F would be obtained for many different spectral bands.
Optimization of the three heuristics (21)-(23) yields a Pareto set, which must then be further ranked to arrive at a single, most suitable solution.

Final decision stage
The outcome of the 3-objective optimization is a set of possible filter combinations forming a Pareto set, and it is unclear which of these combinations should be preferred. A reasonable approach is to find an additional constraint or criterion to narrow down the solution set further -converging to one specific configuration. This stage is done by numerically evaluating temperature recovery performance (as described in Section 2.3) over a broad range of experimentally-obtained emissivity profiles. This step differentiates itself from the heuristics evaluation by obtaining R k (T, ε(λ) ) via evaluation of the integral appearing in equation (5) instead of R k,lin (T, m, n) that uses a locally linear emissivity assumption, appearing in equation (10). Thus, after estimating temperature and finding the recovery error, combinations that exhibit the lowest errors on the real emissivity profiles are preferred.
For this, measured spectral emissivity profiles of aeronautical alloys and thermal-barrier-coated samples (YSZ) were collected from several sources [55][56][57][58], resulting in a dataset of 61 profiles (hereinafter "the default set"). Each material has one or more associated emissivity profiles, differing in oxidation conditions, thermal treatment, surface roughness, and other features. The number of profiles per material is shown in Fig. 4 (left). The nonlinearity of included profiles is summarized graphically in Fig. 4 (right) by plotting the means and the standard deviations of the magnitudes of the x 2 coefficient resulting from fitting a 2nd-degree polynomial to individual emissivity profiles. Spectral emissivity profiles of all materials in the default set are charted in Fig. 5, split into two, where groups (a) and (b) consist of material profiles with negative and positive first derivatives, respectively.
Finally, to assess the temperature recovery performance of a filter combination, recovered temperatures (T rec ) are subtracted from the known ground truth temperatures (T GT ) to obtain the absolute temperature error. In following, for a chosen threshold of temperature error ΔT, a success rate (SR) is defined as SR is a convenient metric for this purpose because it requires no assumptions about the shape of the error distribution.

Pareto set
The minimum condition number, maximum Frobenius norm and minimum spectral width objectives (represented in equations (21)-(23), respectively) along with an initial guess consisting of a random valid filter combination were provided to MATLAB's multiple-objective paretosearch solver. On a system with an Intel i7-6700K, 32 GB of RAM, and an NVIDIA GTX 660 GPU, this computation took around an hour. The solution yielded a Pareto set consisting of 578 unique filter combinations, as shown in Fig. 6. Entries of the Pareto set represent the results of the search stage, consisting of solutions that are optimal to a different extent for each objective. Certain regions of the parameter space are not feasible due to the provided constrains and the bounds of the variables. As a result, the Pareto set contains regions absent of solutions. Since it is unclear which combination is preferable, a decision process to select the most appropriate solution is performed.

Choosing from the Pareto set
In order to select the most appropriate filter combination from the Pareto set, temperature recovery success rates were computed based on simulated noiseless camera signals for all 61 emissivity profiles at 27 temperatures (323 K-973 K in 25 K steps), according to equation (10). As the acceptable error in most practical applications within this temperature range is around 10 K, success rates were calculated according to equation (24) for thresholds of 5 K (SR 5 ), 10 K (SR 10 ), 15 K (SR 15 ) and 20 K (SR 20 ). The highest scoring combinations were returned for each SR threshold, and the filter combination that consistently ranked high was chosen as the optimal. This sort of computations took around a minute per considered material. Fig. 7 shows the filters resulted from this sifting process, hereinafter referred to as the globally-optimal combination 5 Since for a scalar a and a matrix M, κ(a⋅M) = a⋅κ(M), a trivial way to decrease the condition number is finding matrices where the relative element sizes are similar, but the absolute values are as small as possiblewhich coincides with the smallest allowed filter sizes.
(GOC). This combination features filters that are spread out over most of the allowable range, including both sides of the CO 2 band at 4.2-4.4 μm.
The filters are visibly rectangularan allowed degenerate form of a symmetric trapezoid in which the top and bottom bases are of the same width. The success rate for an error threshold of 10 K is SR 10 = 47%. The GOC is marked in Fig. 6 by a square. In terms of heuristics, GOC exhibits a small κ, a large ‖M‖ F and a large Δλ C . This result indicates that the combination of a small condition number and a large Frobenius norm is indeed a useful predictor of temperature recovery capabilities. Furthermore, it confirms the lesser significance of the spectral bandwidth heuristic.
In order to assess the factors influencing the aggregate SR values shown in Fig. 7, recovery errors at individual ground truth temperature T GT and materials were considered, Fig. 8. Large errors were encountered for all materials at low temperatures (320-370 K). Interestingly, materials could be categorized into two distinct groups, according to their temperature recovery errors: the first group consists of DD5, K417G, and K77, which have fairly narrow regions of low error (< |10 K|), roughly between 750 and 850 K; and the second group consists of HastelloyX, Inconel718, TC4, and YSZ, for which the low error region is considerably broader, roughly between 400 and 900 K. After juxtaposing these two groups with the emissivity profiles appearing in Fig. 5, a clear correlation appears between the slopes of the emissivity profiles and the error groups. This analysis indicates that recovery errors are highly dependent on the slope of the emissivity profile and the temperature range.
In order to elucidate these findings, we observe the change of SR 10 with respect to T GT for negative-sloped profiles (m ≤ 0) and positivesloped profiles (m > 0), Fig. 9. Evidently, SR 10 of the GOC is higher in negative-sloped profiles compared to positive ones, reaching a level of 100% for temperatures between roughly 420-500 K. These observations Fig. 4. Composition of the "default set" of normal spectral emissivity profiles. (left) Prevalence of materials withing the dataset, where bar height indicates the number of profiles per material. (right) Parabolicity of materials within the dataset (quantified using the magnitude of the coefficient of x 2 in a 2nd-degree polynomial fit to the emissivity profile), where circles represent mean, error bars represent standard deviation, and the dashed line represents the allowed threshold. "TC4" is another name for the Ti-6Al-4V alloy.

Fig. 5.
Normal spectral emissivity profiles comprising the "default set". The shaded rectangular region represents an absorption band of CO 2 , where filters are disallowed. "Kong-2017" refers to [55]; "Kong-2017a" refers to [56]; "Gonzalez-2012" refers to [57]; and "Gonzalez-2012a" refers to [58]. Fig. 6. Pareto set obtained from the 3-objective optimization. Contour lines are approximate and were created using locally weighted smoothing (Lowess). Optimal filter combinations for specific scenarios are marked with red symbols: GOC (globally optimal combination); POC (optimal combination for positive emissivity slopes); NOC (optimal combination for negative emissivity slopes); and K77, YSZ, DZ125a/b are optimal combinations for specific emissivity profiles. No solutions were found in the lower-left corner, so the extrapolation is not shown.
indicate that narrowing the emissivity and temperature ranges may be a promising strategy for finding highly specialized combinations, potentially resulting in higher SR on the more confined domain.

Effect of segmenting the set of targets
The analysis of the GOC indicated emissivity characteristics have a significant effect on SR. To quantify this effect better, targets were separated into two groups: having either positive or negative emissivity slopes, and optimal filter combinations were obtained for each scenario independently. In these simulations, the temperature range was narrowed to 423-873 K (150-600 • C), a range that avoided the highest errors observed previously. The new emissivity and temperature ranges were provided to the solver so that the optimization bounds could be tightened, and the initial guess updated. The resulting optimal combination for positive slopes (POC) and the optimal combination for negative slopes (NOC) are charted in Fig. 10 and are marked in Fig. 6 by an upward-pointing triangle and a circle, respectively. In terms of heuristics, POC exhibits a smaller (better) κ, a smaller (worse) ‖M‖ F and a larger (worse) Δλ C than GOC. This indicates that a small condition number is the most important criterion for positive slopes. The opposite effect is observed with NOC: κ is smaller (worse), ‖M‖ F is larger (better) and Δλ C is smaller (better) than in GOC. This indicates that a large ‖M‖ F is preferred for negative slopes.
By comparing Fig. 7 with Fig. 10, both filter shapes and aggregate SR appear similar between GOC and POC, which can be explained by the default set containing more positive-than negative-sloped profiles (35 vs. 26), biasing GOC towards the positive majority; secondly, dramatic increases in aggregate SR levels are observed for all ΔT when comparing NOC to GOC for negative targets, best seen in the SR 20 of 100% for NOC.
To better understand the difference between GOC, NOC, and POC, a comparison of SR 5 and SR 10 with respect to T GT , for their respective types of targets, is shown in Fig. 11. Firstly, NOC's improvement in SR 5 is seen across all temperatures, and is especially noticeable for 800 K with an improvement of over 50% over GOC. Improvement is observed in SR 10 as well. Secondly, although POC performs better than GOC if Fig. 7. Optimal filter combination when considering the default set of emissivity profiles for temperatures in the range of 323 K-973 K. This combination is also known as "GOC".  judging by SR 10 , a more complex relation is seen in SR 5 , where POC is only better for 498-698 K and 873 K.
These results show that segmentation of the target set in terms of temperatures and emissivity slopes has a significant impact on the resulting filter combination and the associated temperature recovery success rates. The benefits of segmentation highlight the value of prior information for obtaining highly optimized filter combinations.

Optimization for specific emissivity profiles
Whenever information about relevant spectral emissivity profiles is available to the optimizer, it can be used to determine tighter bounds on and better initial guesses for slope and intercept. When the set of profiles is known, a linear fit can be computed for each profile (as in equation (6)), followed by aggregation of the fit coefficients to obtain mean values for slope and interceptserving as improved initial guesses. Similarly, upper and lower bounds of the fit coefficients, expanded by a safety factor (e.g., 5%), serve as tighter bounds for optimization, reducing the size of the search space and helping avoid some local minima.
The current section demonstrates the above methodology in two scenarios: one that considers a family of profiles, and another that considers individual profiles.
3.3.1.1. Case of a positive-sloped material: K77. In this subsection, the considered set of targets consists only of positive-sloped profiles belonging to the K77 alloy, as shown in Fig. 12. The optimal combination for this set of materials within the narrowed temperature range of 423-873 K is shown in Fig. 13. SR levels for this combination are  comparable to those of NOC, indicating that good performance can be achieved for positive slopes when additional information about the targets is available and provided to the optimizer. In terms of heuristics, interestingly, although K77 profiles exhibit positive slopes, the optimal combination is closer to GOC and NOC than to POC. This combination exhibits a less optimal κ, a less optimal ‖M‖ F and a more optimal Δλ C compared both GOC and NOC. This indicates that for sets of fairly homogenous emissivity profiles the relative impact of Δλ C on the SR is more pronounced.

Non-linear profiles.
In this subsection, the methodology is demonstrated for profiles that were excluded from the default set for exceeding the allowed parabolicity threshold, belonging to the materials DZ125 and YSZ, as shown in Fig. 14. Each of the two profiles is treated individually, and it is pointed out that although knowledge of the exact emissivity profile allows reconstructing temperature mono-spectrally, the analysis still follows the same multispectral procedures as in previous sections. The YSZ profile discussed here differs from the profiles included in the default set by the thickness of the ceramic, previously being 18-24 µm, as opposed to the present 0 µm (ceramic-filled bond coating) [59]. Fig. 15 shows the optimal combination for YSZ, whose most noticeable features are an almost triangular filter (which is not a common shape in practice), and the fact that filters are concentrated in the 3.1-4.1 µm band, a region of relative linearity within the spectral profile of YSZ. In terms of heuristics, this combination features similar trends to the K77 optimum with respect to GOC/NOC, though with a much more noticeable reduction in both κ (worse) and Δλ C (better). Confined to the spectral band of linearity of the considered emissivity profile, this result shows that a large ‖M‖ F is preferred. Fig. 16 shows two possible optimal combinations for DZ125: the first exhibits better performance in SR 5 (42%) while reaching a lower value of SR 20 (53%); whereas the second starts with a lower value of SR 5 (26%) but reaches higher values of SR 20 (79%). In both latter combinations, filters are concentrated within the 3.1-3.5 µm band, where DZ125 exhibits a spectrally linear behavior with a single slope. Another observation is that in the top combination there is a significant overlap between two of the suggested filters, which may indicate that there is little benefit from using 4 channels for this target. The difference in performance between the two combinations is further visualized in Fig. 17, which shows temperatures correctly reconstructed for each SR level: the top combination appears more well-suited to low temperatures, compared to the bottom combination that achieves comparably low errors only for mid-high temperatures. In terms of heuristics, these two combinations feature a somewhat similar (and very low) Δλ C , yet significantly different optimality in κ and ‖M‖ F : the top combination has a lower (better) κ but also a lower (worse) ‖M‖ F . This may indicate that combinations that have a smaller κ tend to stagnate sooner, as can also be observed when comparing GOC and POC. This analysis demonstrates that it is possible to suggest optimal filter combination even for non-linear emissivity profiles; it provides additional evidence that the temperature range of interest plays a vital role in determining combination optimality; and finally, it highlights the tradeoff of choosing an optimal combination based on a specific SR level. For non-linear profiles, the Δλ C heuristic has a significant effect, where it appears to confine filters to linear subranges of the considered profiles. Table 1 presents objective function evaluations and success rates for the three optimized combinations. The first observation that can be made is that SR is much higher for negative profiles. The optimized designs exhibit high SR, per the target-set segmentation: NOC is the best combination for negative slopes, and POC is the best for positive ones, across all SR. Furthermore, the optimized combinations work reasonably well for both types of targets. Lastly, POC appears to outperform GOC for both types of targets, effectively making it "the new GOC" in the       17. Success rates vs. temperature for the two DZ125 combinations. White rectangles indicate that the temperature recovery error is within the acceptable bound for that SR, while black rectangles indicate unacceptably large errors. temperature range of 423-873 K. Based on these observations, it can be said that the algorithm yields successful filter combinations for a broad range of scenarios.

Shape definitions of the optimal filters
For the convenience of the reader, the trapezoidal shape parameters of all filters that make up the various optimal combinations presented in the previous sections are presented in Table 2.

Summary and conclusions
The design space of a multispectral thermography system consists mainly of the choice of spectral filters, whose combination has a significant impact on the system's ability to reconstruct temperatures. This study presented a complete methodology for reconstructing temperature from multispectral data and for numeric evaluation of filter combinations to minimize the temperature recovery error while considering user-supplied information of varying degrees of specificity.
A mathematical model for an IR camera's photo-response when observing linear-emissivity targets was formulated. This model led to the development of an extended optical system calibration technique, able to capture not only the relation between temperature and the detected signal but also between temperature and the drift of the characteristic wavelength. Next, it was shown how to utilize the new calibration technique for temperature recovery in the context of a non-linear least-squares minimization procedure employing a trust-region solver.
The approach taken in this work involves simulating the signals acquired from a multitude of targets having different emissivities for various filter combinations. Later, it is investigated how well these signals can be mapped back to temperatures, assuming linear emissivity bands. The question of choosing optimal filter combinations was approached using an a-posteriori multi-objective optimization procedure based on the simulated data. Band-pass filters were modeled by a trapezoidal shape with 4 DOF (resulting in a 16-DOF representation of a 4-filter combination). Exploration of the design space of filter combinations was performed using a pattern-search solver via the evaluation of three heuristic functions (condition number, Frobenius norm, and combination spectral width) over a grid of perfectly linear simulated targets, resulting in a Pareto set consisting of 578 combinations. The temperature recoverability of combinations in the Pareto set was assessed using a collection of 61 experimentally measured emissivity profiles reported in previous literature. The procedure for determining the optimal combination from the Pareto set was demonstrated for several scenarios differing in emissivity and temperature ranges.
Although complete target-independent thermometry remains out of reach, the present analysis lays the foundations for the design of thermography systems suitable for scenes containing objects with arbitrary temperatures and radiative properties. Considering that filter replaceability is a standard feature in most contemporary thermography systems and that a set of filters costs a fraction of a hyperspectral imager, Table 1 Comparison of objective function values and temperature reconstruction success rates between 423 and 873 K for the optimal filter combinations.  Table 2 Features of the filters in each optimal combination. HWHM stands for "half width at half maximum". the presented approach could be an attractive strategy in many scientific and industrial applications.

Future work
Follow-up experimental work necessitates the production of a lineartransmissivity element (or a suitable alternative) that would allow conducting the extended calibration. Assuming this prerequisite was met, several experimental investigations may be performed. Firstly, the entire measurement and temperature recovery procedure could be tested and compared to other thermometry methods or published experimental results. Secondly, one may assess the predictive ability of the optimization by comparing the predicted temperature recovery errors to those observed in practice. Thirdly, one may test the ability of the heuristic approach to suggest optimal combination by analyzing several filter combinations in the context of a certain target set, then comparing the predicted optimal combination with the experimentally-determined optimum.

Funding
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 853096).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.