Large-scale analysis of neuroimaging data on commercial clouds with content-aware resource allocation strategies

Abstract

The combined use of mice that have genetic mutations (transgenic mouse models) of human pathology and advanced neuroimaging methods (such as magnetic resonance imaging) has the potential to radically change how we approach disease understanding, diagnosis and treatment. Morphological changes occurring in the brain of transgenic animals as a result of the interaction between environment and genotype can be assessed using advanced image analysis methods, an effort described as ‘mouse brain phenotyping’. However, the computational methods involved in the analysis of high-resolution brain images are demanding. While running such analysis on local clusters is possible, not all users have access to such infrastructure and even for those that do, having additional computational capacity can be beneficial (e.g. to meet sudden high throughput demands). In this paper we use a commercial cloud platform for brain neuroimaging and analysis. We achieve a registration-based multi-atlas, multi-template anatomical segmentation, normally a lengthy-in-time effort, within a few hours. Naturally, performing such analyses on the cloud entails a monetary cost, and it is worthwhile identifying strategies that can allocate resources intelligently. In our context a critical aspect is the identification of how long each job will take. We propose a method that estimates the complexity of an image-processing task, a registration, using statistical moments and shape descriptors of the image content. We use this information to learn and predict the completion time of a registration. The proposed approach is easy to deploy, and could serve as an alternative for laboratories that may require instant access to large high-performance-computing infrastructures. To facilitate adoption from the community we publicly release the source code.

Keywords

High-performance computing cloud neuroimaging MRI phenotyping computer vision image analysis resource allocation machine learning

1. Introduction

The introduction of cloud computing has affected how high-performance-computing resources are purchased, allocated, and managed. While previously access to such resources required either the purchase of advanced clusters or the participation in grid services, with commercial clouds, super-computing-level resources are available to everyone. Although initially they were targeted at enterprise users, the use of commercial clouds for research applications is gaining momentum. This transition is expected to facilitate research where large data repositories need to be processed (‘big data’). An application domain that will benefit from this paradigm change is the large-scale analysis of imaging data. Images arising from biological or medical domains (e.g. magnetic resonance imaging (MRI) of the brain) are characterized by high acquisition cost, large resolution, and complex analysis pipelines (i.e. the sequence of image-processing steps). However, since the same pipeline is applied to each dataset they are ideally suited for large data parallelization.

The advent of MRI methods for the study of the brain has significantly advanced our understanding of how this important organ functions in health and disease. In typical studies, subjects are divided in two or more groups, based on experimental design. Each subject is then scanned (possibly several times) with a particular scanning protocol to obtain one or several imaging volumes that reflect the underlying anatomy of the brain or brain tissue status or composition. At the next phase, imaging data must be analyzed to extract several features that improve our understanding of the effects of specific genetic variations or disease on the brain.

An example of an analysis is to identify differences in the volume of specific brain regions among the subjects and groups from brain MRI volumes. While an expert can manually delineate these structures in each subject’s volume, this process is highly time-consuming (six to ten hours depending on imaging resolution and number of structures) and introduces observer variability. Thus, computational solutions that automatically segment the anatomy within the brain are highly desirable.

