1 Introduction

For shape analysis in medicine and biology, outline shapes are fundamental feature for recognition and retrieval in information filtering from large amount of data. Therefore, for fast recognition and retrieval, we are required to extract outline shapes from volumetric medical data. Since the nature of medical and biological data is volumetric, these data are expressed in multiway data array. Using multiway data analysis [1], we extract compressed outline shapes of objects from global texture of volumetric data. For this extraction, we applied three-way tensor principal component analysis [2] to voxel images, although in traditional data analysis, three-way data are embedded in vector space. A small number of major principal components represent the shape of the objects in a voxel image. Applying tensor subspace method (TSM) [3] to these major principal components, we construct tensor-based classification. In the classification, the TSM measures the similarity between a query and a tensor subspace spanned by principal axes of a category.

For numerical computation, we deal with sampled patterns. In traditional pattern recognition, these sampled patterns are embedded in an appropriate-dimensional Euclidean space as vectors. The other way is to deal with sampled patterns as higher-dimensional array data. These array data are expressed by tensor to preserve multilinearity of function in the original pattern space. Tensors allow expressing multidimensional array data in multilinear forms.

For the analysis of three-way array, three-mode factor analysis has been proposed [1] in statistics. This method called Tucker decomposition. The Tucker decomposition with orthogonal constraint is equivalent to the third-order tensor principal component analysis (TPCA) [2]. Applying this three-way TPCA to voxel images, we can directly extract the compressed outline shape of objects from volumetric data. In the decomposition procedure of the TPCA, we can extract multilinear structure of each category from data. Using these structures of each category, we can construct subspace learning method.

For the extraction of the outline of shapes in two-dimensional images, active contour model [4] and active shape model [5] has been proposed. The active contour model, which is known as snakes, extracts the outline of a shape by an energy-minimising spline without specific models for categories. The active shape model extracts the shape of objects using the specific models, which are obtained from learning data of categories. As the extensions of these shape extraction models to points cloud in three-dimensional space, statistical models are proposed [6, 7]. In two- and three-dimensional images, these models rely on the set of points, which represent each boundary of objects on images. These points are manually and semi-automatically extracted in advance as preprocessing, while our method extracts an outline shape without the extraction of points.

Furthermore, we show that the compression of volumetric data by the three-way TPCA can be approximated by the reduction based on the three-dimensional discrete cosine transform (3D-DCT).

2 Extraction Based on Tensor Form

2.1 Tensor Representation for N-way Arrays

A Nth-order tensor \(\mathcal {X}\) defined in \(\mathbb {R}^{I_1 \times I_2 \times \dots I_N}\) is expressed as

$$\begin{aligned} \mathcal {X} = ({x}_{i_1,i_2, \dots i_N}) \end{aligned}$$
(1)

for \({x}_{i_1,i_2, \dots i_N}\in \mathbb {R}\) by Nindices \(i_n\). Each subscript n denotes the n-mode of \(\mathcal {X}\). For the outer products of N vectors, if the tensor \(\mathcal {X}\) satisfies the condition

$$\begin{aligned} \mathcal {X} =\varvec{u}^{(1)} \circ \varvec{u}^{(2)} \circ \dots \circ \varvec{u}^{(N)}, \end{aligned}$$
(2)

where \(\circ \) denotes the outer product, we call this tensor \(\mathcal {X}\) a rank-one tensor. For \(\mathcal {X}\), the n-mode vectors, \(n=1,2,\dots ,N\), are defined as the \(I_n\)-dimensional vectors obtained from \(\mathcal {X}\) by varying this index \(i_n\) while fixing all the other indices. The unfolding of \(\mathcal {X}\) along the n-mode vectors of \(\mathcal {X}\) is defined by

