Comprehensive characterization of human color discrimination thresholds

This section is organized as follows. We begin with a brief overview of the experimental stimuli and task (Task and stimuli), followed by a summary of how our model characterizes the full psychometric field (The WPPM) and a description of the nonparametric adaptive trial-placement procedure used to collect the data (Adaptively sampled trials). Having described these essential methods, we then present our core results (Threshold estimates from the WPPM) and evaluate the validity of our model (Validation of the WPPM). Finally, we compare our findings with previous measurements from the color discrimination literature (Comparison with previous measurements). Additional technical details are provided in Materials and methods and Appendix 1—12.

Participants ($N=8$) performed a 3AFC oddity task. On each trial, three blobby stimuli were shown in a triangular spatial arrangement—two identical reference stimuli and one comparison stimulus with a different surface color (Figure 1A). The comparison stimulus was pseudorandomly assigned to one of the three positions. Participants were asked to identify the odd one out. Stimuli were rendered using the Unity graphics engine, and color was controlled by varying the specified surface reflectance using RGB (red, green, blue) coordinates, with other scene aspects held constant. We used naturalistic stimuli to increase the relevance of our results for understanding color vision in the real world. Hedjar et al., 2025, provide a comparison of color discrimination using stimuli similar to ours versus traditional flat spatially uniform patches.

Figure 1

Download asset Open asset

We made spectral calibration measurements (Brainard et al., 2002) to establish the relationship between RGB and the light emitted from the display. These measurements allowed us to represent the stimuli in terms of the excitations of the human L, M, and S cones elicited by the stimuli, and more generally in any standard color space (Brainard, 1996; Brainard, 2003; Brainard and Stockman, 2010). For this study, stimuli were constrained to lie in the isoluminant plane passing through the monitor’s gray point and bounded by its gamut. This plane was then affine-transformed into a square ranging between −1 and 1 along each axis (Figure 1B; Appendix 1). We refer to the space in which the transformed plane lies as model space because it is directly related to the way we formulated our semiparametric model and also served as a convenient representation for the nonparametric adaptive trial-placement procedure we used.

As an overview of our modeling approach, we fit the color discrimination responses (coded as ‘correct’ or ‘incorrect’) with a novel model, the WPPM—a Bayesian probabilistic model that combines an observer model (specified through a likelihood function) with an expectation of smoothness in the internal noise limiting color discrimination (specified through a prior distribution). Once fit to the data, the WPPM yields a continuously varying field of covariance matrices that characterize the internal noise in the perceptual representation of color stimuli (Figure 1C and D). These covariance matrices, in turn, determine the entire psychometric field.

More specifically, we designed the observer model within the WPPM to formalize the intuition that the stimulus perceived as the most distant from the other two is identified as the ‘odd one out’. The internal representation of each stimulus is assumed to be noisy and modeled as a multivariate Gaussian with the same dimensionality as the stimulus space. We assume the mean of each distribution is given by the corresponding stimulus’ location in model space. In contrast, we allow the covariance matrices to vary across model space to account for differences in the encoding precision of the stimuli. Because discrimination thresholds depend on the relative sizes of signal change and internal noise, an alternative formulation could instead attribute threshold variation across stimulus location to nonlinearities in signal encoding while assuming constant internal noise or to a mixture of nonlinearities and stimulus-varying noise (Zhou et al., 2024). Our formulation should be understood as a characterization of the signal-to-noise properties that limit discrimination and not as a commitment to a particular interpretation of how these properties arise.

On each trial, the observer model has access to one sample from the distribution of each of the three stimuli—two identical reference stimuli and one comparison. The observer model computes the pairwise squared Mahalanobis distance between each pair of noisy samples, using the weighted average of the covariance matrices of the reference and comparison stimuli (Figure 1E). By using Mahalanobis distance to make decisions (instead of, for example, Euclidean distance), the observer accounts for the expected noise structure. The two stimuli whose pairwise distance is smallest are identified as the references, and the remaining stimulus as the comparison (the ‘odd one out’). Because there is no simple closed-form solution for this decision rule (Mullen and Ennis, 1991), we used Monte Carlo simulation to approximate the percent-correct performance (Observer model).

We expect the internal noise that limits color discrimination to vary smoothly across model space—i.e., small changes in the reference stimulus should produce only small changes in the corresponding internal noise. The WPPM reflects this expectation by placing a finite-basis Wishart process prior over the continuous field of covariance matrices (Wilson and Ghahramani, 2011). Intuitively, the Wishart process prior introduces a regularization term into the model: it penalizes rapid variation in the covariance matrix field. The strength of smoothness is controlled by two hyperparameters of the model, ε and γ (Figure 1C and D; Prior over the weight matrix).

To fit the model to each participant’s data, we found the maximum a posteriori estimates of the WPPM parameters, using gradient-based numerical optimization of the log-posterior density, defined as the sum of the log-prior density and log-likelihood function (Model fitting).

The best-fit model parameters, together with the observer model, allow us to read out percent-correct performance for any pair of reference and comparison stimuli. In particular, to read out a one-dimensional psychometric function, we select a reference stimulus and use the observer model to approximate performance as the comparison stimulus varies along a line (Figure 1F, left panels). The threshold distance along the line is defined as the distance that yields 66.7% correct. By repeating this process across many directions, we derive a set of threshold distances around the reference (Figure 1F, right panel). Given our assumption that internal noise follows a multivariate Gaussian distribution, these threshold distances form approximately elliptical contours, which we fit with ellipses for visualization. This approach is consistent with prior work showing that ellipses provide a good approximation of color discrimination thresholds (MacAdam, 1942; Brown and MacAdam, 1949; Noorlander et al., 1981; Noorlander et al., 1983; Poirson and Wandell, 1990; Krauskopf and Gegenfurtner, 1992; Knoblauch and Maloney, 1996; Danilova and Mollon, 2025), despite some reported deviations (Newton and Eskew, 2003; Shepard et al., 2016; Shepard et al., 2017). Notably, while we show threshold contours corresponding to 66.7% correct for visualization, once fit, the WPPM allows us to read out the full psychometric function for any reference and chromatic direction, effectively mapping the entire psychometric field. Given that the psychometric field is derived from the underlying field of covariance matrices that characterize internal noise, the smoothness constraint imposed on the covariance matrices naturally propagates to the threshold contours and the field itself.

Reference and comparison stimuli for each trial were selected using AEPsych (Owen et al., 2021; Letham et al., 2022), an open-source package for adaptive psychophysics. For the adaptive sampling model, we used a probit-Bernoulli Gaussian process (GP) model (Williams and Rasmussen, 2006) with a radial basis function kernel. As with the WPPM, the GP assumes smooth variation in performance across model space, but unlike the WPPM, it does not impose any specific parametric form on the internal noise or thresholds. The semiparametric constraint—multivariate Gaussian-shaped internal noise—was introduced only when fitting the WPPM. For this reason, we describe the adaptive trial-placement procedure as nonparametric (relative to the WPPM), while acknowledging that it incorporates some parametric assumptions that are less restrictive than those of the WPPM. This nonparametric approach ensures that our data collection was not biased by assuming the correctness of the WPPM prior to validation.

Each participant completed 6000 AEPsych-driven trials: the first 900 were generated using quasi-random Sobol' sampling (Sobol’, 1967) to provide an adequate initialization for the GP; for the remaining 5100 trials, the GP was updated continuously based on participants’ responses, and each trial was adaptively selected to be most informative for estimating the thresholds targeted at 66.7% correct (Figure 2A, Appendix 2—figures 1 and 2).

Figure 2

Download asset Open asset

Adaptive sampling with AEPsych requires solving two optimization problems: one for updating the GP model (Williams and Rasmussen, 2006) and another for selecting the next trial using the expected absolute volume change (EAVC) acquisition function (Letham et al., 2022). To reduce computation time, we updated the GP model only every 20 trials. Sometimes, however, either the fitting or the trial selection process did not complete in time for the upcoming stimulus presentation. To avoid perturbing the participants’ rhythm, in these cases we slotted in pregenerated fallback trials (Appendix 3). The number of fallback trials varied across participants, ranging from 0 to 466 (Appendix 2—figure 1). These trials were included along with the 6000 AEPsych-driven trials when fitting the WPPM.

