High-Speed Maximum A Posteriori Birefringence Estimator for
Jones Matrix Optical Coherence Tomography by GPUImplementation
Atsushi Kubota 1)，Shuichi Makita 2), Yoshiaki Yasuno 2)
Sky technology Inc., Tsukuba, Ibaraki, Japan. 1), Computational Optics Group, University of Tsukuba, Tsukuba, Ibaraki, Japan. 2)11228-87
Background and Purpose

GPU architecture

Result

It is known that the birefringence distribution of a sample can be
known by measuring the amount of local phase retardation using
polarization OCT such as Jones-matrix OCT. However, the amountof local phase retardation of a general biological sample is small,and therefore, it generally has a low signal-to-noise ratio (SNR).
Maximum a posterior (MAP) birefringence estimation [1] hasbeen proposed to solve this problem. The maximum likelihood
estimate of birefringence has been obtained from multiple measured
local phase retardation values with a low SNR. This estimator ismathematically exact, however, it requires a long calculation time.
The previous implementation of MAP birefringence estimation isbased on CPU processing, and it took 40 min for the birefringenceestimation of an OCT volume.
In this research, we implemented a MAP birefringenceestimator using GPU and improved the speed of birefringenceestimation. Furthermore, the speed limiting factor was examined bycomparing the speeds using multiple GPUs.

In this implementation, we used GPU of NVIDIA pascal architecture.
In the Pascal architecture, 4 Warp units per Streaming
Multiprocessor (SM), 32 CUDA cores per Warp, and a total of 128
CUDA cores are available. For example, GTX1080Ti (Nvidia, CA,US), one of the Pascal architecture GPUs, has 28 SMs and 3,584
CUDA cores can be used. All CUDA cores in one SM share 96 kB ofShared Memory.
In order to speed up the GPU operation, it is important toefficiently use the Shared Memory and to simultaneously execute
more threads (= CUDA core). It is because the Shared Memory has
a high access speed but a small capacity.

MAP estimation of birefringence was performed using threedimensional JM-OCT data (512 × 640 × 256 pixels) in which
each point was repeatedly measured four times. Table 1
compares the estimation speeds of multiple GPU models andconventional CPU implementations. As a result, it was found thatthe use of the GPU can provide 8-times faster estimation speed
than the CPU.

MethodMAP birefringence estimator
The MAP birefringence estimator performs estimation using a"probability distribution of true birefringence when a certain effective
SNR (ESNR) value and a birefringence (bm) value are measured"(likelihood function).
Fig.1 shows a processing flow of the MAP estimator. Here,prior to the measurement, a likelihood function is calculated for all
combinations of ESNR and birefringence values by the Monte-Carlosimulation and stored in the main memory as a three-dimensional
array L(ESNR, Bm, bt).
Multiple (f-times) measurements at the same location of thesample is performed.
In the subsequent birefringence estimation process, thelikelihood function corresponding to each measured value isreferred to from the three-dimensional array L, and all of them are
multiplied to obtain a true posterior distribution of birefringence.
Finally, the birefringence value at which the posterior
distribution takes the maximum value is used as the MAP estimation
value.
ESNR (x, y, z, f)bm(x, y, z, f)

f=1

L (ESNR, bm, bt)

f=1btbmf=0

Warp Scheduler

Warp Scheduler

Register (16,384 x 32bit)

Register (16,384 x 32bit)

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Warp Scheduler

Warp Scheduler

Register (16,384 x 32bit)

Register (16,384 x 32bit)

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

Core

CoreCoreCoreCoreCore

CoreCoreCoreCoreCoreCoreCoreCoreCoreCoreCoreCore

Core

Core

CoreCore

CoreCoreCoreCore

Warp (=32 CUDA cores)

195.98

240

1280

GTX
1070M

0.47

119.44

384

2048

GTX
1080Ti

0.31

78368

6723584

CPU※

6.30

1,637

(unknown)---4 Warps x 32 CUDA cores= 128 CUDA cores per SM

Shared Memory (96kB)

Table 1 The processing performance of MAP birefringence estimates and GPU specs.

The number of threads that can be processed in parallel differs
among the GPU models. Fig. 4 shows the relationship between
the number of threads that can be processed in parallel and theestimated speed of the OCT slice (512 × 640 pixels). It is foundthat the estimation speed increases in proportion to the numberof threads that can be processed in parallel.
3.5

Fig. 2 Streaming Multiprocessor(SM) configuration of NVIDIA pascal architecture.

GPU implementation
In the process of estimating the MAP of 1 pixel (1 point) in a JMOCT image, it is necessary to store a one-dimensional likelihood
array data of float type 1,024 elements in the same memory and toread and write multiple times. To store this array, 1,024 × 4 bytes =4 kB of memory is required.
In this implementation, this 1-pixel estimation process isperformed by 1-thread (CUDA core), and multiple estimation
processes for multiple pixels are parallelly performed.
The number of parallelly performed estimation is limited by theshared memory size. This GPU has 96 kB shared memory per SM.
And hence, a SM parallelly performs 96 kB / 4 kB = 24 estimationsfor 24 pixels in parallel [Fig 3]. This parallel number is less than thenumber of CUDA cores (128) per SM, so parallelism is not limited by
the number of CUDA cores, but by the size of the shared memory.
CUDACore

CUDACore

24 CUDA Cores

CUDACore

CUDACore

24 Estimation process
Estimation Estimation
process process

Estimation Estimation Estimation
process process process

Shared Memory (96kB)

Π4kB

4kB

4kB

4kB

Fig. 3 MAP estimation process per SM of NVIDIA pascal architecture.
Fig. 1 Schematic diagram of data and their relationship in MAP birefringence estimate.

Total # of
CUDA cores※ Core i7-7700HQ 2.8GHz Processed with python 2.7 with Numpy1.16.4 (OpenBlas)

f=0 f=1 f=0 f=1…

Total # ofparallelly
processable
threadGTX 1080 Ti

ESNR

ESNR(xi yi, zi, 0)

0.77

SM (=4 Warps)

bm (xi, yi, zi, 0)(xi, yi, zi)

GTX 1060

96kB Shared Memory

CUDACore

f=0

CUDA core

Processingtime/slice[sec]4kB

3Estimation speed [Slice/s][1] D. Kasaragod, S. Makita, Y.-J. Hong, and Y. Yasuno, Biomed.
Opt. Express 8, 653–669 (2017).

SM

Processor

Total
processing
time (256slices) [sec]2.5

GTX 1070M
21.5

GTX 1060
10.5

0200

250

300

350

400

450500550

600

650700

Number of parallely processable threads

Fig. 4: The MAP processing speed related to the number of simultaneouslyprocessable threads.

Conclusion
The MAP birefringence estimation processing using GPU wasimplemented. Significantly faster processing speed was achieved
than a CPU implementation. The comparison among GPU models
showed that the processing speed increased linearly with thenumber of threads that could be executed in parallel.
It was found that the processing time of the MAP birefringenceestimation processing is limited by the shared memory size of theGPU SM. Therefore, the following two methods are expected tofurther improve the speed.
(1) Shared memory size per SM increased due to improved GPUarchitecture.
(2) Increase in the number of threads per pixel processing byimproving the operation algorithm.
Owing to this GPU-based MAP birefringence estimation, thebirefringence computation time became a non-significant portion offull JM-OCT image reconstruction, which includes OCTreconstruction, phase stabilization, OCT angiography computation,and degree-of-polarization uniformity computation. In the future,
faster computation of JM-OCT is expected by GPU implementationof other processing than MAP birefringence estimation.