Plant area density (PAD in m2 center dot m- 3) defines the total one-sided total plant surface area within a given volume. It is a key variable in characterizing exchange processes between the atmosphere and land surface. Terrestrial laser scanning (TLS) provides unprecedented detail of the 3D structure of forest canopies. Yet, signal occlusion and uneven sampling density of the TLS point clouds limit our capacity to characterize the 3D distribution of canopy components. Recent studies have made use of statistical estimators of PAD that are applied to TLS point clouds subdivided into three-dimensional (3D) cubes, or voxels. Computation of such metrics under actual field conditions with point clouds containing several millions of returns is challenging. Moreover, rigorous assessment of the estimated PAD and effects of occlusions in forests remain unclear due to laborious, time-consuming, and inaccurate field measurements. In the present study, we present L-Vox, a software that computes PAD per voxel for TLS scans acquired in forest environments, which is based upon recent development of unbiased estimators derived from maximum likelihood. Two applications are presented. First, the software is evaluated for virtual forest plots, which are detailed 3D models of individual trees with corresponding simulated TLS scans, for which reference data are known. Second, L-Vox is applied to actual scans that were acquired in hardwood and coniferous plots in New Brunswick and Newfoundland, Canada. Both test cases were used to investigate the effects of occlusion and the uneven sampling in estimating PAD. The test cases were also used to assess the influence of voxel size and the number of scans per plot on PAD estimates. Our results showed strong correlations between the estimated PAD profile from L-Vox and simulated PAD for virtual forest plots, with a mean R2 = 0.98 and a mean coefficient of variation (CV) = 15.6%. We demonstrated that comparing multi-scan to single scan TLS acquisitions in real forest plots substantially reduced signal occlusion, resulting in an increase up to 50% in PAD values. Effects of voxel size on PAD estimates greatly depended upon the relative size of foliar and woody elements, with an optimal size around 10 cm in coniferous plots. L-Vox proved to be an efficient and accurate tool for computing 3D distributions of PAD from TLS measurements in natural forest environments.