For each participant, we fit the WPPM to the 6000 AEPsych-driven trials, along with any additional fallback trials. To visualize the fits, we read out the elliptical threshold contours around a grid of reference stimuli (Figure 2B for a representative participant). The threshold contours revealed three key regularities: (1) thresholds were lowest for references near the achromatic point defined by the background behind the blobby stimuli, (2) thresholds increased with the distance of the reference from the achromatic point, and (3) the major axes of the elliptical threshold contours tended to be radially oriented with respect to the achromatic point. These regularities are consistent with previous results in the color discrimination literature, as explained further in Comparison with previous measurements.

The data were broadly consistent across participants, in the sense that the three regularities noted above were observed in the individual participant data (Appendix 2.1). In the model-space representation, individual variability was lowest near the achromatic point, where sensitivity was highest, and increased with the distance between the reference and the achromatic point (Figure 2C). Specifically, this variation was quite large in the upper-right quadrant of model space, where ellipse orientations varied considerably. This variation in orientation was also apparent when examining the data in other colorimetric representations (Appendices 7–8). Increases in inter-participant variability with increasing thresholds have been observed in other perceptual discrimination tasks (Girshick et al., 2011; Aguilar et al., 2017; Hong et al., 2021).

To validate the WPPM estimates, we interleaved 6000 validation trials throughout the experiment. These trials were held out from fitting the WPPM. For each participant, we used Sobol' sampling to select 25 reference stimuli and associated chromatic directions, with a unique draw per participant. Along each sampled chromatic direction, we used MOCS to sample 12 comparison levels: 11 were evenly spaced, and one was selected to provide easily discriminable catch trials (Figure 2D). The comparison levels were selected based on a pilot dataset to account for variability in thresholds across different reference stimuli and chromatic directions (see Design for details). Notably, we intentionally avoided densely sampling around a small number of references to minimize differential perceptual learning between the trials used for fitting the WPPM and those reserved for validation (Horiuchi and Nagai, 2024).

For each of the 25 validation references, we fit a Weibull psychometric function to the 240 MOCS trials collected along the sampled chromatic direction and identified the comparison stimulus corresponding to 66.7% correct (see two examples in Figure 2F). We then used the WPPM fit (constrained by nonoverlapping trials) to extract threshold contours for each validation reference stimulus and to read out the threshold along the MOCS chromatic direction (Figure 2E). The 95% bootstrapped confidence intervals for the WPPM estimates overlapped with those from the Weibull fits in all 25 conditions for participant CH (Figure 2E) and in 22–25 conditions across other participants (Appendix 4.1). These results demonstrate a high degree of agreement between thresholds derived from the WPPM psychometric field and those derived from the MOCS validation trials for the 25 discrete conditions. This agreement indicates that the Wishart process prior we imposed did not lead to substantial oversmoothing, as the validation thresholds were estimated independently without any smoothness constraint. Also notable is that the sizes of the 95% bootstrapped confidence intervals for the WPPM and validation thresholds were similar (Figure 2F and G; Appendix 4.1).

To quantify the agreement, we performed a linear regression, constrained to pass through the origin, between the thresholds read out from the WPPM fit and those obtained using the validation trials (Figure 2G). The results further support agreement between the two sets of estimates (mean correlation coefficient = 0.84, range = 0.73–0.96). For seven out of eight participants, the regression slope was not significantly different from 1 (mean slope = 1.04, range = 0.96–1.10) (Figure 2H; see Appendix 4.1 for comparisons for each participant). To assess whether there were more subtle sources of bias not captured by the regression slope, we analyzed the residuals—the differences between the WPPM and validation thresholds. While we found no evidence that residuals depended on the orientation or shape of threshold contours read out from the WPPM fit, we did observe one small but statistically significant relationship: the model slightly overestimated thresholds when validation thresholds were low and underestimated them when validation thresholds were high (slope = −0.176, $t(198)=-5.727$, p<0.001, $R^2=0.142$). However, the magnitude of this bias was small (Appendix 4.2).

As an additional benchmark, we simulated trials and responses from a ground-truth WPPM instance chosen to approximate the CIELAB ΔE94 metric and fit the model to the simulated data (Appendix 5). This allowed us to assess the ability of the WPPM to recover simulated ground truth, which is not possible with human data. The readout threshold ellipses based on the WPPM fit were in good agreement with the ground truth (Appendix 5—figure 2C). We then conducted the same validation analyses on the simulated data as described above. The thresholds read out from the WPPM fit agreed with 23 of the 25 validation thresholds, based on overlapping confidence intervals. A linear regression yielded a correlation coefficient of 0.86 and a slope of 0.92—well within the confidence intervals observed in participants’ data (Figure 2H). Residual analysis revealed a negative correlation with the magnitude of the ground-truth validation thresholds (Appendix 5.6; Appendix 4—figure 9), consistent with trends observed in human participants. With the simulated data, however, we can interpret the magnitude of this bias in the context of the agreement with ground truth and conclude that it is small (Appendix 5—figure 2C). Access to ground truth also provides us with additional ways to visualize patterns in the bias (Appendix 5.7).

Taken together, these results validate the accuracy of the WPPM and highlight the remarkable efficiency of our approach. With 6000 trials, conventional psychophysical methods allowed us to estimate percent-correct performance along only one chromatic direction across 25 references. In contrast, our new approach—combining nonparametric adaptive trial placement with post hoc fitting of the semiparametric WPPM—allowed us to map the entire psychometric field, providing percent-correct performance for any reference–comparison stimulus pair in the isoluminant plane using the same number of trials.

The WPPM is equivariant under affine transformations of color space (Appendix 1.4), allowing threshold contours derived in our model space to be transformed into other colorimetric representations. This flexibility enables direct comparisons with color discrimination thresholds reported in the literature. At the outset, we emphasize that the size and shape of threshold contours depend on ancillary experimental factors, including task design, stimulus spatial and temporal properties, and participants’ state of adaptation. Given these differences, we do not expect quantitative agreement across studies. Nonetheless, such comparisons help set our findings in the context of the literature. To illustrate, we present several such comparisons below, in the colorimetric representations used in the original studies.

We first compared the overall pattern of threshold variation in the isoluminant plane with measurements made by MacAdam, in which color matches were obtained using the method of adjustment (MacAdam, 1942) (Appendix 6). In his seminal work, the ellipses do not represent discrimination thresholds per se, but rather the standard deviation of color matches for each reference stimulus. Nevertheless, we consider his measurements to be comparable to ours, based on the linking assumption that discrimination thresholds are proportional to the internal noise that governs the variability of the appearance-based matches (Crozier and Holway, 1937). We observed a similar global structure in how the orientation and scale of the ellipses vary with reference stimulus. As expected, the absolute sizes of the threshold contours differ between studies (Figure 3A; MacAdam ellipses magnified by 10× and ours magnified by 2×). In addition to differences in stimulus spatial and temporal structure, it is worth noting that in MacAdam’s experiment, participants controlled the stimulus duration themselves (Wandell, 1985) and their state of adaptation differed considerably across reference stimuli (Krauskopf and Gegenfurtner, 1992). Despite these differences, the general correspondence between the datasets is apparent. It is also noteworthy that MacAdam’s results are based on 25,000 adjustments at a limited number of reference locations, whereas our ∼6000 forced-choice responses enabled us to characterize discrimination performance across all in-gamut reference–comparison pairs in the isoluminant plane.

Figure 3

Download asset Open asset

In a more recent study, Danilova and Mollon, 2025, measured threshold contours across a relatively broad region of the isoluminant plane, with sparse sampling of reference stimuli. The experimental paradigm in their study closely resembled ours: both used a fixed adapting point—D65 in their case and the monitor gray point in ours—and employed an oddity task to estimate discrimination thresholds. They used a 4AFC design combined with an adaptive staircase procedure, whereas we used a 3AFC version. To compare our data with theirs, we transformed our discrimination threshold contours read out from the WPPM fit into the same scaled MacLeod–Boynton space used in their study (MacLeod and Boynton, 1979) (Appendix 7). Despite methodological and stimulus differences, our results replicated the overall pattern of variation in ellipse orientation and size across color space. In particular, thresholds were smallest near the adapting point, increased with distance from it, and the ellipses generally pointed toward the adapting point. We observed closer agreement in the absolute sizes of the threshold contours than when comparing with MacAdam’s data (Figure 3B; their ellipses were magnified by 4× and ours by 1.5×). An interesting commonality between our data and those of Danilova and Mollon, 2025, is the rotation of the ellipses at the adapting point relative to the axes of the MacLeod–Boynton space, a rotation seen in all of our participants (Appendix 5—figure 1). See their discussion of this rotation for possible mechanistic interpretations.

