Traditional median smoother for 2D images is insensitive to impulse noise but generates flat areas as unwanted artifacts. The proposed approach to overcome this issue is based on the minimization of the regularized form of total variation functional. At first the continuous functional is defined for n-dimensional signal in an integral form with regularization term. The continuous functional is converted to the discrete form using equidistant spatial sampling in point grid of pixels, voxels or other elements. This approach is suitable for traditional signal and image processing. The total variance is then converted to the sum of absolute intensity differences as a minimization criterion. The functional convexity guarantees the existence of global minimum and absence of local extremes. Resulting non-linear filter iteratively calculates local medians using red-black method of Successive Over/Under Relaxation (SOR) scheme. The optimal value of the relaxation parameter is also subject to our study. The sensitivity to regularization parameter enables to design high-pass and nonlinear band-pass filters as the difference between the image and low-pass smoother or as the difference between two different low-pass smoothers, respectively. Various median based approaches are compared in the paper.