We develop a novel stochastic derivative estimation framework for sample performance functions that are discontinuous in the parameter of interest, based on the multidimensional Leibniz integral rule. When discontinuities arise from indicator functions, we embed the indicator functions into the sample space, yielding a continuous performance function over a parameter-dependent domain. Applying the Leibniz integral rule in this case produces a single-run, unbiased derivative estimator. For general discontinuous functions, we apply a change of variables to shift parameter dependence into the sample space and the underlying probability measure. Applying the Leibniz integral rule leads to two terms: a standard likelihood ratio (LR) term from differentiating the underlying probability measure and a surface integral from differentiating the boundary of the domain. Evaluating the surface integral may require simulating multiple sample paths. Our proposed Leibniz integration framework generalizes the generalized LR (GLR) method and provides intuition as to when the surface integral vanishes, thereby enabling single-run, easily implementable estimators. Numerical experiments demonstrate the effectiveness and robustness of our methods.