In the next comparison, we turned to the study by Krauskopf and Gegenfurtner, 1992, whose measurements were concentrated within a small region near the achromatic point. Their experiment used a fixed adapting point and a 4AFC oddity task, with individual thresholds estimated using a three-down-one-up staircase procedure. To enable direct comparison, we read out threshold contours from our model at the same set of reference stimuli they used. Their results revealed two key features: (1) the threshold contour was smallest at the adapting point, and (2) as the reference moved away from it, the contours generally became increasingly elongated along the axis pointing toward the adapting point. Both features were observed in our data, albeit with some inter-participant variability (Figure 3C; Appendix 8). Our measurements differ from theirs, however, in the orientation of the ellipse at the adapting chromaticity: in our data, the ellipses are rotated with respect to the DKL axes for all our participants.

Lastly, we compared our results with iso-distance contours obtained with different versions of the CIELAB ΔE color-difference metrics (CIE, 2004; CIE, 1995; Luo et al., 2001; CIE, 2001). Although ΔE metrics were developed to describe suprathreshold perceptual color differences for a stimulus configuration that differs from ours, comparisons with threshold-level measurements are of interest—particularly because of the widespread use of ΔE to equate perceptual differences in studies of cognitive processes (Winawer and Witthoft, 2023; Garside et al., 2025). To derive a threshold contour for any given reference stimulus, we identified the comparison stimuli corresponding to $\mathrm{\Delta }E=2.5$ across multiple chromatic directions and fit an ellipse. While the choice of $\mathrm{\Delta }E=2.5$ is arbitrary, it primarily affects the overall size of the contour rather than its shape. The comparison reveals that the iso-distance contours of the original CIELAB ΔE76, which remain widely used, bear little resemblance to our threshold contours (Figure 3D, left panel; Appendix 9). The large deviations we observed between ΔE76 and our data provide further caution against the practice of using ΔE76 to predict perceptual color differences. In contrast, the more recent ΔE94 and ΔE00 metrics provided a much closer match (Figure 3D, middle and right panels), with only modest deviations from our measurements. These deviations may arise from differences between threshold and suprathreshold perceptual judgments, as well as from discrepancies in experimental conditions between our study and those used to constrain the parameters of the CIELAB ΔE metrics. An important feature of our data is that it enables such comparison with any perceptual metric across the isoluminant plane.

In this study, we demonstrated a data-efficient approach for achieving a comprehensive characterization of human color discrimination thresholds. Participants performed a 3AFC oddity task and completed 6000 trials that were specifically targeted near threshold using a nonparametric adaptive trial-placement procedure (Owen et al., 2021; Letham et al., 2022). We then developed and fit a novel WPPM to these adaptively sampled trials, along with a small number of fallback trials. The WPPM defines a continuous mapping from each reference stimulus to a covariance matrix that characterizes the associated internal noise. This mapping, in turn, enables predictions of discrimination performance for any pair of reference and comparison stimuli, effectively mapping out the full four-dimensional psychometric field. To evaluate model validity, we interleaved 6000 additional validation trials to estimate 25 probe psychometric functions. The results revealed that thresholds read out from the WPPM closely matched those derived from the validation trials, supporting the model’s accuracy. Thus, by combining nonparametric adaptive trial placement with post hoc fitting of the semiparametric WPPM, we achieved an unprecedentedly comprehensive characterization of color discrimination in the isoluminant plane.

Our measurements align qualitatively with previous studies that either used sparse sampling or targeted a small region of color space (MacAdam, 1942; Krauskopf and Gegenfurtner, 1992; Danilova and Mollon, 2025). Moreover, our measurements provide a more comprehensive characterization, in that the WPPM allows direct readout of a threshold contour at any reference stimulus without the need for additional measurement. Additionally, for studies examining how thresholds vary with factors such as stimulus size, presentation duration, or adaptation state, our approach offers a scalable and data-efficient way to measure how these factors affect the psychometric field. Finally, we have performed simulations and collected preliminary data that indicate it will be feasible to fully characterize the color discrimination psychometric field across the three-dimensional gamut of our display (Hong et al., 2026), a goal that has previously been described as ‘hopelessly difficult’ (Schrödinger, 1920).

A key assumption of the WPPM is that internal noise varies smoothly across stimulus space. This smoothness assumption is implemented through a prior on the variance of the weights applied to the model’s basis functions (Prior over the weight matrix). The smoothness prior plays a nontrivial role in the final WPPM estimates and therefore requires careful selection.

Our cross-validation analyses indicate that the smoothness hyperparameters used in the main analyses fall within a regime that balances oversmoothing against excessive uncertainty in the estimates (Appendix 10.1). When the smoothness imposed by the prior is too strong, the model produces overly uniform threshold estimates that fail to capture structure in the data. When the prior smoothness is too weak, the estimates become more variable. Consistent with this bias–variance tradeoff, agreement between WPPM and validation thresholds starts off low under strong smoothness, increases as the constraint is relaxed, and declines once it becomes too weak (Appendix 10.2). These two analyses narrow the range of sensible hyperparameter values and support our choice of $ϵ=0.4$ and $\gamma =0.0003$ for the main analyses.

As a general matter, determining appropriate prior hyperparameter values can be challenging when interpreting data with Bayesian models. This problem is at the heart of empirical Bayesian approaches, in which prior hyperparameters are estimated from the data, and then the data are analyzed with these estimates (Efron, 2012). A full empirical Bayesian approach of this sort currently exceeds what we can compute under reasonable time constraints. In our experience, evaluating model performance across a range of hyperparameters using cross-validation helps identify the region that balances oversmoothing against excessive uncertainty in the estimates, while the inclusion of validation trials helps identify the regime that maximizes agreement between WPPM predictions and validation thresholds (Appendix 10.3).

The number of basis functions included in the model is another modeling choice we had to make. To evaluate this choice, we examined the fitted weights as a function of Chebyshev polynomial order and found that they decayed to near zero at the highest polynomial order used in the model. This indicates that, given our choice of hyperparameters, including additional basis functions would not materially affect the inferred psychometric field (Appendix 2—figure 4).

Consistent with a well-established body of evidence, we found that thresholds were smallest near the achromatic reference, reflecting heightened sensitivity at the adapting point (Craik, 1938; Brown, 1952; Hurvich and Hurvich-Jameson, 1961; Pointer, 1974; Loomis and Berger, 1979; Krauskopf and Gegenfurtner, 1992). In addition, threshold contours were oriented toward the achromatic center, in agreement with previous findings (Krauskopf and Gegenfurtner, 1992; Gegenfurtner, 2025; Danilova and Mollon, 2025). Although the WPPM characterizes the data in terms of stimulus-dependent noise, the fitted psychometric field can be used to evaluate mechanistic models that posit specific transformations between stimuli and their internal representations.

The observation that the size and orientation of the elliptical threshold contours vary with the reference stimulus rules out mechanistic models that posit a linear transformation of cone excitations into three post-receptoral channels followed by fixed additive noise. Such models predict identical ellipses across the stimulus space. Moreover, the observation that the orientation of the elliptical threshold contours changes across reference stimuli also rules out mechanistic models in which a linear transformation of cone excitations is followed by limiting noise applied independently to each of the three channels. These models allow variation in the lengths of the major and minor ellipse axes, but predict that the orientation of these axes will be the same for all reference stimuli.

Cone-opponent models that posit noise and nonlinearities at multiple stages of processing, possibly with an overcomplete cone-opponent representation to capture parallel channels along the visual pathways, may be able to account for the observed data, as may models that invoke higher-order mechanisms (e.g. mechanisms narrowly tuned for hue). For more on relevant ideas, see Wyszecki, 1982; Wandell, 1995; Chen et al., 2000; Eskew, 2009; Stockman and Brainard, 2010; Hansen and Gegenfurtner, 2013; Shevell and Martin, 2017. Notably, mechanistic models are often tested using additional manipulations such as adaptation and noise masking; our approach can be extended to incorporate manipulations of such factors (Zhang et al., 2026), as well as of stimulus spatial and temporal structure and retinal location.

