I'll answer my own question and hope this information is useful to someone. I'll take $\hbar = 1$ and will deal with systems of one degree of freedom (the generalisation should be obvious). The normalisation factors are very confusing, so I'll omit most of the reasons for them to be what they are.
First, let's take a look at the good old continuous Wigner function. There are many ways to define it. I'll describe two.
Wigner function as coefficients in an expansion in terms of Weyl-Heisenberg operators
Consider the operator $\hat{T}_{(a,b)} = e^{-i(a\hat{p} - b\hat{q})}$. This operator is the generator of translations in position by $b$ and in momentum by $a$. We are, thus, dealing with a formalism that puts positions and momenta on the same ground (we say we're dealing with quantum mechanics in phase space). Associated to the translation operators we can define another operator as its symplectic Fourier transform, i.e, $\hat{R}_{(q,p)} = \frac{1}{2\pi} \int e^{-i(ap - bq)} \hat{T}_{(a,b)}$. It can be shown that these operators form a representation of the group of reflections and translations in a 2-dimensional phase space. We can also show that these operators are orthonormal with respect to the Hilbert-Schmidt metric,
$$ Tr \left(\hat{T}^\dagger_{(a,b)} \hat{T}_{(a',b')} \right) = 2 \pi \delta(a' - a) \delta(b'-b) \\
Tr \left(\hat{R}^\dagger_{(q,p)} \hat{R}_{(q',p')} \right) = 2 \pi \delta(q' - q) \delta(p'-p) .$$
This means that we can expand any operator acting in our Hilbert space in a reflection (or translation basis). In the expansion of the density matrix as a continuos linear combination of reflections,
$$\hat{\rho} = \frac{1}{2\pi} \int W(q,p) \hat{R}_{(q,p)} dq dp ,$$
the coefficients are defined as the Wigner function associated to the system. Such equation can be inverted to explicitly give us
$$W(q,p) = Tr(\hat{\rho} \hat{R}) .$$
Wigner function from projection properties
Say we create a function from the density matrix
$$W(q,p) = Tr(\hat{\rho} \hat{R}) ,$$
and that we wish this operator to have the following property: the integral of $W$ over the strip of phase space bounded by the parallel lines $aq + bp = c_1$ and $aq + bp = c_2$ is the probability that the operator $a\hat{q} + b\hat{p}$ will take a value between $c_1$ and $c_2$. This, by itself, it enough to define the operator $\hat{R}$ and fix the form of the Wigner function (you only need to invert a Radon transform to get it).
Now, let's go back to discrete phase space. After reading a lot of articles on the subject I noticed what they do is to define the Wigner function either as coefficients, or following the projection property. But, since they're trying to define a quasi-probability distribution to a system with discrete degrees of freedom (mostly spin), the "phase space" needs to be discrete. Actually more than discrete, it needs to be a torus, because the system has only a finite number of accessible states (contrary to the continuous case, where there were infinite accessible states). The thing here is that what you need to do is to either create reflection operators on this discrete set or associate an eigenstate of some observable with lines in phase space and use it as a guide to create a Wigner function-like object, to mimic the continuous case. Nevertheless, you are trying to associate a net of points with states. This is very different from the continuous case, where you where starting from a phase space in the sense of classical mechanics ($q$ and $p$ and etc).
So, my conclusion is the following: the "discrete phase spaces" that appear in the Wigner function formulation of finite state quantum mechanics have nothing to do with quantum mechanics in phase space. These strange, non-unique Wigner functions are nothing more than an efficient way of performing state tomography. The question regarding the possibility of discretisation of phase space lost their meaning, and should be forgotten. I'll, nevertheless, still focus in the problem of compactifying cotangent fibres.
Now, two points of view were very interesting. The first one deals with really creating a group-quotient structure in a classical phase space projecting the Wigner representation from the plane to the torus. By a main result in the spectral theory of compact operators, compactifying the domain of the elements of the Banach space the operator is defined upon results in the countability of its spectrum. This means that projecting the Wigner function from the plane to any compact domain (the torus is the most important case) will cause the operators that act in the Hilbert space to discretise their spectra. This means that we can in fact find the discrete analog of the Wigner function through phase space formalism, but this only makes sense when we're measuring systems where position and momentum are discrete. Spins and angular momenta still need the artificial Wigner-like functions commented above. (The article that develops this is Annals of Physics 276, 223-256 (1999).)
The second point of view that was really interesting deals with problems concerning the dimension of our state-space. Since the dimension is finite, we need a finite field in order to make mathematical sense out of, for example, foliations of finite-dimensional nets (if we were not on a field, two non parallel lines would intersect in more than one point). Now, finite-fields are only possible when the quotient set is taken to be $mod p$, where $p$ is a prime. In order to deal with systems whose degrees of freedom do not follow the prime condition, the authors read a little about Galois theory and added points inside the state space without changing the modular group basis. This was one of the cutest applications of Galois theory in applied physics I've ever seen. (The article that deals with this is Physical Review A 70, 062101 (2004).)
I'm sorry for the long answer. I hope it's useful to someone.
