Horizon picking is a fundamental and crucial step for seismic interpretation, but it remains a time-consuming task. Although various automatic methods have been developed to extract horizons in seismic images, most of them may fail to pick horizons across discontinuities such as faults and noise. To obtain more accurate horizons, we have developed a dynamic programming algorithm to efficiently refine manually or automatically extracted horizons so that they can more accurately track reflectors across discontinuities, follow consistent phases, and reveal more geologic details. In this method, we first compute an initial horizon using an automatic method, manual picking, or interpolation with several control points. The initial horizon may not be accurate, and it only needs to follow the general trend of the target horizon. Then, we extract a subvolume of amplitudes centered at the initial horizon and meanwhile flatten the subvolume according to the initial horizon. Finally, we use the dynamic programming to efficiently pick the globally optimal path that passes through the global maximum or minimum amplitudes in the subvolume. As a result, we are able to refine the initial horizon to a more accurate horizon that follows consistent amplitude peaks, troughs, or zero crossings. Because our method does not strictly depend on the initial horizon, we prefer to directly interpolate an initial horizon from a limited number of control points, which is computationally more efficient than automatically or manually picking an initial horizon. In addition, our method is convenient to interactively implement to update the horizon while editing or moving the control points. More importantly, these control points are not required to be exactly placed on the target horizon, which makes the human interaction highly convenient and efficient. We develop our method with multiple 2D and 3D field examples that are complicated by noise, faults, and salt bodies.