We studied a relatively young cohort of eight participants and found broadly consistent patterns across individuals, while also observing individual differences. Such variation has provided valuable insights into the mechanisms of color vision (Bosten, 2022) and is also of interest for understanding how much any given individual is likely to differ from an average characterization. A successful mechanistic model should allow investigation of whether our observed individual differences can be attributed to individual variation in biological factors known to influence color vision, such as preretinal absorption, photopigment spectral sensitivity, and the ratio of L to M cones in the mosaic (Neitz and Jacobs, 1986; Brainard et al., 2000; Kremers et al., 2000; Carroll et al., 2002; Hofer et al., 2005; Bosten, 2022; Rezeanu et al., 2023).

To enable studies involving larger and more diverse populations, further improvements in the efficiency of our approach are likely achievable. In the present study, we used a nonparametric adaptive trial-placement procedure to avoid biasing data collection by assuming the correctness of the WPPM. Given the validation of the WPPM presented here, future studies could instead incorporate adaptive trial-placement strategies tailored to the model, thereby improving efficiency. Alternatively, one could leverage the current dataset to develop stronger priors that capture the regularities we observed and use these priors to guide more efficient trial placement. Stronger priors could also increase the quality of estimates available from a fixed set of trials, although care should be taken to ensure that the prior does not overly constrain the estimates. Another approach is to develop mechanistic models with relatively few parameters, which could be estimated efficiently using parametric adaptive sampling (Watson, 2017). Finally, a complementary strategy is to increase the rate at which participants provide information about thresholds through more efficient psychophysical experimental paradigms (Agosti et al., 2026; Barnett et al., 2025; Burge and Cormack, 2024).

Another aspect of the WPPM framework that can be adapted for different applications is the mapping between stimulus space and model space, as well as the choice of basis functions. In the present work, we defined model space to be bounded by [−1, 1] for mathematical convenience, as this is the domain on which our chosen basis functions—the two-dimensional Chebyshev polynomials—are defined. These basis functions could in principle be replaced by alternatives, such as Zernike polynomials (Zernike, 1934; Thibos et al., 2000) or Fourier basis functions (Stein and Shakarchi, 2011), which may be better suited for stimulus domains that are disk-shaped or lack clear boundaries. The number of basis functions can also be adjusted based on prior knowledge about the expected smoothness of the psychometric field for a given stimulus domain. More generally, other approaches that leverage the smoothness of psychophysical performance and physiological responses have been developed (Rad and Paninski, 2010; Gravesen, 2015; Savin and Tkacik, 2016; Waz et al., 2025), including recent work on color discrimination (Koenderink et al., 2026).

A longstanding and fundamental question in vision science is whether it is possible to develop a perceptual metric that accurately predicts both threshold-level and suprathreshold judgments of color difference. For example, considerable effort has gone into attempts to find color representations in which the perceptual difference between two color stimuli is predicted by the Euclidean distance between their coordinates, as in the original 1976 CIELAB and CIELUV ΔE metrics (Brainard, 2003; Robertson et al., 1977). Our measurements directly establish a locally Euclidean metric for threshold-level differences. While threshold behavior is well described as locally Euclidean, suprathreshold judgments have been shown to violate the assumptions of a globally Euclidean geometry (Wuerger et al., 1995; Ennis and Zaidi, 2019). In particular, perceptual similarity judgments at larger distances often fail to satisfy key Euclidean properties. For example, judgment variability does not necessarily increase with Euclidean distance (Wuerger et al., 1995), and a stimulus that is equidistant from two endpoints is not necessarily perceived as equally similar to both (Ennis and Zaidi, 2019).

An alternative framework, originally proposed by Fechner, 1860, and explored subsequently (Schrödinger, 1920; MacAdam, 1979; Wyszecki, 1982; Zaidi, 2001; Koenderink, 2010; Bujack et al., 2022; Helmholtz, 2024; Stark et al., 2025), suggests that suprathreshold differences may be understood as the accumulation of small threshold-level differences along a path between stimuli. In this framework, color space is taken to be a Riemannian manifold—a space that is locally Euclidean but may be globally curved. The perceptual distance between two colors is hypothesized to correspond to the geodesic—the shortest path between them in terms of accumulated thresholds. This distance is computed by integrating local thresholds along all possible paths between the two points and selecting the path with the smallest total. In our observer model, this integration is effectively equivalent (up to a constant) to integrating internal noise along the path.

Testing this geodesic hypothesis requires knowledge of how internal noise (or thresholds) varies across color space, as this determines the geodesics. Our measurements provide the necessary knowledge for the isoluminant plane, enabling direct empirical tests of the geodesic hypothesis within this slice of color space, as well as elaborations of this hypothesis (Bujack et al., 2022; Stark et al., 2025). The results of such tests may depend on the particular experimental paradigms used to assess suprathreshold perceptual differences.

Because there is no guarantee that the geodesics between two stimuli in the isoluminant plane are themselves confined to this plane within full three-dimensional color space, testing the geodesic hypothesis in this plane based on our current data would be considered provisional. Nonetheless, such tests would provide valuable exploration of the perceptual geometry revealed by our measurements. As noted above, our approach makes it feasible to extend the measurements to full three-dimensional color space (Hong et al., 2026), which, when completed, will allow subsequent investigations to overcome this limitation.

It is possible that the geodesic hypothesis—and more generally the idea that threshold-level judgments can predict suprathreshold judgments—will fail. Nonetheless, we view understanding whether, when, and how such failures occur as central to guiding the development of a successful account of suprathreshold color-difference perception.

Our approach is generalizable to a wide range of perceptual tasks. A key insight that makes comprehensive characterization of human color discrimination thresholds feasible is the assumption—shared by both our model and the models implemented in AEPsych—that internal noise, and thus thresholds, vary smoothly across stimulus space. This smoothness assumption is not unique to color perception; it applies broadly to other domains. Indeed, smoothly varying elliptical or ellipsoidal thresholds have been reported in studies of motion perception (Reisbeck and Gegenfurtner, 1999; Champion and Freeman, 2010), auditory speed discrimination (Freeman et al., 2014; Carlile and Leung, 2016; Bertonati et al., 2021), motion-in-depth (Wardle and Alais, 2013), and numerosity perception (Cicchini et al., 2016; Cicchini et al., 2019; Cicchini et al., 2023). These parallels highlight the broader relevance of our framework and suggest that combining nonparametric adaptive trial placement with the WPPM could be a powerful strategy for characterizing perceptual thresholds across diverse domains.

This study was preregistered at a public repository. As described in the preregistration document, exploratory analyses were conducted on data from one participant (CH) prior to preregistering an initial hyperparameter choice of $ϵ=0.5$ for the main analysis. After data collection was completed, we performed the hyperparameter sweeps (Appendix 10). These led to our final choice of $ϵ=0.4$ and $\gamma =0.0003$.

Eight participants (six female, aged 22–30 years; seven right-handed) were recruited for the study. Six were paid volunteers who were naive to the purpose of the study. The remaining two were experimenters and participated without additional compensation. All participants had normal or corrected-to-normal vision (20/40 or better in each eye), assessed using a Snellen eye chart, and normal color vision, assessed using Ishihara plates. The study was approved by the Institutional Review Board at the University of Pennsylvania, and written informed consent was obtained from all participants prior to the experiment.

Stimuli were presented using an Alienware computer (Aurora R11) running Windows 10 Enterprise, equipped with an Intel Core i7-10700K processor and an NVIDIA GeForce RTX 3080 GPU. The display was a DELL U2723QE monitor (59.8 cm width, 33.6 cm height, 3840×2160 resolution, 60 Hz refresh rate). The monitor was positioned 189 cm from the chinrest, subtending 18.0 × 10.2 degrees of visual angle (dva). Monitor color and luminance measurements were obtained with a Klein K-10A colorimeter and a SpectraScan PR-670 radiometer. The display resolution was approximately 200 pixels/dva, above the typical human foveal resolution limit.

