Abstract. A systematic and comprehensive framework for finite impulse response (FIR) lowpass/fullband derivative kernels is introduced in this paper. Closed form solutions of a number of derivative filters are obtained using the maximally flat technique to regulate the Fourier response of undetermined coefficients. The framework includes arbitrary parameter control methods that afford solutions for numerous differential orders, variable polynomial accuracy, centralized/staggered schemes, and arbitrary side-shift nodes for boundary formulation. Using the proposed framework four different derivative matrix operators are introduced and their numerical stability is analyzed by studying their eigenvalues distribution in the complex plane. Their utility is studied by considering two important image processing problems, namely gradient surface reconstruction and image stitching. Experimentation indicates that the new derivative matrices not only outperform commonly used method but provide useful insights to the numerical issues in these two applications.Key words. numerical differentiation, boundary condition, derivative matrix, fullband/lowpass FIR differentiation, maxflat technique, inverse Vandermonde, gradient surface recovery, image stitching 1. Introduction. Finite difference (FD) methods are used extensively in variety of image processing tasks. For example, they are used in digital approximation of order moments such as gradients, Hessian, and high-order tensor. They find numerous applications in either forward imaging problems to encode certain features such as edge and corner [11,44,60,61,64,69,78], or discretizing the differential operators used in inverse problems for image restoration [1,8,23,33,59,67,76,77]. Despite extensive applications, the problem of numerical differentiation suffers from several deficiencies such as the lack of estimation accuracy, sensitivity over perturbing artifacts, and improper formulation of boundary condition (BC). Such deficiencies can easily make the overall imaging applications highly complicated (ill-posed) and sometimes turn the whole exercise to be useless. Recent works highlight the need for a systematic and comprehensive framework of numerical differentiation to overcome these issues by reducing the complexities of error propagation in the system design and furthermore prevent computational burdens. Applications can be found in gradient surface reconstruction [30,33], edge detection [2,3,5,10], edge synthesis [22,35,81], feature extraction [16,17,30], and many more. This article concerns with the main question of what is the best numerical approach to approximate derivatives in an image processing related task? With this question in mind we present a comprehensive framework for differentiation applied in forward and inverse imaging problems. The "forward" here is referred to find a solution of direct approximation in a measurement system A(X) = Y where A is an operator and X is a given signal/image to find Y . While, the "inverse" is referred to inferring X from a simi...