Gaussian process (GP) regression is widely used for flexible and non-parametric Bayesian modeling of data arising from underlying smooth functions. This paper introduces a solution to GP regression when the observations are subject to value-based censoring. We derive exact and closed-form expressions for the conditional posterior distributions of the underlying functions in both the single-curve fitting case and in the case of a hierarchical model where multiple functions are modeled simultaneously. Our method can accommodate left, right, and interval censoring, and is directly applicable as an empirical Bayes method or integrated in a Markov–Chain Monte Carlo sampler for full posterior inference. The method is validated through extensive simulations, where it substantially outperforms naive approaches that either exclude censored observations or treat them as fully observed values. We give an application to a real-world dataset of longitudinal HIV-1 RNA measurements, where the observations are subject to left censoring due to a detection limit.