A Tensor Factorization Method for 3D Super-Resolution with Application to Dental CT
Janka Hatvani Adrian Basarab Jean-Yves Tourneret Miklós Gyöngy Denis Kouamé








Hatvani et al.:
A Tensor Factorization Method for 3D Super-Resolution with application to Dental CT

A Tensor Factorization Method for 3D
 Super-Resolution
 with Application to Dental CT

    Janka Hatvani, Adrian Basarab, Jean-Yves Tourneret, Miklós Gyöngy and Denis Kouamé 

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.
This research has been partially supported by the European Union, co-financed by the European Social Fund (EFOP-3.6.3-VEKOP- 16-2017-00002 and by Pázmány University KAP17-19.
This work was supported by the thematic trimester on image processing of the CIMI Labex, Toulouse, France, with the Program ANR-11-IDEX-0002-02 under Grant
ANR-11-LABX-0040-CIMI.

J. Hatvani, A. Basarab, and D. Kouamé are with IRIT, CNRS UMR 5505, Université Paul Sabatier Toulouse 3, University of Toulouse, France.
J. Hatvani, and M. Gyöngy are with the Faculty of Information Technology and Bionics,
Pazmany Peter Catholic Unviersity, Budapest Hungary
J. Y. Tourneret is with  IRIT/INP-ENSEEIHT/Tésa, University of Toulouse, France.
    December 30, 2023
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





Available super-resolution techniques for 3D images are either computationally inefficient prior-knowledge-based iterative techniques or deep learning methods which require a large database of known low- and high-resolution image pairs. A recently introduced  tensor-factorization-based approach offers a fast solution without the use of known image pairs or strict prior assumptions. 
In this article this factorization framework is investigated for single image resolution enhancement  with an off-line estimate of the system point spread function. The technique is applied to 3D cone beam computed tomography for dental image resolution enhancement. To demonstrate the efficiency of our  method, it is compared to a  recent state-of-the-art  iterative technique using low-rank and total variation regularizations. 
In contrast to this comparative technique, the proposed reconstruction technique gives a 2-order-of-magnitude improvement in running time – 2 minutes compared to 2 hours for a dental volume of 282×266×392 voxels. Furthermore, it also offers slightly improved quantitative results (peak signal-to-noise ratio, segmentation quality). Another advantage of the presented technique is the low number of hyperparameters. As demonstrated in this paper, the framework is not sensitive to small changes of its parameters, proposing an ease of use.





3D super-resolution, single image super-resolution, tensor factorization, cone beam computed tomography, dental application





§ INTRODUCTION


Root canal treatment is carried out on a regular basis in dental centers in order to save decayed and infected teeth. In spite of their popularity, the success rate of the treatment is only 60-85% <cit.>.
For an improvement of the process the dentists need a better visualization of the canal, as its length, diameter and curvature are all important factors for planning the therapy <cit.>. Therefore, further research on the visualization possibilities of the pulp cavity is necessary, as stated by the European Commission on Radiation Protection in 2012 <cit.>.

Dental offices use cone beam computed tomography (CBCT) for determining the 3D structure of the teeth. Apart from the detector size, the spatial resolution of such imaging devices is also affected by partial volume effect, noise, and beam hardening, resulting in a typical value of 500 µm. This resolution is not sufficient in endodonty since the diameter of the canal is usually in the range of 0.16-1.60 mm and the apical, narrower segment is more important for planning the treatment <cit.>.
On the other hand, the resolution of micro-CT (µCT) is sufficient for precise measurements on the cavity, but the physical dimensions of the system only permit the imaging of extracted teeth. The long acquisition time and high radiation dose also prohibit in vivo measurements. 

The post-processing application of image super-resolution (SR) algorithms is an intensely investigated field in the image processing community. Classical SR techniques can combine information from a sequence of measurements <cit.>, from different modalities <cit.>, or in the simplest case they try to improve the resolution of a single image <cit.>. Many SR methods assume that the low resolution (LR) image of interest is obtained from the high resolution (HR) image by blurring and decimation with a residual additive noise. The SR problem can then be formulated as an inverse problem, which is ill-posed and thus requires an appropriate regularization to provide suitable solutions.A standard regularization often employed is total variation (TV) leading to piecewise smooth solutions <cit.>. Low-rank <cit.>, or wavelet representations <cit.> have  also proved to be efficient tools for SR. A method based on a sparse representation was applied to 3D MRI images with a patch-based structural similarity constraint in <cit.>.  Convolutional neural networks have also shown interesting properties for SR, where the network is trained to map an LR image to its HR counterpart <cit.>. However, this technique requires a large training dataset, which is not always available. 
Furthermore, only a few of the above-mentioned techniques is available for 3D volumes (e.g., <cit.>), and they all suffer from heavy computational costs, preventing them from application in practice.

A new hyperspectral-multispectral image fusion technique using tensor factorization (TF) was introduced in <cit.>. This technique combines a multispectral image (with high spatial and low spectral resolutions) and a hyperspectral image (with low spatial and high spectral resolution) to obtain an SR image (with high spatial and high spectral resolutions). 

One advantage of the tensor-based method of <cit.> is that it does not need to unfold the image of interest into a 2D matrix as in many existing SR methods <cit.>,<cit.>. As a consequence, this method avoids any loss of information about the locality of the image pixels and does not require to introduce spatial regularization (such as the TV of the image).

This paper investigates a 3D single image SR (SISR) method based on TF for CBCT images of teeth as an attempt to approximate their HR µCT pairs. The idea is to decompose the image of interest using its canonical polyadic decomposition (CPD). The CPD of a tensor is a representation based on a sum of an appropriate number of rank-1 tensors, which number depends on the structure of the image. It will be shown in this paper that this representation leads to a notably fast and efficient reconstruction method. The described method is compared to a state-of-the-art iterative deconvolution technique with low-rank and TV regularization (LRTV). For validation the peak signal-to-noise ratio (PSNR) is calculated and the canal is segmented, permitting volumetric and diametric comparison.  

The rest of this paper is organized as follows. First, the tensor operations used in this paper are defined and a connection between image complexity and tensor decomposition is drawn. In Section <ref> the proposed TF-SISR method is first defined for the CBCT resolution enhancement problem, followed by the data acquisition and the estimation of the blurring point spread function (PSF), ending with the introduction of the evaluation metrics. Section <ref> compares the images obtained by the two different SISR methods and discusses the possibilities and limits of tensor factorization. Finally a conclusion is drawn about the applicability of the introduced SR technique to dental imaging with some possible future work.




§ TENSORS AND IMAGE COMPLEXITY
 



 §.§ Notations

For easier distinction, 2D matrices are denoted using uppercase letters (e.g., A) and 3D tensors by bold uppercase letters (e.g., A). The uppercase letter with an overline (eg., A) denotes a set of 2D matrices. 



 §.§ Factorization, mode product, matricization


In this section, operations from tensor algebra necessary for the proposed method are summarized. Readers may refer to <cit.> and <cit.> for further details.

A tensor is a generalization of vectors and matrices, where the order of the tensor indicates the dimensionality. A 3D CT image volume is a third-order tensor X∈ℝ^I× J× K from which one dimensional fibers can be extracted. Depending on which indices are fixed, there are mode-1 fibers denoted as X(:,j,k) vectors (columns), mode-2 fibers denoted as X(i,:,k) vectors (rows) and mode-3 fibers denoted as X(i,j,:) vectors. The outer product (denoted by ∘)  between one mode-1, one mode-2 and one mode-3 array forms a rank-1 third order tensor, written as

    X = u∘ v∘ w, 
       u∈ℝ^I,v∈ℝ^J , w∈ℝ^K, X∈ℝ^I× J× K
       where
       X(i,j,k) = u(i) v(j) w(k).


The smallest number of rank-1 tensors that can sum up to form the tensor X is called the tensor rank of X, denoted by F. The resulting factorization of X is called the CPD of X expressed as

    X = ∑_f=1^F U^1(:,f)∘ U^2(:,f)∘ U^3(:,f)
       where
       X(i,j,k) =  ∑_f=1^F U^1(i,f)  U^2(j,f)  U^3(k,f).


U = {U^1,U^2,U^3}
is a set of three 2D matrices, {U^1∈ℝ^I× F, U^2∈ℝ^J× F, U^3∈ℝ^K× F}, known as the decomposition of the tensor X.
For illustration, the reader may refer to Fig. <ref>. In the following, the operation in (<ref>) will be denoted as

    [[U^1, U^2, U^3]]  = ∑_f=1^F U^1(:,f)∘ U^2(:,f)∘ U^3(:,f).



An important property of the CPD is that this decomposition is essentially unique (allowing permutations within U^1, U^2, U^3). Thus U can be identified almost surely if its tensor rank F is smaller than an upper bound. Chiantini et al. <cit.> proved that if I ≥ J ≥ K, with F ≤ 2^⌊ log_2J ⌋ +⌊ log_2K ⌋ - 2, the CPD of the rank-F tensor X∈ℝ^I× J× K is essentially unique. This condition allows identifiability of the CPD even for tensors with high rank. For example, in the application addressed herein, a typical CBCT volume with 260× 260× 300 pixels can be decomposed uniquely even if the tensor rank of the image is as high as 2^14 = 16384.

Next, the multiplication between a 2D matrix and a 3D tensor referred to as the mode product is defined. This multiplication can be performed along all three dimensions, and in each case the mode-n fibers of the tensor are extracted and are pre-multiplied by the matrix one-by-one. The mode-n products (n ∈{1,2,3}) of X∈ℝ^I× J× K with P_1∈ℝ^I^*× I , P_2∈ℝ^J^*× J, P_3∈ℝ^K^*× K are denoted as ×_n, and are defined as

    X×_1 P_1 =    X_1∈ℝ^I^*× J× K
    where    X_1(:,j,k) =  P_1X(:,j,k) 
    X×_2 P_2 =    X_2∈ℝ^I× J^*× K
    where    X_2(i,:,k)= P_2X(i,:,k) 
    X×_3 P_3 =     X_3∈ℝ^I× J× K^*
    where    X_3(i,j,:) =P_3X(i,j,:)


where I^*,J^*,K^* are arbitrary integer numbers. In Fig. <ref>, the principle of the mode-1 product, X×_1P_1 = X_1 is illustrated, where the columns of the tensor are premultiplied by P_1, leading to a shrinkage along the first dimension.  




Using the factorization of X in (<ref>) and (<ref>) the mode-n products can also be written as

    X×_1 P_1×_2 P_2×_3 P_3 = [[ P_1 U^1,P_2U^2,P_3U^3 ]].


Finally, the matricization or unfolding of the tensor X∈ℝ^I× J× K
from 3D to 2D is defined. Note that this matricization can be realized again along the three dimensions. For a mode-n matricization the mode-n fibers are extracted and form the columns of X^(n) in lexicographical order expressed as

    X^(1)=[   X(:,1,1),X(:,2,1),...
       X(:,J,1), X(:,1,2)...,X(:,J,K)]
    X^(1)   ∈ℝ^I× JK
    X^(2)=[   X(1,:,1),X(2,:,1),...
       X(I,:,1), X(1,:,2)...,X(I,:,K)]
    X^(2)   ∈ℝ^J× IK
    X^(3)=[   X(1,1,:),X(2,1,:),...
       X(I,1,:),X(1,2,:)...,X(I,J,:)]
    X^(3)   ∈ℝ^K× IJ.


The same operation can be realized using the decomposition U of X. For this the Khatri-Rao product – denoted as ⊙ – is necessary. It operates on matrices having the same number of columns, and calculates their column-wise Kronecker-product (A⊙ B = C, where A∈ℝ^I× F, B∈ℝ^J× F, C∈ℝ^IJ× F). Using this notation the matricization can be written as

    X^(1) = U^1(U^3⊙ U^2)^T
    X^(2) = U^2(U^3⊙ U^1)^T
    X^(3) = U^3(U^2⊙ U^1)^T.




 §.§ Connection between the tensor rank and the image complexity
 

The tensor rank is related to the complexity of the image. To illustrate this claim, an illustrative set of examples is provided in Fig. <ref>. In these examples the notations and dimensionality of  (<ref>) are used with I=J=K=2. Fig. <ref> a) shows that the image with a single dark pixel (representing 1) in the white volume (representing 0) has a tensor rank of F = 1. More precisely

    u_a ∘ v_a ∘ w_a = [1,0] ∘ [0,1] ∘[1,0] = X_a
       X^(3)_a = 
    [
    0    1    0    0  
     0    0    0    0
    ].


