This paper presents a regularization scheme for the nearly singular integrals used for 3D elastostatic boundary element analysis. For the regularization process, the local projection coordinates of the source point are first located via an iteration procedure. For planar elements, the boundary integrals are analytically integrated by parts to smooth the drastic fluctuations of their integrands so that the regularized forms can be numerically integrated by any conventional schemes in an usual manner. The validity of the formulations is numerically tested using the Gauss Quadrature scheme. The test shows the accuracy is satisfactory for the distance ratio (distance: Element characteristic length) falling below micro-scale. To further demonstrate our successful implementation, a numerical example is studied with verifications compared with ANSYS analysis.