One approach that has been particularly popular in the literature is a registration-based atlas propagation (Figure 1). Assuming the existence of an atlas $(A_{1}^{(I)}, A_{1}^{(L)})$ , where $A_{1}^{(I)}$ is an intensity volume and $A_{1}^{(L)}$ is the corresponding labeling (i.e. its anatomical segmentation) previously delineated by an expert, the labels can be propagated to another intensity volume $S_{1}^{(I)}$ , if a mapping between the two intensity volumes is known. The result of such a process is the labeled subject $S_{1}^{(L)}$ . An intensity-based registration finds a linear or, better, a non-linear mapping between two volumes using only voxel intensities. Generally, for one atlas and $N$ subjects ( $S_{n}^{(I)}$ , for $n = 1, \dots, N$ ) this process requires $N$ registrations. Several processing platforms implement this approach, for example ANTs (http://stnava.github.io/ANTs/), Freesurfer (http://surfer.nmr.mgh.harvard.edu/), and LONI (http://www.loni.ucla.edu/Software/).

Recently, it has been suggested that it is best to use multiple atlases when performing a registration-based label propagation (Heckemann et al., 2006) to account for variability between the subject population and the atlas. In this context several registrations are performed and the best atlas is identified, usually the one most similar to the subject, to propagate the labels. Alternatively, all atlases are used in a fusion-based segmentation approach (Heckemann et al., 2006; Jia et al., 2012). In recognition that obtaining multiple atlases is demanding and that atlases may not always reflect population variability, several intermediate templates are also used. These templates are either previous subjects, or even synthetic MRI volumes created by learning the variability between previously processed subject populations (Jia et al., 2012). As Figure 2 illustrates, in this case all subjects $S_{n}^{(I)}$ are registered using multiple templates $T_{p} = (T_{p}^{(I)}, T_{p}^{(L)})$ , for $p = 1, \dots, P$ , which are registered to multiple atlases $A_{m} = (A_{m}^{(I)}, A_{m}^{(L)})$ , for $m = 1, \dots, M$ , using the intensity component. The final anatomical labels ${\tilde{S}}_{n}^{(L)}$ for each subject are obtained by fusing the results of all registration-based propagations.

As the number of subjects ( $N$ ), templates ( $P$ ), and atlases ( $M$ ) increases, the complexity increases (in total $P \times (M + N)$ non-linear registrations are required). Depending on the non-linear registration algorithm, the complexity of this process increases even further. Current state-of-the-art algorithms that provide highly accurate, diffeomorphic, and invertible non-linear registrations between volume pairs are increasingly complex since they optimize models with tens of millions of degrees of freedom while they reach cubic theoretical complexity (Klein et al., 2009).

It is clear that this analysis is computationally demanding and is just an example of the many complex neuroimaging analysis pipelines used throughout the community (Kim et al., 2009; Antoniu et al., 2010). Traditionally, the majority of such analyses have occurred in single powerful workstations with limited throughput. Recently, since non-linear registration of volume pairs (as well as the majority of image analysis pipelines) is suitable for efficient data parallelism, the use of computing clusters to achieve higher throughput has been proposed (see the reviews of Dinov et al., 2010; Frisoni et al., 2011).

In recognition that some researchers may not have access to these options, the scientific community reacted by establishing consortia (e.g. LONI, neuGRID, outGRID, CBRAIN) that aim at making distributed (grid) computing available for neuroscience and neuroimaging research (Dinov et al., 2010; Frisoni et al., 2011). These grid-level services usually provide a user-friendly graphical interface for designing and putting together any image analysis pipeline, with respect to both tools and data. However, for executing this process seamlessly, access to the computer grid associated with the consortia underwriting these efforts is assumed.

Unfortunately, users that have limited or no access to powerful computer grids or clusters are left with limited computational options. Even users that do have access to such consortia or local clusters can in some cases benefit from having additional computational capacity (for example, increasing throughput to meet a deadline). The users can configure and install the underlying computational framework of LONI on another grid environment or the cloud (e.g. they can install the NITRC environment on Amazon EC2, http://www.nitrc.org/plugins/mwiki/index.php/nitrc:User_Guide_-_NITRC_Computational_Environment), but they have to export the image analysis pipeline as a virtual machine and take care of the details of its execution and the management of the cloud resources. To manage the workflow there exist several systems (Soma-workflow (Laguitton et al., 2011), Megha Workflow Management System (Pandey et al., 2009), Gridbus (Buyya and Venugopal, 2004), Cloudbus (Buyya et al., 2009), Aneka (Vecchiola et al., 2009), CometCloud (Kim and Parashar, 2011)), and several demonstrations of those for imaging applications are available (Kim et al., 2009; Antoniu et al., 2010). However, managing such systems (e.g. installation on the cloud environment, and configuration of services) and deployment of jobs may prove challenging and time-consuming for most users. Admittedly, these limitations might explain why cloud computing has not been adopted by the broad neuroscience and neuroimaging community for compute-intensive analyses.

In this paper, we use a commercial cloud computing platform (PiCloud, http://www.multyvac.com/) to perform one such intensive analysis procedure: a registration-based multi-atlas, multi-template anatomical segmentation for estimating the size of anatomical brain regions on the basis of magnetic resonance images. Here we show that when PiCloud, which is being developed to facilitate transparent access to commercial cloud infrastructures, is combined with sophisticated analysis and resource optimization methods, it offers an appealing alternative for several laboratories and types of users. Our approach:

relies on a general-purpose cloud environment and not on dedicated grids and clusters;

utilizes the PiCloud platform, which offers flexible, simple, and transparent allocation of computational resources via a Python application programming interface;

is based on Python, which is simpler for building complicated image analysis pipelines than other solutions (e.g. using C or C++ programming languages) as others have also found (Millman and Brett, 2007);

is sufficiently easy to be adopted and extended by non-computer-experts;

allocates cloud resources efficiently to minimize monetary cost; and

incorporates a method to learn and predict execution time by taking into account the image content to facilitate resource allocation.

We demonstrate our findings using data from mouse phenotyping experiments, with the goal of identifying morphological differences between two different population groups. Our cloud-based neuroimaging pipeline identifies such differences in several anatomical parts of the mouse brain. Furthermore, using the same imaging data and execution times we show that it is possible to reduce monetary cost by intelligently allocating resources. This paper builds upon our previous work (Minervini et al., 2012), where we first propose the use of commercial clouds (PiCloud) for neuroimaging.

The rest of this paper is organized as follows: in Section 2.1 we outline the proposed image analysis pipeline, and in Section 2.2 we describe the proposed cloud implementation. Section 2.3 details the proposed framework for intelligent resource allocation on the cloud. Section 3 outlines implementation details and presents and discusses our experimental findings, while, finally, Section 4 offers conclusions and directions for future work.

2. Proposed methods

2.1. Registration-based brain anatomical segmentation pipeline

To demonstrate the potential of high throughput processing on the cloud we implement a typical example of a compute-intensive analysis: a multi-atlas, multi-template anatomical segmentation such as the one illustrated in Figure 2. We rely on a highly accurate non-linear registration algorithm to match a moving volume (e.g. a subject $S_{n}^{(I)}$ ) and a fixed volume (e.g. a template $T_{p}^{(I)}$ ). This registration finds a diffeomorphic transformation that is symmetric, invertible, topology-preserving, and differentiable (Avants et al., 2008; Klein et al., 2009). This allows us to propagate information from the fixed to the moving volume and vice versa in an accurate manner. Briefly described, the goal of non-linear registration is to find the optimal transformation, $ϕ$ , within a specified transformation space which maps each location $x$ of volume $S (x)$ to a location in volume $T (z)$ such that a specified cost function relating the similarity between $S$ and $T$ is minimized. Recently, diffeomorphisms ( $ϕ$ ) have been proposed that are differentiable, invertible, and topology-preserving. With diffeomorphic registration, the following variational function is optimized:

{v_{1}^{*}, v_{2}^{*}} = \underset{v_{1, 2}}{argmin} {\int_{0}^{0.5} | | L v_{1} (x, t) | |^{2} d t + \int_{0}^{0.5} | | L v_{2} (x, t) | |^{2} d t + λ \int_{Ω} Π_{~} (S ° ϕ_{1} (x, 0.5), T ° ϕ_{2} (x, 0.5)) d Ω}

(1)

where $v (x, t) = v_{1} (x, t)$ when $t \in [0, 0.5]$ and $v (x, t) = v_{2} (x, 1 - t)$ when $t \in [0.5, 1]$ , $Π_{~}$ is a volume similarity metric, $L$ is an appropriate norm, $λ$ is a regularizer, and $Ω$ is the domain of the diffeomorphism. While finding the optimal transform is complex and ill-defined, under certain assumptions, symmetric diffeomorphic solutions are obtained numerically using greedy optimization strategies (Avants et al., 2008). To increase the probability that the greedy strategy finds a global minimum and to reduce the number of necessary iterations, a multi-resolution strategy using up to $LV$ levels of the Gaussian pyramid is adopted. A transformation is found first for an under-sampled version of the original data, and then is refined at each step while reaching the final full resolution.

The results of each registration are linear (affine) registration parameters as well as forward (subject-to-template, or template-to-atlas) and inverse (template-to-subject, or atlas-to-template) warp fields (essentially voxel-to-voxel mappings) (Avants et al., 2008). Using the inverse warp field and affine parameters obtained we estimate $S_{n}^{(L)}$ using the given labeled template $T_{p}^{(L)}$ , denoted by $S_{n}^{(L)} | T_{p}^{(L)}$ . The same process can be applied to obtain $T_{p}^{(L)} | A_{m}^{(L)}$ . Technically speaking we do not perform these steps independently, that is, first propagating from atlas to template and then from template to subject, but we concatenate the transforms to minimize interpolation errors.

Finally, we combine all intermediate segmentations ( $S_{n}^{(L)} | T_{1}^{(L)}, \dots, S_{n}^{(L)} | T_{P}^{(L)}$ ) by a fusion process to obtain the final result ${\tilde{S}}_{n}^{(L)} = ⋃_{p = 1}^{P} S_{n}^{(L)} | T_{p}^{(L)}$ , implemented as a majority vote algorithm that maintains for each voxel the most populous label across the $P$ candidates (Heckemann et al., 2006).

2.2. Proposed cloud implementation

To translate our image analysis pipeline to the cloud we use the PiCloud platform. This platform provides the user with transparent access to cloud resources currently supplied by Amazon Web Services (AWS). A Python library and middleware (implemented as an API) enables the use of cloud storage and computational power directly via the Python interpreter. It relieves the user from dealing with several issues such as virtual servers, connections, resource scaling, job scheduling, and billing, which may prevent non-computer-experts from using cloud computing services.

Among the types of service offered by PiCloud (as of July 2013), we consider the following:

c1 core type, offering 1 compute unit, 300 MB of memory, and low I/O performance at US$0.05/hour;

c2 core type, offering 2.5 compute units, 800 MB of memory, and moderate I/O performance at US$0.13/hour;

f2 core type, offering 5.5 compute units, 3.7 GB of memory, and moderate I/O performance at US$0.22/hour;

m1 core type, offering 3.25 compute units, 8 GB of memory, and high I/O performance at US$0.30/hour,

where core types are similar to what Amazon refers to as instances, a compute unit refers to the AWS definition of equivalent CPU capacity of a 1.0–1.2 GHz 2007 Opteron or 2007 Xeon processor, and I/O performance refers to internal communication within the cloud infrastructure. Each core type can be used in two modes: on-demand or real-time core provisioning. PiCloud charges on a per millisecond rate and not for the full hour when using on-demand cores.

The PiCloud service can provide a virtually unlimited number of cores, hence not imposing any strict limits on the total number of potential jobs parallelized simultaneously. This, coupled with the fact that the jobs are self-contained, translates into a high speedup of the overall process.

Subject to resource availability and internal proprietary resource management, any job executed on demand on PiCloud is placed in a queue, incurring a variable (unknown to the user) waiting time. Allocating a priori a job for real-time operation guarantees immediate execution, with a fixed waiting time. If the usage meets a certain per hour minimum ( $< 60 %$ of the hour), the real-time allocation is free of charge. Data storage and download incur a cost, whereas data uploading and intra-cloud communications (between storage and node) are free of charge.

We configure a Ubuntu GNU/Linux virtual machine (actually a Linux Container) following the PiCloud instructions for setting an environment to include customized libraries and applications required for our operation. We optimize (by enabling multi-threading) and compile the ANTs toolkit, which offers registration (ANTS), warping (WarpImageMultiTransform), and label fusion (ImageMath), and other stand-alone applications. Furthermore, NeuroDebian ( http://neuro.debian.net/ ) repositories have also been included to facilitate installation of other neuroimaging applications (Frisoni et al., 2011).

We use Python scripts to execute the workflow illustrated in Figure 3 for a given set of input volumes (subjects, templates, and atlases). More specifically:

We use the PiCloud file storage interface, in particular the cloud.files.put function, to upload the input data (i.e. subjects, templates, and atlases) to the cloud from our local workstation.

We use Python wrapper functions, an abstract example of which is shown in Listing 1, that execute the components of the image analysis pipeline for each input dataset: they retrieve the input data (subjects, templates, or atlases) from cloud storage to a node, run the applications (registration or fusion) composing the pipeline as subprocesses, and finally save the output files on cloud storage.

A single Python instruction of the PiCloud library (cloud.map) executes (maps) the functions for each dataset pair (e.g. an atlas and a template, or a template and a subject), specifying the type of computational resource (core type, on demand or real-time) and custom environment to use. This function essentially implements the core of the data parallelism that we employ here.

We retrieve the output data (e.g. labeled subjects ${\tilde{S}}_{n}^{(L)}$ ) from the cloud, using the cloud.files.get function of the PiCloud file storage interface.

Once all jobs have completed (registrations and fusions), we download to our local workstation all final labeled volumes ( ${\tilde{S}}_{n}^{(L)}$ ) for subsequent statistical analysis. The web interface and the APIs available by PiCloud provide several statistics (execution times), metadata, and output logs on each job, which are collected for analysis and comparisons.

To execute jobs on the PiCloud system (via the cloud.map arguments) we must identify the type of service (real-time or on demand) and the core type (c1, c2, f2, or m1) among the function’s arguments. Each option involves different charges and maintains different guarantees for service. One can assign all jobs to the same service type for simplicity, however, this is far from being an optimal allocation. While for certain cases assigning the type of core is trivial (e.g. the c1 does not have enough memory to execute a registration), when several core types satisfy the requirements a resource allocation methodology is necessary. Several approaches exist to assign resources (and schedule them if needed; see Lin and Lu, 2011; Mao and Humphrey, 2011) by satisfying various performance criteria, such as (1) cumulative, makespan, or minimum execution time or cost or deadline (Garg et al., 2010; Liu et al., 2010; Ramakrishnan et al., 2011; Reynolds et al., 2011; Tabak et al., 2013); (2) constrained formulations such as minimizing cost under time constraints (Van den Bossche et al., 2010; Kumar et al., 2011); and (3) allocating resources under budget constraints (Shi and Hong, 2010). Nevertheless, all options essentially care for the proper assignment of each job to a core and service type. Note that we do not refer here to cost optimization of renting services from a mix of providers (Tordsson et al., 2012).

More critically, most approaches do require an estimate of the execution time for each job and application to make efficient and accurate allocation and scheduling decisions (Kim et al., 2009; Garg et al., 2010). Typically, these estimates are obtained by analytical modeling of the underlying source code (which is not always available), or using empirical and historical data of each job and application type (Albers et al., 2009; Kim et al., 2009; Garg et al., 2010; Matsunaga and Fortes, 2010). However, when the job to be executed (in our case a registration) is homogeneous (i.e. the same job and parameters are executed on different data of the same size), we would have to rely on an average estimate of execution time (Smith et al., 1998; Kapadia et al., 1999; Li et al., 2009). Unfortunately, this option underperforms when the execution time depends on the actual data content (not only their size) and involves non-linear iterative components. Image registration is a perfect example of a highly non-linear process that iterates until a numerical minimum is found and thus is highly dependent on the input volumes. In fact, the expected convergence time varies largely among volume pairs (Kim et al., 2009).

Thus, in order to intelligently allocate resources in an environment that is composed of similar jobs but where the data content changes, an approach that utilizes this data dependence to more accurately estimate execution time is necessary.

2.3. Predicting the execution time of image analysis applications

Here we propose for the first time a learning-based method, which utilizes prior history of execution times and actual imaging data to learn a model that relates image content and execution time for each individual registration. Although we focus on registration, our formulation is general enough to accommodate several image analysis applications that involve iterative optimization components. It may be utilized wherever estimates of execution time are needed (e.g. in any resource management framework; see Kim et al., 2009; Garg et al., 2010).

To showcase the effect of accurate prediction of execution time from a resource management perspective we consider a proof-of-concept problem: minimizing the total cost of performing a multi-atlas, multi-template segmentation for a set of subjects, templates, and atlases. Under certain assumptions (e.g. ignoring total completion time, or makespan, and not incorporating the unknown to the user wait time for on-demand cores) this optimization reduces to minimizing the individual cost of each registration. The main innovation is the method to accurately predict the execution time of non-linear processes. Once a good estimate of time is available it can be incorporated into any complex optimization and resource-scheduling process.

For convenience of presentation and without loss of generality, we refer to a registration between a subject $S_{n}^{(I)}$ (i.e. the moving volume) and a template $T_{p}^{(I)}$ (i.e. the fixed volume), and assume that all volumes have the same dimensions and data size (as such we ignore time for internal-to-the-cloud data transfers, i.e. $application time \approx execution time$ ). We denote as $f_{$}$ the function that maps the execution time $t_{n, p}^{c}$ of a pair of volumes $S_{n}^{(I)}$ , $T_{p}^{(I)}$ to monetary cost for a given core type $c$ (c1, c2, f2, or m1). Under our assumptions, the minimum total cost is achieved by minimizing the cost of each individual registration and as such we are interested in which core type $c$ minimizes

c_{n, p}^{*} = \underset{c}{argmin} f_{$} (t_{n, p}^{c})

(2)

where $c_{n, p}^{*}$ identifies the minimum cost allocation for a particular volume pair.

Unfortunately, the execution time of an application is not known a priori and a model to estimate it is necessary. We model the execution time $t_{n, p}^{c}$ of a pair of volumes $S_{n}^{(I)}$ , $T_{p}^{(I)}$ on core type $c$ as

t_{n, p}^{c} = \frac{f_{app}}{f_{CPU}^{c}}

(3)

where $f_{CPU}^{c}$ reflects the underlying CPU capabilities of core type $c$ , while $f_{app}$ is a function that reflects how many processing cycles the application requires. Although in equation (3) a simple linear relation is adopted, more complex models that relate application cycles to processing time exist; however, the following analysis still holds and would only require that another equation (3) be adopted.

Most works in the literature related to scheduling or resource-provisioning would assume either the worst-case scenario for the number of cycles needed for an application (by source code profiling; see Li et al., 2009), or they would rely on a historical average estimate (Smith et al., 1998; Kapadia et al., 1999), or learned models of the usage of computational resources for an application (Albers et al., 2009; Matsunaga and Fortes, 2010). For the same application (e.g. a registration) with data size equal for each case, this means that the execution time is estimated either using a constant (a worst-case value is assumed) or as the historical mean (or maximum) of previous executions for the same core type. However, in this context this will not lead to good allocation, because it does not take into account the effect of the actual data content on the execution time of the application (e.g. how many iterations it takes to converge).

We propose that in the context of non-linear image registration, we can actually model the application component more accurately by incorporating the image content within the formulation as

f_{app} = f_{param} (τ) \cdot f_{cont} (S_{n}^{(I)}, T_{p}^{(I)}, τ)

(4)

where $f_{param}$ is a multiplicative component that relates how different registration parameter options (represented as a set τ) of the application affect the required number of cycles per iteration, while $f_{cont}$ represents how many iterations are necessary to achieve the minimum in the non-linear optimization function of the transformation model given parameters τ. Substituting equation (4) into equation (3) we obtain

t_{n, p}^{c} = \frac{f_{param} (τ)}{f_{CPU}^{c}} \cdot f_{cont} (S_{n}^{(I)}, T_{p}^{(I)}, τ)

(5)

For simplicity but without loss of generality, we can assume that all jobs are executed with the same parameters:

t_{n, p}^{c} = {const}^{c} \cdot f_{cont} (S_{n}^{(I)}, T_{p}^{(I)}, τ)

(6)

It is known that intensity-based registration algorithms utilize local similarity metrics, such as normalized cross-correlation, mean squared error, and normalized mutual information, to drive the registration process (Avants et al., 2008). Thus, the similarity of two volumes should affect how fast and easy it will be to register them. Figure 4(top) provides a motivating example, showing color-coded local correlation values (using a $w \times w \times w$ sliding window, with $w = 5$ ) for two different volume pairs at full resolution: one with a slow execution time, and one faster. We see that the slower pair has more areas with low similarity (colored in red hues) when compared to the faster one. As commonly done in neuroimaging, here we assume that the two volumes are first registered linearly to the same common space and as such, gross volumetric differences are non-existent, and only local variation must be accounted for.

Based on this motivation, we proceed by modeling $f_{cont} (S_{n}^{(I)}, T_{p}^{(I)}, τ)$ to reflect the similarity between $S_{n}^{(I)}$ and $T_{p}^{(I)}$ . We must pick a function that is computationally easy to estimate (in other words, it does not take long to compute), but should reflect the diversity of volume pairs a registration algorithm encounters. It should also relate to the similarity criteria the registration algorithm used (hence our inclusion of the parameter set τ to reflect the type of similarity metric and the neighborhood size utilized in the registration). Using, for example, the average across all slices of local similarities should be a good indicator of similarity, but it is not the only one. In fact, the statistical distribution and the topology of local similarities convey important information, as Figure 4(middle) illustrates. Therefore, we propose to use statistical features that characterize the distribution of similarity values and shape descriptors to approximate the topological distribution of high-similarity values between two volumes.

In order to reflect the contribution to time of processing the original volumes hierarchically, using $LV$ levels of the Gaussian pyramid obtained by subsampling the original resolution, we represent as $V_{n, p}^{g}$ the 3D matrix of the collected values of similarity between two volumes $S_{n}^{(I)}$ , $T_{p}^{(I)}$ at the $g$ ^th level of the Gaussian pyramid obtained using a local similarity criterion in a $w \times w \times w$ neighborhood as defined in the set τ.

On the basis of this volume $V_{n, p}^{g}$ we extract the mean, standard deviation, kurtosis, skewness, and median of the vector $ρ_{n, p}^{g}$ , which is the collection of local correlation values of the volume $V_{n, p}^{g}$ .

The above statistical features do not take into account the topological distribution of correlated areas. Thus, we also adopt several easy-to-compute topological 3D shape descriptors introduced in the discrete space (a voxel representation) of binary volumes. We obtain the binary volume $| V_{n, p}^{g} |_{h}$ after thresholding with value $h$ the absolute value volume $| V_{n, p}^{g} |$ , such that an empty voxel (0) means low similarity and a filled voxel (1) means high similarity. We adopt several descriptors, such as sphericity ( $Sp$ ) (Shen and Makedon, 2006), discrete compactness (Cd) (Bribiesca, 2008), and fill ( $Fl$ ) that are computationally efficient. Using these topological features we aim to capture the compactness and density of the similarity volume, since it is known that recovering several small deformations that are spaced far apart requires more effort from a registration standpoint than a larger but co-located and compact deformation.

Sphericity $Sp$ Sphericity is defined as

Sp = \frac{A^{3}}{36 π V^{2}}

(7)

where $A$ is the area of the enclosing surface and $V$ is the volume of the space. Considering every measurement to be unit length, computing $V$ reduces to counting the filled voxels in the space while $A$ corresponds to the sum of the areas of the external polygons (1 in each case) of the voxels which are on the visible faces of $V_{n, p}^{g}$ . Sphericity achieves its minimal value of $1$ for a perfect sphere shape. When comparing the sphericity of multiple objects, the values are scaled to produce a meaningful, comparable result in the range $[0, 1]$ .

Discrete compactness $Cd$ (Bribiesca, 2008) is computed as

Cd = \frac{A_{c}}{A_{c_{\max}}}

(8)

where $A_{c}$ is the contact surface area which corresponds to the sum of the areas of the contact surfaces (1 in each case) common to each two voxels, while $A_{c_{\max}} = 3 (n - n^{2 / 3})$ for an $n$ -voxel representation. $Cd$ is in the range $[0, 1]$ , with $0$ achieved when all voxels are disconnected and $1$ for a perfect cube shape.

Fill $Fl$ is a basic measure defined as

Fl = \frac{V}{V_{total}}

(9)

where $V_{total}$ counts the total voxels in the object representation. $Fl$ is again in the range $[0, 1]$ .

For each new registration we compute the local correlation volume $V_{n, p}^{g}$ (and $| V_{n, p}^{g} |_{h})$ and extract the previously defined features (a process which is extremely fast compared to registration). Based on such features and a model of execution time, which we learn as we will describe subsequently, we estimate the time to perform the registration, in order to identify the core allocation that minimizes monetary cost in our PiCloud implementation. For a collection of new registrations we split the jobs accordingly, by passing the respective volume pairs and core types as cloud.map arguments.

We use supervised learning approaches to obtain models of the execution time $t_{n, p}^{c}$ of intensity-based volume registration from historical data. Given a collection of $ℓ$ previously processed data samples (input volume pairs, and time) on a particular computational set-up $c$ , we represent the local similarity features (extracted at each of the $LV$ levels of the Gaussian pyramid) as a matrix $X \in R^{ℓ \times d}$ , where $d = 8 \cdot LV$ due to the five statistical and the three topological features at each of the $LV$ levels.

We linearize equation (6) by taking the logarithm on each side; therefore we proceed by modeling $\log (t_{n, p}^{c})$ as the dependent variable in all of the regression formulations that follow. Accordingly, a predicted execution time $\log ({\hat{t}}_{n, p}^{c})$ is transformed back to the original scale (seconds), prior to using its value to calculate the monetary cost of executing the corresponding registration on core type $c$ . To predict execution time of non-linear registration, we adopt a variety of linear and non-linear regression models from the machine-learning literature (i.e. linear regression, random forest (Breiman, 2001), and support vector regression (SVR; see Drucker et al., 1996; Smola and Schölkopf, 2004)), which we briefly review in the following paragraphs.

The first type of regression analysis we consider is a simple linear model, which expresses execution time (actually, the logarithm of time) as a linear combination of the input features:

y_{i} = β_{0} + β_{1} x_{i, 1} + \dots + β_{d} x_{i, d}

(10)

where $y_{i} \in R$ is the response (i.e. the collection of logarithmic scaled execution times), $β_{0} \in R$ is the intercept, $β_{1}, \dots, β_{d} \in R$ are the regression coefficients, and $x_{i, j} \in R$ are the predictors, with $i = 1, \dots, ℓ$ and $j = 1, \dots, d$ . We estimate the model parameters that identify the regression hyperplane using a least-squares formulation. With the coefficients ${\hat{β}}_{0}, \dots, {\hat{β}}_{d}$ now known, we predict the computational time for an unseen volume pair using equation (10). We should note that in equation (10), if we ignore the contribution of the content (in other words assume $β_{1} = \dots = β_{d} = 0$ ), solving the least-squares problem to find the optimal constant is equivalent to using the historical mean of execution time at each core type as a predictor of future execution time.

To see if some features (in their extracted level) are more important than others, we optimize the number of features employed by the linear model using a feature selection strategy widely adopted for sparse approximations, in order to discard unnecessary predictors and obtain a possibly more interpretable model (Andrle et al., 2004). We explicitly impose the maximum allowed representation error $ε \in R$ , with respect to the model obtained using the full support of predictors (i.e. all features are used). We adopt a greedy method that gradually reduces the support of the solution by iteratively removing the predictors that offer the smallest model error increase, while remaining below the target imposed representation error $ε$ (Minervini et al., 2013).

In addition to the linear model, we also investigate possible non-linear relations between execution time and features, using random forest regression (Breiman, 2001) and SVR (Drucker et al., 1996; Smola and Schölkopf, 2004). Random forest regression (Breiman, 2001) uses an ensemble of $n_{tree}$ unpruned regression trees ${h (x, Θ_{k}), k = 1, \dots, n_{tree}}$ , where $x$ is a $d$ -dimensional feature vector, and $Θ_{k}$ are random vectors determining bootstrap samples of the original dataset. When growing a tree, the best split at each node is selected among a random subset of $m_{try} << d$ features. A new data sample $x'$ is evaluated at each tree, whose output is the numerical value corresponding to the terminal node reached by applying the set of rules for that tree. The prediction of the random forest for $x'$ is the average across all individual tree predictions. Since features are selected randomly when building the trees, random forest regression performs feature selection implicitly.

On the other hand, from the area of maximum margin classifiers, SVR (Drucker et al., 1996) is a statistical learning method for the computation of a regression function from labeled training data. We first use the basic linear formulation of SVR, that operates in the original $d$ -dimensional input space, thus restricting the regression function to be a hyperplane. We also investigate the non-linear fit to the training data, using a function $φ : R^{d} \to H$ to map the input patterns $x_{i} \in R^{d}$ into a Hilbert space $H$ , possibly higher-dimensional, where the best approximating function could be linear. SVR is formulated as a convex optimization problem (Drucker et al., 1996):

\begin{matrix} \underset{w, ξ, ξ^{*}, b}{minimize} \frac{1}{2} ∥ w ∥^{2} + C \sum_{i = 1}^{ℓ} (ξ_{i} + ξ_{i}^{*}) \\ subject to w^{⊤} φ (x_{i}) + b - y_{i} \leq ε + ξ_{i}, \\ y_{i} - w^{⊤} φ (x_{i}) - b \leq ε + ξ_{i}^{*}, \\ ξ_{i}, ξ_{i}^{*} \geq 0, \end{matrix}

(11)

where $ℓ$ is the number of training samples, $w$ and $b$ are the coefficients of the hyperplane, $C > 0$ is a regularization parameter controlling the flatness of the regression function, $ε$ determines the width of the insensitivity zone (i.e. the margin), $ξ_{i}, ξ_{i}^{*}$ are the slack variables, $y_{i}$ is the target value, and $i = 1, \dots, ℓ$ . The optimization problem of equation (11) is solved in its dual form using the Lagrange multipliers method, to identify the support vectors that define the regression function in the input space. In general, evaluating $φ$ explicitly could be costly (or even unfeasible), hence the mapping is computed efficiently using a kernel function $k : R^{d} \times R^{d} \to R$ , such that $k (x_{i}, x_{j}) = 〈 φ (x_{i}), φ (x_{j}) 〉_{H}$ , where $〈 \cdot, \cdot 〉_{H}$ denotes the inner product in the high-dimensional space. We include in our study the following kernels widely adopted in statistical learning theory:

Linear kernel: $k (x_{i}, x_{j}) = x_{i}^{⊤} x_{j}$ ,

Gaussian kernel: $k (x_{i}, x_{j}) = \exp (- γ ∥ x_{i} - x_{j} ∥^{2})$ ,

Sigmoid kernel: $k (x_{i}, x_{j}) = \tanh (γ x_{i}^{⊤} x_{j} + r)$ , and

Polynomial kernel: $k (x_{i}, x_{j}) = (γ x_{i}^{⊤} x_{j} + r)^{a}$ .

The parameters of the SVR learning algorithm ( $C, \in$ ) and kernel functions ( $γ$ , $r$ , $a$ ) are selected using a grid search strategy.

We use different regression models in order to maximally cover the advantages offered by each one in recovering/interpreting the underlying data structure. For example, the linear model albeit simple is generally able to produce good, interpretable results (especially with the addition of a feature selection step) and can be easily extended to process data in an online fashion, thus updating the model on the fly with every new addition of data. On the other hand, SVR methods are able to accommodate for non-linear dependencies, even though the generated models rarely have a real-world interpretation and a considerable amount of time has to be allocated to the delicate task of parameter tuning. Finally, random forest also allows for non-linearity and additionally can accommodate categorical variables in the model (e.g. if we were to use one model to describe all possible parameter and environment setups) while also providing a natural feature selection mechanism (variable importance estimation). Again, the selection of the model may be time-consuming since many realizations of the random forest are computed for the same parameters and the best model is kept.

In the following section we present the performance of the proposed prediction approaches compared to the classical one of using the historical mean of execution time.

3. Results and discussion

3.1. Experimental setup

To demonstrate the potential of our approach, in this paper we use imaging data from small animal (mouse) phenotyping experiments. We use a total of $N = 20$ subject volumes, $P = 6$ template volumes (from other mouse experiments), and M = 1 atlas (from Dorr et al., 2008). While all the mice are of the same age, 10 mice are common wild types and 10 are transgenic individuals carrying a mutation thought to mimic genetic variation related to a rare human disease. The subjects are scanned ex vivo with a 7T Bruker Scanner, using a T2-weighted RARE sequence at a 0.09 mm isotropic resolution, resulting in volumes of $177 \times 200 \times 177$ voxels. All volumes are first processed to isolate the brain (skull stripping) and to correct for intensity inhomogeneity using standard approaches. Our aim is to identify and quantitate the presence of brain morphological abnormalities in transgenic individuals.

Parameters for all ANTs-based registrations are normalized cross-correlation as similarity criterion in a window of size $w = 5$ voxels (i.e. $5 \times 5 \times 5$ neighborhood), $LV = 4$ , a maximum of 100 non-linear iterations at each level, and Gaussian regularization and symmetric normalization parameters as recommended in Avants et al. (2008).

We execute the image analysis pipeline using the proposed PiCloud implementation only on ‘real-time’ cores. Thus, six jobs (registration between atlas and templates) are required for the first part of the analysis, and 120 jobs (registration between templates and subjects), for the second part. Also, 20 fusion jobs are necessary. Each job is assigned to a core instance on the PiCloud, with locks to complete the registration processes (126 jobs) before launching the fusion jobs. The benefit of using ‘real-time’ cores is that each instance (and job) is launched simultaneously after a fixed provisioning time.

To compare the gain in efficiency by parallelizing this process using the PiCloud platform, we also execute the same process on a local workstation (eight core 3.4 GHz Intel Core i7 and 16 GB RAM) using one (serially) or eight cores (parallel) without any optimization.

3.2. Phenotyping findings

Figure 5 shows a 3D rendering of four structures of interest for one of our subjects. Collectively Table 1 shows average (and standard deviation) volume sizes for the four areas of interest, which are obtained by measuring the number of voxels that have a given label number in the labeled subject volume ${\tilde{S}}_{n}^{(L)}$ . Observe that the transgenic mice exhibit larger cerebellar and thalamic volumes, with respect to the wild type group. However, total intracranial volume (TIV) is shown to be statistically different among the two groups, where the transgenic mice exhibit larger brain sizes. To identify differences in relative volume, we obtain normalized structure volumes, where the size of each structure is normalized by the TIV of that subject. It is evident now that overall the size for all of these four structures is altered in the transgenic group, indicating gross morphometric differences in the brain structure of this transgenic mouse line. (All statistical findings are obtained using Mann–Whitney–Wilcoxon non-parametric tests.)

Table 1.

Morphological phenotyping results for the two mouse lines used where structure volumes are shown as mean (standard deviation), and *, *** indicate significant difference at the 0.05, 0.005 levels, respectively.

	Volume (mm³)		Normalized volume (%)
Structure	Wild type	Transgenic	Wild type	Transgenic
Cortex	67.2 (2.85)	67.44 (5.31)	13.73 (0.18)	12.85 (0.5)***
Cerebellum	57.7 (4.22)	44.83 (5.08)***	11.78 (0.54)	8.53 (0.58)***
Thalamus	14.99 (0.54)	17.09 (1.27)***	3.06 (0.09)	3.26 (0.19)*
Hypothalamus	8.27 (0.3)	8.43 (0.47)	1.7 (0.06)	1.61 (0.05)***
TIV	489.41 (17.6)	524.53 (31.06)*	–	–

3.3. Computational results

In this section we present our findings related to the computational aspects of the process as well as the results of the learning-based modeling of the execution time. All models use the previously defined probability densities and the topological features of the local correlation volumes at all $LV = 4$ levels.

As depicted in Figure 3, the workflow assumes that the user uploads the input data to the cloud and, upon completion, downloads the output data. We measure the upload and download rate from our site in Italy to PiCloud (hosted in East Coast AWS) at 3.7 and 7 Mbit/s respectively. Assuming on average 3 MB for each volume, the total times spent uploading and downloading data are 2.9 min and 69 s respectively, since we upload 27 volumes and download only 20 volumes (the final labeled subjects). Average times (for intra-cloud data transfer and processing time) for each non-linear registration as well as monetary cost are shown in Table 2. Observe that with respect to processing time, data transfer time is negligible, demonstrating the compute-intensive characteristics of this analysis, making it ideal for data parallelism, which is known to scale efficiently.

Table 2.

Average times and cost per non-linear registration on cloud and local workstation shown as mean (standard deviation).

	Average time (min)
Resource type	Data transfer^a	Processing	Average cost^b (US$)
PiCloud c2	0.27 (0.17)	58.28 (22.83)	0.13 (0.05)
PiCloud m1	0.19 (0.13)	43.83 (15.16)	0.22 (0.08)
PiCloud f2	0.19 (0.02)	37.23 (10.44)	0.14 (0.04)
Workstation	–	21.47 (7.07)	–

This accounts for time spent moving around data within the PiCloud infrastructure, reflecting both network-level I/O latency and also latency to the local (node) storage. For reference, input data size = $5040 \pm 136$ kB, output data size = $44, 605 \pm 362$ kB.

Costs according to PiCloud pricing of July 2013.

The makespan of the overall process (in our context the execution time of the longest registration job for the first part of the analysis, followed by the longest registration for the second part, and finally the longest fusion job) for our study on PiCloud requires 8, 5.1, and 4.2 hours for c2, m1, and f2 core types, respectively. On the other hand, using all eight cores on the local workstation requires a total of 18 hours. The total time at the workstation without any parallelization (one core) requires over 140 hours (approximately six days). As the number of registrations or the time to perform a registration increase (higher resolution, different parameters) the local workstation is not able to handle the load.

Considering a hypothetical experiment (a similar one is described in Jia et al., 2012) with $M = 5$ atlases, $P = 200$ synthetic templates, and $N = 45$ subjects will require $P \times (M + N) = 10, 000$ registrations. A workstation with eight cores would require approximately 57 days. Even assuming the existence of a cluster/grid service that would guarantee 32 simultaneous jobs, this hypothetical experiment would need 14 days. On the other hand, PiCloud offers high scalability by allowing the user to provision the required resources for real-time service; as such from a computation view point the makespan of this experiment will be similar to doing a small experiment. Naturally the only limiting factors may be the time spent uploading and downloading data to and from the cloud (insignificant in our case), the guarantee of service up-time, and the accumulated total cost. It is therefore important to identify strategies that may help mitigate the cost involved.

To showcase the potential of the proposed learning framework for predicting execution time of non-linear registration (for minimizing cost), we use real execution times and the 120 volume pairs of the second stage of the registration process (the right-hand part of the graph in Figure 2). These registrations are executed in all possible core services (c2, f2, m1) and thus we collect 360 execution times in total. For reference, the total cost of executing all 120 registrations on c2, f2, and m1 is US$17.2, US$18.1, and US$30.1, respectively. In the following we refer only to c2 and f2, because m1 is more expensive than f2 and slower in all cases. Since for each registration we know the actual execution time at c2 or f2 core service, using equation (2) we identify the optimal scenario of assigning a registration to the core type that minimizes cost. We label this assignment as $c_{n, p}^{*}$ and the actual time as $t_{n, p}^{c}$ . This mix of c2 and f2 core services is referred to as the best-case scenario in the following.

Subsequently we use the volume pairs and execution times for each core service to train all the models described in Section 2.3. We use the same similarity criterion (correlation) and the same neighborhood size (cube of five voxels) as those used for ANTs registrations to obtain $V_{n, p}^{g}$ and extract the statistical and topological features from this correlation map. We use a Matlab implementation of random forest (http://code.google.com/p/randomforest-matlab/) and built-in functions for linear regression. For SVR we use LIBSVM (version 3.17; see Chang and Lin, 2011).

To test the accuracy of predicting the execution time of an unseen volume pair $S_{\bar{n}}^{(I)}$ , $T_{\bar{p}}^{(I)}$ we adopt an unbiased cross-validation strategy, resembling the leave-one-out validation approach. We take the original dataset consisting of 120 registrations (volume pairs) and their execution times at each core type, and for each volume pair we perform a round of cross-validation, using that pair for testing and a subset of the remaining for training. To ensure that the models are trained without encountering any of the testing volumes, and hence are blinded to the testing data, we identify which volumes are involved (for example, $S_{1}^{(I)}$ and $T_{2}^{(I)}$ ). To ensure that no bias exists, we eliminate from the dataset all pairs that involve either of the two volumes in a registration (for our example all registrations that involve either $S_{1}^{(I)}$ or $T_{2}^{(I)}$ are excluded). We use the remaining volume and time pairs (called here the unbiased training set) to train the models described in Section 2.3, using parameters reported in Table 3. Subsequently, we apply each model to the unseen volume pair $S_{\bar{n}}^{(I)}$ , $T_{\bar{p}}^{(I)}$ , obtaining the corresponding predicted execution time ${\hat{t}}_{\bar{n}, \bar{p}}^{c}$ .

Table 3.

Regression model parameters used in the experiments.

Model	Parameters
Linear model	–
Linear model + feature selection	$ε = 16 %$
Random forest	$n_{tree} = 90$ , $m_{try} = 6$
SVR + linear kernel	$C = 6$ , $ε = 0.1$
SVR + Gaussian kernel	$C = 2$ , $ε = 0.1$ , $γ = 0.1$

We repeat this process first on execution times of the c2 core and then for the f2 core, and using equation (2) we estimate the core assignment $c_{\bar{n}, \bar{p}}^{*}$ that minimizes cost. For each volume pair, we estimate the time difference between the actual and predicted times $t_{\bar{n}, \bar{p}}^{c} - {\hat{t}}_{\bar{n}, \bar{p}}^{c}$ , and also record the agreement of assignment, in other words, if $c_{\bar{n}, \bar{p}}^{*} = {\hat{c}}_{\bar{n}, \bar{p}}^{*}$ . We repeat this process for all 120 volume pairs and calculate the mean absolute percentage error (MAPE) of predicting the execution time for each core service, and the percentage of the 120 pairs for which the assignment to a core service is the same as for the best case.

Using such evaluation criteria allows us to compare the performance of the different models when predicting the execution time of a new registration. We use MAPE since this indicator is widely adopted to determine the quality of estimation, especially in the case of time series forecasts, and allows us to evaluate directly any differences in time prediction between core types. (If we were to use actual times this comparison would not be possible since one core type is faster than the other.)

Given the estimated core assignment ${\hat{c}}_{\bar{n}, \bar{p}}^{*}$ , we calculate the average execution time (across all jobs) and the total cost for that allocation. For comparison we use the historical mean of execution times of the training set at each core service as an alternative to the proposed approaches.

Estimating the execution times of non-linear registration processes as accurately as possible is an important aspect of resource provisioning when considering allocation for monetary cost minimization and also for job scheduling. Table 4 outlines the main results obtained regarding the execution time estimation and provides a comparison with the historical mean. (SVR with sigmoid and polynomial kernel report low performance in all experiments, and are excluded from the comparison.) Overall, the MAPE in predicting the execution time of our learning-based approach using a linear model is 17.6% and 15.5% for the c2 and f2 core services respectively. On the other hand, using simply the mean of historical data, the error is 30.1% and 28.7% respectively. Thus, our approach provides almost a 50% improvement over current practice. It appears that a linear model describes the relationship between the features used and execution time well, and offers comparable performance in both training and testing sets, indicating good generalization. The reduced linear model with feature selection (using on average six features) does not improve performance, demonstrating the importance of all features (statistical and topological extracted at all resolution levels). The other non-linear models do not seem to provide any benefit in terms of testing performance. While random forest does offer the lowest training error ( $\approx 7 %$ ) it does not generalize well, indicating probable over-fitting.

Table 4.

Cross-validated results of regression analysis for c2 and f2 core type execution times, reported as MAPE.

	Training		Testing
Model	c2 (%)	f2 (%)	c2 (%)	f2 (%)
Proposed
Linear model	13.0	11.2	17.6	15.5
Linear model + feature selection	14.3	12.3	24.9	17.9
Random forest	7.3	6.2	17.7	16.0
SVR + linear kernel	13.0	11.7	20.0	16.9
SVR + Gaussian kernel	13.5	12.4	18.2	16.2
Baseline
Historical mean	27.3	26.2	30.1	28.7

The first direct application of estimating the job execution times is job allocation. When using the PiCloud platform, the central question is to which type of core to allocate each individual job. We are interested in the overall allocation performance (percentage-wise, how many jobs are allocated to the same type of core as in the best case) and the specific allocation performance (percentage-wise, for each core type, the allocation accuracy as compared to the best case). The complete results are depicted in Table 5. The linear model agrees with the best case in 77.5% of the cases, whereas using the historical mean to drive decisions results in an overall agreement of only 50%. Remarkably, the historical mean agrees less than 3% in f2 allocations with respect to the best case. Such biased behavior is not observed with the linear model. Overall, all of the proposed models offer superior performance compared to the historical mean, and the random forest offers the highest accuracy among all (88.3%). It appears that there is slightly higher accuracy when allocating to c2 rather than f2 core type. However, this is not due to less accuracy in the model (no such difference in execution time prediction was observed before) but because there are fewer allocations to f2 in the best case (32% jobs are allocated to f2, while 68% are allocated to c2). Finally, it appears here that all models perform similarly and sometimes better than the linear model. This is due to the fact that the allocation accuracy is a byproduct of a cost decision. That is, after execution time prediction, the allocation is decided according to cost. Most of the difference in the execution time (estimation error) is masked by the non-linear pricing strategy of real-time cores of PiCloud. On average, as shown in Table 2, a registration takes 37 min on a f2 core. This is remarkably close to the 60% of the hour (36 min) decision threshold for real-time core charge. Thus, what matters for allocation purposes, at least for f2 core type, is whether both real and predicted execution times are below or above this threshold, and not exactly how accurate the time prediction is. However, if the pricing strategy changes one would expect these accuracy figures to change as well, and this is the compelling reason for reporting prediction time as well, since execution time is the input to any decision problem.

Table 5.

Cross-validated results for resource allocation accuracy.

	Training accuracy	Testing accuracy
Model	Overall^a (%)	Overall (%)	c2 Allocation^b (%)	f2 Allocation^b (%)
Proposed
Linear model	83.1	77.5	81.5	69.2
Linear model + feature selection	81.8	73.3	70.4	79.5
Random forest	92.1	88.3	88.9	87.2
SVR + linear kernel	83.8	80.0	82.7	74.4
SVR + Gaussian kernel	86.1	83.3	84.0	82.1
Baseline
Historical mean	62.5	49.2	71.6	2.6

Overall accuracy denotes the percentage of correct allocations with respect to the best case.

c2 (resp. f2) allocation accuracy denotes the percentage of correct c2 (resp. f2) allocations with respect to the best case.

Encouraged by the performance indicators obtained when estimating execution times and allocations, we apply our model to the estimation of the monetary cost of using the PiCloud platform and illustrate the potential monetary gains for intelligently allocating resources. We use the proposed approach to drive core allocations which in turn reduce total monetary cost as seen in the bar plot of Figure 6 (for clarity not all variations are shown). The proposed approaches offer a cost close to the best-case scenario. In the same figure the blue line shows the average execution time of each case. Overall, driving decisions using any of the proposed approaches results in an average execution time of 47.1 min on par with the 46.4 min of the best case. Thus, adopting any of the proposed approaches one can achieve both lower cost and better performance compared to using the c2 service (longer time), the f2 service (more costly), or the historical mean (more costly and longer time). Naturally, the prediction accuracy will improve by adding more data and a more advanced similarity model (one that is still computationally efficient to estimate).

Several conclusions can be drawn from these findings. The linear, random forest and SVR + Gaussian kernel models perform very well in terms of time and cost estimation, leading to an approximate 12–15% cost reduction compared with ‘All f2’: the most expensive solution. We emphasize that all proposed methods offer quantitatively better results than the historical mean. This leads to the conclusion that inquiry into the content of the input image data has a great importance for execution time estimation and, ultimately, cost prediction and optimization. This approach seems to provide more insight into the correct allocation patterns than only using data dimension and/or historical information.

In this paper we use a few data points but we can assess how the methods will behave when applied to larger datasets. Assuming the availability of a large number of independent and homogeneous jobs we expect the result of the historical mean to improve significantly. Of course, this setting also implies that the results obtained by the proposed methods will improve as well. Although the difference between the historical mean solution and the proposed solutions may decrease in this case, still we expect the proposed ones to always outperform the naive historical mean approach and provide a more reliable time estimation mechanism, for example for scheduling applications. Conversely, if we consider a more complicated setup (such as more core types, more complex imposed cost function), a more heterogeneous set of imaging volumes, or an application where accurate time prediction is critical (such as the scheduling of large time-dependent workflows of imaging applications), the results of the proposed solutions will perform well while we expect a degradation in the results of the historical mean.

4. Conclusions

In this paper we demonstrate that offloading intensive analyses to the cloud is possible and beneficial, since they are ideally suited to data parallelism (i.e. they are compute-intensive and not limited by inter-process or intra-cloud communication). While the PiCloud platform is convenient and facilitates the use of cloud infrastructures, it is possible to implement similar pipelines using other libraries if the user so chooses. Nevertheless, no matter what environment is used (cloud, cluster, or multi-core workstations) optimizing resource allocation is necessary and the proposed framework for incorporating data dependency does improve the efficiency of resource allocation methods. While here we demonstrate this potential using simple similarity-driven features (statistics of their distribution and topological shape descriptors) for estimating the computational burden of non-linear registrations, more sophisticated methods can be developed. It is clear that the more volumes are registered the more accurate our estimation will be. We are currently developing methods to integrate this prediction process seamlessly, and to extend it to predicting execution time of other core neuroimaging analysis steps (e.g. segmentation). We will release in the public domain all source code material to facilitate adoption by users that are interested in performing such an analysis on the PiCloud platform. Although here we rely on mouse data, the same code and approach can be used with any brain MRI data (human or other), and can be easily extended to accommodate different pipelines as well. With this paper our goal is to provide a compelling example of how the cloud paradigm is used to create powerful and comprehensive image analysis pipelines that are simple to use and available to a large community of scientists, thus increasing the democratization of science. The possibility to perform hundreds of parallel executions at a fraction of the cost of owning a workstation opens the possibility to many laboratories around the world to develop state-of-the-art image analysis pipelines, without having to rely only or at all on (shared) distributed computing facilities.

Footnotes

Acknowledgements

The authors would like to thank Alberto Galbusera for assistance with data collection, and PiCloud Inc., particularly Ken Elkabany, for providing technical assistance.

Funding

This work was partially supported by a Marie Curie Action: ‘Reintegration Grant’ (grant number 256534) of the EU’s Seventh Framework Programme (FP7).

Author biographies

Massimo Minervini received BSc and MSc degrees in Computer Science from the University of Bari, Italy, in 2008 and 2010, respectively. He is currently pursuing a PhD in Computer Science and Engineering at the IMT Institute for Advanced Studies, Lucca, Italy, investigating application-aware image compression for plant phenotyping. His research interests include data compression, computer vision, machine learning, and applications to life science.

Cristian Rusu received a PhD in Signal Processing from the University Politehnica of Bucharest, Romania, in 2012. He is currently a post-doctoral researcher with the Pattern Recognition and Image Analysis group at the IMT Institute for Advanced Studies, Lucca, Italy. His current research interests include signal processing, numerical linear algebra, sparse representations, dictionary learning, and applications to image processing.

Mario Damiano obtained a BSc in Biomedical Engineering from the University of Pisa, Italy, in 2008 with a thesis on algebraic reconstruction techniques in CT images. In 2011 he earned an MSc in Clinical Engineering from the University of Trieste and became a research assistant at the Italian Institute of Technology, developing pipelines for the morphological analysis of CNS rodent models in MRI. He is now pursuing a specialist Master in Management of Clinical Engineering at the University of Trieste, focusing on the maintenance of medical devices and health technology assessment.

Valter Tucci graduated in Psychology in 2000 at the University of Padua, Italy, studying the cardiovascular changes associated with NREM and REM sleep states in humans. He received a PhD in Sleep Medicine from the University of Genoa, Italy, in 2004, investigating the physiological and cognitive traits of narcoleptic patients. Between 2003 and 2008 he was a post-doctoral researcher and then investigator scientist at the Medical Research Council Harwell, Oxford, UK, where he focused on behavioral and cognitive phenotyping of mouse models. In 2008 he was awarded the medal ‘Mente & Cervello’ and the medal by the President of the Italian Republic as ‘Young Investigator’. He is currently team leader of the Neurobehavioural Genetics Group at the Italian Institute of Technology in Genoa, Italy. His research is focused on the analysis of the effects that genetic and epigenetic mechanisms exert on sleep and cognition.

Angelo Bifone is the Director of the Center for Neuroscience and Cognitive Systems of the Italian Institute of Technology in Rovereto, Italy. He received a PhD in Physics from the Scuola Normale Superiore in Pisa, Italy, and a Master in Business Administration from the University of Bologna. He trained and worked at the University of Berkeley, CA, at the University of Leiden, The Netherlands, and at the Institute of Cancer Research, London, UK. Before joining the Italian Institute of Technology in 2010, he headed the Neuroimaging Department of the GlaxoSmithKline Medicines Research Center in Verona, Italy, for eight years. His main research focus is MRI, and its application to the study of brain structure and function.

Alessandro Gozzi received an MSc degree in Biotechnology and a PhD in Biomedical Imaging from the University of Verona, Italy. He later joined GlaxoSmithKline, a research-based pharmaceutical company, where he developed and implemented functional MRI methods to probe pharmacological mechanisms in preclinical models of CNS disorders. During his research, he has identified and described the circuital substrates modulated by centrally active substances belonging to different pharmacological classes, and within multiple areas of experimental neuroscience (schizophrenia, depression, drug addiction, etc.). He is currently group leader at the Istituto Italiano di Tecnologia in Rovereto, Italy, where he leads an MRI laboratory equipped with a preclinical 7 Tesla scanner for in vivo mouse neuroimaging. His research interests focus on the application of advanced functional and structural MRI methods to describe healthy and pathological neurofunctional states, and the investigation of circuital basis of behavior.

Sotirios A Tsaftaris is with the IMT Institute for Advanced Studies, Lucca, Italy, where he is the Director of the Pattern Recognition and Image Analysis Unit. Previously, he held a joint research assistant professor appointment with the Departments of Electrical Engineering and Computer Science and Radiology at Northwestern University, IL, USA. Currently, he maintains an adjunct appointment with Northwestern. He has published extensively, particularly in interdisciplinary fields, with almost 90 journal and conference papers in his active record. His research interests are medical image analysis, computational biology, machine learning, and large-scale analysis of imaging data. Professor Tsaftaris is a Murphy Fellow, a Marie Curie Fellow, and a Fellow of the Alexander S Onassis Public Benefit Foundation.

References

Albers

Suijs

de With

(2009) Triple-C: Resource-usage prediction for semi-automatic parallelization of groups of dynamic image-processing tasks. In: IEEE international symposium on parallel distributed processing (IPDPS), pp. 1–8.

Andrle

Rebollo-Neira

Sagianos

(2004) Backward-optimized orthogonal matching pursuit approach. IEEE Signal Processing Letters 11(9): 705–708.

Antoniu

Bougé

Thirion

. (2010) AzureBrain: Large-scale joint genetic and neuroimaging data analysis on Azure clouds. Available at: http://www.irisa.fr/kerdata/lib/exe/fetch.php?media=pdf:inria-microsoft.pdf . (accessed 14 December 2013)

Avants

Epstein

Grossman

. (2008) Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis 12(1): 26–41.

Breiman

(2001) Random forests. Machine Learning 45(1): 5–32.

Bribiesca

(2008) An easy measure of compactness for 2D and 3D shapes. Pattern Recognition 41(2): 543–554.

Buyya

Venugopal

(2004) The Gridbus toolkit for service oriented grid and utility computing: An overview and status report. In: IEEE international workshop on grid economics and business models (GECON), pp. 19–66.

Buyya

Pandey

Vecchiola

(2009) Cloudbus toolkit for market-oriented cloud computing. In: Martin

Gansen

Chunming

(eds.), Cloud Computing (Lecture Notes in Computer Science, vol. 5931). Berlin/Heidelberg: Springer, pp. 24–44.

Chang

Lin

(2011) LIBSVM: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology 2(3).

10.

Dinov

Lozev

Petrosyan

. (2010) Neuroimaging study designs, computational analyses and data provenance using the LONI pipeline. PLoS ONE 5(9): e13070.

11.

Dorr

Lerch

Spring

. (2008) High resolution three-dimensional brain atlas using an average magnetic resonance image of 40 adult C57Bl/6J mice. NeuroImage 42(1): 60–69.

12.

Drucker

Burges

CJC

Kaufman

. (1996) Support vector regression machines. In: Advances in neural information processing systems 9 (NIPS), pp. 155–161.

13.

Frisoni

Redolfi

Manset

. (2011) Virtual imaging laboratories for marker discovery in neurodegenerative diseases. Nature Reviews Neurology 7(8): 429–438.

14.

Garg

Buyya

Siegel

(2010) Time and cost trade-off management for scheduling parallel applications on utility grids. Future Generation Computer Systems 26(8): 1344–1355.

15.

Heckemann

Hajnal

Aljabar

. (2006) Automatic anatomical brain MRI segmentation combining label propagation and decision fusion. NeuroImage 33(1): 115–126.

16.

Jia

Yap

Shen

(2012) Iterative multi-atlas-based multi-image segmentation with tree-based registration. NeuroImage 59(1): 422–430.

17.

Kapadia

Fortes

JAB

Brodley

(1999) Predictive application-performance modeling in a computational grid environment. In: 8th international symposium on high performance distributed computing, pp. 47–54.

18.

Kim

Parashar

(2011) CometCloud: An autonomic cloud engine. In: Rajkumar

James

Andrzej

(eds.), Cloud Computing: Principles and Paradigms. New York, NY: Wiley, pp. 275–297.

19.

Kim

Parashar

Foran

. (2009) Investigating the use of autonomic cloudbursts for high-throughput medical image registration. In: 10th IEEE/ACM international conference on grid computing, pp. 34–41.

20.

Klein

Andersson

Ardekani

. (2009) Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration. NeuroImage 46(3): 786–802.

21.

Kumar

Feng

Nimmagadda

. (2011) Resource allocation for real-time tasks using cloud computing. In: 20th international conference on computer communications and networks (ICCCN), pp. 1–7.

22.

Laguitton

Rivière

Vincent

. (2011) Soma-workflow: A unified and simple interface to parallel computing resources. In: MICCAI workshop on high performance and distributed computing for medical imaging.

23.

Singh

. (2009) Machine learning based online performance prediction for runtime parallelization and task scheduling. In: International symposium on performance analysis of systems and software (ISPASS), pp. 89–100.

24.

Lin

(2011) Scheduling scientific workflows elastically for cloud computing. In: International conference on cloud computing (CLOUD), pp. 746–747.

25.

Liu

Jin

Chen

. (2010) A compromised-time-cost scheduling algorithm in SwinDeW-C for instance-intensive cost-constrained workflows on a cloud computing platform. International Journal of High Performance Computing Applications 24(4): 445–456.

26.

Mao

Humphrey

(2011) Auto-scaling to minimize cost and meet application deadlines in cloud workflows. In: International conference for high performance computing, networking, storage and analysis, pp. 1–12.

27.

Matsunaga

Fortes

(2010) On the use of machine learning to predict the time and resources consumed by applications. In: 10th IEEE/ACM international conference on cluster, cloud and grid computing (CCGrid), pp. 495–504.

28.

Millman

Brett

(2007) Analysis of functional magnetic resonance imaging in Python. Computing in Science & Engineering 9(3): 52–55.

29.

Minervini

Damiano

Tucci

. (2012) Mouse neuroimaging phenotyping in the cloud. In: 3rd international conference on image processing theory, tools and applications (IPTA), Istanbul, Turkey, pp. 55–60.

30.

Minervini

Rusu

Tsaftaris

(2013) Learning computationally efficient approximations of complex image segmentation metrics. In: 8th international symposium on image and signal processing and analysis (ISPA), Trieste, Italy, pp. 60–65.

31.

Pandey

Karunamoorthy

Gupta

. (2009) Megha workflow management system for application workflows. In: IEEE science & engineering graduate research expo, Melbourne, Australia.

32.

Ramakrishnan

Chase

Gannon

. (2011) Deadline-sensitive workflow orchestration without explicit resource control. Journal of Parallel and Distributed Computing 71(3): 343–353.

33.

Reynolds

Winter

Terstyanszky

. (2011) Scientific workflow makespan reduction through cloud augmented desktop grids. In: 3rd international conference on cloud computing technology and science (CloudCom), pp. 18–23.

34.

Shen

Makedon

(2006) Spherical mapping for processing of 3D closed surfaces. Image and Vision Computing 24(7): 743–761.

35.

Shi

Hong

(2010) Resource allocation with a budget constraint for computing independent tasks in the cloud. In: 2nd international conference on cloud computing technology and science (CloudCom), pp. 327–334.

36.

Smith

Foster

Taylor

(1998) Predicting application run times using historical information. In: Dror

Larry

(eds.), Job Scheduling Strategies for Parallel Processing (Lecture Notes in Computer Science, vol. 1459). Berlin/Heidelberg: Springer, pp. 122–142.

37.

Smola

Schölkopf

(2004) A tutorial on support vector regression. Statistics and Computing 14(3): 199–222.

38.

Tabak

Cambazoglu

Aykanat

(2013) Improving the performance of independent task assignment heuristics MinMin, MaxMin and Sufferage. IEEE Transactions on Parallel and Distributed Systems.

39.

Tordsson

Montero

Moreno-Vozmediano

. (2012) Cloud brokering mechanisms for optimized placement of virtual machines across multiple providers. Future Generation Computer Systems 28(2): 358–367.

40.

Van den Bossche

Vanmechelen

Broeckhove

(2010) Cost-optimal scheduling in hybrid IaaS clouds for deadline constrained workloads. In: 3rd international conference on cloud computing (CLOUD), pp. 228–235.

41.

Vecchiola

Chu

Buyya

(2009) Aneka: A software platform for .NET-based cloud computing. High Speed and Large Scale Scientific Computing 18: 267–295.