In Fig. <ref> b), two neighboring pixels are dark. This does not change the complexity of the image, since one outer product can still describe this volume. Indeed, we have

    u_b ∘ v_b ∘ w_b = [1,0] ∘ [0,1] ∘[1,1] =X_b
       X^(3)_b  = 
    [
    0    1    0    1  
     0    0    0    0
    ].

In Fig. <ref> c), two fibers are linearly dependent with 2×[21,35]=[42,70], which also makes one outer product sufficient for decomposing the tensor since

    u_c ∘ v_c ∘ w_c = [5,3] ∘ [1,2] ∘[7,0] = X_c
       X^(3)_c = 
    [
    35    70    0    0  
     21    42    0    0
    ].


This set of illustrative examples shows that for images with piecewise constant regions (like the neighboring cells in Fig. <ref> b)) or with low matrix rank (as the linearly dependent fibers in Fig. <ref> c)) a smaller tensor rank can be expected. More generally, the tensor decomposition (<ref>) tends to promote solutions with small tensor ranks. This property is useful in the case of denoising, when independent outlier pixels have to be eliminated.  A degraded image may contain larger constant areas, with higher dependency between neighboring rows and columns. It means that describing these images will also be more efficient with a tensor of small rank. Thus these simple examples allow us to understand why CBCT images can be represented by a reduced number of rank-1 tensors.
















