Discontinuities in data often provide vital information, and representing these discontinuities sparsely is an important goal for approximation and compression algorithms. Little work has been done on efficient representations for higher dimensional functions containing arbitrarily smooth discontinuities. We consider the N-dimensional Horizon class -- N-dimensional functions containing a C^K smooth (N-1)-dimensional singularity separating two constant regions. We derive the optimal rate-distortion function for this class and introduce the multiscale surflet representation for sparse piecewise approximation of these functions. We propose a compression algorithm using surflets that achieves the optimal asymptotic rate-distortion performance for Horizon functions. This algorithm can be implemented using knowledge of only the N-dimensional function, without explicitly estimating the (N-1)-dimensional discontinuity.