A new full-field digital image correlation method is proposed for identification of three-dimensional displacement field within a soft body. Two sets of three-dimensional image data, which consist of intensity values with discretely sampled positions, are obtained by X-ray CT images for undeformed and deformed states of the body. A displacement field is tentatively represented by B-spline basis functions with unknown coefficients on the overall identification domain. Following the conventional digital image correlation method, the error function of displacement field is defined by the sum of error squares of intensity values on the corresponding material points in the two sets of image data, where the relationship between deformed and undeformed positions of the material point is described by the tentative displacement field. The iterative solution search is employed for modification of the unknown coefficients to identify the actual displacement field. Use is made of the penalty function method with weak constraint of incompressibility to cope with the indefiniteness of solution caused by the intensity flatness and measurement errors of CT images. An efficient parallel algorithm for the solution search is devised by utilizing the locality of B-spline function. The decomposed coefficient matrix for each decomposed domain, derived from the linearized equation of the Levenberg-Marquardt method, is constituted parallel at computational nodes. The modified conjugate gradient method, where the sequence of a vector variable calculation is modified, efficiently solves the linear equation, using the decomposed coefficient matrices. A displacement field within a specimen under compressive load is identified for the validation of the proposed method.