§ METHODS
 


 §.§ Problem formulation


The image degradation model considered herein is the one classically used in SISR methods. It relates the LR image (CBCT in the case of the current dental application) to an HR image (considered to be close to the µCT). The HR image X∈ℝ^I× J× K is corrupted by a decimation operator D with rate r, a blurring kernel H, and some added noise N,  resulting in the LR image Y∈ℝ^I/r× J/r× K/r such that

    vec(Y) = DHvec(X) + vec(N)


where vec(·) vectorizes the elements of the 3D tensor in lexicographical order. We assume that H∈ℝ^IJK× IJK is the block-circulant version of the 3D Gaussian kernel h  to avoid circular convolution. A 3D Gaussian kernel h is separable along the three dimensions to 1D kernels as h = h_1∘ h_2∘ h_3 and is usually assumed for a blurring PSF <cit.>. For Gaussian kernels h_1,h_2,h_3 with standard deviations σ_1,σ_2,σ_3 the corresponding block-circulant matrices are H_1 ∈ℝ^I× I,H_2∈ℝ^J× J,H_3∈ℝ^K× K.
The decimation operator downsamples the image by an integer number, by averaging blocks of r neighboring pixels in each direction. In matrix form the downsampling operators for the three dimensions are D_1 ∈ℝ^I/r× I, D_2∈ℝ^J/r× J, D_3∈ℝ^K/r× K.
This formulation of D corresponds for instance to the physical process of a large CBCT detector element collecting rays over a larger area, than µCT does. This matrix also has better inversion properties compared to the regular decimation operator which discards pixels at a rate r.