The Alienware computer was used solely for stimulus presentation, whereas adaptive sampling of the stimuli was performed on a separate custom-built PC with a high-performance Gigabyte motherboard (X299X Aorus Master), an NVIDIA GeForce RTX 3070 GPU and a 12-core Intel i9-10920X processor. This computer also ran Windows 10 Enterprise. The two computers communicated via a shared network disk, using a custom protocol based on text files that both computers could read and write.

A USB speaker (3 W output power, 20 kHz frequency response) was used for playing auditory feedback, and a gamepad controller (Logitech Gamepad F310) was used for registering trial-by-trial responses.

The visual scene (Appendix 11—figure 1A) was constructed in Unity (v2022.3.24f1) and rendered using its standard shader. The scene consisted of three identical blobby three-dimensional objects, each created in Blender (v4.0) with a matte, nonreflective surface. On each trial, the surface color of the blobby objects was varied by adjusting their RGB values in Unity. The three blobby objects (2.5×2.5 dva; 203,900 pixels each) were arranged in a triangular configuration (Figure 1A). Each blobby object was centered and floating inside its own cubic room (3.3×3.3 dva; $x=0.302$, $y=0.322$, $Y=66.1\mathrm{c}\mathrm{d}/{\mathrm{m}}^2$). Each room, along with the blobby stimulus inside it, was illuminated exclusively by an achromatic spotlight positioned in front of the object and set to maximum intensity ($R=G=B=1$). The three rooms were presented against a spatially uniform gray background (18.0×10.2 dva; $x=0.306$, $y=0.326$, $Y=116.8\mathrm{c}\mathrm{d}/{\mathrm{m}}^2$). The centers of the blobby objects were 3.7 dva apart.

We used a SpectraScan PR-670 to measure the monitor’s primaries and gamma function as rendered through Unity (Appendix 11.1). These measurements directly characterized the relationship between the specified RGB values for the blobby stimuli and the light emitted from the display. The same calibration was repeated for all three blobby stimuli, confirming consistent color behavior across screen locations. Based on these results, a single gamma correction—derived from the bottom-right stimulus—was applied to all three objects during the experiment. This correction was validated by remeasuring the output with gamma correction applied, showing good alignment with the predicted identity line. To confirm stability over time, we repeated the calibration 1 month into data collection and observed negligible changes.

Additionally, we used a Klein K-10A colorimeter to verify that the system achieved sufficient color depth. For this check, a single blobby stimulus was presented at the center of the screen, rather than in the full triangular arrangement. Measurements confirmed that Unity and our video chain were able to produce at least 12-bit color precision via the native 8-bit output and implicit spatial dithering that occurred across the surface of the blobby object during rendering (Appendix 11.2).

We restricted our stimuli to lie within the isoluminant plane that passes through the monitor’s gray point (i.e. $R=G=B=0.5$). To define the boundaries of this plane, we identified the limits of RGB values that remained within the monitor’s gamut. These boundary points formed a parallelogram in RGB space. We then computed an affine transformation that maps this parallelogram onto a square bounded by [−1, 1] (Appendix 1). The forward and inverse transformations enabled conversion between RGB and model space: stimuli were rendered in RGB space, while trial placement and model fitting were performed in model space.

We used AEPsych (v7.3) to sample a total of 6000 reference–comparison stimulus pairs. The first 900 trials were generated using Sobol' sampling (Sobol’, 1967), a ‘space-filling’ design based on a low-discrepancy quasi-random sequence. The remaining 5100 trials were adaptively selected to efficiently estimate thresholds across the entire psychometric field. Each stimulus pair was defined in the two-dimensional model space. As a result, the psychometric field comprised four variables: two specifying the reference stimulus, $x_0\in {\mathbb{R}}^2$, and two specifying a difference vector, $\mathrm{\Delta }\in {\mathbb{R}}^2$, which was added to the reference to define the comparison stimulus $x_1=x_0+\mathrm{\Delta }$. Reference values were constrained to [−0.75, 0.75] along each model dimension. Each element of Δ was constrained to [−0.25, 0.25] to ensure that all comparison stimuli remained within the $[-1,1{]}^2$ bounds of model space. During the initial 900 Sobol'-sampled trials, the difference vector Δ was scaled by one of three factors (1/4, 2/4, or 3/4) before being added to the reference stimulus. This controlled the distance between the reference and comparison stimuli, effectively modulating task difficulty. These scaling factors were evenly distributed and pseudorandomized across trials. For the remaining 5100 trials, all four variables were adaptively selected using AEPsych’s optimization procedure. Specifically, the underlying GP model was updated every 20 trials, and new trials were selected using the EAVC acquisition function (Letham et al., 2022), targeting 66.7% correct across the entire psychometric field.

In addition to the 6000 AEPsych-driven trials, we interleaved 6000 validation trials sampled using MOCS. Each participant was tested on 25 reference stimuli: one was fixed at the achromatic point and the remaining 24 were Sobol'-sampled within the isoluminant plane bounded by [−0.6, 0.6] along each model dimension. For each reference, a chromatic direction was Sobol'-sampled between 0° and 360°. Each validation condition consisted of 12 stimulus levels: 11 equally spaced along the sampled direction and one easily discriminable level, with each level repeated 20 times. These levels were determined based on a pilot dataset described in the preregistration documents.

The validation trials were pregenerated for each participant, pseudorandomized so that every 300 validation trials contained all unique trials (25 conditions × 12 levels). To minimize differential learning effects between AEPsych-driven and validation trials, we pregenerated a randomized sequence in which the two trial types were arranged in alternating pairs, with the order within each pair shuffled. However, because AEPsych occasionally required more time to determine the next trial placement, this sequence could not always be followed in real time. For this reason, we implemented a fallback trial strategy (Appendix 3): if, for any trial, AEPsych had not computed trial placement in time, the next validation trial was inserted to keep the experiment moving. If necessary, subsequent validation trials were queued, but this was capped at a lead of four validation trials ahead of AEPsych trials. Once the cap was reached and AEPsych was still not ready, pregenerated fallback trials were presented instead. These fallback trials were Sobol'-sampled with the difference vector Δ scaled by one of three factors (2/8, 3/8, or 4/8) to manipulate task difficulty. Validation trials resumed once AEPsych caught up. Notably, the fallback trials (range: 0–466) were included alongside the 6000 AEPsych trials when fitting the WPPM.

Participants performed a 3AFC oddity task. Each trial began with a fixation cross presented at the center of the screen for 0.5 s, followed by a blank screen for 0.2 s. Then, three blobby stimuli appeared inside the cubic rooms for 1 s. After participants responded, a blank screen was shown for 0.2 s, followed by auditory and visual feedback indicating accuracy (‘correct’ with a beep or ‘incorrect’ with a buzz). Consecutive trials were separated by a 1.5 s inter-trial interval (ITI). Participants were instructed that they could move their eyes freely during the stimulus presentation, but should maintain fixation while the fixation cross was on the screen.

Most participants (seven out of eight) completed 12 sessions. Each session began with 40 practice trials to familiarize participants with the task. This was followed by 1000 experimental trials–consisting of 500 AEPsych-driven trials and 500 predetermined validation trials–plus a small number of fallback trials. The validation trials were randomized, and the two trial types were fully intermixed. Participants took a break every 200 trials. Each session took approximately 1.5 hr to complete. In total, these seven participants completed between 12,256 and 12,466 trials, depending on the number of fallback trials inserted. Participant CH completed 12,000 trials across 10 sessions, without any fallback trials. As a result, the ITI was slightly longer for this participant, but we expect this to have had a negligible effect on performance.

Our implementation of the WPPM relies on two core assumptions about color perception: (1) the internal noise that limits color discrimination follows a multivariate Gaussian distribution, centered on the corresponding stimulus, with a covariance matrix that captures both the size and orientation of the noise, and (2) the covariance matrix varies smoothly across model space, without rapid local variations. In the following subsections, we describe the WPPM in five parts. First, we define the observer model, which predicts percent-correct performance for a given pair of reference and comparison stimuli by modeling both the noisy internal representations and the decision rule. Second, we describe how we use a finite-basis Wishart process to parameterize the entire field of covariance matrices across model space, along with the factors that control its smoothness. Third, we describe the weak prior imposed on the covariance matrix field to favor smooth variation. Fourth, we explain how, given a specification of the covariance matrix field, we compute the likelihood and thereby the posterior probability of the model parameters. Finally, we show how, once the model is fit, the covariance matrix for any reference–comparison stimulus pair can be read out and combined with the observer model to predict percent-correct performance, including threshold contours around any reference stimulus.

