title: Orientation-Disentangled Unsupervised Representation Learning for Computational Pathology
authors: Maxime W. Lafarge Josien P. W. Pluim Mitko Veta
text: 
  
  [
    \begin@twocolumnfalse
  
  
  
  
  
  
  
  
  
  
  
  
  Orientation-Disentangled Unsupervised Representation Learningfor Computational Pathology
  
      
  Maxime W. Lafarge,
  Josien P.W. Pluim,
  Mitko Veta
      
  ==========================================================================================
  
  
  
  
  
  0.9
  
  
  Department of Biomedical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands 
  
  
  
  
  
  
  
  
  Unsupervised learning enables modeling complex images without the need for annotations.
  The representation learned by such models can facilitate any subsequent analysis of large image datasets.
  
  However, some generative factors that cause irrelevant variations in images can potentially get entangled in such a learned representation causing the risk of negatively affecting any subsequent use.
  The orientation of imaged objects, for instance, is often arbitrary/irrelevant, thus it can be desired to learn a representation in which the orientation information is disentangled from all other factors.
  
  Here, we propose to extend the Variational Auto-Encoder framework by leveraging the group structure of rotation-equivariant convolutional networks to learn orientation-wise disentangled generative factors of histopathology images.
  This way, we enforce a novel partitioning of the latent space, such that oriented and isotropic components get separated.
  
  We evaluated this structured representation on a dataset that consists of tissue regions for which nuclear pleomorphism and mitotic activity was assessed by expert pathologists.
  We show that the trained models efficiently disentangle the inherent orientation information of single-cell images.
  In comparison to classical approaches, the resulting aggregated representation of sub-populations of cells produces higher performances in subsequent tasks.
  
  
  
  
    \end@twocolumnfalse
  ]
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  
  Dimensionality reduction is an efficient strategy to facilitate the analysis of large image datasets by representing individual images by a small set of informative variables, which can be used in place of the original images.
  Unsupervised learning methods can be used to obtain such an informative latent representation from a given dataset without the need for expert annotations.
  For this purpose, popular unsupervised learning frameworks such as the Variational Auto-Encoder <cit.> or flow-based approaches <cit.> can be used to model a joint distribution between an image dataset and a set of latent generative factors.
  As these frameworks provide a posterior distribution over a space of latent variables, they enable the estimation of the latent factors of new previously unseen images, that can then be used for any subsequent task.
  
  
  However, irrelevant factors that affect the appearance of images but are independent of the factors of interest can get entangled in the learned representation <cit.>.
  These irrelevant factors can be treated as nuisance variables that affect the learned representation in an unpredictable way <cit.>, consequently perturbing any downstream analysis[We refer to downstream task/analysis to express any task/analysis performed on an image dataset for which a learned representation is used in place of original images.] performed on a distribution of generative factors.
  Therefore, there is a need for disentangling such nuisance variables from the informative generative factors of interest.
  
  
  In computational pathology, these nuisance variables are known to affect the generalization power of machine learning models.
  They affect the appearance of the images across slides, scanners and hospitals and can be associated with the inevitable variations in tissue slide preparation and scanner-dependent digitization protocols.
  
  In a supervised learning context, strategies were developed to filter-out such irrelevant factors from the learned representation; popular methods applied in computational pathology include: staining normalization <cit.>, staining/style transfer <cit.>, data augmentation <cit.>, domain-adversarial training <cit.> and rotation-equivariant modeling <cit.>.
  
  In this paper, we focus on a specific generative factor that can be considered as a nuisance variable in some specific tasks: the orientation of individual image patches.
  In digital pathology, the orientation of localized image patches in a dataset of WSIs is arbitrary in the sense that tissue structures are likely to be observed in any orientation, as opposed to natural images or organ-level medical images for which the orientation of the imaged objects is typically not uniformly distributed.
  
  We propose an unsupervised learning framework to model a partitioned latent space of generative factors in which specific independent latent variables either code for oriented or non-oriented (isotropic) morphological components of histopathology images.
  
  
  
  
    
  Motivation
  We identified several points that motivate the development of methods to handle nuisance variables for computational pathology in an unsupervised learning context:
  
  	
  	
    - Using an informative representation in place of large and complex images can reduce the computational cost and facilitate training of subsequent task-specific models.
  In particular, such a representation can be used to directly process Whole Slide Images (WSIs) via patch-based compression <cit.> or to represent bags that consist of a high number of image patches in two-stage multiple-instance-learning frameworks <cit.>.
  By removing irrelevant factors from the representation, such existing frameworks can be further improved.
  
  	
    - A representation learned without supervision can better conserve the extent of the morphological information of tissue images making it suitable for a wide range of potential downstream tasks.
  	This is opposed to using the representation of a supervised model that potentially discards information that is irrelevant for the task for which it was trained, but that might be relevant for other downstream tasks.
  
  	
    - Latent variable models equipped with a generative component enable visual inspection of the individual learned factors.
  This can support the interpretation of a model, as a tool to gain insights into the morphological factors that are predictive for a given downstream task.
  
  
  
  
  
  
  
  
  
  
    
  The proposed method
  To enable such a partitioning of the latent space, we leveraged the structure of 2-group convolutional networks <cit.> to build the components of a new VAE that produces a pair of rotation-equivariant and rotation-invariant embeddings of image patches.
  This means that instead of representing an image by a vector of scalar latent variables (z_0 in Figure. <ref> (a)), the proposed framework learns in parallel a vector of real-valued isotropic variables (z^ISO) and a vector of angular orientation variables (z^ORI, see Figure. <ref> (c)) that enable the disentanglement of the orientation information in images.
  
  In Figure. <ref>, we illustrated the orientation disentanglement property obtained with the proposed framework, in which the effect of varying each type of generative factor can be observed in the generated examples.
  The resulting structure of the learned representation is also illustrated in Figure. <ref>, in which we show that the translation between existing datapoints in the latent space along the oriented or isotropic dimenisons causes distinctive generative effects.
  
  The independence between the two types of variables is guaranteed by the structure of the 2-CNN-based auto-encoder.
  In the spirit of the original VAE framework, we propose an extension of the objective function so as to encourage the mutual independence between all the introduced latent variables.
  
  
  
  
  
  We made a comparative analysis of the proposed framework with a baseline VAE and well-established hand-crafted measurements in the context of histopathology image analysis.
  We trained and evaluated the models using a dataset of nucleus-centered images extracted from histological cases of 148 breast cancer patients, exposing a large variability of nuclear morphology to be learned. 
  We evaluated the quality of the unsupervised learned representation by training simple logistic regression models on downstream classification tasks. We compared the ability of the learned embedding to predict the pleomorphism grade and tumor proliferation grade associated to each case of a hold-out test set using multi-class ROC-AUC metrics.
  
  
  
    
  Contributions
  
  
  	
  	
    - To our knowledge, this is the first time that an auto-encoder is proposed to explicitly disentangle the orientation information in images by learning a 2-part structured representation consisting of rotation-invariant real-valued variables and rotation-equivariant angle variables.
  	
    - We propose to use 2-structured CNNs to generate latent variables with guaranteed equivariance/invariance properties.
  	
    - We show that such an unsupervised model can quantify nuclear phenotypical variation in histopathology images and that the learned representation can be used in dowstream analysis to predict slide-level target values.
  
  
  
  
  
  
  
  
  § RELATED WORK
  
  
  
  
  
    
  Representation Learning for Microscopy Image Analysis
  Automated quantification of morphological features of single-cell images has been a paradigm for comparing populations of cells in high-throughput studies across microscopy modalities and applications <cit.>.
  
  In particular, machine learning methods were proposed to learn representation directly from image data.
  Transfer learning methods use the internal representation of deep learning models that were trained with an independent dataset and task <cit.>.
  Weakly-supervised and self-supervised models use the representation of deep learning models that are trained with the data at hand but optimized to solve an auxiliary pretext task <cit.>.
  
  We argue that methods that rely on training a generative adversarial network and that exploit the inner feature maps of the discriminator network as a representation fall into this category as these feature maps do not necessarily correspond to generative factors.
  These methods rely on the hypothesis that deep learning models can learn generic features that will generalize to an external task without further assumptions.
  However, such generalization is not guaranteed for a difficult medical-related task whose domain is too different from the task and dataset that were used for the original training of the models.
  Also, these methods often rely on subsequent fine-tuning using the data at hand, and such a transferred learned representation cannot give any direct visual insight into the nature of the individual features.
  
  Generative models and in particular latent variable models, were developed to learn the generative factors of cell images as a representation to be used in downstream tasks. 
  These methods are based on (sparse) dictionary learning <cit.>, sparse auto-encoding <cit.>, variational auto-encoding <cit.>, conditional auto-encoding <cit.> or other auto-encoding frameworks <cit.>.
  
  
  
  
    
  Other Applications of Representation Learning for Computational Pathology
  We consider unsupervised representation learning of random patches of digital slides related work: such representations facilitated the achievement of downstream tasks <cit.>.
  However, when supervision is possible, patch-based representation can be achieved in a end-to-end fashion in a multiple-instance-learning framework for a given task <cit.>.
  In <cit.>, the authors compared different latent variable models in order to compress WSIs and enable their processing in a single run, and investigated their potential on downstream tasks.
  Studies on realistic generation of histopathology images <cit.> showed that decoder-based generative models can embed the fine-grained morphological structures of tissues in a low-dimensional latent space, which is in line with the motivation of our work.
  
  
  
  
    
  Structured Latent Variable Models
  Although unsupervised latent variable modeling is appealing, <cit.> and <cit.>, supported by the work of <cit.> showed that learning disentangled generative factors is not possible without constraining the learned representation.
  This argument justifies the limited performances of the learned representation of baseline VAE models in downstream tasks.
  As a solution, methods were proposed to structure the VAE latent space as a form of inductive prior: hyperspherical latent structure <cit.>, supervised nuisance variable separation <cit.> or domain-wise latent space factorization <cit.> for example.
  
  In the context of single-cell representation learning, <cit.> proposed a structured latent space via a conditional VAE model that encourages separation of cell/nuclear shape information from sub-cellular component localization. 
  
  The framework proposed in this paper is in the direction of research of these prior works but specifically address the spatial orientation of the generative factors of cell images.
  
  
  
  
    
  Rotation-Equivariance in Convolutional Networks
  Deep Learning methods were proposed to learn representations that are equivariant to the orientation of images.
  These methods enable learning a representation that changes in a deterministic way when the input image is rotated.
  In particular, group convolutional networks <cit.> extend standard CNNs by replacing the convolution operation.
  Advantages of using group-structured convolutional networks were shown on computational pathology tasks in a supervised training context <cit.>.
  
  Here we leverage the structure of 2-CNNs in an unsupervised context as a new way to structure the latent space of VAE-based models.
  
  
  
  
  
  
  
  § DATASETS
  
  
  To train and compare the models investigated in this study, we used one dataset for training purposes and to assess slide-level classification performances and another patch-based benchmark dataset to assess cell-level classification performances.
  
  
  
  
    
  TUPAC-ROI
  We used a dataset of 148 WSIs of Hematoxylin-Eosin-stained (H&E) tissue slices of breast cancer patients.
  These WSIs are part of the training set of the TUPAC16 challenge <cit.> and are originally provided by The Cancer Genome Atlas Network <cit.>.
  Three Regions of Interest (ROIs) were annotated by a pathologist to indicate tumor regions with high cellularity, that pathologists would typically select for cancer grading.
  Note that we used this subset of WSIs as it was the only subset for which ROIs were provided.
  
  
  <cit.> provided several patient-level metrics on these cases, including molecular and genetic information as well as expert-based visual morphological assessment (mitosis grading, tubular formation, pleomorphism grading).
  We used the pleomorphism grade and tumor proliferation grade (discrete grades in {1, 2, 3}) associated to each WSI as a target value to evaluate the quality of the learned representations investigated in this study.
  In the case of tumor proliferation grading, we make the assumption that this value can be associated to cell-level patterns and their local distributions.
  
  
  To reduce the inter-case staining variability, we pre-processed all images by applying the well-established staining normalization method described in <cit.>.
  We applied our internal nuclei segmentation deep learning model <cit.> within each ROI, and used the center of mass of the segmented instances as an estimate of the nuclei center locations.
  Image patches of size 68×68px^2 at a resolution of ∼0.25μm/px centered on these locations were extracted and constitute the effective dataset of cell-centered images we used to train and test our models.
  We made a training-validation-test split of this dataset (including respectively 104, 22 and 22 cases). See the supplementary material for details about the class distributions across the splits.
  We will refer to this refined dataset as TUPAC-ROI.
  	
  
  
  
  
    
  CRCHistoPhenotypes
  In order to assess the single cell-level quality of the trained models and investigate the transfer ability of the learned representation to other tissue types, we used the classification subset of the CRCHistoPhenotypes - Labeled Cell Nuclei Data (CRCHP) dataset provided by <cit.>.
  This dataset consists of 22,444 localized nuclei extracted from 100 ROIs, themselves originating from WSIs of H&E stained histology images of colorectal adenocarcinomas.
  Cell-type labels (epithelial, inflammatory, fibroblast, miscellaneous) were provided for each nucleus.
  We made a training-validation-test split of this dataset free of any ROI-overlap (including respectively 11,090, 3,133 and 8,221 nuclei).
  We resampled and cropped image patches centered at the nuclei locations so that the resolution and dimensions matched the ones of the TUPAC-ROI dataset.
  Finally, we applied the same staining normalization protocol as for the TUPAC-ROI dataset.
  
  
  
  
  
  
  
  
  § METHODS
  
  
  
  This section describes the proposed framework, first summarizing the baseline VAE framework, then presenting its expansion using 2,N-group convolutions and finally developing how we enable disentanglement of the orientation information.
  
  We formalize the representation learning problem in this generative modeling setting as the problem of learning the joint distribution p(x, z) of the observed images x with their latent generative factors z.
  Typically, we want to estimate the distribution that maximizes the marginal likelihood p(x) of this model for a given dataset.
  
  
  
  
   §.§ Variational Auto-Encoder
  
  
  In the VAE framework <cit.>, the likelihood of the observed images given a latent embedding p_ψ(x | z) is modeled by a decoder CNN with parameters ψ.
  It is assumed that the latent z are drawn from a given prior distribution p(z), typically a multivariate normal distribution.
  By introducing an approximation of the posterior on the latent q_ϕ(z | x) modeled by a CNN encoder (parameterized by ϕ), <cit.> propose to optimize ψ and ϕ by maximizing a tractable lower bound on the marginal log likelihood, as written in Eq.<ref>.
  
  
      .9!
      ℒ_ VAE(x, z; ψ, ϕ) =
      𝔼_q_ϕ(z | x) [ log p_ψ(x | z) ]
      -β·D_KL[q_ϕ(z | x)  ||  p(z)) ]
  
  The Kullback-Leibler divergence term D_KL, encourages the distribution of the sampled latents to be close to the prior distribution.
  The β hyper-pararameter controls the strength of this constraint as introduced by <cit.>.
  
  
  
  
    
  Orientation Encoding
  The encoder and decoder CNNs of conventional VAE models are built as a series of alternating trainable 2D convolution operations, non-linearity activation functions and down/up-pooling operations.
  For a given image x∈𝕃_2[ ℝ^2 ], the encoder CNN q_ϕ generates low-dimensional embedding samples z∼ q_ϕ(z | x) (of M elements) with z= [ z_i ]_i=1^M  and z_i∈ℝ.
  
  In the context of tissue imaging, we make the hypothesis that every image can be decomposed as a pair x = (x_0, θ_0) of independent variables such that p(x) = p(x_0) · p(θ_0).
  With this formulation, we assume the existence of a reference distribution of images x_0∼ p(x_0) that get rotated by an angle drawn from a uniform distribution θ_0∼ U(0, 2π).
  This uniformity assumption is logical for histopathology since tissue slices are prepared independently of their orientation and thus tissue images can be acquired in any possible orientation.
  
  We write x = ℒ_θ_0[ x_0 ] as the relationship between these variables, with ℒ_θ_0: 𝕃_2[ ℝ^2 ] →𝕃_2[ ℝ^2 ] the left-regular representation on 2D images of the rotation group 2, parameterized by θ_0.
  In this notation θ_0 indicates the action of a planar rotation R_θ_0∈2, such that ℒ_θ_0[ x_0 ](u) = x_0(R^-1_θ_0·u), given a vector location u∈ℝ^2.
  
  
  From this point of view, θ_0 is a generative factor of the observed images, and so it would be expected that the model learns to isolate this factor in the learned latent z.
  Without loss of generality, we can assume that an optimal model would learn to decompose the latent variables such that z = [z_0, θ_0] with z_0 the subset of latent variables that are independent of θ_0.
  
  This decomposition would imply that the posterior distribution of θ_0 is equivariant under the deterministic action of the rotation group on the image domain, and via the concatenation of rotations on the domain of orientations (written as the addition of angles modulo 2π in Eq. <ref>).
  The desired (conditional) independence relationship between z_0 and θ_0 is equivalent to the posterior being invariant under the same group actions (see relationship of Eq. <ref>).
  
  .9!
    
  
  	
  	
      p( θ_0 + θx = ℒ_θ[ x_0 ]) = p( θ_0x = x_0)  ; ∀θ∈ [0, 2π) 
      
      p( z_0x = ℒ_θ[ x_0 ]) = p( z_0x = x_0)  ; ∀θ∈ [0, 2π)
  
  
  
  
  However empirical experiments have shown that perfect independence of the latent variables is hard to achieve in a generative model and that such disentanglement of the generative factors is typically not obtained without further constraints on the models <cit.>.
  At that, the assumed uniform distribution of θ_0 is not encouraged by the Gaussian prior distribution of the formulation of the standard VAE.
  Therefore, we propose to consider a network architecture that explicitly encodes the orientation information of the images, guarantees the relationships of Eq. <ref>-b and enables the modeling of the Bayesian network illustrated in Figure. <ref>-c.
  
  
  
  
  
  
  
   §.§ 2 Variational Auto-Encoder
  
  
  Although conventional CNNs are equivariant to translations (since 2D convolutions are equivariant to planar translations), they are not guaranteed to be equivariant with respect to rotations of the input images. 
  
  
  
    
  Group Structured CNNs
  Group convolution operations were proposed to give CNNs the structure of the roto-translation group 2 := ℝ^2⋊2 <cit.>.
  The internal feature maps of CNNs with such group structure can be treated as 2-images F ∈𝕃_2[2] and the application of convolutional operations with 2-image kernels are equivariant under the action of the elements of 2.
  This architecture provides an end-to-end roto-translation equivariance property to CNNs.
  
  
  The architecture of a 2-CNN can be implemented by discretizing the sub-group 2,N := ℝ^2⋊2,N by sampling 2 with the elements that correspond to the N rotation angles of {2πn / N n=0…N-1} <cit.>.
  This way, the internal feature maps of the network can be implemented as tensors of shape H × W × N × M with H and W the size of the spatial dimensions, N the size of the discretized orientation-axis and M the number of channels in the layer.
  
  
  
    
  Application to the VAE framework
  We propose to replace the 2D-convolution operations of the conventional CNNs in the VAE framework by 2-group convolutions to yield rotation-equivariant mappings from the input images to the sampled latent variables and from the latent variables to the reconstructed images.
  This change of architecture also relies on the replacement of the first layer of the encoder by a lifting layer to produce 2-image representation maps; and the replacement of the penultimate layer of the decoder by a projection layer to output 2D images <cit.>.
  
  
  The bottleneck of a conventional CNN-based encoder corresponds to feature vectors of size M.
  Likewise, in a 2-CNN-based VAE, the feature vectors at the bottleneck of the encoder can be defined in terms of the rotation elements of the circle group 2 solely.
  We treat these feature vectors as 𝒰_g[ f ] where 𝒰_g is the left-regular group-representation of 2 on functions f ∈𝕃_2[2] with g = R_θ∈2 (we simplify this notation to 𝒰_θf).
  
  
  In practice, we consider the sub-group 2,N, so that the 2,N-vectors 𝒰_θf can be implemented as tensors of shape N×M with N the size of the discretized orientation-axis and M the number of latent variables.
  The difference in embedding structure between conventional VAEs and 2-CNN-based VAEs is illustrated in Figure. <ref>.
  
  
  
    
  2,N-Structured Latent Variables
  Instead of considering real-valued latent variables, we propose to model the latent variables as 2,N-vector-valued random variables Z.
  The group structure of the 2-modified encoder enables to model the posterior distribution q_ϕ(Zx) with the property of being equivariant under the action of 2,N.
  Thus, the modeled distribution verifies the relationship of Eq. <ref>.
  
  
      .9!
      q_ϕ(𝒰_θ[ Z ] x = ℒ_θ[ x_0 ]) = q_ϕ(Zx = x_0)  ; ∀θ∈SO(2,N)
  
  
  The 2-CNN decoder takes the samples Z as input and models the likelihood p_ψ(x | Z) as a multivariate Gaussian with identity covariance.
  Note that this 2-CNN-based VAE can be trained with the same objective as conventional VAEs (see Eq.<ref>) after adjusting the prior on the latent to a multivariate normal distribution that matches the dimensions of Z.
  
  
  
  
    
  Consequences for Downstream Analysis
  
  By construction, the equivariance property of the variational posterior q_ϕ guarantees that rotating the input images by an angle measure θ will cause a cycle-shift on the posterior distribution as expressed by the relationship of Eq. <ref>.
  Likewise, the rotational equivariance of p_ψ implies that cycle-shifting the values of a latent sample will cause a rotation of the reconstructed images.
  
  
  As a result, the generative process of the images does not depend directly on θ_0 anymore: this variable becomes encoded in the 2,N latent variables as a hidden shift on the orientation-axis.
  Still, the dependence of each variable Z_i,j to θ_0 makes downstream analysis of these generative factors subject to the variability of this arbitrary orientation within a dataset.
  
  
  
  
   §.§ Separation of Isotropic and Oriented Latents
  
  
  In order to achieve disentanglement of the orientation information in the latent variables, we make the hypothesis that the set of generative factors can be split in two sets of independent variables: 
  a set of real-valued variables z^ISO = [ z^ISO_i ]_i=1^M that codes for non-oriented/isotropic features in the images, and a set of angle variables z^ORI = [ z^ORI_i ]_i=1^M with values in [0,2π] that code for oriented structures in the images.
  
  
  To achieve such partitioning of the latent space, we design the 2,N-CNN encoder to approximate two posterior distributions q_ϕ(z^ISO | x) and q_ϕ(z^ORI | x) by producing three output components that parameterize these distributions (as illustrated in Figure. <ref>):
  
  	
  	
    - Two sets of 2,N-vectors that are projected via the mean operator along the orientation-axis resulting in a mean vector μ^ISO∈ℝ^M and a variance vector (σ^ISO)^2∈ (ℝ^+)^M. 
  	
    - One set of softmax-activated 2,N-vectors Q^ORI_i that correspond to discretized approximations of q_ϕ(z^ORI_i | x) as defined in Eq. <ref> with i=1 … M.
  	Here the softmax function is used to ensure that each vector Q^ORI_i represents a probability mass function.
   
  
      Q^ORI_i,j = ∫_(j-1)2π/N^j2π/N q_ϕ(z^ORI_i | x)  dz^ORI_i  ;   j=1 … N
  
  
  Finally, the z^ISO_i can be drawn from 𝒩(μ^ISO_i, (σ^ISO_i)^2) and the z^ORI_i can be directly drawn from q_ϕ(z^ORI_i | x) using the approximations Q^ORI_i (the implementation of these sampling procedures are detailed in the next paragraph).
  
  
  By construction, the μ^ISO_i and σ^ISO_i are rotation-invariant, and thus ensure that z^ISO verifies the invariance relationship of Eq. <ref>.
  The modeled posteriors q_ϕ(z^ORI_i | x) follow the equivariance relationship of Eq. <ref>, as θ_0 becomes encoded as a shared hidden shift (modulo 2π) across the variables z^ORI_i.
  
  
  
  
  
  
    
  Implementation Details on Latent Sampling
  During training, the stochastic sampling process of the z^ISO and z^ORI variables requires implementation via differentiable operations in the computational graph, to enable gradient back-propagation through the encoder.
  
  As proposed by <cit.>, we implemented the sampling process of z^ISO via the reparameterization trick to obtain sampled latents z^ISO = μ^ISO + ϵ^ISO·σ^ISO (as illustrated in Figure. <ref>).
  Here, ϵ^ISO∼𝒩(0,1) is an auxiliary noise variable to simulate sampling from a Gaussian distribution.
  
  
  To approximate the sampling process of z^ORI we implemented an Inverse Transform Sampling layer that returns a set of angle measures drawn from the distributions coded by Q^ORI.
  This layer computes a continuous inverse cumulative distribution Q^ -1_i for each vector Q^ORI_i and calculates each angle measure as z^ORI = Q^ -1_i(ϵ^ORI) with ϵ^ORI∼𝒰(0, 1) a uniformly distributed auxiliary noise variable.
  An example of this two-step procedure is shown in Figure. <ref>.
  
  
  So as to conserve the end-to-end equivariance property of the framework, the latent variables are reshaped as 2,N-vectors, such that the decoding procedure can occur in the same conditions as described in Section. <ref>.
  This is done by expanding and repeating the values of z^ISO along the orientation-axis,
  and via label smoothing to encode the sampled angles z^ORI_i into 2,N-vectors.
  
  
  
  
    
  Extended Objective
  Since the images are generated from two independent sources of generative factors we simply further developed the lower bound in the VAE formulation such that the proposed model can be trained by maximizing the loss written in Eq. <ref>.
  The constraints are computed by fixing the prior p(z^ISO) as a multivariate normal distribution and the prior p(z^ORI) as a multivariate uniform distribution on [0, 2π].
  
  .9!
    
  	
      ℒ_θ-VAE(x, z^ISO, z^ORI; ψ, ϕ)  =    𝔼_q_ϕ(z^ISO, z^ORI | x)[ log p_ψ(xz^ISO, z^ORI) ] 
         -β^ISO·D_KL[q_ϕ(z^ISO | x)  ||  p(z^ISO) ] 
         -β^ORI·D_KL[q_ϕ(z^ORI | x)  ||  p(z^ORI) ]
  
     
  
  
  
    
  Consequences for Downstream Analysis
  The isotropic generative factors z^ISO are guaranteed to be independent to θ_0 by construction, so they can be compared and aggregated within/across populations of tissue image patches independently of their spatial orientation.
  Likewise, the distribution of angles z^ORI in a given population characterizes the variability of oriented features independently of isotropic factors.
  
  
  
  
    §.§.§ Complementary Reconstruction Loss
  
  Conventional VAE models are known for generating/reconstructing blurry images of a lower quality than original images.
  Poor reconstructions might imply that the high-frequency details in the images do not get encoded in the latent representation and might entail poor performances in downstream tasks.
  This limited quality of generated images is often associated to the pixel-wise reconstruction term of the VAE objective (see Eq. <ref>).
  
  To ensure the reconstruction term enables generation of realistic images, we used the extension of the VAE objective that was described in <cit.>.
  This method is inspired by the initial work of <cit.> that introduces a discriminator CNN D with parameters χ in the VAE framework.
  
  This discriminator is trained to classify batches balanced between real images x∼ p(x) and reconstructed images x∼ p_ψ(x | z).
  The internal feature maps of D for a given input image x are defined as D_i(x) with i a layer index.
  The D_i(x) can be considered as additional Gaussian observations with identity covariance drawn from p_χ(D_i(x) | z).
  We thus define extra reconstruction losses ℒ^D_i that we use to complement the model objective as defined in Eq. <ref> where γ is a weighting hyper-parameter.
  
  
  .9!
    
  	
      ℒ_ VAE+(x, z; ψ, ϕ, χ)  = ℒ_ VAE(x, z; ψ, ϕ) + γ∑_iℒ^D_i(x, z; ψ, ϕ, χ) 
         ℒ^D_i(x, z; ψ, ϕ, χ) = 𝔼_q_ϕ(z|x)[ log p_χ( D_i(x) z ) ]
  
     
  
    
  The training procedure of the full model pipeline consists of alternating training steps between updating ψ and ϕ by maximizing ℒ_ VAE+ and updating χ by minimizing the cross-entropy loss on the predictions of D.
  
  This extension of the VAE framework has the benefit to keep the rest of the model formulation intact and was shown to be beneficial in autoencoder-based applications involving cell images <cit.>.
  The discriminator is used at training time only and is discarded at inference time, which thus does not cause any computational overhead for downstream analysis.
  All the models investigated in this paper were extended by this method.
  
  
  
  
  
  
  
  § EXPERIMENTS
  
  
  In this section we detail the network architectures we designed to implement the proposed orientation-disentangled VAE (Section. 
  <ref>) and a comparable baseline VAE model (Section. <ref>).
  We also describe their training procedures and the evaluation protocols we applied to compare the quality of the learned representation and gain insights into the effect of the disentanglement and newly introduced hyper-parameters.
  
  
  
  
   §.§ Model Architectures
  
  We designed the conventional CNN encoders and SE(2,N)-CNNs encoders of the VAEs as straight-forward sequences of four blocks that each consists of a convolutional layer, a batch normalization layer (BN), a leaky reLU non-linearity and a max-pooling layer.
  
  In the case of the 2,N-CNNs, conventional ℝ^2-convolutional layers were replaced by 2,N convolutional layers and BN layers were adjusted to include the orientation-axis in the computation of the batch statistics.
  We also introduced intermediate projection layers similar to the locally rotation-invariant CNNs proposed by <cit.> within all the hidden layers in order to reduce the computational cost of the network.
  
  In order to generate embedding samples at the bottleneck of the models during training, we implemented the computational sampling procedures presented in Section. <ref> for the conventional VAEs and in Section.<ref> for the orientation-disentangled VAEs.
  
  The decoder networks correspond to mirrored versions of their encoder counterparts via the use of up-sampling layers and transposed convolutions.
  We included a mean-projection layer at the end of the 2,N-CNN decoder to project the features maps on ℝ^2.
  Finally we added an extra 1×1-convolutional layer to the decoders to output images whose dimensions match the input dimensions.
  
  For all the networks, we used kernels with spatial size 5×5 so as to enable proper rotation of the 2,N-kernels as described in <cit.>.
  We fixed the angular resolution of the 2,N layers to N=8 as it was previously shown to give optimal performances <cit.>.
  To ensure fair comparison of the models, we balanced the number of channels in each layer such that the total number of weights between the two types of VAEs is approximately the same.
  In order to have a fair comparison,, we also fixed the total number of variables in the bottleneck of all the models (size of z is 64 in the conventional VAEs and z^ISO and z^ORI are both of size 32 in the orientation-disentangled VAEs).
  
  More details about the model architectures can be found in the Appendix.
  
  
  
   §.§ Training Procedures
  
  All the models were trained using the training set of the TUPAC-ROI dataset described in Section. <ref>.
  We used mini-batches that consist of 35 image patches, such that the distribution of the WSIs of origin within each batch was approximately uniform.
  We used the Adam optimizer to update the weights of the encoders and decoders (learning rate 0.001, β_1=0.9, β_2=0.999), and Stochastic Gradient Descent with momentum to update the weights of the discriminator (learning rate 0.001, momentum 0.9).
  All convolutional kernels were regularized via decoupled weight decay with coefficient 1×10^-4.
  We stopped the training process after convergence of the loss on the validation set (approximately 20×10^3 iterations).
  The weighting coefficient of the disciminator-based reconstruction loss was fixed to γ=0.01 across all experiments.
  
  In order to assess the effect of the weighting coefficient of the prior constraint as it was evidenced by <cit.> we trained the models with varying values of β, β^ISO and β^ORI in { 0.1, 0.5, 1.0, 2.0, 4.0 }.
  
  
  
   §.§ Downstream Analysis
  
  In order to compare the quality of the learned representation at the bottleneck of the different trained VAE models, we trained subsequent logistic regression models that take the learned representation as input to solve downstream tasks.
  
  We investigated two types of downstream tasks to evaluate the value of the learned representation in different contexts: patient-level classification tasks (predicting the pleomorphism grade and tumor proliferation grade of a given WSI) using the TUPAC-ROI dataset and single-cell-level classification tasks (predicting the cell-type of a given nucleus) using the CRCHP dataset.
  
  Each downstream task was investigated with different types of representation of the data points based on the different latent variables estimated in the VAEs (see z in Section. <ref>, z^ISO and z^ORI in Section. <ref>.
  For comparison, we also considered well-established nuclear morphometric measurements for comparison (mean nuclear area and standard deviation of the nuclear area) as well as combinations of representations.
  
  
  
    
  Patient-level classification tasks: we first aggregated the representation of all nucleus-centered image patches within every ROI of the dataset (see Section. <ref>).
  For the representations corresponding to z and z^ISO, we considered the mean μ and μ^ISO of the estimated posterior distributions as cell-level embeddings, and we further computed ROI-level embeddings by computing the mean of all the embeddings obtained within a ROI.
  We then trained a single-layered multi-class logistic regression model to minimize the cross-entropy loss given the WSI-level ground-truth labels.
  
  For the representation corresponding to z^ORI, we used their discrete distribution at the ROI-level as an aggregated representation (see Q^ORI in Section. <ref>).
  Likewise, we trained a two-layer multi-class logistic regression model with an intermediate maximum projection layer to ensure rotational-invariance of the predictions.
  
  
  
    
  Single-cell-level classification task: we directly used the mean μ and μ^ISO of the posterior distributions on z and z^ISO as representation of nucleus-centered image patches and trained a single-layered multi-class logistic regression model to minimize the cross-entropy loss given cell-type ground-truth labels.
  
  All logistic regression models were regularized with L2-weight decay and the associated coefficient was fine-tuned on the validation set.
  
  
  
  
  
  
  
  
  
  
  
  
  § RESULTS
  
  
  This section details the evaluation protocols and metrics we used to assess the quality of the learned representation we obtained with the methods presented in Section. <ref>.
  We present the downstream performances for each type of representation on the tasks described in Section. <ref>.
  
  
  
  
  
  
   §.§ Slide-Level Downstream Performances
  
  We evaluated the trained multi-class logistic regression models on the hold-out test set described in Section. <ref>.
  For each model, we considered a set of binary classifiers that correspond to pairwise comparisons of each class against the other classes.
  For each set of such binary classifiers, we computed a set of Receiving Operating Curves (ROCs) based on the model predictions.
  Within each set of ROCs, we computed the corresponding Areas Under the Curve (AUC) and the mean of AUCs (mAUC) across that set to summarize the AUC metric given the multi-class setting at hand.
  
  To assess the robustness of the learned representations given perturbations of the training data, we resampled the training set ten times, re-trained the models and reported the mean and standard deviation of the mAUCs across these repeats.
  
  
  
  
    
  Pleomorphism Grading
  The results for the pleomorphism grade prediction task are summarized in Figure. <ref>.
  We report an improvement of the mAUC of 0.038 from the isotropic learned representation in comparison to the learned representation of the baseline VAE.
  The isotropic learned representation performed better than the oriented learned representation for all the tested values of β^ISO.
  We obtained a consistent additional improvement of performance when combining the isotropic learned representation with segmentation-based features.
  
  
  
  
  
  
    
  Tumor Proliferation Grade Prediction
  The results for the tumor proliferation grade prediction task are summarized in Figure. <ref>.
  We do not report any significant improvement of the mAUC from using isotropic or oriented learned representation in comparison to the learned representation of the baseline VAE and to segmentation-based features.
  We report a consistent additional improvement of performance when combining the isotropic learned representation with segmentation-based features.
  
  
  
  
   §.§ Cell-Level Downstream Performances
  
  We used the same protocol to evaluate the cell-type classification models.
  The results for the cell-type classification task are summarized in Figure. <ref>.
  We report a consistent improvement of the mAUC from using the isotropic representation in comparison to the learned representation of the baseline VAE.
  The oriented learned representation performed worse than the representation of the baseline VAE or the isotropic learned representation for all the tested values of β and β^ISO.
  
  
  
  
  
  
  
  
  
  § DISCUSSION AND CONCLUSIONS
  
  
  In this study, we proposed a novel rotation-equivariant variational auto-encoder framework that learns two types of generative factors: isotropic real-valued components and oriented angular components.
  We showed that this two-fold low-dimensional structure can efficiently represent histopathology images.
  We investigated in a controlled experimental setup, the predictive power of the learned representation using the proposed frameworks and show its advantage in comparison to unsupervised baseline counterparts.
  The difference of generative action of each type of variable was qualitatively demonstrated via smooth transitions of generated examples given interpolated embeddings (see Figures. <ref>-<ref>).
  
  
  Qualitatively, we observed that the isotropic learned representation captures morphological factors such as stain ratios, nuclear sparsity, thickness of the nuclear boundary (see Figure. <ref>(b)) whereas the oriented learned representation codes for the radial location of the surrounding objects (non-centered neighboring nuclei) and asymmetric structures (see Figure. <ref>(a)).
  
  
  Quantitatively, using isotropic representation was always better or as good as the representation learned by conventional VAE or segmentation-based nuclear area features (see Section. <ref>).
  This is in agreement with our hypothesis that the orientation information that is entangled in the representation learned by the baseline VAE affects the quality of the aggregated representation, and subsequent downstream performances.
  This was both observed for patient-level grading tasks based on aggregated representation and for cell-level classification tasks.
  This is also in line with previous results reported in studies on supervised rotation-invariant CNN models trained to solve computational pathology tasks <cit.>.
  
  
  Although our results suggest that unsupervised learned representation is an efficient alternative to hand-crafted feature-based representations, the fact that the combination of hand-crafted nuclear feature representation with unsupervised learned representation gave a consistent improvement of performances also reveals the limitations of the proposed framework.
  Indeed, complex knowledge-based quantities such as the mean nuclear area that were relevant for the task at hand could not be extracted from the learned representation as they do not necessarily correspond to independent generative factors that the models learned.
  However we conjecture that this limitation might be due to the restricted architectural design of the model that we chose for the comparative analysis.
  We thus believe that this limitation can be potentially overcome by using more complex and sophisticated architectures for the latent variable models and downstream classification models.
  
  The performances achieved on the pleomorphism grading task using the isotropic learned representation indicate that these features were more predictive than the oriented features to solve this specific task.
  This is expected as by definition the pleomorphism grade was labeled based on rotation-invariant nuclear morphological factors.
  However on the tumor proliferation grading task, oriented features were slightly more predictive than isotropic features.
  This suggests that these two types of features were equally informative of this specific patient-level grade.
  
  Besides these applications, the proposed framework can be used as a generic tool to quickly gain insights, and with minimal training effort, into the slide-level predictive value of fixed-scale image patches.
  Indeed, we showed that logistic regression models could be trained using aggregated representation of cell populations to predict patient-level target values such as the pleomorphism grade and tumor proliferation grade.
  But the same approach could potentially be applied to estimate any other slide-level value.
  Also, extensions to other convolutional latent variable models are possible including more complex architectures (such as flow-based generative models), other families of variational distributions (such as structured posterior distributions of the oriented variables) and other training paradigms (such as in a semi-supervised framework).
  For computational convenience, we proposed implementing the sampling of the angle variables by means of a straight-forward inverse sampling layer, however, other end-to-end sampling strategies could be investigated for further improvement.
  This framework is also transferable to other problems in which one wants to model a posterior distribution on the rotation group.
  Other interesting applications for future work include pre-training for patch-based classification tasks, and compression of WSIs.
  
  
  
  
  
  
  
  § REFERENCES
  
  
  
  unsrtnat
  
  
  
  
  
  
  
  
  
  [pages=1-4]./appendix.pdf
  
  
  