Let U = { U^1 ∈ℝ^I× F, U^2 ∈ℝ^J× F, U^3∈ℝ^K× F} be the CPD  of X. 
The image degradation problem can be rewritten following (<ref>) using the separated kernels

    Y    = X×_1 D_1H_1×_2 D_2H_2×_3 D_3H_3 + N
        = [[ D_1H_1U^1,D_2H_2U^2,D_3H_3U^3 ]] + N.


The SR task can be defined as finding the set of matrices U which is the solution of the following minimization problem

    min_UY - [[ D_1H_1U^1,D_2H_2U^2,D_3H_3U^3 ]]_F^2


where ·_F denotes the Frobenius norm of a tensor defined as the square root of the sum of its squared elements. This cost-function is different from the minimization problem of <cit.> in the following aspects. First, only one measured datavolume is used here in contrast to the two measurements in the fusion problem of <cit.>. Second, here the degradation happens in all three dimensions between the HR and LR image, while in <cit.> the hyperspectral measurement is degraded in the first two dimensions, the multispectral volume in the third dimension. 

As problem (<ref>) is NP-hard an alternating optimization method is investigated, minimizing the cost function in (<ref>) sequentially for U^1, U^2, U^3.
Building a tensor from its decomposition  (<ref>) consists of a summation of F outer products. Taking the square of this sum of F terms in (<ref>) would result in the sum of
F+2-12 terms,

