When biological specimens are cut into physical sections for three-dimensional (3D) imaging by confocal laser scanning microscopy, the slices may get distorted or ruptured. For subsequent 3D reconstruction, images from different physical sections need to be spatially aligned by optimization of a function composed of a data fidelity term evaluating similarity between the reference and target images, and a regularization term enforcing transformation smoothness. A regularization term evaluating the total variation (TV), which enables the registration algorithm to account for discontinuities in slice deformation (ruptures), while enforcing smoothness on continuously deformed regions, was proposed previously. The function with TV regularization was optimized using a graph-cut (GC) based iterative solution. However, GC may generate visible registration artifacts, which impair the 3D reconstruction. We present an alternative, multilabel TV optimization algorithm, which in the examined samples prevents the artifacts produced by GC. The algorithm is slower than GC but can be sped up several times when implemented in a multiprocessor computing environment. For image pairs with uneven brightness distribution, we introduce a reformulation of the TV-based registration, in which intensity-based data terms are replaced by comparison of salient features in the reference and target images quantified by local image entropies.