$$\begin{aligned} \mathcal {X}_{(n)} \in \mathbb {R}^{I_n \times I_{n'}}, \end{aligned}$$
(3)

where \(I_{n'} = I_1 \times I_2 \times \dots I_{n-1} \times I_{n+1} \times \dots \times I_N\), and the column vectors of \(\mathcal {X}_{(n)}\) are the n-mode vectors of \(\mathcal {X}\). The n-mode product \(\mathcal {X}\times _n \varvec{U}\) of a matrix \(\varvec{U} \in \mathbb {R}^{J_n\times I_n}\) and a tensor \(\mathcal {X}\) is a tensor \(\mathcal {G} \in \mathbb {R}^{I_1 \times I_2 \times \dots \times I_{n-1} \times J_n \times I_{n+1}\times \dots \times I_N}\) with elements

$$\begin{aligned} {\begin{matrix} g_{i_1,i_2,\dots , i_{n-1},j_n,i_{n+1},\dots ,i_N} = \sum _{i_n=1}^{I_n} x_{i_1,i_2,\dots , I_N} u_{j_n,i_n}, \end{matrix}} \end{aligned}$$
(4)

by the manner in ref. [8]. For the m- and n-mode product by matrices \(\varvec{U}\) and \(\varvec{V}\), we have

$$\begin{aligned} \mathcal {X}\times _n \varvec{U} \times _m \varvec{V} = \mathcal {X} \times _m \varvec{V} \times _n \varvec{U} \end{aligned}$$
(5)

since n-mode projections are commutative [8]. We define the inner product of two tensors \(\mathcal {X}=(x_{i_1,i_2,\dots ,i_N})\) and \(\mathcal {Y}=(y_{i_1,i_2,\dots ,i_N})\) in \(\mathbb {R}^{I_1 \times I_2 \times \dots \times I_N}\) by

$$\begin{aligned} \langle \mathcal {X}, \mathcal {Y} \rangle \ \ \ = \sum _{i_1} \sum _{i_2} \dots \sum _{i_N} {x}_{i_1,i_2,\dots ,i_N} {y}_{i_1,i_2,\dots ,i_N}. \end{aligned}$$
(6)

Using this inner product, the Frobenius norm of a tensor \(\mathcal {X}\) is

$$\begin{aligned} \Vert \mathcal {X} \Vert _{\mathrm {F}} = \sqrt{\langle \mathcal {X}, \mathcal {X} \rangle } = \Vert \mathcal {X} \Vert _{\mathrm {F}} = \Vert \mathrm {vec} \ \mathcal {X} \Vert _2, \end{aligned}$$
(7)

where vec and \(\Vert \cdot \Vert _2\) are the vectorisation operator and Euclidean norm of a tensor, respectively. For the two tensors \(\mathcal {X}_1\) and \(\mathcal {X}_2\), we define the distance between them as

$$\begin{aligned} d( \mathcal {X}_1, \mathcal {X}_2 ) = \Vert \mathcal {X}_1 - \mathcal {X}_2\Vert _{\mathrm {F}}. \end{aligned}$$
(8)

Although this definition is a tensor-based measure, this distance is equivalent to the Euclidean distance between the vectorised tensors \(\mathcal {X}_1\) and \(\mathcal {X}_2\) from Eq. (7).

2.2 Projection to Tensor Subspace

As the tensor \(\mathcal {X}\) is in the tensor space \(\mathbb {R}^{I_1} \otimes \mathbb {R}^{I_2} \otimes \dots \otimes \mathbb {R}^{I_N}\), the tensor space can be interpreted as the Kronecker product of N vector spaces \(\mathbb {R}^{I_1},\mathbb {R}^{I_2}, \dots , \mathbb {R}^{I_N}\). To project \(\mathcal {X} \in \mathbb {R}^{I_1} \otimes \mathbb {R}^{I_2} \otimes \dots \otimes \mathbb {R}^{I_N}\) to another tensor \(\mathcal {Y}\) in a lower-dimensional tensor space \(\mathbb {R}^{P_1} \otimes \mathbb {R}^{P_2} \otimes \dots \otimes \mathbb {R}^{P_N}\), where \(P_n \le I_n \) for \(n=1,2, \dots , N\), we need N orthogonal matrices \(\{ \varvec{U}^{(n)} \in \mathbb {R}^{I_n \times P_n} \}_{n=1}^N\). Using the N projection matrices, the tensor-to-tensor projection (TTP) is given by

$$\begin{aligned} \hat{\mathcal {X}} = \mathcal {X} \times _{1} \varvec{U}^{(1) \top } \times _{2} \varvec{U}^{(2) \top } \dots \times _{N} \varvec{U}^{(N) \top }. \end{aligned}$$
(9)

This projection is established in N steps, where at the nth step, each n-mode vector is projected to a \(P_n\)-dimensional space by \(\varvec{U}^{(n)}\). The reconstruction from a projected tensor \(\hat{\mathcal {X}}\) is achieved by

$$\begin{aligned} {\mathcal {X}} = \hat{\mathcal {X}} \times _1 \varvec{U}^{(1)} \times _2 \varvec{U}^{(2)} \dots \times _N \varvec{U}^{(N)}. \end{aligned}$$
(10)

2.3 Principal Component Analysis for Third-Order Tensors

A third-order tensor \(\mathcal {X}\in \mathbb {R}^{I_1 \times I_2 \times I_3}\), which is the array \(\varvec{X} \in \mathbb {R}^{I_1 \times I_2 \times I_3} \), is denoted as a triplet of indices \((i_1, i_2, i_3)\). We set the identity matrices \(\varvec{I}_j, \ j=1, 2, 3\) in \(\mathbb {R}^{I_j\times I_j}\). Here we summarise higher-order singular value decomposition (HOSVD) [9] for third-order tensors. For a collection of tensors \(\{\mathcal {X}_i\}_{i=1}^N \in \mathbb {R}^{I_1 \times I_2 \times I_3}\) satisfying the zero expectation condition \(\mathrm {E} (\mathcal {X}_i )\) \(= 0\), we compute

$$\begin{aligned} \hat{\mathcal {X}_i} = \mathcal {X}_i \times _1 \varvec{U}^{(1)\top } \times _2 \varvec{U}^{(2)\top } \times _3 \varvec{U}^{(3)\top }, \end{aligned}$$
(11)

where \(\varvec{U}^{(j)}=[\varvec{u}^{(j)}_1,\dots ,\varvec{u}^{(j)}_{I_j}]\), that minimises the criterion

$$\begin{aligned} J_- = \mathrm {E} \left( \Vert \mathcal {X}_i - \hat{\mathcal {X}}_i \times _1 \varvec{U}^{(1)} \times _2 \varvec{U}^{(2)} \times _3 \varvec{U}^{(3)} \Vert _{\mathrm {F}}^2 \right) \end{aligned}$$
(12)

and maximises the criterion

$$\begin{aligned} J_{+} = E(\Vert \hat{\varvec{X}}_i \Vert _{\mathrm {F}}^2 ), \end{aligned}$$
(13)

with respect to the conditions \(\varvec{U}^{(j)\top }\varvec{U}^{(j)} =\varvec{I}_j\). By fixing \(\{ \varvec{U}^{j} \}_{j=1}^N\) except \(\varvec{U}^{(j')}\), \(j'=\{ 1, 2, 3 \}\), we have

$$\begin{aligned} J_j = \mathrm {E} \left( \Vert \varvec{U}^{(j)\top } \mathcal {X}_{i,(j)} \mathcal {X}_{i,(j)}^{\top } \varvec{U}^{(j)} \Vert ^2_{\mathrm {F}} \right) . \end{aligned}$$
(14)

Eigendecomposition problems are derived by computing the extremes of

$$\begin{aligned} E_j= & {} J_j+ tr((\varvec{I}_j - \varvec{U}^{(j)\top }\varvec{U}^{(j)} ) \varvec{\varSigma }^{(j)}), \, j=1, 2, 3. \end{aligned}$$
(15)

For matrices \(\varvec{M}^{(j)} = \frac{1}{N} \sum _{i=1}^N \mathcal {X}_{i,(j)}\mathcal {X}_{i,(j)}^{\top }\), \(j=1,2,3\) of \(\mathrm {rank} \, \varvec{M}^{(j)}=K\), the optimisation of \(J_-\) derives the eigenvalue decomposition

$$\begin{aligned} \varvec{M}^{(j)} \varvec{U}^{(j)}=\varvec{U}^{(j)}\varvec{\varSigma }^{(j)}, \end{aligned}$$
(16)

where \(\varvec{\varSigma }^{(j)} \in \mathbb {R}^{I_j \times I_j}, \ j=1, 2, 3\), are diagonal matrices satisfying the relationships \(\lambda ^{(j)}_k =\lambda ^{(j')}_k, \ k \in \{ 1, 2, \dots , K\}\) for

$$\begin{aligned} \varvec{\varSigma }^{(j)}= & {} \mathrm {diag} (\lambda ^{(j)}_1,\lambda ^{(j)}_2\cdots , \lambda ^{(j)}_K,0\cdots ,0). \end{aligned}$$
(17)

For the optimisation of \( \{ J_j \}_{j=1}^3\), there is no closed-form solution to this maximisation problem [9]. Algorithm 1 is the iterative procedure of multilinear TPCA [2]. We adopt Algorithm 1 [2] to optimise \( \{ J_j \}_{j=1}^3\). For \(\varvec{p}_{k} \in \{ \varvec{e}_k\}_{k=1}^K\), we set orthogonal projection matrices \(\varvec{P}^{(j)} = \sum _{k=1}^{k_j} \varvec{p}_{k} \varvec{p}_{k}^{\top }\) for \(j=1, 2, 3\). Using these \(\{ \varvec{P}^{(j)} \}_{j=1}^3\), the low-rank tensor approximation [9] is given by

$$\begin{aligned} \mathcal {Y} = \mathcal {X} \times _1 (\varvec{P}^{(1)}\varvec{U}^{(1)})^{\top } \times _2 (\varvec{P}^{(2)}\varvec{U}^{(2)})^{\top } \times _3 (\varvec{P}^{(3)}\varvec{U}^{(3)})^{\top }, \end{aligned}$$
(18)

where \(\varvec{P}^{(j)}\) selects \(k_j\) bases of projection matrices \(\varvec{U}^{(j)}\). The low-rank approximation using Eq. (18) is used for compression in TPCA.

figure a

For HOSVD for third-order tensors, we have the following theorems.

Theorem 1

The compression computed by HOSVD is equivalent to the compression computed by TPCA.

(Proof) The projection that selects \(K=k_1 k_2 k_3\) bases of the tensor space spanned by \(u^{(1)}_{i_1} \circ u^{(2)}_{i_2} \circ u^{(3)}_{i_3}, \ i_j = 1, 2, \dots , k_j\) for \(j= 1, 2, 3\), is

$$\begin{aligned} \begin{aligned}&(\varvec{P}^{(3)}\varvec{U}^{(3)} \otimes \varvec{P}^{(2)}\varvec{U}^{(2)} \otimes \varvec{P}^{(1)}\varvec{U}^{(1)}) \\&\ \ \ = (\varvec{P}^{(3)} \otimes \varvec{P}^{(2)} \otimes \varvec{P}^{(1)}) (\varvec{U}^{(3)} \otimes \varvec{U}^{(2)} \otimes \varvec{U}^{(1)}) = \varvec{P}\varvec{W}, \end{aligned} \end{aligned}$$
(19)

where \(\varvec{W}\) and \(\varvec{P}\) are the projection matrix and a unitary matrix, respectively. Therefore, HOSVD is equivalent to TPCA for third-order tensors. ( Q.E.D.)

Furthermore, we have the following theorem.

Theorem 2

The HOSVD method is equivalent to the vector PCA method.

(Proof) The equation

$$\begin{aligned} \mathcal {Y} = \mathcal {X} \times _1 (\varvec{P}^{(1)}\varvec{U}^{(1)})^{\top } \times _2 (\varvec{P}^{(2)}\varvec{U}^{(2)})^{\top } \times _3 (\varvec{P}^{(3)}\varvec{U}^{(3)})^{\top } \end{aligned}$$
(20)

is equivalent to

$$\begin{aligned} \mathrm {vec} \mathcal {Y} = (\varvec{P}^{(3)}\varvec{U}^{(3)}\otimes \varvec{P}^{(2)}\varvec{U}^{(2)} \otimes \varvec{P}^{(1)}\varvec{U}^{(1)})^{\top } \mathrm {vec}\mathcal {X} = (\varvec{PW})^{\top }\mathrm {vec}\mathcal {X}. \end{aligned}$$
(21)

( Q.E.D.)

This theorem implies that the 3D-DCT-based reduction is an acceptable approximation of HOSVD for third-order tensors since this is an analogy of the approximation of PCA for two-dimensional images by the reduction based on the two-dimensional discrete cosine transform [10].

2.4 Reduction by Three-Dimensional Discrete Cosine Transform

For sampled one-dimensional signal \(\varvec{x}=(x_1,x_2,\dots ,x_n)^{\top }\), we have transformed signal \(\tilde{\varvec{x}} = (\tilde{x}_1,\tilde{x}_2,\dots ,\tilde{x}_n)^{\top }\) by using discrete cosine transform (DCT)-II [11, 12]. Matrix representation of the DCT transform is given by

$$\begin{aligned} \tilde{\varvec{x}} = \varvec{D}\varvec{x}, \ \ \varvec{D} = ((d_{ij})), \ \ d_{ij} = \cos \left( \frac{\pi }{n} \left( (j-1)+\frac{1}{2} \right) (i-1) \right) , \end{aligned}$$
(22)

where \(i, j = 1, 2, \dots n\).

For a set of a third-order tensor \( \{ \varvec{X}_i \}_{i=1}^N\) such that \( \varvec{X}_i \in \mathbb {R}^{n \times n \times n}\), setting a DCT matrix \(\varvec{D} \in \mathbb {R}^{n \times n}\) and projection matrices \(\varvec{P}^{(1)} \in \mathbb {R}^{k_1 \times n}\), \(\varvec{P}^{(2)} \in \mathbb {R}^{k_2 \times n}\) and \(\varvec{P}^{(3)} \in \mathbb {R}^{k_3 \times n}\), we define the 3D-DCT-based reduction by

$$\begin{aligned} \hat{\varvec{X}}_i = \varvec{X}_i \times _1 (\varvec{P}^{(1)}\varvec{D}) \times _2 (\varvec{P}^{(2)} \varvec{D}) \times _3 (\varvec{P}^{(3)}\varvec{D}), \end{aligned}$$
(23)

where \(k_1, k_2, k_3 < n\). The 3D-DCT-based reduction is an acceptable approximation of the compression by the PCA, TPCA and HOSVD. If we apply the fast Fourier transform to the computation of the 3D-DCT for tensor projections of each mode, the computational complexity is \(\mathcal {O}(n \log n)\).

3 Classification Based on Tensor Form

We adopt multilinear tensor subspace method for third-order tensors. This method is a three-dimensional version of the 2DTSM [3]. For a third-order tensor \(\mathcal {X}\), setting \(\varvec{U}^{(j)}, \ j=1, 2, 3\), to be orthogonal matrices, we call the operation

$$\begin{aligned} \mathcal {Y} = \mathcal {X} \times _1 \varvec{U}^{(1)\top } \times _2 \varvec{U}^{(2)\top } \times _3 \varvec{U}^{(3)\top } \end{aligned}$$
(24)

the orthogonal projection of \(\mathcal {X}\) to \(\mathcal {Y}\). Therefore, using this expression for a collection of matrices \(\{\mathcal {X}_i \}_{i=1}^M\), such that \(\mathcal {X}_i \in {\mathbb {R}^{I_1 \times I_2 \times I_3}}\) and \(\mathrm {E}(\mathcal {X}_i)=0\), the solutions of

$$\begin{aligned} \{ U^{(j)} \}_{j=1}^3 = \mathrm {arg} \ \mathrm {max} \ \mathrm {E} \left( \frac{ \Vert \mathcal {X} \times _1 \varvec{U}^{(1)\top } \times _2 \varvec{U}^{(2)\top } \times _3 \varvec{U}^{(3)\top } \Vert _{\mathrm {F}} }{\Vert \mathcal {X}_i \Vert _{\mathrm {F}} }\right) \end{aligned}$$
(25)

with respect to \(\varvec{U}^{(j)\top }\varvec{U}^{(j)} =\varvec{I}\) for \(j=1, 2, 3\) define a trilinear subspace that approximates \(\{\mathcal {X}_i \}_{i=1}^M\). Therefore, using orthogonal matrices \(\{ \varvec{U}_{k}^{(j)} \}_{j=1}^3\) obtained as the solutions of Eq. (25) for kth category, if a query tensor \(\mathcal {G}\) satisfies the condition

$$\begin{aligned} \text{ arg } \left( \underset{l}{\max } \frac{\Vert \mathcal {G} \times _1 \varvec{U}_{l}^{(1)\top } \times _{2} \varvec{U}_{l}^{(2)\top }\times _3 \varvec{U}_{l}^{(3)\top } \Vert _{\mathrm {F}}}{\Vert \mathcal {G} \Vert _{\mathrm {F}}} \right) =\{ \varvec{U}_{k}^{(j)} \}_{j=1}^3, \end{aligned}$$
(26)

we conclude that \(\mathcal {G} \in \mathcal {C}_k\), \(k,l=1,2, \dots , N_{\mathcal {C}}\), where \(\mathcal {C}_k\) and \(N_{\mathcal {C}}\) are the tensor subspace of kth category and the number of categories, respectively. For the practical computation of projection matrices \(\{ \varvec{U}_{k}^{(j)} \}_{j=1}^3\), we adopt the iterative method described in Algorithm 1.

4 Numerical Examples

We present two examples for extraction of outline shapes of volume data, and abilities of our method for classification of volumetric data. For experiments, we use the voxel images of human livers obtained as CT images. This image set contains 25 male-livers and seven female-livers. Note that these voxel images are aligned to their centre of gravity. In the experiments, we project these voxel images to small size tensors. For the projections, we adopt TPCA and the 3D-DCT. In the iterative method of TPCA, setting the number of bases to the size of the original tensors in Algorithm 1, we call the method full projection (FP). If we set the number of bases to smaller than the size of the original tensors in Algorithm 1, we call the method full projection truncation (FPT). Table 1 summarises the sizes and numbers of original and dimension-reduced voxel images.

Firstly, we show the approximation of a voxel image of a liver by three methods. The FP, FPT and 3D-DCT reduce the size of the data from \(89 \times 97 \times 76\) voxels to \(32 \times 32 \times 32\) voxels. Figure 1 illustrates volume rendering of original data and reconstructed data by these compressed tensors. Compared to Figs. 1(a) and 1(e), in Figs. 1(b)–(c) and (f)–(h), the FP, FPT and 3D-DCT preserve outline shapes of liver. In Fig. 1, the reconstructed data by the 3D-DCT gives a closer outline shape and more similar interior texture to those of the original than the FP and FPT. In Fig. 1, these results show that projections to small-size tensors extract outline shapes.

For the analysis of projected data by the FP, FPT and 3D-DCT, we decompose these projected tensors by Algorithm 1. Here, we set the size of bases in Algorithm 1 to \(32 \times 32 \times 32\) and use 35 projected tensors of livers for each reduction methods. In decompositions, we reordered eigenvalues \(\lambda ^{(j)}_i\), \(j=1, 2, 3\), \(i=1, 2, \dots , 32\) of the three modes to \(\lambda _i\), \(i=1, 2, \dots , 96\) in the descending order. Figure 2 shows the cumulative contribution ratios of reordered eigenvalues for the projected tensors obtained by the FP, FPT and 3D-DCT. Figure 3 illustrates reconstructed data obtained by using the 20 major principal components.

In Fig. 2, profiles of curves for three methods are almost coincident while the CCR of the 3D-DCT is a little bit higher than the others. In three methods, the CCRs become higher than 0.8 if we select more than 19 major principal components. In Fig. 3, shapes and interior texture for three methods are almost the same. In Figs. 3(d)–(f), the interior texture of a liver is not preserved and the outer shape is burred. In these results for three methods, major principal components represent outline shapes.

Secondly, we show results of the classification of voxel images of livers by the TSM. For the classification, we use 25 male-livers and seven female-livers since the sizes and shapes of livers between male and female are statistically different. Figures 4(a) and 4(b) illustrates the examples of livers of male and female, respectively. We use the voxel images of livers of 13 males and 4 females as training data. The residual voxel images are used as test data. In the recognition, we estimate the gender of livers. The recognition rate is defined as the successful estimation ratio for 1000 gender estimations. In each estimation of a gender for a query, queries are randomly chosen from the test dataset. For the 1-, 2- and 3-modes, we evaluate the results for multilinear subspaces with sizes from one to the dimension of the rejected tensors. Figure 4(c) shows the results of the classification. The TSM give 90 % recognition rate at the best with tensor subspace spanned by every two major principal axis of the three modes.

Table 1. Sizes and number of volumetric data of livers. \(\sharp \)data represents the number of livers obtained from different patients. The data size is the original size of the volumetric data. The reduced data size is the size of the volume data after tensor-based reduction.
Fig. 1.
figure 1

Original and reconstructed volumetric data of liver data. (a) shows the rendering of original data. (b)–(d) show the rendering of reconstructed data after the FP, FPT and 3DDCT, respectively. (e)–(f) illustrate axial slice images of these volumetric data in (a)–(d), respectively. The sizes of reduced tensors are shown in Table 1.

Fig. 2.
figure 2

Cumulative contribution ratios for three compressed tensors. For compression, we adopt FP, FPT and 3D-DCT. For the computation of the cumulative contribution ratio of eigenvalues obtained by the FP, we used all eigenvalues of modes 1, 2 and 3 after sorting them into descending order.

Fig. 3.
figure 3

Reconstruction by using only major principal components of the decomposition by the FP. Top and bottom rows illustrate volume rendering and axial slice of reconstructed data, respectively. For reconstruction, we use the 20 major principal components. Left, middle and right columns illustrate the results for the tensors projected by the FP, FPT and 3D-DCT.

Fig. 4.
figure 4

Recognition rates of gait patterns and liver data for original and compressed tensors. (a) and (b) illustrate the examples of livers of male and female, respectively. (c) shows the recognition rates for three compression methods. For compression, we use the HOSVD, FP, FPT and 3D-DCT. In (c), the horizontal and vertical axes represent the compression ratio and recognition ratio [%], respectively. For the original size \(D = 89 \times 97 \times 76\) and the reduced size \(K = k \times k' \times k''\), the compression ratio is given as D / K for reduced size k.

5 Conclusions

We applied the three-way TPCA to the extraction of outline shapes of volumetric data and classification of them. In the numerical examples, we demonstrated that the three-way TPCA extracts the outline shape of human livers. Furthermore, the TSM accurately classified the extracted outline shapes. Moreover, we showed that the 3D-DCT-based reduction approximated both the outline shape and texture of livers.

This research was supported by the “Multidisciplinary Computational Anatomy and Its Application to Highly Intelligent Diagnosis and Therapy” project funded by a Grant-in-Aid for Scientific Research on Innovative Areas from MEXT, Japan, and by Grants-in-Aid for Scientific Research funded by the Japan Society for the Promotion of Science.