where ·· denotes the binomial coefficient.
Instead when minimizing over U_n, the tensors are mode-n-matricized using (<ref>) and (<ref>) leading to

    min_U^11/2Y^(1) - D_1H_1U^1(D_3H_3U^3⊙ D_2H_2U^2)^T_F^2
    min_U^21/2Y^(2) - D_2H_2U^2(D_3H_3U^3⊙ D_1H_1U^1)^T_F^2
    min_U^31/2Y^(3) - D_3H_3U^3(D_2H_2U^2⊙ D_1H_1U^1)^T_F^2.


Note that the unfolding is performed in each direction sequentially, conserving the 3D local information. These minimizations can now be carried out easily, leading to

    U^1 = (D_1H_1)^+Y^(1)(D_3H_3U^3⊙ D_2H_2U^2)^+T
    
    U^2 = (D_2H_2)^+Y^(2)(D_3H_3U^3⊙ D_1H_1U^1)^+T
    
    U^3 = (D_3H_3)^+Y^(3)(D_2H_2U^2⊙ D_1H_1U^1)^+T


where ^+ is the regularized Moore-Penrose pseudo-inverse defined as 

    A^+ = A^TA + ϵ^2I)^-1A^T


where ϵ is a hyperparameter used to provide a stable inverse (this procedure is classically referred to as diagonal loading).

To implement the proposed TF-SISR method the TensorLab toolbox <cit.> was used in Matlab 2017b. In the algorithm U
was initialized with elements from the standard normal distribution and U^1,U^2,U^3 were updated iteratively several times as described in Algo. <ref>.








 §.§ Data acquisition


The dataset used for testing contains images of 13 teeth which were extracted for health reasons and donated anonymously for research. This set consists of all different tooth types, including incisors, canines, premolars and molars.
A Carestream 81003D

system was used for CBCT imaging.
The linewidth resolution of the CBCT machine was 500 µm and the volumes had a voxel size of 80× 80× 80 µm^3. 

For evaluation purposes the reconstructed HR images were compared to µCT images acquired from the same samples.
The µCT acquisitions were obtained with a Quantum FX system from Perkin Elmer, with a voxel size of 40× 40× 40 µm^3 and linewidth resolution of 50 µm.








 §.§ PSF estimation


