The three-dimensional time-resolved Lagrangian particle tracking (3D TR-LPT) technique has recently advanced flow diagnostics by providing high spatiotemporal resolution measurements under the Lagrangian framework. To fully exploit its potential, accurate and robust data processing algorithms are needed. These algorithms are responsible for reconstructing particle trajectories, velocities, and differential quantities (e.g., pressure gradients, strain- and rotation-rate tensors, and coherent structures) from raw LPT data. In this paper, we propose a 3D divergence-free Lagrangian reconstruction method, where three foundation algorithms—constrained least squares (CLS), stable radial basis function (RBF-QR), and partition-of-unity method (PUM)—are integrated into one comprehensive reconstruction strategy. Our method, named CLS-RBF PUM, is able to (1) directly reconstruct flow fields at scattered data points, avoiding Lagrangian-to-Eulerian data conversions; (2) assimilate the flow diagnostics in Lagrangian and Eulerian descriptions to achieve high-accuracy reconstruction; (3) process large-scale LPT datasets with more than hundreds of thousand particles in two dimensions (2D) or 3D; (4) enable inter-frame and inter-particle interpolation while imposing physical constraints (e.g., divergence-free for incompressible flows) at arbitrary time and locations. Validation based on synthetic and experimental LPT data confirm that our method can achieve the above advantages with accuracy and robustness. Graphical abstract