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.