The proposed method in Algo. <ref> assumes that the PSF is known. In practice the blurring kernel has to be measured or estimated, which is usually carried out empirically in many existing works <cit.>. Here the blurring kernel was assumed to be Gaussian and its standard deviation was estimated from the observed data. Employing direct inverse filtering on each sample image, the µCT volume was divided by the CBCT volume in the frequency domain to obtain the Fourier transform of the PSF. The high-frequency band was suppressed by a Hanning-window before computing the inverse Fourier-transform and averaging the 13 estimated PSFs. For further details see <cit.>. Finally a 3D Gaussian function was fitted to the averaged PSF to estimate the standard deviations σ_1, σ_2 and σ_3.



















 §.§ Metrics
 




The comparison of the image volumes was carried out through two metrics.
The first one measures the PSNR between the estimated SR CBCT and µCT volumes. It is calculated by dividing the square of the dynamic range with the mean square error between the enhanced image and the µCT image, expressed in dB.

The second, more application oriented metric  consists of comparing volumes segmented from µCT and SR CBCT. The canal root was segmented with a dedicated adaptive local thresholding method (see <cit.> for details). The canal area and Feret's diameter (the longest distance between two parallel straight lines that are tangent to the shape) were calculated for each radicular axial slice. The measured values are compared through the mean of absolute differences. The differences of the  canal volumes were also measured  using the Dice coefficient <cit.>.
Finally MeVisLab <cit.> was used for visualizing the segmentation results.








§ RESULTS
 



 §.§ Comparison to an existing 3D SR method




The state-of-the-art LRTV introduced in 2013 in <cit.> was used as a benchmark to compare the performance of the proposed method[The Matlab code associated with LRTV is available at https://bitbucket.org/fengshi421/superresolutiontoolkithttps://bitbucket.org/fengshi421/superresolutiontoolkit].
 Among the relatively small collection of 3D SISR techniques, LRTV provided competitive results compared to other popular methods (cubic interpolation, non-local means, TV-based up-sampling) <cit.>.  It uses low-rank and total-variation regularizers, minimizing the following 
cost-function

    X̂ =    _XDHX-Y^2 
       + λ_RRank(X) + λ_TVTV(X),

where λ_R and λ_TV are hyperparameters.
The minimization problem <ref> is solved by the alternating direction method of multipliers (ADMM), which requires to adjust two additional hyperparameters, the penalty term ρ and an iteration number n_ADMM. One of the subproblems within the ADMM scheme is solved by gradient descent with an additional iteration number n_grad and an update rate denoted as dt.



The parameters used for testing can be seen in Table <ref>. They were tuned manually to get the highest possible improvement of the PSNR. The tests were run on a standard PC with an Intel(R) Core(TM) i7 2×2.5GHz processor and 16 GB of RAM.





The two methods were tested for three samples from the dataset, including an incisor, a premolar and a molar. The sizes of the sample volumes, the PSNR calculated against the µCT images and the execution times  are provided in Table <ref>. Compared to the CBCT images the PSNR improves similarly for the LRTV (average of 1.2 dB) and the TF-SISR (average of 1.5 dB) methods  with the chosen parameters. However, this enhancement is achieved at a much lower computational cost: 10 iterations of TF-SISR run 100 times faster than 5 iterations of LRTV. This faster execution time is important since it permits a wider range of applications, including those requiring a rapid diagnosis during a medical examination.
In Fig. <ref>, the quality of the enhanced volumes is visualized, showing that the canal is better defined and contrasted compared to the CBCT image, suggesting better segmentation properties.



    < g r a p h i c s >

[]Results on Sample #1. In the rows the CBCT, LRTV output, TF-SISR output and µCT images can be seen, whereas the columns correspond to one axial, a coronal and a sagittal slice. The CBCT image  is shown on the same scale as the HR images, for better comparison. The location of the slices within the volume is illustrated on the CBCT images in colored lines.





For further analysis the root canal was segmented from each volume, using the segmentation method described in Section <ref>. Qualitative and quantitative results are provided in Fig. <ref> and in Table <ref>. In particular, Table <ref> shows differences between the estimates and the values obtained using the µCT image for three parameters (Feret diameter, area of the canal and Dice coefficient). The estimated Feret diameter improves similarly with both SR techniques compared to the CBCT images with an averaged improvement of 63 µm for LRTV and 81 µm for TF-SISR. The second line of the table shows how the area of the canal on the axial slices is changing from one method to another. Note that the LRTV method shows a higher difference compared to the original CBCT (by 0.0256 mm^2), suggesting that the TV regularization overestimates the canal. This observation is also confirmed in Fig. <ref>, as the LRTV volumes have a more blueish color corresponding to positive differences. The TF-SISR provides the best overall performance with an improvement of 0.0152 mm^2 on average. The last metric in the table is the Dice-coefficient, also showing some improvement in the overlap of the canals, by 1% using the LRTV and 2% with the TF-SISR. Fig. <ref> displays zoomed-in sections of the apical part of the canal, as this part is the most important during the treatment. Considering these results, the TF-SISR method shows slightly better segmentation properties than the LRTV technique,  while offering a great reduction in running time.








 §.§ Adjusting the parameters of the TF-SISR method





The impact of the tensor rank and the iteration number was investigated using Sample #1.

Fig. <ref> a) shows that the runtime increases linearly with the number of iterations, as expected.

