This is a post-peer-review, pre-copyedit version of an article published in IEEE Transactions on Medical Imaging. The final authenticated version is available online at: http://dx.doi.org/10.1109/TMI.2018.2852685. © 2019 IEEE. Personal use of this material is permitted. Permission from IEEE must 1 be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. Mosaic-Based Color-Transform Optimization for Lossy and Lossy-to-Lossless Compression of Pathology Whole-Slide Images Miguel Hernández-Cabronero∗ , Victor Sanchez, Member, IEEE, Ian Blanes, Senior Member, IEEE, Francesc Aulı́-Llinàs, Senior Member, IEEE, Michael W. Marcellin, Fellow, IEEE, and Joan Serra-Sagristà, Senior Member, IEEE. Abstract—The use of whole-slide images (WSIs) in pathology life-critical in isolated, impoverished or otherwise underserved entails stringent storage and transmission requirements because regions where experts are not available. Moreover, virtual of their huge dimensions. Therefore, image compression is an pathology can increase pathologist availability [2] and reduce essential tool to enable efficient access to these data. In particular, color transforms are needed to exploit the very high degree of costs in pathology laboratories [3]. Recently, computer-aided inter-component correlation and obtain competitive compression diagnosis (CAD) algorithms are being developed [4]–[12] to performance. Even though state-of-the-art color transforms re- assist the pathologists in the diagnosis, especially in tedious move some redundancy, they disregard important details of the and error-prone tasks. WSIs can also enhance performance to compression algorithm applied after the transform. Therefore, help pathologists detect and classify clinical markers, and to their coding performance is not optimal. We propose an optimiza- tion method called Mosaic Optimization for designing irreversible create annotated databases for training, research and reference and reversible color transforms simultaneously optimized for any purposes [13]. given WSI and the subsequent compression algorithm. Mosaic An important difficulty when dealing with WSIs is their Optimization is designed to attain reasonable computational huge size. A 15×15 mm specimen scanned at 40× magnifica- complexity and enable continuous scanner operation. Exhaustive tion requires a 60000×60000, 24-bit RGB color image, result- experimental results indicate that, for JPEG 2000 at identical compression ratios, the optimized transforms yield images more ing in over 10 GB of raw data. Even at a lower, more common similar to the original than other state-of-the-art transforms. 0.5 µm resolution, images attain sizes of 30000 × 30000 and Specifically, irreversible optimized transforms outperform the 2.5 GB of uncompressed data [1], [14]. Karhunen-Loève Transform (KLT) in terms of PSNR (up to Furthermore, multiple WSIs –sometimes exceeding 35 [15]– 1.1 dB), the HDR-VDP-2 visual distortion metric (up to 3.8 dB) are commonly produced per case. Thus, hundreds of Terabytes and accuracy of computer-aided nuclei detection tasks (F1 score up to 0.04 higher). Additionally, reversible optimized transforms are produced yearly by a single WSI scanner [15], [16]. These achieve PSNR, HDR-VDP-2 and nuclei detection accuracy gains massive data volumes hinder the adoption of virtual pathol- of up to 0.9 dB, 7.1 dB and 0.025, respectively, when compared ogy. Excessive bandwidth demands can handicap the remote to the reversible color transform (RCT) in lossy-to-lossless visualization of WSIs, while stringent storage necessities – compression regimes. including archiving and backup systems– can diminish the Index Terms—Color-Transform Optimization, Whole-Slide Im- cost-effectiveness of digital pathology workflows. In such a ages, Image Compression. scenario, image compression is paramount to alleviate these storage and transmission problems. I. I NTRODUCTION Lossless coding provides perfect reconstruction of the data produced by the WSI scanner. To the best of our knowledge, In the traditional pathology workflow, pathologists examine only one work has proposed a lossless compression algorithm specimens under an optical microscope to detect biological specific for these images. In [17], High Efficiency Video structures relevant to the diagnosis. In virtual pathology, tra- Coding (HEVC) [18] intra coding mode is adapted to the ditional glass slides are being replaced by digital whole-slide multidirectional patterns that are typically present in these images (WSIs). WSIs allow concurrent examination of a case images. In particular, sample-by-sample differential pulse- by several pathologists, possibly in different geographical lo- code modulation is combined with uniform angular modes cations, to reduce medical errors [1]. Remote diagnosis can be to improve prediction across strong edges. However, perfect reconstruction comes at the cost of relatively low compression This work has been funded by the EU Marie Curie CIG Programme under Grant PIMCO, the Engineering and Physical Sciences Research Council ratios, commonly between 2:1 and 4:1, which may not suffice (EPSRC), UK, and by the Spanish Ministry of Economy and Competitiveness for interactive and remote WSI visualization purposes. and the European Regional Development Fund under Grant TIN2015-71126- On the other hand, lossy compression can offer almost R (MINECO/FEDER, UE) and by the Catalan Government under grant 2014SGR-691. arbitrary ratios at the cost of introducing some distortion in the ∗ M. Hernández-Cabronero is with the University of Warwick and the images. A number of studies –involving pathologists and real University of Arizona, USA (email: miguelh@email.arizona.edu). V. Sanchez specimens– have shown that the images reconstructed with is with the University of Warwick, UK. M. W. Marcellin is with the University of Arizona, USA. F. Aulı́-Llinàs, I. Blanes and J. Serra-Sagristà are with the lossy coding methods are suitable for clinical and research Universitat Autònoma de Barcelona, Spain. purposes. Compression ratios up to 20:1 with JPEG 2000 [19] 2 yield statistically indistinguishable diagnostic results from follows, this method is referred to as Mosaic Optimization. those obtained using traditional glass slides [20]–[22], and do For any input image, Mosaic Optimization produces a color not affect the pathologists’ confidence on the diagnosis [23]. transform that is optimized simultaneously for the image and Average compression ratios of approximately 13:1 produce for the compression algorithm applied after the transform. As images that pathologists deem indiscernible from the original a result, the produced transform attains compression perfor- WSIs [24]. Compression ratios up to 200:1 can hinder human- mance results that improve upon the state of the art. made diagnoses [21], but not some carcinoma-detection CAD To allow reasonable computational complexity, the pro- algorithms [23]. Traditional JPEG [25] has also been shown to posed Mosaic Optimization employs a novel mosaic-based produce reconstructed images that can be employed for diag- sampling technique. Representative blocks of the image are nosis [26]. Produced results are on par with the ones obtained first extracted using a fast heuristic, and then an optimization using glass slides [27], [28], although exact compression ratios process is applied to a mosaic of the extracted blocks. Mosaic for JPEG have not been reported. An accurate rate-allocation Optimization builds upon the Full Optimization method, orig- method specific for WSIs has also been proposed [29] for inally presented in [31]. Even though both methods produce lossy regimes with HEVC. In most commercial WSI scanners, transforms that yield very similar compression performance, standard JPEG 2000 or JPEG compression is employed [14]. Full Optimization has a prohibitive computational complexity, Color transforms improve compression by removing inter- i.e., exceeding 2 hours per image on average–, which makes component redundancy [30]. Due to the especially strong sim- it unusable in practical scenarios. The proposed Mosaic Opti- ilarities among color components in WSIs, applying adequate mization method enables execution times close to 2.4 minutes transforms can greatly improve compression performance. on the workstation employed for this research, i.e., 54 times WSI-specific color transforms are scarce [31], [32], although faster than Full Optimization and comparable to the time it other transforms can also be applied to exploit inter-component takes to scan the WSIs [14], [27]. This allows continuous similarities. In HEVC, fixed color-space translations can be pipelining of the scanning and the compression of WSIs using applied to the prediction residuals using the Adaptive Color reasonable computing power and requiring additional storage Transform [33]. JPEG 2000 Part 1 supports the irreversible space for only one scanned image. Mosaic Optimization also color transform (ICT) and the reversible color transform builds upon the Single Patch Optimization method intro- (RCT) to decorrelate color components [19]. Arbitrary linear duced in [32]. Although both exhibit similar computational transforms can be used in Part 2 of the standard [34]. In complexity, Mosaic Optimization yields better compression commercial WSI scanners [14] and WSI compression works performance than Single Patch Optimization, as discussed in the literature [20]–[24], the ICT is the most common color later in Section IV. We also propose an extension to design transform. The Karhunen-Loève Transform (KLT) typically optimized reversible color transforms that can be used for produces highly competitive compression results, which nor- lossless and progressive lossy-to-lossless compression. This mally improve upon the ICT. For Gaussian sources, the KLT can be useful in scenarios where the original data must be is known to provide perfect color component decorrelation. stored losslessly, e.g., due to legal restrictions, and bandwidth- Notwithstanding, WSIs do not follow a Gaussian distribution. efficient lossy representations are necessary, e.g., for remote Moreover, the KLT does not take into account the details of the visualization of the WSIs. compression algorithm applied after the transform. Hence, its Exhaustive experimental results are provided to evaluate the compression performance is not optimal either. The Optimal color transforms produced with Mosaic Optimization when Spectral Transform algorithm (OST [35]) produces linear applied to WSIs depicting different tissue types. The PSNR– transforms that minimize the mutual information among color based on the mean squared error–, the HDR-VDP-2 [39] met- components when they are considered as random variables. ric –calibrated to emulate a human observer– and the accuracy Mutual information of two random variables X and Y can be of a CAD nuclei detection algorithm [10] are employed for expressed as I(X; Y ) = H(X) − H(X|Y ), where H(·) de- this purpose. Results indicate that the Mosaic Optimization notes the marginal entropy and H(·|·) the conditional entropy. produces irreversible transforms that outperform the ICT and Therefore, lower mutual information implies less redundant the KLT in terms of PSNR, HDR-VDP-2 and nuclei detection information in X and Y , which often improves compression performance. The optimized reversible transforms improve performance. The OST considers the spatial wavelet transform upon the RCT with respect to the PSNR, HDR-VDP-2 and and the quantization of JPEG 2000, but does not consider other CAD metrics in lossy-to-lossless regimes, i.e., when only part important aspects of JPEG 2000 such as its entropy coder. of a losslessly compressed codestream is used for reconstruc- Therefore, its actual compression performance is also subop- tion. Both the optimized transforms and the RCT produce very timal. Color deconvolution methods are successfully applied in similar bitrates in lossless compression regimes. CAD algorithms for WSIs [36]–[38]. These methods separate To recapitulate, we propose a color transform optimization image colors by using so-called deconvolution matrices, which method based on [31]. The main novelties of this work are can be considered color transforms as well. However, these (a) an original mosaic-based method to reduce computational transforms are not designed for compression and do not yield complexity that outperforms the one described in [32]; (b) competitive coding performance in practice. an extension of [31] to allow progressive lossy-to-lossless This work introduces a fast mosaic-based optimization compression; and (c) an exhaustive performance analysis that method for designing color transforms that overcomes the extends those of [31], [32] with more test images, new tissue inherent drawbacks of the approaches described above. In what types and comparison with more color transform design algo- 3 discussion for all types of WSI, a total of 16 types of tissue are hereinafter considered. The average inter-component correlation for these images exceeds 0.95, whereas for the ISO-12640-1 [41] and ISO-12640-2 [42] sets of natural images the average correlation is 0.81 and 0.66, respectively. Due to the high inter-component correlation present in WSIs, employing adequate color transforms is crucial to ob- tain competitive compression performance. Traditional static transforms such as YCb Cr [19] and YIQ [43] translate each RGB pixel into a different color space using a fixed 3 × 3 matrix M = (mi,j ):  u   m  r  (a) (b) (c) 1,1 m1,2 m1,3 v = m2,1 m2,2 m2,3 g , (1) Fig. 1: Sample WSIs depicting (a) lymphatic tissue; (b) w m3,1 m3,2 m3,3 b pancreatic tissue; (c) renal tissue. where (u, v, w) represents the pixel in the new color space and mi,j are real-valued entries. Typical transforms are not rithms. With the proposed method, the coding performance of designed for the particular characteristics of WSIs. Thus, they existing compression algorithms can be significantly enhanced fail to properly decorrelate color components, which in turn and, improving upon [31], only reasonable computational may negatively affect compression. complexity is required. Moreover, the proposed lossy-to- Data-dependent approaches adapt the coefficients of the lossless extension broadens the applicability of this work to transform matrix M to the image to be compressed. For situations where a lossless copy of the image must be stored. instance, the well-known Karhunen-Loève Transform (KLT) The rest of the paper is structured as follows. Section II removes any inter-component correlation in the transformed describes key properties of WSIs that motivate the use of color image components. Other methods based on independent transformations for compression. Section III describes the component analysis such as the Optimal Spectral Transform proposed color-transform optimization framework, including (OST) [35] minimize the mutual information. In Strutz’s the mosaic-based approach. Section IV provides extensive work [44], an adaptive transform is implemented by choosing experimental results and Section V concludes the paper. an optimal mode based on the entropy of the image trans- formed with that mode. Another reversible transform was proposed by Kim et al., which applies an additional lifting II. W HOLE -S LIDE I MAGES AND C OLOR T RANSFORMS step to decorrelate the U and V components produced by the A key step to generate WSIs is to stain the biological ICT [45]. These data-dependent approaches yield significant sample to highlight structures relevant to the diagnostic task at compression performance gains as compared to non-adaptive hand. The most commonly employed stain for many types of methods. Although data-dependent methods in the literature WSIs is Hematoxylin and Eosin (H&E) [40], which exposes are designed to optimize criteria such as inter-component cell nuclei and cytoplasm. The stained slide is then inserted correlation or mutual information, none of these methods in a special scanner with objective lenses for magnification, considers all aspects of the coding algorithm employed after which registers a color image depicting the complete glass the color transform. For instance, the KLT does not take slide. This takes less than 3 minutes per slide [14], [27]. into account any spatial decorrelation applied after the color Optical magnification factors of 20×, 40× or even transform. The OST includes the 2D discrete wavelet trans- 100× [1], [10], [27] are necessary to provide a microscope-like form (DWT) in its model, but not the probability model viewing experience. These magnification factors correspond used to drive the subsequent entropy coding stage. In turn, to resolutions of 0.5 µm, 0.25 µm and 0.1 µm per pixel, Strutz’s method evaluates each candidate color transform by respectively, and are often used with numerical aperture (NA) applying LOCO-I’s [46] spatial predictor after the spectral of 0.75 [14]. A comprehensive description of the hardware decorrelation. The entropy of the produced residuals is then and methods employed to scan WSIs is available in [14]. used as a surrogate for the actual bitrate of the compressed As a result of the high magnification levels, the 24-bit RGB bitstream. Hence, existing data-dependent color transforms do WSIs output by the scanners have very high spatial resolutions, not attain optimal compression performance either. which motivates compression, but also limits the complexity of the algorithms that can be used under reasonable time and III. C OLOR T RANSFORM O PTIMIZATION memory constraints. Due to the staining process, each image exhibits a limited range of color gamut. H&E-stained images typically contain This section describes a color transform optimization algo- violet, blue and red hues only, as can be seen in Fig. 1. rithm that attains higher compression efficiency for WSIs than Therefore, WSIs exhibit inter-component redundancy signif- the KLT and other state-of-the-art component transforms. Sec- icantly higher than regular photographic imagery. Table I tion III-A describes a numerical transform-optimization frame- shows average values if these key characteristics for 50 WSIs work and Section III-B describes the Mosaic Optimization considered in this work. In order to provide a meaningful algorithm, which is a fast approximation to this framework. 4 TABLE I: Properties of 50 H&E-stained WSIs considered in this work. Set #Images #Tissue Dimensions Average Inter-Component Correlation Types Min Size Max Size ρ(R, G) ρ(R, B) ρ(G, B) Average Lymphatic 6 1 12288×17408 17408×25600 0.964 0.976 0.983 0.975 Pancreatic 10 1 12470×12938 31655×32619 0.918 0.919 0.990 0.942 Renal 7 1 9352×8139 34018×16327 0.903 0.908 0.991 0.934 Mixed 27 13 34018×16327 32001×21186 0.955 0.954 0.992 0.967 All images 50 16 9352×8139 31655×32619 0.942 0.943 0.990 0.958 An extension of Mosaic Optimization that allows progressive lossy-to-lossless transmission is introduced in Section III-C. Note that, even though other parts of the compression pipeline can also be optimized to provide improved compres- sion results [45], the scope of this work is limited to studying and optimizing the color transform stage for WSIs. A. Optimization Framework The design of efficient color transforms can be regarded as Fig. 2: Diagram of the E VALUATE D ISTORTION cost function. an optimization problem. For instance, the KLT is defined to remove inter-component correlation, while the OST minimizes inter-component mutual information. These optimization cri- teria generally have a positive effect on compression, although squared error (MSE) is employed due to its low computational they are not specifically devised to do so except under certain complexity and its efficient integration with the JPEG 2000 assumptions such as high-rate scalar quantizers and Gaussian standard. As discussed later in Section IV, optimizing for MSE sources for the KLT, which may not hold in practice and ultimately improves results for the HDR-VDP-2 metric and for therefore lead to suboptimal performance [47]. Contrarily, the the accuracy of CAD nuclei detection tasks. proposed framework is conceived to optimize compression performance directly. Given an input image I, our framework Deriving an analytical expression for (2) is not practical aims at finding a transform matrix M that minimizes the cost for coding systems such as JPEG 2000. This is due to the function E VALUATE D ISTORTION, complex non-linear nature of reversible transforms, spatial decorrelation transforms, quantizers, entropy coders and rate- M = argmin {E VALUATE D ISTORTION(M 0 , I)} . (2) allocation algorithms. Therefore, unlike for the KLT or the M0 OST, which assume simpler models, closed-form solutions This function appraises a color transform M 0 by comparing cannot be easily obtained for (2). I with the reconstructed image obtained after compression and decompression with M 0 . This scheme is depicted in To produce competitive transforms, we employ our previ- Fig. 2. Note that, for each I, the same bitrate R must be ously proposed iterative numerical optimization framework, used for compression in E VALUATE D ISTORTION for all color T RANS O PT [31], summarized in Fig. 3. In the Evaluate Candi- transforms M 0 so that their relative performance can be dates box, several candidate transform matrices are compared compared. A single M that works well for a wide range of using the E VALUATE D ISTORTION cost function described R can be obtained by using R = RICT (I)/2, where RICT (I) above. Here, the identity transform (i.e., not applying any color is the compression rate attained for I by lossy JPEG 2000 transform), the KLT, the YCb Cr and the YIQ (a rotation of using Kakadu v7.8 [48] with its default parameters with the YCb Cr ) transforms are employed as candidates. Note that the irreversible color transform, without limiting the target bitrate. set of initial transforms was empirically chosen to maximize The value of R for each test image considered in this work compression performance while minimizing the number of is available as supplementary content at http://ieeexplore.ieee. such transforms. Only the best of these, i.e., the one that yields org. the smallest cost for the input image, is further considered. By construction, the resulting transform M is not only opti- This transform, denoted “Best candidate M 0 ” in Fig. 3, mized for I, but also for the compression algorithm employed undergoes iterative optimization based on Powell’s conjugate in the cost function. As discussed in Section I, JPEG 2000 direction algorithm [49]. In the i-th iteration of this algorithm, has been successfully employed for compression of WSIs. the coefficients of M 0 are modified using a conjugate gradient Part 2 of the standard provides the syntax for user-defined method to produce Mi0 , to which E VALUATE D ISTORTION is color transforms [34]. Therefore, JPEG 2000 is considered applied to obtain Costi . Powell’s algorithm employs Mj0 and 0 exclusively hereafter. Nevertheless, any other compression Costj , j ≤ i, to generate Mi+1 . This process, depicted with algorithm that admits arbitrary linear color transforms could blue dashed arrows in the figure, is repeated Niter times. The also be used within this framework. output optimized transform Mopt is then defined as the Mi0 Any distortion metric can be used to compare the input with the smallest associated Costi among all iterations. No and the reconstructed images in Fig. 2. In this work, the mean improvements are observed for Niter > 250. 5 Algorithm 1 Mosaic-based generation of a smaller image. 1: function G ENERATE M OSAIC(inputImage, blockSize, ϕ) 2: gridBlocks ← C REATE G RID(inputImage, blockSize) 3: for all gridBlocks as block do 4: block.colors ← C OUNT D IFFERENT C OLORS(block) 5: end for 6: sortedBlocks ← S ORT B Y C OLORS(gridBlocks) 7: candidateBlocks ← sortedBlocks[|sortedBlocks| · ϕ: 8: |sortedBlocks|] 9: selectedBlocks ← U NIFORM S AMPLE(candidateBlocks) 10: mosaicImage ← C OMBINE C ELLS(selectedBlocks) Fig. 3: Diagram of the T RANS O PT transform optimization 11: return mosaicImage framework. 12: end function Fig. 4: Diagram of the mosaic-based fast approximation to T RANS O PT. B. Mosaic Optimization When the T RANS O PT framework is applied to WSIs, it generates color transforms that improve upon the KLT and other state-of-the art methods in terms of compression per- Fig. 5: Crop of a mosaic image produced with Algorithm 1. formance. However, this approach –hereafter referred to as Full Optimization– exhibits a run time exceeding 130 minutes averaged for the 50 images referenced in Table I, as discussed For this work, random uniform sampling is used, denoted by later in Section IV. A fast mosaic-based approximation to Full U NIFORM S AMPLE in the Algorithm. Optimization, denoted Mosaic Optimization, is proposed next. The resulting mosaic image is then created by concatenating Over 95% of the run time of T RANS O PT is spent in the the selected blocks (line 10), as depicted in Fig. 5. Note that by E VALUATE D ISTORTION function in compression, decompres- using adequate values of ϕ, blocks depicting only background sion and comparison tasks. To accelerate the execution, a small regions, which carry no useful information for the optimization version S of the input image is generated first, and then the process, are rarely included in the mosaic images. optimization framework is applied to S, as depicted in Fig. 4. When the produced mosaic image is used in the It is important to generate a representative version of the full- E VALUATE D ISTORTION function, the JPEG 2000 compressor size image so that the resulting color-transform MS yields is invoked with tiles of size blockSize, so that the 2D good coding performance for the original image. The proposed spatial wavelet transform is never performed across the (ar- method for generating S, referred to as G ENERATE M OSAIC, tificial) block boundaries. Additionally, identical bit budgets is provided in Algorithm 1. The image is first divided in a are assigned to each mosaic block to avoid individual blocks rectangular grid with blocks of size blockSize×blockSize dominating the rate allocation and the cost function. Since (line 2). Some of the blocks are then selected to build the the selected blocks are chosen in a spatially uniform way, mosaic image as follows. The number of different pixel colors the mosaic image is representative of all regions depicting is computed for each grid block (lines 3–5). Here, a color stained tissue. Moreover, the probability of a given biological is defined as the 3D vector (r, g, b), and two colors (r, g, b) structure appearing in a block is comparable to the fraction and (r0 , g 0 , b0 ) are considered identical if and only if r = r0 , of tissue area depicting this structure. Therefore, this mosaic- g = g 0 and b = b0 . Blocks that contain very few different based approach produces color transforms MS optimized to colors typically depict portions of the background, which is produce good overall performance across the original WSI. normally bright and homogeneous, as can be observed in This is in contrast with the fast method presented in [32], Fig. 1. Conversely, blocks with many different colors tend to where the small image version contains pixels from a single depict stained tissue and be more representative of the relevant area of the image. parts of the image. Consistent with this, grid blocks are sorted The maximum area of the mosaic image, hereinafter referred by the number of colors (line 6) and only the top fraction to as maximum mosaic area (MMA), determines the time ϕ is kept as candidate blocks for the mosaic (line 7), where complexity of the proposed algorithm. Note that it is the MMA ϕ is a real parameter in (0, 1). Finally, blocks are selected and not ϕ that has a measurable impact on the size of the with uniform probability from the list of candidates until the produced mosaic images, and therefore on execution times. In total number of pixels reaches a certain target value (line 9). general, smaller MMAs need shorter execution times, although 6 coding performance may be affected when the area is too small. We found that an MMA equal to 0.25% of the original image is adequate to achieve competitive coding performance and computational complexity. This value is employed in the experimental results. Notwithstanding, this parameter can be adjusted to attain faster execution times with a small impact in coding performance. As discussed later in Section IV, using this value results in execution times comparable to the scan- ning time of many commercial WSI scanners [14], with coding performance similar to Full Optimization. The size of the blocks, determined by blockSize, must be selected carefully. Too small blocks yield non-representative performance from Fig. 6: Diagram of R EVERSIBLE E VALUATE D ISTORTION. JPEG 2000 because of boundary effects at tile edges when the DWT is applied. Too large blocks reduce the number of selected grid blocks, thus limiting the diversity of image (TCGA [53]) data portal. Images were selected to depict 16 regions and features that can be represented in the mosaic different tissue types and a broad range of hues so that results image. Block sizes between 128 × 128 and 256 × 256 are can be representative of WSIs in general. The Mosaic Opti- found to produce the best results. mization method and its lossy-to-lossless extension are applied to each image to generate optimized irreversible and reversible C. Progressive Lossy-to-Lossless Optimization transforms, respectively. The obtained transform matrices are The optimized color-transform matrices produced by Mo- available as supplementary content at http://ieeexplore.ieee. saic Optimization have real-valued entries. They are irre- org. The Kakadu v7.8 implementation of JPEG 2000 Part 2 versible and cannot be directly employed for efficient lossless is employed for compression and decompression [48]. In all or progressive lossy-to-lossless compression of WSIs. In some cases, Niter = 250. Our evaluations shown that blockSize = situations it is desirable or even necessary to retain the original 128 and ϕ = 0.5 produce the best overall results, hence image data without loss, e.g., due to archiving policies or legal only these parameters are hereinafter considered. To evaluate constraints. Mosaic Optimization can be extended to produce complexity, total execution times have been measured on a optimized, reversible integer-to-integer color transforms. WSIs dedicated workstation with a 4-core i5-4570 CPU at 3.2 GHz can be losslessly compressed with these transforms, achieving and 16 GB of RAM. All reported time results are relative to compression ratios similar to other state-of-the-art lossless this platform. coding algorithms. The losslessly compressed image files Section IV-B compares Mosaic Optimization to the identity can be transmitted partially or progressively, yielding rate- transform (i.e., no color transform), the ICT, the KLT, the distortion results that improve upon the state of the art. OST [35], the YCoCg [54], the color transforms proposed The proposed extension of the T RANS O PT framework to by Strutz et al. [44] and Kim et al. [45], the Full Opti- lossless and lossy-to-lossless scenarios is based on the work mization method described in Section III-A, and the Single by Galli et al. [50]. In that work, transform matrices are Patch Optimization method introduced in [32]. The 50 images factorized as a product of permutation matrices and triangular summarized in Table I are compressed and reconstructed using elementary reversible matrices (TERMs), similar to a well- each color transform at 15 bitrates ranging from 0.04 to 1.6 known lower-upper LU matrix decomposition [51], where bits per pixel per component (bpppc). These correspond to TERMs are defined as triangular matrices containing only compression ratios from 200:1 to 5:1, typically considered ±1 in the diagonal. The resulting factorization is equivalent in the literature [20]–[24]. Results for HEVC [18] using the to a reversible integer-to-integer transform that approximates Adaptive Color Transform1 are also provided for reference. the original lossy transform [50]. In [52], another algorithm The optimized reversible transforms are compared in for producing reversible transforms based on the matrix fac- Section IV-C to the default reversible color transform of torization proposed in [51] is described. Both [50] and [52] JPEG 2000 (RCT [19]), the OST [35], the YCoCg-R trans- produce accurate approximations. Galli’s method is employed form [55], the reversible KLT (RKLT [50]) of each image, exclusively hereafter, and referred to as L OSSLESS A PPROX. and the Single Patch Optimization method. Images are loss- The E VALUATE D ISTORTION function is modified to as- lessly transformed, compressed and reconstructed at 13 bitrates sess the compression performance of the reversible ap- linearly spaced between 0.04 bpppc and 1.2 bpppc. Five proximation of M 0 in each iteration of the design pro- of the 50 images are not included in the comparison since cedure. The resulting evaluation function, denoted by their lossless compression bitrate is lower than 1.2 bpppc. R EVERSIBLE E VALUATE D ISTORTIONand depicted in Fig. 6, Reversible transforms are not compared to HEVC since it substitutes E VALUATE D ISTORTION in T RANS O PT. does not provide progressive lossy-to-lossless compression for single images. IV. E XPERIMENTAL R ESULTS Several image distortion metrics have been employed to A. Materials and Methods assess the severity of the distortion introduced in WSIs [23], To evaluate the performance of the proposed framework, 1 Produced with the HM reference software v16.5, available at https://hevc. we use 50 H&E images from the Cancer Genome Atlas hhi.fraunhofer.de/ 7 TABLE II: Average BD∗ -PSNR and BD∗ -HDR-VDP-2 results in dB for irreversible color transforms. Color transform Lymphatic Pancreatic Renal Mixed All Images PSNR VDP PSNR VDP PSNR VDP PSNR VDP PSNR VDP No transform 42.56 19.13 51.22 38.86 48.31 32.09 47.31 27.07 47.67 29.20 YCoCg 48.33 41.15 51.57 42.06 47.71 30.53 47.82 31.47 48.63 34.93 ICT 47.85 38.72 51.26 42.37 47.01 29.43 47.31 31.38 48.14 34.22 KLT 46.89 38.01 52.38 42.65 48.96 32.45 48.49 30.82 49.14 34.30 OST 46.79 38.82 52.53 41.26 49.50 33.79 47.65 30.78 48.79 34.29 Full Optimization 48.37 42.36 52.70 43.00 49.43 33.41 49.18 31.74 49.83 35.54 Single Patch Optimization 48.02 41.64 50.74 39.69 49.54 32.85 48.56 31.14 49.07 34.37 Mosaic Optimization blockSize = 128, ϕ = 0.5 48.33 42.82 52.63 43.92 49.27 32.52 49.02 31.77 49.70 35.65 [24], [28], [29], [56]. In [56], the ability of several metrics to TABLE III: Average execution times in minutes for obtaining detect synthetic noise is evaluated. PSNR is employed in [29] irreversible color transforms. to compare different rate-allocation algorithms for WSIs. Algorithm Lym. Pan. Ren. Mix. All In [24], the PSNR, the SSIM and a Visual Discrimination KLT 0.2 0.3 0.2 0.3 0.3 Model (VDM) based on the works by Daly [57] and Lubin [58] OST 19.7 15.1 17.6 19.0 18.0 are investigated as tools to guide compression of WSIs. In [23], Full 65.6 113.3 90.5 162.3 130.8 Single Patch 1.2 2.2 1.4 2.7 2.3 [28], the accuracy of CAD algorithms is employed to evaluate Mosaic imaging and compression systems, respectively. Consistent blockSize = 128, 1.3 2.5 1.6 2.9 2.4 ϕ = 0.5 with this, the coding performance of the different transforms is appraised using the PSNR, the High Dynamic Range Vi- sual Difference Predictor v2.2.1 [39] (HDR-VDP-2)2 and the using the F1 score defined as effect on a recent CAD nuclei detection algorithm [10]. The HDR-VDP-2 metric extends [57] and [58] to predict visible 2 · TP F1 Score = , (4) differences between two images, and has been shown [39] 2 · TP + FP + FN to provide higher accuracy than other visual metrics such as where T P (true positives) are the nuclei detected in the SSIM, MS-SSIM [60] and mDCT-PSNR [61]. original and the reconstructed image, F P (false positives) are To compare the average PSNR performance of each trans- the nuclei detected only in the reconstructed image and F N form across all tested bitrates, a modified version of Bjonte- (false negatives) are the nuclei detected in the original image gaard’s BD-PSNR [62] metric is defined as but not in the reconstructed image. F1 scores lie in [0, 1], and 1 Z rmax higher values indicate better detection performance. ∗ BD -PSNR = PSNR(r) dr. (3) rmax − rmin rmin B. Irreversible Transforms Here, PSNR(r) is the PSNR between the original and recon- structed image at rate r, and rmin and rmax are the minimum Average BD∗ -PSNR and BD∗ -HDR-VDP-2 results for lossy and maximum reconstruction bitrates. This integral is numeri- color transforms are provided in Table II for each test set cally approximated with Simpson’s rule at the aforementioned individually3 . Consistent with [31], [32], and as expected due reconstruction target bitrates. To compare the HDR-VDP-2 to the inter-component correlation described in Table I, all metric across the range of tested bitrates, the BD∗ -HDR-VDP- of the transforms (optimized or not) outperform the case of 2 metric is defined analogously. The BD∗ -PSNR and BD∗ - not using any color transform. As can be observed, the ICT – HDR-VDP-2 metrics provide, respectively, the average PSNR commonly used in WSI scanners [14] and in the WSI compres- and HDR-VDP-2 values in dB across the tested bitrate ranges, sion literature [20]–[24]– yields average BD∗ -PSNR results hence higher values correspond to smaller distortion. 1 dB lower than the KLT, and very similar BD∗ -HDR-VDP- To evaluate the effect of the transforms on the CAD nuclei 2 results to the KLT. The YCoCg transform attains similar detection algorithm described in [10], each of the 100 H&E performance to the ICT, with some improvements in terms images employed in that study is compressed and recon- of PSNR. In turn, the optimized transforms produced by the structed at the same range of bitrates as described above. Full, Single Patch and Mosaic modes yield average results that The CAD algorithm is then applied to the original and the are more competitive than the ICT and the KLT. In particular, reconstructed images using the most competitive irreversible Full Optimization and Mosaic Optimization improve upon the and reversible transforms. The nuclei coordinates detected ICT and the KLT for each set and for both distortion metrics. in the original images are used as ground truth and the Average results for the OST are slightly lower than those of performance of each transform at each bitrate is calculated the KLT in terms of PSNR, and very similar in terms of the HDR-VDP-2 metric. 2 The probability map produced by the HDR-VDP-2 method is summarized 3 BD∗ -PSNR and BD∗ -HDR-VDP-2 results for HEVC are not included into a single value in dB using the Minkowsky mean as described in [59]. in Table II because the maximum bitrates are below 1 bpppc for most The map is calculated using the default observer model, assuming an RGB images, even for the lowest QP value of 1. Hence, a fair comparison for display with 100 cd/m2 peak luminance, 1 cd/m2 black level, and a resolution rates r ∈ [rmin , rmax ] is not possible. Results for HEVC are included in of 30 pixels per degree. Figs. 7a and 7b. 8 (a) (b) (c) Fig. 7: Average rate-distortion comparison of irreversible transforms for all images (a) using the PSNR metric; (b) using the HDR-VDP-2 metric; (c) using the F1 score based on [10]. Fig. 8: Magnified sample crop (50×90) of image NCI02 of the Lymphatic set. (left) Original; (center) reconstructed at 0.1 bpppc using Mosaic Optimization; (right) reconstructed at 0.1 bpppc using the KLT. Areas of interest are highlighted in yellow. Full Optimization yields better BD∗ -PSNR results than i.e., 20% data volume reduction. As can also be observed in Mosaic Optimization for all test sets. This is because all Figs. 7a and 7b, the proposed method is better than HEVC image samples are available for optimization in the Full mode. with the Adaptive Color Transform for the PSNR and HDR- Notwithstanding, Mosaic Optimization provides very similar VDP-2 metrics. results for each set, between 0.04 dB and 0.16 dB lower Average F1 scores for the KLT, for the ICT and for Mosaic than the Full mode. In terms of BD∗ -HDR-VDP-2, Mosaic Optimization obtained for the CAD nuclei detection algorithm Optimization exhibits gains over Full Optimization of up to described in [10] are plotted in Fig. 7c. As can be observed, 0.92 dB for the Pancreatic set, and 0.11 dB on average for all transforms produce F1 scores above 0.9 at bitrates as all images. This can be explained by the fact that Mosaic low as 0.4 bpppc. Mosaic Optimization produces F1 scores Optimization prioritizes candidate grid blocks depicting tissue. consistently higher than the ICT and the KLT for all tested Hence, the resulting color transforms enable more accurate bitrates. These results suggest that, for any compression ratio, reconstruction of visually relevant areas, e.g., edges between the Mosaic Optimization method enables more accurate CAD nuclei and cytoplasm. analysis results than the ICT and the KLT. Average PSNR and HDR-VDP-2 rate-distortion results are Average execution times for obtaining different data- plotted in Figs. 7a and 7b, respectively. Results indicate that dependent irreversible color transforms are provided in Ta- PSNR gains between 0.6 dB and 1.1 dB and HDR-VDP- ble III. Note that these times refer to the time it takes to 2 gains up to 3.8 dB over the KLT are obtained with the generate the color transform matrices, not the time it takes proposed method at the same bitrates. Therefore, at identical to compress an image using that transform matrix. For all compressed file sizes, images reconstructed using the Mosaic tested image sets, the KLT is the fastest algorithm. However, Optimization method exhibit smaller distortion –e.g., fewer as discussed above, its coding performance is lower than that morphological and chromatic artifacts– than those recon- of the optimization methods. Full Optimization has an average structed using the KLT, as illustrated in Fig. 8. Note that a low execution time exceeding 130 minutes, rendering it impractical bitrate (0.1 bpppc) is used in Fig. 8 so that artifacts are more for most use cases. On the other hand, Mosaic Optimization easily appreciable. These artifacts disappear at higher bitrates. requires under 2 minutes for the Lymphatic and Renal image Mosaic Optimization attains similar PSNR and HDR-VDP-2 sets, and under 3 minutes for the Pancreatic and Mixed sets. levels –i.e., similar quality– at bitrates between 0.05 bpppc These execution times are comparable to the scanning time of and 0.10 bpppc lower than the KLT, i.e., smaller compressed many commercial scanners, i.e., up to 3 minutes depending files. Thus, the transmission of images at similar quality on the objective magnification [14], [27]. levels can be completed faster. For instance, PSNR values of To recapitulate, the optimized transforms produced by the approximately 47 dB are produced on average at 0.5 bpppc proposed method outperform the ICT –the most usual trans- for the optimized transforms and at 0.6 bpppc for the KLT, form in commercial scanners and in the WSI compression 9 TABLE IV: Average BD∗ -PSNR and BD∗ -HDR-VDP-2 results in dB for reversible color transforms. Color transform Lymphatic Pancreatic Renal Mixed All Images PSNR VDP PSNR VDP PSNR VDP PSNR VDP PSNR VDP No Transform 35.37 11.83 44.07 18.48 41.03 12.28 41.21 12.95 40.59 13.84 RCT 42.87 18.64 45.43 20.04 41.13 13.13 42.26 15.34 42.82 16.85 YCoCg-R 39.88 12.02 42.25 18.13 41.10 12.83 41.77 12.56 41.79 13.74 Strutz’s [44] 42.91 18.97 45.46 20.05 41.15 13.17 42.44 15.48 42.87 16.93 Kim’s [45] 41.56 11.82 42.12 18.61 41.38 12.72 41.70 12.65 42.51 13.79 OST 41.58 20.08 44.84 19.15 41.69 13.21 42.28 14.53 42.62 16.15 RKLT 42.06 21.19 44.84 18.71 41.35 13.04 42.43 15.62 42.72 16.64 Single Patch Optimization 42.88 22.29 45.37 20.12 41.75 13.26 42.91 15.13 43.24 17.36 Mosaic Optimization blockSize = 128, ϕ = 0.5 43.11 25.07 45.47 21.09 41.76 13.33 42.86 15.51 43.27 17.61 (a) (b) (c) Fig. 9: Average rate-distortion comparison of reversible transforms for all images (a) using the PSNR metric; (b) using the HDR-VDP-2 metric; (c) using the F1 score based on [10]. literature– as well as transforms such as the KLT and the OST about 0.6 bpppc (i.e., a compression ratio of about 13:1) have for traditional, visual and CAD-based distortion metrics. been described to produce reconstructed images indiscernible from the original [24]. These ranges are well aligned with the bitrate ranges for which the proposed Mosaic method yields C. Reversible Transforms compression gains. Average BD∗ -PSNR and BD∗ -HDR-VDP-2 results for the Lossless compression bitrates are provided in Table V. On reversible transforms are provided in Table IV. As can be average for all images, all tested methods yield similar results. observed, the proposed Mosaic Optimization method improves This indicates that the progressive lossy-to-lossless coding upon the RCT, the OST, Single Patch and the RKLT in almost gains produced by the Mosaic Optimization method are not all cases. In general, larger improvements are observed for at the expense of hindering lossless compression efficiency, as the BD∗ -HDR-VDP-2 than for the BD∗ -PSNR metrics. Again, compared to other works in the literature. this is explained by the fact that optimization is designed to Average F1 scores for the RCT and for Mosaic Optimization focus on visually relevant areas of the WSIs. Note that the obtained for the CAD nuclei detection algorithm described compression performance differences between the irreversible in [10] are plotted in Fig. 9c. As can be observed, both transforms described in Section IV-B can be largely explained the RCT and Mosaic Optimization allow accurate CAD-based by the different DWTs used for 2D decorrelation, i.e., the CDF analysis of the images, with F1 scores above 0.9 for bitrates 9/7 DWT for irreversible transforms and the integer-to-integer of 0.27 bpppc or higher. Identical F1 scores can be obtained CDF 5/3 DWT for reversible ones. at bitrates up to 0.1 bpppc lower with Mosaic Optimization As shown in Figs. 9a and 9b, significant PSNR improve- than with the RCT. ments over the RKLT are observed for bitrates of 0.2 bpppc Average execution times for producing reversible trans- and higher, with gains up to 0.9 dB. HDR-VDP-2 improve- forms using different methods are shown in Table VI. Note ments become significant from 0.4 bpppc onwards, with max- that compression time is not included in this table. As can imum gains of 3.7 dB. Results suggest that similar distortion be observed, generating the reversible optimized transforms can be attained for smaller data volumes with the proposed requires about 19 minutes on average for all images. The method, e.g., an average PSNR of 44 dB can be obtained higher time complexity, compared to the irreversible transform transmitting 14% less data. Note that bitrates of 0.4 bpppc optimization methods discussed in the previous section is (i.e., compression ratios of 20:1) are often described in the expected, due to the additional approximation stage required in WSI compression literature [21]–[24] and a average bitrates of each invocation of the R EVERSIBLE E VALUATE D ISTORTION 10 TABLE V: Average lossless compression rates in bpppc. scores 0.025 higher than the RCT at the same bitrates. Similar Color transform Lym. Pan. Ren. Mix. All distortion levels are obtained at bitrates up to 14% lower than No Transform 3.41 1.98 1.96 2.07 2.14 for the RCT. Moreover, Mosaic Optimization is designed to RCT 2.44 1.90 1.98 1.97 2.00 attain execution times comparable to that of a typical WSI YCoCg-R 2.36 1.87 1.99 2.02 1.99 scanning process. This enables continuous operation of the Strutz’s [44] 2.35 1.86 1.94 1.99 1.96 Kim’s [45] 2.42 1.84 2.01 2.03 1.99 WSI scanner with only reasonable computational resources OST 2.96 2.01 1.99 2.06 2.15 dedicated to compression. In summary, the proposed color- RKLT 2.87 2.09 1.99 1.94 2.09 Single Patch 2.71 2.01 1.93 1.95 2.02 transform optimization method can produce tangible benefits Mosaic for WSI-centered laboratories in practical scenarios. blockSize = 128, 2.69 1.96 1.87 1.92 1.99 ϕ = 0.5 ACKNOWLEDGMENTS TABLE VI: Average execution times in minutes for obtaining We thank M. Barret for kindly providing an implementation reversible color transforms. of the OST algorithm from [35], and the authors of [10] Algorithm Lym. Pan. Ren. Mix. All for sharing an implementation of their CAD nuclei detection RKLT 0.44 0.52 0.41 0.49 0.47 algorithm. OST 19.94 15.33 17.84 19.21 18.48 Strutz’s [44] 0.65 1.74 0.81 1.58 1.39 Kim’s [45] 0.51 1.20 0.54 1.14 0.99 R EFERENCES Single Patch 8.74 23.64 10.50 22.61 19.46 Mosaic [1] M. May, “A better lens on disease,” Scientific American, vol. 302, no. 5, blockSize = 128 pp. 74–7, May 2010. ϕ = 0.5 9.40 21.26 10.91 22.72 19.18 [2] J. Ho, S. M. Ahlers, C. Stratman, O. Aridor, L. Pantanowitz, J. L. Fine, J. A. Kuzmishin, M. C. Montalto, and A. V. Parwani, “Can digital pathology result in cost savings? A financial projection for digital pathology implementation at a large integrated health care organization,” function, depicted in Fig. 6. Journal of Pathology Informatics, vol. 5, no. 33, 2014. [3] M. Isaacs, J. K. Lennerz, S. Yates, W. Clermont, J. Rossi, and J. D. Pfeifer, “Implementation of whole slide imaging in surgical pathology: A V. C ONCLUSIONS value added approach,” Journal of Pathology Informatics, vol. 2, no. 39, 2011. WSIs can help pathologists improve their diagnostic accu- [4] A. Madabhushi, “Digital pathology image analysis: opportunities and racy, reach underserved regions and enhance their productivity, challenges,” Imaging in Medicine, vol. 1, no. 1, pp. 7–10, 2009. thus reducing the operational costs of pathology laboratories. [5] O. Sertel, J. Kong, H. Shimada, U. V. Catalyurek, J. H. Saltz, and M. N. Gurcan, “Computer-aided Prognosis of Neuroblastoma on Whole-slide However, the enormous data volumes associated with WSIs Images: Classification of Stromal Development.” Pattern recognition, make image compression necessary to enable efficient stor- vol. 42, no. 6, pp. 1093–1103, Jun. 2009. age and agile remote access. In particular, lossy compres- [6] M. D. DiFranco, G. O’Hurley, E. W. Kay, R. W. G. Watson, and P. Cunningham, “Ensemble based system for whole-slide prostate cancer sion can greatly reduce compressed data volumes without probability mapping using color texture features.” Computerized Medi- affecting pathologists’ diagnoses or their confidence therein. cal Imaging and Graphics, vol. 35, no. 7-8, pp. 629–45, Jan. 2011. WSIs exhibit exceptionally high inter-component correlation, [7] S. Samsi, A. K. Krishnamurthy, and M. N. Gurcan, “An efficient compu- tational framework for the analysis of whole slide images: Application to which can be exploited using color transforms that enhance follicular lymphoma immunohistochemistry,” Journal of Computational coding performance. Color transforms in the literature are Science, vol. 3, no. 5, pp. 269–279, Sep. 2012. designed to minimize theoretical criteria such as correlation, [8] J. Xu, L. Xiang, Q. Liu, H. Gilmore, J. Wu, J. Tang, and A. Madabhushi, “Stacked Sparse Autoencoder (SSAE) for Nuclei Detection on Breast but do not fully take into account all details of the coding Cancer Histopathology Images,” IEEE Trans. Med. Imag., vol. 35, no. 1, algorithm applied after the transform. Therefore, in practice, pp. 119–130, Jan. 2016. their compression performance is suboptimal. The proposed [9] S. B. Ginsburg, G. Lee, S. Ali, and A. Madabhushi, “Feature Importance in Nonlinear Embeddings (FINE): Applications in Digital Pathology,” Mosaic Optimization method is tightly integrated with the IEEE Trans. Med. Imag., vol. 35, no. 1, pp. 76–88, Jan. 2016. implementation of the coding algorithm used in the Kakadu [10] K. Sirinukunwattana, S. E. A. Raza, Y. W. Tsang, D. R. J. Snead, implementation and the produced color transforms improve I. A. Cree, and N. M. Rajpoot, “Locality Sensitive Deep Learning for Detection and Classification of Nuclei in Routine Colon Cancer upon the ICT, the KLT and other methods in terms of Histology Images,” IEEE Trans. Med. Imag., vol. 35, no. 5, pp. 1196– mean squared error (PSNR), visual distortion metrics (HDR- 1206, May 2016. VDP-2) and F1 score for a recent CAD nuclei detection [11] A. Vahadane, T. Peng, A. Sethi, S. Albarqouni, L. Wang, M. Baust, K. Steiger, A. M. Schlitter, I. Esposito, and N. Navab, “Structure- algorithm. At identical rates, PSNR values up to 1.1 dB Preserving Color Normalization and Sparse Stain Separation for Histo- higher, HDR-VDP-2 results up to 3.8 dB higher and average logical Images,” IEEE Trans. Med. Imag., vol. 35, no. 8, pp. 1962–1971, F1 scores up to 0.04 higher than the KLT are obtained. Aug. 2016. [12] B. E. Bejnordi, M. Balkenhol, G. Litjens, R. Holland, P. Bult, N. Karsse- Conversely, equal distortion levels are produced for bitrates meijer, and J. A. W. M. van der Laak, “Automated Detection of DCIS in between 0.05 bpppc and 0.17 bpppc lower, i.e., file sizes up Whole-Slide H&E Stained Breast Histopathology Images,” IEEE Trans. to 20% smaller than for the KLT. An extension to design Med. Imag., vol. 35, no. 9, pp. 2141–2150, Sep. 2016. [13] The Digital Pathology Association. (2017) Whole Slide Imaging reversible transforms is proposed as well to allow lossless and Repository. [Online]. Available: https://digitalpathologyassociation.org/ progressive lossy-to-lossless compression. Without hindering whole-slide-imaging-repository lossless compression performance, the optimized transforms [14] N. Farahani, A. Parwani, and L. Pantanowitz, “Whole slide imaging in pathology: advantages, limitations, and emerging perspectives,” Pathol- yield average PSNR values up to 0.9 dB higher, average ogy and Laboratory Medicine International, vol. Volume 7, p. 23, Jun. HDR-VDP-2 results up to 7.1 dB higher, and average F1 2015. 11 [15] D. R. J. Snead, Y.-W. Tsang, A. Meskiri, P. K. Kimani, R. Crossman, [36] A. C. Ruifrok and D. A. Johnston, “Quantification of histochemical N. M. Rajpoot, E. Blessing, K. Chen, K. Gopalakrishnan, P. Matthews, staining by color deconvolution,” Analytical and quantitative cytology N. Momtahan, S. Read-Jones, S. Sah, E. Simmons, B. Sinha, S. Suor- and histology, vol. 23, no. 4, pp. 291–299, 2001. tamo, Y. Yeo, H. El Daly, and I. A. Cree, “Validation of digital pathol- [37] M. K. K. Niazi, K. Yao, D. Zynger, S. Clinton, J. Chen, M. Koyuturk, ogy imaging for primary histopathological diagnosis,” Histopathology, T. LaFramboise, and M. Gurcan, “Visually meaningful histopathological vol. 68, no. 7, pp. 1063–1072, 2016. features for automatic grading of prostate cancer,” IEEE Journal of [16] A. Huisman, A. Looijen, S. M. van den Brink, and P. J. van Diest, Biomedical and Health Informatics, 2016, In press. “Creation of a fully digital pathology slide archive by high-volume tissue [38] M. K. K. Niazi, A. V. Parwani, and M. N. Gurcan, “Computer-assisted slide scanning,” Human pathology, vol. 41, no. 5, pp. 751–757, 2010. bladder cancer grading: α-shapes for color space decomposition,” in [17] V. Sanchez, F. Aulı́-Llinàs, J. Bartrina-Rapesta, and J. Serra-Sagristà, SPIE Medical Imaging, 2016, pp. 979 107–1 – 979 107–8. “HEVC-based Lossless Compression of Whole Slide Pathology Images,” [39] R. Mantiuk, K. J. Kim, A. G. Rempel, and W. Heidrich, “HDR-VDP- in IEEE Conference on Signal and Information Processing (GlobalSIP). 2: a calibrated visual metric for visibility and quality predictions in IEEE, Dec. 2014, pp. 297–301. all luminance conditions,” in ACM Transactions on Graphics (TOG), [18] G. J. Sullivan, J. R. Ohm, W. J. Han, and T. Wiegand, “Overview of the vol. 30, no. 4. ACM, 2011, p. 40. High Efficiency Video Coding (HEVC) Standard,” IEEE Trans. Circuits [40] A. Khan, N. Rajpoot, D. Treanor, and D. Magee, “A Non-Linear Map- Syst. Video Technol., vol. 22, no. 12, pp. 1649–1668, Dec. 2012. ping Approach to Stain Normalisation in Digital Histopathology Images [19] D. S. Taubman and M. W. Marcellin, JPEG 2000: Image Compression using Image-Specific Colour Deconvolution,” IEEE Trans. Biomed. Eng., Fundamentals, Standards and Practice. Kluwer Academic Publishers, vol. 61, no. 6, pp. 1729–1738, 2014. Nov. 2001. [41] Graphic technology - Prepress digital data exchange - Part [20] T. Kalinski, R. Zwönitzer, S. Sel, M. Evert, T. Guenther, H. Hofmann, 1: CMYK standard colour image data, International Standard J. Bernarding, and A. Roessner, “Virtual 3D Microscopy Using Multi- Organization (ISO) Std. 12 640-1, 1997. [Online]. Available: http: plane Whole Slide Images in Diagnostic Pathology,” American Journal //www.iso.org/iso/catalogue detail.htm?csnumber=40928 of Clinical Pathology, vol. 130, no. 2, pp. 259–264, 2008. [42] Graphic technology - Prepress digital data exchange - Part 2: RGB en- [21] T. Kalinski, R. Zwönitzer, F. Grabellus, S. Y. Sheu, S. Sel, H. Hofmann, coded standard colour image data, International Standard Organization J. Bernarding, and A. Roessner, “Lossy compression in diagnostic virtual (ISO) Std. 12 640-2, 2004. 3-dimensional microscopy – where is the limit?” Human Pathology, vol. 40, no. 7, pp. 998–1005, 2009. [43] W. H. Buchsbaum, Color TV Servicing, E. Cliffs, Ed. Prentice Hall, [22] T. Kalinski, R. Zwönitzer, F. Grabellus, S. Y. Sheu, S. Sel, H. Hofmann, 1975. and A. Roessner, “Lossless compression of JPEG2000 whole slide [44] T. Strutz, “Multiplierless reversible color transforms and their automatic images is not required for diagnostic virtual microscopy,” American selection for image data compression,” IEEE Trans. Circuits Syst. Video Journal of Clinical Pathology, vol. 136, no. 6, pp. 889–895, 2011. Technol., vol. 23, no. 7, pp. 1249–1259, Jul. 2013. [23] S. Doyle, J. Monaco, A. Madabhushi, S. Lindholm, P. Ljung, L. Ladic, [45] S. Kim and N. I. Cho, “Hierarchical prediction and context adaptive J. Tomaszewski, and M. Feldman, “Evaluation of effects of JPEG2000 coding for lossless color image compression,” IEEE Trans. Image compression on a computer-aided detection system for prostate cancer Process., vol. 23, no. 1, pp. 445–449, Jan 2014. on digitized histopathology,” in IEEE International Symposium on [46] M. J. Weinberger, G. Seroussi, and G. Sapiro, “The loco-i lossless Biomedical Imaging: From Nano to Macro (ISBI), Apr. 2010, pp. 1313– image compression algorithm: Principles and standardization into jpeg- 1316. ls,” IEEE Transactions on Image processing, vol. 9, no. 8, pp. 1309– [24] J. P. Johnson, E. A. Krupinski, M. Yan, H. Roehrig, A. R. Graham, and 1324, 2000. R. S. Weinstein, “Using a visual discrimination model for the detection [47] M. Effros, H. Feng, and K. Zeger, “Suboptimality of the Karhunen- of compression artifacts in virtual pathology images,” IEEE Trans. Med. Loeve transform for transform coding,” IEEE Trans. Inf. Theory, vol. 50, Imag., vol. 30, no. 2, pp. 306–14, Feb. 2011. no. 8, pp. 1605–1619, Aug 2004. [25] JPEG Requirements and guidelines, ISO/IEC IS Std. 10 918-1, 1994. [48] Kakadu JPEG2000 implementation, available [Online]: http://www. [26] A. Sharma, P. Bautista, and Y. Yagi, “Balancing image quality and kakadusoftware.com, 2017. compression factor for special stains whole slide images,” Analytical [49] J. Nocedal and S. Wright, Numerical Optimization. Springer, 2006. Cellular Pathology, vol. 35, no. 2, pp. 101–106, 2012. [Online]. Available: https://www.springer.com/us/book/9780387303031 [27] D. D. Rhoads, N. F. Habib-Bein, R. S. Hariri, D. J. Hartman, S. E. [50] L. Galli and S. Salzo, “Lossless hyperspectral compression using KLT,” Monaco, A. Lesniak, J. Duboy, M. E.-S. Salama, and L. Pantanowitz, in IEEE International Geoscience and Remote Sensing Symposium, “Comparison of the diagnostic utility of digital pathology systems for IGRSS, vol. 1, Sep. 2004, p. 316. telemicrobiology,” Journal of Pathology Informatics, vol. 7, p. 10, 2016. [51] P. Hao and Q. Shi, “Matrix factorizations for reversible integer map- [28] L. Barisoni, J. P. Troost, C. Nast et al., “Reproducibility of the ping,” IEEE Trans. Signal Process., vol. 49, no. 10, pp. 2314–2324, NEPTUNE descriptor-based scoring system on whole-slide images and Oct. 2001. histologic and ultrastructural digital images,” Modern Pathology, vol. 29, [52] S. C. Pei and J. J. Ding, “Improved reversible integer-to-integer color no. 7, pp. 671–684, Jul. 2016. transforms,” in Proceedings of the IEEE International Conference on [29] V. Sanchez, F. Auli-Llinas, R. Vanam, and J. Bartrina-Rapesta, “Rate Image Processing, ICIP, Nov. 2009, pp. 473–476. control for lossless region of interest coding in HEVC intra-coding [53] “The Cancer Genome Atlas, National Cancer Institute, National Institute with applications to digital pathology images,” in IEEE International of Health,” Available [Online]: https://cancergenome.nih.gov/. Conference on Acoustics, Speech and Signal Processing (ICASSP), Apr. 2015, pp. 1250–1254. [54] H. Malvar and G. Sullivan, “Transform, scaling & color space impact of [30] V. Goyal, “Theoretical foundations of transform coding,” IEEE Trans. professional extensions,” ISO/IEC JTC1/SC29/WG11 and ITU-T SG16 Signal Process., vol. 18, no. 5, pp. 9–21, 2001. Q, vol. 6, 2003. [31] M. Hernández-Cabronero, F. Aulı́-Llinàs, V. Sanchez, and J. Serra- [55] S. S. Henrique S. Malvar, Gary J. Sullivan, “Lifting-based reversible Sagristà, “Transform Optimization for the Lossy Coding of Pathology color transformations for image compression,” pp. 7073.1 – 7073.10, Whole-Slide Images,” in Proceedings of the Data Compression Confer- 2008. ence, DCC, Mar. 2016. [56] A. Jiménez, G. Bueno, G. Cristóbal, O. Déniz, D. Toomey, and [32] M. Hernández-Cabronero, V. Sanchez, F. Aulı́-Llinàs, and J. Serra- C. Conway, “Image quality metrics applied to digital pathology,” SPIE Sagristà, “Fast MCT Optimization for the Compression of Whole-Slide Photonics Europe, vol. 9896, p. 98960S, 2016. Images,” in Proceedings of the IEEE International Conference on Image [57] S. Daly, Digital Images and Human Vision. MIT Press, 1993, ch. Processing, ICIP, Sep. 2016, pp. 2370–2374. The Visible Differences Predictor: An Algorithm for the Assessment of [33] L. Zhang, J. Chen, J. Sole, M. Karczewicz, X. Xiu, and J. Z. Xu, Image Fidelity, pp. 179–206. “Adaptive Color-Space Transform for HEVC Screen Content Coding,” [58] J. Lubin, Vision Models for Target Detection and Recognition. World in Proceedings of the Data Compression Conference, DCC, April 2015, Scientific, 1995, vol. 2, ch. A visual discrimination model for imaging pp. 233–242. system design and evaluation, pp. 245–357. [34] JPEG 2000 image coding system: Core coding system, ISO/IEC IS Std. [59] K. J. Kim, B. Kim, R. Mantiuk, T. Richter, H. Lee, H. S. Kang, J. Seo, 15 444-2, 2004. and K. H. Lee, “A Comparison of Three Image Fidelity Metrics of Dif- [35] I. P. Akam Bita, M. Barret, and D. T. Pham, “On optimal transforms in ferent Computational Principles for JPEG2000 Compressed Abdomen lossy compression of multicomponent images with JPEG2000,” Signal CT Images,” IEEE Trans. Med. Imag., vol. 29, no. 8, pp. 1496–1503, Processing, vol. 90, no. 3, pp. 759–773, 2010. Aug. 2010. 12 [60] Z. Wang, E. P. Simoncelli, and A. C. Bovik, “Multiscale structural Workshop on Quality of Multimedia Experience (QoMEx, Jul. 2009, pp. similarity for image quality assessment,” in Asilomar Conference on 53–58. Signals, Systems and Computers, vol. 2, Nov. 2003, pp. 1398–1402. [62] G. Bjontegaard, “Calculation of average PSNR differences between RD- [61] T. Richter, “On the mDCT-PSNR image quality index,” in International curves,” Tech. Rep. VCEG-M33 ITU-T Q6/16, Apr. 2001.