Accurate and efficient finite element (FE) analysis of problems with strongly heterogeneous material distributions depends on the quadrature scheme used to compute element-level quantities. When describing material variations through a voxelized geometry, for instance obtained from computed tomography (CT), resolving this information using standard Gauss quadrature can lead to significant integration errors, particularly when multiple voxels with a high contrast are embedded within a single FE. Although distributing quadrature points on each voxel is possible, this approach quickly becomes computationally unfeasible. Pre-integration is a common remedy that exploits the piece-wise constant nature of voxel-based materials by pre-computing element stiffness matrices for each voxel. However, the cost of this approach scales with the number of voxels per element and it does not directly transfer to hierarchical mesh refinement algorithms. To address this challenge, we propose an efficient and discretization-independent voxel integration scheme based on a moment-fitting approach, as originally presented in. The method constructs a custom quadrature rule for each finite element, including voxel-level material information directly into the integration weights. This allows exact integration of the element stiffness matrix for polynomial FE bases while requiring significantly fewer quadrature points than other voxel-resolving approaches, thus decreasing the effective runtime for the integration routine. The proposed scheme is compatible with refined FE meshes and accommodates heterogeneous material distributions. This talk explains the methodology, presents a comparison with existing quadrature approaches, and showcases the strength with a numerical example on real CT scan data.