Most seismic horizon extraction methods are based on seismic local reflection slopes that locally follow seismic structural features. However, these methods often fail to correctly track horizons across discontinuities such as faults and noise because the local slopes can only correctly follow laterally continuous reflections. In addition, seismic amplitude or phase information is not used in these methods to compute horizons that follow a consistent phase (e.g., peaks or troughs). To solve these problems, we have developed a novel method to compute horizons that globally fit the local slopes and multigrid correlations of seismic traces. In this method, we first estimate local reflection slopes by using structure tensors and compute laterally multigrid slopes by using dynamic time warping (DTW) to correlate seismic traces within multiple laterally coarse grids. These coarse-grid slopes can correctly correlate reflections that may be significantly dislocated by faults or other discontinuous structures. Then, we compute a horizon by fitting, in the least-squares sense, the slopes of the horizon with the local reflection slopes and multigrid slopes or correlations computed by DTW. In this least-squares system, the local slopes on the fine grid and the multiple coarse-grid slopes will fit a consistent horizon in areas without lateral discontinuities. Across laterally discontinuous areas where the local slopes fail to correctly correlate reflections and mislead the horizon extraction, the coarse-grid slopes will help to find the corresponding reflections and correct the horizon extraction. In addition, the multigrid correlations or slopes computed by dynamic warping can also assist in computing phase-consistent horizons. We apply the proposed horizon extraction method to multiple 2D and 3D examples and obtain accurate horizons that follow consistent phases and correctly track reflections across faults.