Covariance matrix field

Request a detailed protocol

The WPPM specifies a covariance matrix at any selected reference stimulus across the entire isoluminant plane. Each matrix specifies the perceptual noise in terms of the variance along the two model dimensions (${\sigma }_{\text{dim1}}^2$, ${\sigma }_{\text{dim2}}^2$) and their covariance (${\sigma }_{\text{dim1},\text{dim2}}$) (Figure 1C and D).

The covariance matrix field is constructed using one-dimensional Chebyshev polynomial basis functions (Chebyshev, 1853). We chose Chebyshev polynomials because they allow for the expression of smoothness over a bounded interval without imposing periodic boundary conditions. Let $x=[x_{\text{dim1}},x_{\text{dim2}}]$ denote a location in the two-dimensional model space. The basis functions are defined recursively for each model space dimension as given here for $x_{\text{dim1}}$:

(9) ${\displaystyle {\displaystyle T_0(x_{\text{dim1}})=1,}}$

(10) ${\displaystyle {\displaystyle T_1(x_{\text{dim1}})=x_{\text{dim1}},}}$

(11) ${\displaystyle {\displaystyle T_{i+1}(x_{\text{dim1}})=2x_{\text{dim1}}\cdot T_i(x_{\text{dim1}})-T_{i-1}(x_{\text{dim1}}),}}$

where $x_{\text{dim1}},T_i(x_{\text{dim1}})\in {\mathbb{R}}^n$, and $n$ is the number of discretized points along that stimulus dimension, which can be chosen flexibly to achieve any desired resolution. We construct two-dimensional basis functions by taking the outer product:

(12) ${\displaystyle {\displaystyle {\varphi }_{i,j}(x)=T_i(x_{\text{dim1}})\cdot T_j(x_{\text{dim2}}),}}$

where ${\varphi }_{i,j}\in {\mathbb{R}}^{n\times m}$, with $n=m$ representing the number of discretized points along each dimension of model space. We limited the number of basis functions to five per dimension, i.e., $i,j\in \{0,1,\dots ,4\}$, resulting in a total of 5 × 5 = 25 two-dimensional basis functions (Figure 4, first panel). The polynomial order of each two-dimensional basis function is given by $i+j$, with higher-order basis functions describing more rapidly varying patterns.

Figure 4

Download asset Open asset

The basis functions were weighted by a learnable parameter matrix, $\mathbf{W}\in {\mathbb{R}}^{5\times 5\times 2\times 3}$, where the first two dimensions index the Chebyshev basis functions along each model space dimension ($i,j\in \{0,1,\dots ,4\}$), and the last two dimensions index the output components ($k\in \{1,2\}$ and $l\in \{1,2,3\}$). The weighted basis functions are expanded into an overcomplete representation ${\mathbf{U}}_{k,l}\in {\mathbb{R}}^{n\times m}$ (Figure 4, second panel) as follows:

(13) ${\displaystyle {\displaystyle {\displaystyle {\mathbf{U}}_{k,l}(x)=\sum _{i=0}^4\sum _{j=0}^4{W}_{i,j,k,l}\cdot {\varphi }_{i,j}(x).}}}$

This overcomplete representation was then combined with its own transpose to yield a positive semi-definite covariance matrix (Figure 4, third panel), $\mathrm{\Sigma }(x)\in {\mathbb{R}}^{2\times 2}$ for $x$ at any discretized point in model space, i.e.,

(14) ${\displaystyle {\displaystyle \mathrm{\Sigma }(x)=\left[\begin{array}{cc}{\sigma }_{\text{dim1}}^2& {\sigma }_{\text{dim1},\text{dim2}}\\ {\sigma }_{\text{dim1},\text{dim2}}& {\sigma }_{\text{dim2}}^2\end{array}\right]={\mathbf{U}}_{k,l}(x)\cdot {\mathbf{U}}_{k,l}(x{)}^{\mathrm{\top }}.}}$

Notably, in our implementation, rather than matching the dimensionality of the intermediate representation U to that of Σ, we adopt an overcomplete parameterization motivated primarily by practical considerations. When we restricted the dimensionality indexed by l to 2, the optimization occasionally became ill-conditioned, leading to singular or unstable solutions. Expanding l to 3 substantially improved numerical stability and made the fitting procedure more robust. Increasing l beyond 3, however, would introduce additional degrees of freedom. We therefore selected $l=3$ as a compromise between model flexibility and numerical stability. Regardless of whether the representation is square or overcomplete, the resulting matrices are symmetric and positive semi-definite.

The weights are the free parameters of the model, determining the output covariance matrix field. The model is highly flexible, capable of generating a wide range of covariance matrix fields, from smooth to rapidly varying (Figure 1C and D).

Color calibration analyses were performed using MATLAB R2023b. We computed inverse gamma lookup tables from the measured gamma functions (Appendix 11) and derived transformation matrices to convert values from model space to RGB space (Appendix 1). Stimulus presentation, including gamma correction, was implemented in Unity, coded in C#.

All analyses were conducted in Python 3.11 using a variety of open-source packages. Model fitting was implemented primarily using JAX (Bradbury et al., 2018). Behavioral data were separated into AEPsych-driven plus fallback trials and validation trials. The WPPM was fit exclusively to the AEPsych and fallback trials. To assess variability in model estimates, we performed 120 bootstrap resamplings of the AEPsych-driven trials, preserving the original ratio between Sobol', adaptively sampled, and fallback trials in each resampled dataset. The WPPM was then refit to each of the 120 bootstrapped datasets.

To compute a 95% bootstrap confidence interval for the threshold contours, we first computed the summed normalized Bures similarity (NBS) score (Muzellec and Cuturi, 2018) between the predictions from the model fit to each bootstrapped dataset and those from the model fit to the original dataset, evaluated on a finely sampled grid of reference stimuli (−0.85 to 0.85 with 103 uniformly spaced points). Higher scores indicate greater similarity to the predictions from the original dataset. We then sorted the model fits by their summed NBS scores and retained the top 114 (95% of 120) fits. The confidence interval bounds were defined by the union and intersection of the threshold contours, subsequently computed for any reference stimulus using this fixed set of retained model fits.

For the held-out validation trials, we computed the Euclidean distance between each reference stimulus and its paired comparison stimulus. For each of the 25 conditions, a Weibull psychometric function was fit to the binary color discrimination responses, with the guess rate fixed at 33.3% correct. Threshold was defined as the distance to the comparison stimulus corresponding to 66.7% correct. To estimate variability, we bootstrapped each condition 120 times and computed 95% confidence intervals for the threshold estimates.

To assess the agreement between the thresholds predicted by the WPPM and those estimated from the validation trials, we performed linear regression (constrained to pass through the origin) between the two sets of predictions using the original dataset, as well as for each of the 120 paired bootstrapped datasets. We then sorted the resulting slopes and correlation coefficients and computed 95% confidence intervals separately for each. Additionally, we computed the number of conditions for which the 95% bootstrap confidence intervals of the WPPM-predicted thresholds and the validation thresholds overlapped as an additional measure of agreement.

This section summarizes the colorimetric transformations between RGB space of our monitor, the 2D model space, and standard color spaces.

DKL space provides a representation of the isoluminant plane with the adapting point at the origin (Derrington et al., 1984; Brainard, 1996). We defined our DKL space with respect to the CIE 2° physiologically relevant cone fundamentals and the corresponding photopic luminosity function. The adapting point was specified by the cone excitations elicited by the monitor’s gray point ($R=G=B=0.5$), such that the isoluminant plane passed through this point. Within the isoluminant plane, we densely sampled chromatic directions spanning 360° around the origin. For each direction, we searched outward from the origin until reaching the edge of the monitor gamut (details are provided in Appendix 1.1). Repeating this procedure across all directions yielded a set of gamut boundary points that formed a parallelogram in RGB space. From these boundary points, we identified the four vertices and recorded their coordinates in both DKL and RGB spaces (Appendix 1—table 1). We then used the corresponding vertex coordinates to derive a projective transformation matrix (homography) that maps coordinates from DKL space to model space (see Appendix 1.2 and Appendix 1—table 2). Although a homography provides a general solution for mapping between arbitrary quadrilaterals, the resulting matrix revealed that the transformation in our case was affine. We therefore used the affine transformation and its inverse to convert between linear RGB values and model space (see Appendix 1.3 and Appendix 1—table 2).

