  Abridged: Debris disks are valuable systems to study dust properties. Because
they are optically thin at all wavelengths, we have direct access to the
properties of dust grains. One very promising technique to study them is to
measure their phase function. Disks that are highly inclined are promising
targets as a wider range of scattering angles can be probed. The phase function
is usually either inferred by comparing the observations to synthetic disk
models assuming a parametrized phase function, or estimating it from the
surface brightness of the disk. We argue here that the latter approach can be
biased due to projection effects leading to an increase in column density along
the major axis of a non flat disk. We present a novel approach to account for
those column density effects. The method remains model dependent, as one still
requires a disk model to estimate the density variations as a function of the
scattering angle. This method allows us however to estimate the shape of the
phase function without having to invoke any parametrized form. We apply our
method to SPHERE/ZIMPOL observations of HR\,4796 and highlight the differences
with previous measurements. Our modelling results suggest that the disk is not
vertically flat at optical wavelengths. We discuss some of the caveats of the
approach, mostly that our method remains blind to real local increase of the
dust density, and that it cannot yet be readily applied to angular differential
imaging observations. Similarly to previous studies on HR\,4796, we still
cannot reconcile the full picture using a given scattering theory to explain
the shape of the phase function, a long lasting problem for debris disks.
Nonetheless, we argue that similar effects as the ones highlighted in this
study can also bias the determination of the phase function in total intensity.
