In [Turek, S.: A Generalized Mean Intensity Approach for the Numerical Solution of the Radiative Transfer Equation, Preprint, 1994], the concept of the generalized mean intensity has been proposed as a special numerical approach to the (linear) radiative transfer equations which can result in a significant reduction of the dimension of the discretized system, without eliminating any information for the specific intensities. Moreover, in combination with Krylov-space methods (CG, Bi-CGSTAB, etc.), robust and very efficient solvers as extensions of the classical approximate Lambda-iteration have been developed. In this paper, the key tool is the combination of special renumbering techniques together with finite difference-like discretization strategies for the arising transport operators which are based on short-characteristic upwinding techniques of variable order and which can be applied to highly unstructured meshes with locally varying mesh widths. We demonstrate how such special upwinding schemes can be constructed of first order, and particularly of second order accuracy, always leading to lower triangular system matrices. As a consequence, the global matrix assembling can be avoided (`on-the-fly`), so that the storage cost are almost optimal, and the solution of the corresponding convection-reaction subproblems for each direction can be obtained very efficiently. As a further consequence, this approach results in a direct solver in the case of no scattering, while in the case of non-vanishing absorption and scattering coefficients the resulting convergence rates for the full systems depend only on their ratio and the absolute size of these physical quantities, but not on the grid size or mesh topology. We demonstrate these results via prototypical configurations and we examine the resulting accuracy and efficiency for different computational domains, meshes and problem parameters.