Appendix 3—figure 1

Download asset Open asset

For each participant, we ran 6000 AEPsych-driven trials interleaved with an additional 6000 validation trials. Although our initial plan was to present these 12,000 trials in a predetermined randomized sequence, we quickly realized that this approach was impractical. Under this design, AEPsych would only have the ITI to compute the next trial placement—a window that is difficult to optimize. A long ITI risks participant fatigue or loss of attention, whereas a short ITI does not give AEPsych enough time to complete its computations. To achieve both a smooth experimental flow and adequate computation time for AEPsych, we implemented a fallback trial strategy using a dual-computer setup.

In this setup, stimulus presentation was handled by a display computer, while adaptive trial placement using AEPsych ran on a control computer. The two computers communicated via a shared network disk using text files that both computers could read and write. This decoupled design provided modular separation between code specialized for stimulus presentation and code specialized for trial placement. It should also make the trial placement code easier to port to different stimulus display systems.

With the dual-computer design, AEPsych had at least 2.9 s to compute the next trial after the participant’s response (Appendix 3—figure 1). This window spanned both the post-stimulus period of the current trial (0.2 s blank, 0.5 s feedback, 1.5 s ITI) and the pre-stimulus period of the upcoming trial (0.5 s fixation and 0.2 s blank). Importantly, this computation window began only after the participant responded because AEPsych requires the participant’s response to update its model. The computation window ended just before the stimulus presentation of the next trial, when AEPsych had to deliver the values for the upcoming reference and comparison stimuli.

The fallback trial strategy ensured continuous stimulus presentation. If AEPsych failed to return a new trial within the 2.9 s window, we defaulted to presenting the next available MOCS trial from the pre-determined randomized sequence. In such cases, AEPsych’s computation continued in a separate thread and attempted to meet the following decision deadline, which is approximately 7 s after the previous one, including an additional 1 s stimulus presentation and an approximate 0.2 s response time. If AEPsych again missed this deadline, the next opportunity came at around 11.1 s. This staggered scheduling ensured that trials continued smoothly while allowing AEPsych sufficient time to compute adaptive placements when possible.

A potential drawback of the fallback trial strategy is that it could disrupt the intended interleaving of adaptive and validation trials, potentially introducing differential learning effects. To mitigate this, we capped how far MOCS trials could advance relative to AEPsych trials. This cap was set at four trials. If this limit was reached and no AEPsych trial was ready, we inserted pregenerated Sobol'-sampled trials instead. These trials were created in advance using participant- and session-specific random seeds and were separate from those selected by AEPsych.

To evaluate how well thresholds read out from the WPPM fits aligned with those estimated via Weibull fits, we simulated a dataset with a known ground truth. The following subsections outline the key steps in this process.

Appendix 5—figure 1

Download asset Open asset

Appendix 5.2: Simulation of trials near threshold contours

To introduce some variability, we added bivariate Gaussian noise to each comparison stimulus at the iso-distance contour in model space. The noise standard deviation was proportional to the Euclidean distance between the reference stimulus ${\mathbf{x}}_0$ and the comparison stimulus ${\mathbf{x}}_1$ in model space. The jittered comparison stimulus ${\mathbf{x}}_1^{\prime }$ was computed as:

(S9) ${\displaystyle {\displaystyle {\mathbf{x}}_1^{\prime }={\mathbf{x}}_1+\mathit{ϵ},\mathit{ϵ}\sim \mathcal{N}\left(\mathbf{0},0.3\cdot |{\mathbf{x}}_1-{\mathbf{x}}_0{|}^2\cdot \mathbf{I}\right).}}$

We modeled performance using a Weibull psychometric function with fixed shape parameters. Specifically, for each pair of reference and jittered comparison stimuli, we computed ΔE and used the Weibull curve to obtain the predicted percent correct:

(S10) ${\displaystyle {\displaystyle \mathrm{\Psi }(\mathrm{\Delta }E)=\gamma +(1-\gamma )\left(1-e^{-(\mathrm{\Delta }E/\alpha {)}^{\beta }}\right)=\frac{1}{3}+\frac{2}{3}\left(1-e^{-(\mathrm{\Delta }E/3.189{)}^{1.505}}\right).}}$

The values of α and β were selected such that the psychometric curve returns 66.7% correct when $\mathrm{\Delta }E=2.5$ (Appendix 5—figure 1C). A binary (correct or incorrect) response was sampled from a Bernoulli distribution using this predicted probability:

(S11) ${\displaystyle {\displaystyle r\sim \mathrm{B}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{o}\mathrm{u}\mathrm{l}\mathrm{l}\mathrm{i}(\mathrm{\Psi }(\mathrm{\Delta }E)).}}$

The comparison stimuli were selected to be near threshold, while the reference stimuli and chromatic directions were Sobol'-sampled to ensure uniform coverage of model space without repeated trials (Appendix 5—figure 1D). In total, we simulated 18,000 trials.

Appendix 5—figure 2

Download asset Open asset

Appendix 5.3: Fitting the WPPM and treating the model fit as ground truth

We fit the WPPM to the 18,000 simulated trials in model space and treated the resulting fit as the ground truth for simulating responses for both AEPsych and validation trials (Appendix 5—figure 1E, color lines). We chose to use the WPPM fit as the ground truth—rather than percent-correct performance derived from CIELAB ΔE94 with the Weibull psychometric function—because our goal here was to evaluate how well the WPPM can recover ground truth that is itself described by the WPPM.

Appendix 5—figure 3

Download asset Open asset

Appendix 5.6: Statistical analysis of residuals between WPPM readouts and validation thresholds

Appendix 5—figure 4

Download asset Open asset

We applied the same statistical analysis to the simulated data as we did for the human data (Appendix 4.2). Consistent with the human results (Appendix 4—figure 9), we found no strong evidence that residuals systematically varied with the orientation or shape of the elliptical threshold contours read out from the WPPM fits. However, we did observe a significant negative correlation between the residuals and the magnitude of the validation thresholds (slope = $-0.393,t(23)=-5.273$, p<0.001, $R^2=0.547$; Appendix 5—figure 4; Appendix 5—table 1). As noted earlier, the size of this bias is small compared to the overall range of validation thresholds.

Appendix 5.7: Comparison between WPPM estimates and simulation ground truth

To evaluate whether the WPPM readouts systematically deviated from the ground truth, we sampled thresholds over a fine grid of reference locations (15 × 15 points evenly spaced between −0.7 and 0.7 in model space) and compared them with the corresponding ground-truth thresholds. As a comparison metric, we used the Bures–Wasserstein (BW) distance (Bhatia et al., 2019), which quantifies the dissimilarity between two positive semi-definite covariance matrices, ${\mathrm{\Sigma }}_1$ and ${\mathrm{\Sigma }}_2$. Intuitively, it captures the ‘effort’ required to morph one ellipse into another. Mathematically, the BW distance is defined as:

(S12) ${\displaystyle {\displaystyle d({\mathrm{\Sigma }}_1,{\mathrm{\Sigma }}_2)=[tr({\mathrm{\Sigma }}_1)+tr({\mathrm{\Sigma }}_2)-2\cdot tr{\left({\mathrm{\Sigma }}_1^{1/2}\cdot {\mathrm{\Sigma }}_2\cdot {\mathrm{\Sigma }}_1^{1/2}\right)}^{1/2}{]}^{1/2}.}}$

The BW distance is non-negative and equals zero only when the two matrices being compared are identical. Smaller distances indicate greater similarity between the threshold ellipses.

Appendix 5—table 1

Predictor	Term	Coef	Std Err	t	p	[0.025, 0.975] CI	R2
Absolute difference of angles	Intercept	–0.007	0.005	–1.532	0.139	[–0.017, 0.003]	0.028
	Slope	0.000	0.000	0.812	0.425	[–0.000, 0.000]
Aspect ratio	Intercept	–0.007	0.012	–0.579	0.568	[–0.033, 0.019]	0.003
	Slope	0.002	0.008	0.263	0.795	[–0.014, 0.018]
