This thesis establishes a mathematical toolbox for numerical difference (ND) approximation that can be applied in numerous imaging problems such as restoration and feature detection. The proposed framework is based on integrating the interpolating power polynomials with maximally-linear-flat design methods to regulate the desired ND transfer functions in Fourier domain. Variety of filters are extracted that are capable of addressing different orders of differentiation, accuracies, and cut-off frequency parametrization. The numerical stability of the designed coefficients are highly robust through different filter orders. The second phase of this thesis studies a wide range of imaging applications that utilize the proposed ND filters to avoid unnecessary complexities in practical designs. One of the main applications is studied over the total variation (TV) regularization, which employs functional derivatives to regulate problems in image restoration. We propose to approximate the TV norm space by means of the latter ND filters. Proper mathematical foundations are presented to understand the performance of the TV models in different signal transitions. In particular, we focus on de-noising and compressed video sensing applications and illustrate with a number of experiments to study the impact of alternative ND methods. A proper convex minimization algorithm is proposed to utilize the latter filters to recover images from their degraded measurements. Furthermore, the problem of optical flow in computer vision is considered in this thesis. A comparative study is represented over the well-known taxonomy that links lowpass ND design to the gradient constancy constraints and avoid pre-processing steps for image smoothing. Curvature estimation and edge detection problems in imaging are also studied to modify the pertinent feature estimation from noisy sampling. The experimental evaluations on selected imaging problems outperform the existing state-of-the-art methods.