Fig. <ref> a) also shows that the PSNR converges rapidly to its maximum value (close to 24.5), which is an interesting property of the proposed method.
Fig. <ref> a) shows how the solutions qualitatively evolve with the iteration number. For improved visibility the difference from the n_TF = 10 case is shown in the figure. In the case of n_TF = 10 a second test was run, and the difference was calculated compared to this result, as the random initialization of U results in slightly different outputs. It can be seen that as the iteration number increases, the difference becomes less structured: in the third column the shape of the tooth is almost invisible and it is lost in random noise.


According to the upper limit of F for a unique CPD, in the case of Sample #1, F≤ 2^14 should be efficient. However, numbers higher than 2000 caused memory problems, and were therefore not tested. Fig. <ref> b) shows that the computational time increases exponentially with the rank F since the algorithm requires the inversion of larger matrices in U. It can also be seen that the PSNR stabilizes for ranks larger than F=500. Some sample images can be seen in Fig. <ref> b) showing that low values of the rank F lead to large constant blocks in the image, which is characteristic of a low-rank or TV regularization. For higher numbers finer details become visible.

Note that our results indicate that neither of the above parameters have to be estimated precisely. After a small number of iterations the result converges. F can be considered as a prior information on the complexity of the image. Using higher values a more natural result can be obtained, but above a threshold the method will not give more precise outputs.


















§ CONCLUSION



In this work, a tensor-factorization-based method was proposed for the 3D single image super-resolution problem. This method showed interesting computational advantages compared to currently available regularization-based methods, with slightly improved image quality compared to the investigated LRTV technique. The runtime of this method was about 100 times faster than with LRTV, allowing a wider field of applications. The method also uses significantly less parameters (tensor rank and iteration number) that can be easily adjusted by visual inspection of the reconstruction results.
Dental CBCT volumes used as experimental data showed improved PSNR and canal-segmentation properties, with moderately better results than the LRTV method.  Considering these results, the method was found to be promising for 3D single image super-resolution.

The prior information in regularization-based techniques is often empirical and guides the solution. Future work can investigate if such classical priors could be included in this framework, and whether they would improve the result. As <cit.> proposed a solution also for embedded PSF estimation, its application to TF-SISR  could be a potential direction of further research.





§ ACKNOWLEDGMENT



The authors would like to thank Jérôme Michetti for the data acquisition and the work on the visualization of the teeth.


bib/IEEEtran