Validation thresholds	Intercept	0.028	0.006	4.429	0.000	[0.015, 0.041]	0.547
	Slope	–0.393	0.075	–5.273	0.000	[–0.547, –0.239]

The results showed that BW distance generally increased as the reference color moved farther from the achromatic point (Appendix 5—figure 5A), suggesting that the WPPM has more difficulty accurately capturing large threshold contours in regions with higher internal noise. To provide a benchmark for what constitutes a substantial mismatch, we computed the BW distance between each ground-truth ellipse and a circle with radius equal to the largest major axis length among all ground-truth ellipses. The maximum of these values served as a reference point (shown as the upper limit of the color bar in Appendix 5—figure 5A). Overall, the mismatches observed in our simulations were modest—well below the level expected if the model were fundamentally mischaracterizing the threshold shapes.

We also examined differences in the estimated major axis lengths. The WPPM showed slight underestimation in the upper region and overestimation in the lower region of the space (Appendix 5—figure 5B; also visible, though small, in Appendix 5—figure 2C). Similar to the BW analysis, these deviations were relatively small compared to the overall range of ground-truth values. Together, these results indicate that the WPPM provides a close and robust approximation of the true threshold contours, with only minor local deviations.

Appendix 5—figure 5

Download asset Open asset

In this section, we compare our measurements with those from Danilova and Mollon, 2025, by transforming our results into the chromaticity space used in their study, a scaled version of the MacLeod–Boynton space (MacLeod and Boynton, 1979). While a direct transformation path exists from our model space to theirs (model space → RGB → LMS → MacLeod–Boynton → scaled MacLeod–Boynton), it assumes that the adaptation point and isoluminant plane are identical between the two studies, which is not the case. To account for these differences, we instead took a detour through the DKL space (Derrington et al., 1984), where cone-opponent mechanisms are explicitly defined, and adaptation is more easily controlled. Specifically, we followed the transformation chain: model space → RGB_us → LMS_us → ΔLMS_us → DKL → ΔLMS_dm → LMS_dm → MacLeod–Boynton → scaled MacLeod–Boynton. Here, the subscript ‘us’ refers to values computed using our study’s cone fundamentals, luminosity function, and adaptation point, while ‘dm’ denotes those used in Danilova and Mollon, 2025. This approach allowed us to approximate how our stimuli would be represented in their perceptual framework, enabling visual comparison of the threshold contours. The comparison reveals a general qualitative agreement between their measurements and ours (Appendix 7—figure 1).

Appendix 7—figure 1

Download asset Open asset

We compared our threshold estimates with those reported by Krauskopf and Gegenfurtner, 1992. To do so, we first read out, for each participant, the threshold contour corresponding to 66.7% correct at the achromatic reference color in model space and transformed it to DKL space (Derrington et al., 1984). We then normalized the DKL cardinal axes so that the threshold contour at the achromatic reference had unit length along both axes. This normalized space—referred to here as stretched DKL space—is the coordinate system in which Krauskopf and Gegenfurtner, 1992, conducted their measurements. Finally, we converted the 16 reference stimuli used in their study into our model space, read out the corresponding threshold ellipses, and transformed them into stretched DKL space to enable direct comparison (Appendix 8—figure 1). We observed generally good agreement between their measurements and those from some of our participants, such as CH and FM. Notably, however, individual differences are evident, particularly in the upper-right quadrant of stretched DKL space (Appendix 8—figure 2). In addition, at the adapting chromaticity, our ellipses are consistently rotated relative to the DKL axes, as noted in the main text.

Appendix 8—figure 1

Download asset Open asset

Appendix 8—figure 2

Download asset Open asset

Appendix 9—figure 1

Download asset Open asset

We obtained elliptical thresholds on a grid of reference stimuli in model space from each participant’s WPPM fits (Appendix 2—figure 1). We then transformed these ellipses from model space to L*a*b* space. Specifically, values in model space were first converted to RGB values via an affine transformation (Appendix 1.3). The resulting RGB values were converted to XYZ using the CIE 1931 2° standard-observer color matching functions (ISO/CIE, 2019). Finally, the XYZ values were converted to L*a*b* using an open-source Python package (Taylor, 2017). The adapting background ($R=G=B=0.5$) was used as the reference white in the XYZ-to-L*a*b* transformation. In addition, we computed threshold contours directly in L*a*b* space, defining them as iso-distance contours at a fixed perceptual distance of $\mathrm{\Delta }E=2.5$.

Comparisons revealed that the iso-distance contours from ΔE94 and ΔE00 provided reasonable approximations to our model-predicted thresholds (Appendix 9—figures 2 and 3), with only modest deviations. In contrast, the ΔE76 contours—despite their continued widespread use—diverged substantially from our measurements (Appendix 9—figure 1).

Appendix 9—figure 2

Download asset Open asset

Appendix 9—figure 3

Download asset Open asset

Recall from the main text that the log-likelihood function implied by the WPPM observer model can be written in terms of

which has no simple closed-form solution and must be estimated by Monte Carlo simulation. This quantity of interest takes the form of a cumulative distribution function $g(u)=\mathrm{Pr}[v\le u]$ for some random variable v and scalar constant u. In this section, we describe how to approximate the log-likelihood in a manner that is compatible with automatic differentiation libraries, which enables gradient-based optimization of the log-posterior density.

Let $P_{\theta }$ denote some probability distribution parameterized by θ. Given n independent and identically distributed random variables, $v_1,\dots ,v_n\sim P_{\theta }$, we would like to form an estimate of the cumulative distribution function, $g(u)=\mathrm{Pr}[v\le u]$, where $v\sim P_{\theta }$. A simple and well-known estimate is the empirical cumulative distribution function:

where 1[·] is the indicator function—i.e., 1[A] evaluates to one if the event A occurs and evaluates to zero otherwise. In many respects, this is a perfectly fine estimator. For example, the celebrated Dvoretzky–Kiefer–Wolfowitz inequality (Dvoretzky et al., 1956) states that this estimate converges exponentially fast to the true cumulative distribution function as $n\to \mathrm{\infty }$.

In our setting, we would like to not only evaluate $g(u)$ for any given u, but to also evaluate $\partial g(u)/\partial {\theta }_j$ for all parameters ${\theta }_1,\dots ,{\theta }_J$ that define the underlying distribution $P_{\theta }$. Equation S14 provides an estimate of $g(u)$ but it is unfortunately not differentiable with respect to $v_1,\dots ,v_n$ because 1$\left[v_i\le u\right]$ is a discontinuous step as a function of u. A straightforward and intuitive solution is to replace this step function with a smooth sigmoid function. We formalize this approach below, showing that it can be motivated by forming a smoothed estimate of the underlying density function.

Specifically, let $K(v)$ denote a smooth, non-negative function that integrates to one and satisfies $K(v)=K(-v)$. Suppose that $P_{\theta }$ has a density function $f(v)$. Then, given $v_1,\dots ,v_n\sim P_{\theta }$, we can estimate the density function as

where $h>0$ is a user-specified hyperparameter called the bandwidth. Equation S15 is known as a kernel density estimate (Wasserman, 2006). Asymptotically, $\hat{f}$ approaches the true density function f as $n\to \mathrm{\infty }$ and $h\to 0$. Intuitively, larger values of h lead to smoother density estimates, which is preferable in sample-limited (i.e. small n) regimes.

Now define:

which is a smooth sigmoid function centered at $v_i$, and consider the following estimate of the cumulative distribution function:

Notice that in the limit of $h\to 0$, we recover the empirical cumulative distribution estimator $\hat{g}={\hat{g}}_{\text{emp}}$ because, in this limit, we have that $I_{v_i}(u)=\mathbf{1}\left[v_i\le u\right]$. We can further justify Equation S17 as a reasonable estimator of g by recognizing it as the integral of the density estimate in Equation S15. That is,

Our refined estimator $\hat{g}$ is clearly differentiable whenever we choose $K(v)$ to be a smoothly differentiable function. In our model fitting routine, we chose $K(v)$ to be the density of a standard logistic distribution:

This smoothing kernel has heavy tails, which we reasoned would enable numerically stable auto-differentiation routines even when h is chosen to be small. Another feature is that the integrated density is the well-known logistic function:

which is familiar and easy to compute.