A two-point flux like this is not convergent if the mesh is not "orthogonal", in the sense that the edge/face between two cells is orthogonal to the line segment joining the cell centroids. If your mesh is orthogonal, you would use the distance between centroids for $h_{ij}$ above.
If you would like a method to work on more general meshes within the cell-centered finite volume framework, you have to do a bit of extra work. A classical approach is to reconstruct inside cells, as outline in these notes. This is not very robust, especially in the presence of irregular coefficients.
Another alternative is to consider the mixed finite element method using a BDM-1 space (constant pressure, linear velocity on faces/edges) and choose a reduced quadrature consisting of one point per vertex. That special quadrature causes the velocity mass matrix to be block diagonal, with one block per vertex, so the velocity degrees of freedom can be eliminated, returning you to a cell-centered non-mixed formulation. This is illustrated in Figure 2.3 from Wheeler and Yotov (2006).
After velocity is eliminated, the cell-centered operator is SPD and cells are coupled only to other cells they share a vertex with, which is sparser than the gradient reconstruction method mentioned earlier. With this method, pressure is second order superconvergent and cell centers and velocity is first order convergent, even on irregular meshes. If the mesh is within $O(h^2)$ of affine, velocity is also second order superconvergent at face centers. This method naturally handles general tensor-valued coefficients. This method is a multipoint flux approximation.
It is known (see page 10 of Edwards) that 9-point schemes for regular quadrilateral meshes cannot be monotone for general tensor-valued coefficients. One approach to maintain monotonicity is to use a nonlinear discretization.
