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 propose 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 only needs to follow the general trend of the target horizon. Then, we extract a sub-volume of amplitudes centered at the initial horizon and meanwhile flatten the sub-volume according to the initial horizon. We finally use the dynamic programming to efficiently pick the globally optimal path that passes through global maximum or minimum amplitudes in the sub-volume. 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. As 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 be interactively implemented 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 demonstrate our method with multiple 2D and 3D field examples that are complicated by noise, faults, and salt bodies.