Philippe Monmousseau
email: philippe.monmousseau@enac.fr

Aude Marzuoli
email: amarzuoli3@gatech.edu

Christabelle Bosson
email: cbosson@alumni.purdue.edu

Eric Feron
email: eric.feron@aerospace.gatech.edu

Daniel Delahaye
email: delahaye@recherche.enac.fr

Doorway to the United States: An Exploration of Customs and Border Protection Data

Keywords: big data, passenger-centric, wait time prediction

This paper presents a data-driven study of wait time patterns for international arriving passengers across all 61 terminals from the 44 airports of entry of the United States. Each airport is an independent entity which operates with various airlines and handles demand volumes differently. This induces seasonal variation in service quality from one airport to another. Exploring six years worth of data, this paper investigates the current and long-term performance trends -an increasing number of flights versus a decreasing number of customs booths -of all airports of entry from a passenger perspective. A performance analysis is then conducted that compares average wait times of incoming passengers, considering incoming traffic ratios and allocated resources. Leveraging machine learning algorithms, six regression algorithms are trained and tested to accurately predict passenger wait times through customs at selected airports. An analysis of the performance of these models shows that the best approach -using a Gradient Boosting regressor for each terminal of entry -can capture the daily and seasonal variations of traffic patterns and immigration booth availabilities with a mean absolute error of less or equal to 5 minutes for twentyeight terminals of entry and less than 10 minutes for all terminals. Observations show significant disparities across airports that may be explained by the foreign/US passenger ratio and the quality of booth management.

I. INTRODUCTION

The U.S. Department of Transportation's Bureau of Transportation Statistics (BTS) reported that U.S. airlines carried an all-time high number of passengers in 2016 -823.0 million systemwide, 719.0 million domestic and 103.9 million international [1]. Passengers are at the core of the Air Transportation System and, yet, obtaining quantitative information about passenger movements remains difficult, particularly on the ground. NextGen [START_REF]NextGen Integration and Implementation Office[END_REF] in the United States and SESAR [START_REF] Ky | SESAR: towards the new generation of air traffic management systems in europe[END_REF] in Europe have been advocating a shift from flight-centric metrics to passenger-centric metrics to evaluate the performance of the Air Transportation System. Both NextGen and SESAR intend to not only improve the predictability and resilience of the Air Transportation System, but also to reduce door-to-door travel time for passengers. While passengers whereabouts are known once they board an airplane, their behavior inside the airport is less understood, and passenger-related information is siloed between different stakeholders, such as airlines, airport security, customs and border protection. Foreign international air travelers arriving in the U.S. spend billions of dollars while visiting. Torres et al. showed that consumption of goods and services grows with the time spent by passengers in the leisure areas. [START_REF] Torres | Passenger waiting time in an airport and expenditure carried out in the commercial area[END_REF]. Thus less time spent queueing at various control points may result in time and financial benefits for everyone.

Airports are known to be the main bottlenecks of the Air Transportation System, yet passenger satisfaction is largely affected by their experience at the airport [START_REF] Pruyn | Effects of waiting on the satisfaction with the service: Beyond objective time measures[END_REF]. Airport management is distributed across various stakeholders, including airlines, federal agents, local airport authorities and third parties, whose decisions may, unbeknownst to them, impact each other's performance [START_REF] Nikoue | Passenger Flow Predictions at Sydney International Airport: A Data-Driven Queuing Approach[END_REF]. Wait times at security and border patrol play a big role in assessing passenger satisfaction : capacity constraints and inefficiencies at airport entry roads, parking, security, immigration, customs, gates, ramp areas, runways are the primary causes of congestion and of the ensuing delays. Since September 11, 2001, airport screening procedures in the U.S. have been continuously evolving. For example, the passenger screening process is now trying to strike a balance between security and customer service (i.e. minimizing wait times). Using data from 2002 and 2003, Gkritza et al. [START_REF] Gkritza | Airport security screening and changing passenger satisfaction: An exploratory assessment[END_REF] showed that, while wait times at security screening points are significant determinants of passenger satisfaction, many other factors come into play.

Moreover, delays on the ground have a disproportionate impact for passengers who often experience lengthy delays before being re-booked if they miss their connecting flights [START_REF] Wang | Methods for analysis of passenger trip performance in a complex networked transportation system[END_REF]. For international passengers arriving from overseas to the United States, immigration checks are mandatory, whether they are U.S. citizens or foreigners.

This paper focuses on the passenger experience while going through U.S. customs and border protection. Roberts et al. studied the evolution of wait times at airports over a few years [START_REF] Roberts | Analysis of primary inspection wait time at U.S. ports of entry[END_REF]. They showed that average passport inspection wait time at 24 U.S. airports rose by 25% during 2010-2013. They focused on JFK airport. At JFK, nearly 3 million passengers (25 percent of total arrivals) experienced a delay of more than 1 hour, putting them at risk of missing a connecting flight, with the 11 percent who had a total delay of more than 2 hours, possibly missing connections at a higher rate. Extended passport inspection waits were the sole source of missed connection risk for 13 percent of passengers and one of the reasons behind missed connections for many more passengers.

There has been little research in the systematic analysis of passenger wait times at customs across airports. Sankaranarayanan et al. [START_REF] Sankaranarayanan | An exploratory data analysis of airport wait times using big data visualisation techniques[END_REF] performed an exploratory analysis of airport wait times on customs, border protection data taken from top 3 busiest airports (Atlanta, Chicago, and Los Angeles) from the United States, highlighting the effects of seasonality. Johnstone et al. presented a dynamic queue controller to generate realistic queue formation and behaviour within a discrete event environment at airports in Australia [START_REF] Johnstone | A dynamic architecture for increased passenger queue model fidelity[END_REF]. Robertson et al. [START_REF] Robertson | The role of modeling demand in process re-engineering[END_REF] modeled passenger arrivals to estimate how many passengers arrive at the airport per day and time of day. They also estimated passenger arrival patterns at different processing points (check-in, baggage security, security checkpoint, etc). Nikoue et al. [START_REF] Nikoue | Passenger Flow Predictions at Sydney International Airport: A Data-Driven Queuing Approach[END_REF] focused on Sydney airport, where they studied flight information, cell phone geolocations from MAC device time stamps, and immigration records from the Australian Department of Immigration and Indigenous Affairs. They modeled arrivals of passengers at Sydney International airport, simulated walk times to immigration and service at immigration queues, to study the effects of a variation in staffing levels at immigration on congestion. However, to the best of the authors' knowledge, no work is available on comparing performance across airports or predicting performance at any airport.

The present paper leverages publicly available data from the United States Customs and Border Protection (CBP) [START_REF] Customs | Airport wait times[END_REF]. As stated on their website, "CBP closely monitors the flight processing times, commonly referred to as wait times, for arriving flights at the busiest international airports." The data provided in the online reports show the number of passengers processed on flights arriving in each hour based on how long it took for those passengers to clear Passport Control.

This paper tackles the following research questions:

• Which are the best airports to enter the U.S. for U.S. citizens and foreigners, in terms of wait times? • Which airports best manage their customs area?

• Can wait times per hour at airports be reliably predicted from historical data?

The paper is organized as follows. Section II explores the main trends at different scales visible in the CBP data, from average wait times, to passenger volumes and flight volumes. Section III compares the performances of the different airports of arrival. Section IV proposes a machine learning approach to predict passenger wait time per hour at any airport and examines the performance of this approach. Section V details the conclusion of the paper and offers future research perspectives.

II. EXPLORATION

A. Dataset contents

The data from CBP [START_REF] Customs | Airport wait times[END_REF] contains the following fields once a time period of interest is selected:

B. Long term evolution

Looking at the overall evolution from 2013 to 2019, some clear trends appear as illustrated in Figures 1 &2. Figure 1 shows the evolution of the total number of arriving international flights per day versus the total number of open booths per day. While the number of flights is steadily increasing over the years, the number of open booths is slowly decreasing. Figure 4 shows the box-and-whisker plot variations of the average wait time across all terminals per month, i.e. for each month it shows the median average wait time and the first and third quartile along with whiskers for a better visualization of the range of the data. This figure highlights the four-month periodic variation discovered in Figure 3. There are three highs during a year -around the winter holidays from December to January, April and August -interlaced with periods with shorter wait times. This four-month periodic behavior does not seem to be induced by a similar behavior in the number of arriving passengers, as shown in Figure 5.

Figure 6 shows the evolution of average wait time for all passengers per hour of the day from 2013 to 2019. As previously observed, the average time for 2013 is higher than all the others, which are tightly packed. The average wait time for 2019 is lower than the previous years since only data from January are available at the time of this study.

C. Wait time differences

The trends presented so far were for any passenger. Yet, U.S. citizens and foreigners enter the U.S. through separate lines, and statistics on wait times are available for each category specifically. For the same number of passengers, processing foreigners at a booth typically takes longer, since it requires checking more paperwork, such as supporting entrance documents and visas, whereas U.S. citizens only need to present their passports. The average wait time for all passengers over the year 2017 is 16.7 minutes and 16.2 minutes over the year 2018. Figure 8 shows the break down by U.S. passengers and non-U.S. passengers for these two years. Non-U.S. passengers on average spend twice as much time in line at immigration, and experience higher volatility in wait times. Figure 9 shows the different wait times distribution per year for US citizens and non-US citizens from 2013 to 2019. The average wait time for non-US citizens does not vary much over these years, while being twice as more important as for US citizens throughout these years. Though 2019 seems to be better, observations from Section II-B indicates that the month of January is not representative of the yearly distribution.

A. Overall comparison

Figure 10a shows the yearly box-plots per airports of the average wait time for all passengers for the year 2018. Interestingly both Orlando airports in Florida (SFB and MCO) show particularly high wait times with high volatility. Second in line for overall longest wait times is Hawai (SPN), before airports from the San Francisco area (SMF and SFO). Ontario International airport (ONT) is not considered in this top five due to its low volatility compared to SFO.

This high wait time is not necessarily due to lack of means. Figure 10b shows the distribution of the number of open booths per airport for the year 2018. SFO has one of the highest median of open booths, which contrasts with its wait time performance.

Figure 10c shows the distribution of the number of arriving passengers per airport for the year 2018. A first observation is that airports with high passenger volume volatility have in general a high flexibility regarding the number of open booths.

A comparison of the ratio of arriving passengers and open booths will appear in a follow-up publication.

B. Worst and Best case scenarii

From Figure 4 & 5, one can infer that the worst month for entering the United States is August: the highest median average wait time with high volatility combined with one of the largest distribution of arriving passengers. Figure 11a shows the distribution of the average wait times for all passengers per airports in August over the last six years. Airports ONT and CVG were removed since there is only two years of data for these two airports. The same two airports from Orlando Florida appear in the top five worst performing airports. On the West coast, the previously spotted airports are joined by OAK and FAT join Orlando in the top six worst performing airports.

The performance of the first four airports are even worse when considering the number of arriving passengers shown in Figure 11b. They do not have the highest median number of arriving passenger per hour nor the highest volatility.

On the other hand, from Figure 4 & 5, February appears like the best month to enter the US: lowest average wait time and volatility as well as the lowest volume of arriving passengers. When combining these observations with the number of arriving passengers distribution shown in Figure 12b, CLT has more merit seeing how few passengers enter through PBI during that month of the year. 13b for a better comparison between the two categories. A first observation is that airports with short wait times for US citizens do not necessarily shorter wait times for non-US citizens than airports with high wait times for US citizens. For example, SPN has the second worst wait time for non US passengers while having the third shortest wait time for US passengers. Only two airports (STL and PBI) have a lower median average wait time for non-US passengers than for US passengers. Figure 13c, which shows the distribution of the average wait time ratio of non-US citizens over US citizens, indicates that only five airports have a median ratio outside the range [1,2.5]: STL and PBI with a ratio lower than one as noted previously, along with MCO, Guam airport (GUM) and SPN with ratios greater than 2.5.

IV. HOURLY WAIT TIME PREDICTION ACROSS AIRPORTS

Typical modeling of queues at airports relies on queuing theory studied in Operations Research to evaluate queue length and service time [START_REF] Kleinrock | Computer applications[END_REF]. In this paper, we choose to adopt a different approach. Leveraging machine learning techniques, our goal is to predict the average wait time per hour at any

A. Features and regressors

The problem of interest falls in the category of regression techniques [START_REF] Bishop | Pattern recognition and machine learning[END_REF], more specifically under time series forecasting. The data set is partitioned into a train set and a test set. Each row in the data set corresponds to a particular hour, and the corresponding label is the average wait time for all passengers for the next hour. The following base set of features is created: The different regression models are implemented using Python as the programming language and using the Scikitlearn library [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF].

To obtain the best possible performance, we experiment with various algorithms, each having its specific advantages and drawbacks. Below is a brief overview of each algorithm tested:

1) Linear Regression assumes that the relationship between the input variables and the measured variable is linear with some noise, and estimates the parameter vector with ordinary least squares minimization. It is the simplest regression method and is easily interpretable. 2) Ridge Regression, also known as Tikhonov regularization, extends the ordinary linear regression with a penalty term in the objective function proportional to the error norm. This improves the conditioning of the problem and reduces overfitting. 3) Lasso Regression performs a linear regression with regularization as well as a variable selection. 4) Random Forest Regression is an ensemble technique that fits a number of classifying decision trees on various sub-samples of the dataset and uses averaging to improve the predictive accuracy and control over-fitting 5) Gradient Boosted Tree Regression is an ensemble technique relying on decision trees as weak learners and resampling the training samples while assigning weights to the samples. It optimizes a cost function over function space by iteratively choosing a function that points in the negative gradient direction. It is typically more robust than basis learners. 6) AdaBoost is also an ensemble of decision trees relying on boosting. However, AdaBoost and Gradient Boosting optimize different loss functions.

B. Performance measures

In order to measure the quality of the predictions, two different performance measures were computed: the R 2 score and the mean absolute error (MAE).

The R 2 score, also known as the coefficient of determination, is defined as the unity minus the ratio of the residual sum of squares over the total sum of squares: R 2 = 1i (yi-fi) 2 i (yi-ȳ) 2 , where y is the value to be predicted, ȳ its mean and f is the predicted value. It ranges from -∞ to 1, 1 being a perfect prediction and 0 meaning that the prediction does as well as constantly predicting the mean value for each occurence. In the case of a negative R 2 , then the model has a worse prediction than if it were predicting the mean value for each occurence and therefore yields no useful predictions. The closer the R 2 score is to 1, the more accurate the prediction is.

Regarding the mean absolute error, the smaller its value is, the more accurate the prediction is. It is calculated using the following formula: MAE = 1 n i |f i -y i | where n is the number of values being predicted.

C. Simple prediction benchmarks

Four different simple prediction benchmarks were tested in order to have a better understanding of the ease of prediction. For the first three, the prediction consisted in taking the value from the previous year, the previous day or the previous hour. The last one consists in constantly predicting the mean value of the the training set.

As a more elaborate comparison benchmark, Facebook's time-series forecasting tool Prophet [START_REF] Taylor | Forecasting at Scale[END_REF] was trained on the actual training set and its performance was measured on the data from the year 2018. The Prophet tool is based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality [START_REF] Facebook | Prophet -Forecasting at scale[END_REF]. It is described as robust to outliers and missing data with no parameter tuning necessary, therefore the default parameters of the Prophet tool were used for this forecasting benchmark.

D. Performance analysis

1) One-hot encoding analysis: A single model for all terminals was created by adding one-hot encoding features for the different terminals, i.e. 59 binary features were added, each one indicating whether a specific terminal is considered or not. This single model was trained on the data from 2013 to 2017 and tested on the data from 2018. 59 other models, one for each terminal, were trained and tested on the same data fitlered by terminal. Figure 14 plots the MAE per terminal for each considered regressor. It shows that for each regressor the single terminal method outperforms the one-hot encoding method for a majority of terminals. For example, in the case of the Gradient Boosting regressor, having a different regressor per arriving terminal is better than the one-hot encoding method for forty terminals out of fifty-nine. 2) Benchmark comparison: The comparison with the chosen benchmarks was done using the single-terminal models, i.e. for each regressor type, one model was trained per arrival terminal. The R 2 score and mean absolute error of these models are aggregated in boxplots presented in Figure 15 alongside the boxplot performance of the five benchmarks.

Figure 15a shows that he Gradient Boosting regressors have the best R 2 performance, its R 2 scores being greater than 0 for 43 terminals out of 59, which is to be compared to 28 out of 59 for the first benchmark, the Prophet tool. Gradient Boosting has also the smallest R 2 performance deviation of all tested models. Though Ada Boost has the highest R 2 score of all models for one terminal (SFB Terminal A), of the six chosen regressors, it is the only model with a median R 2 score less than 0, along with the five benchmarks. The three linear models (Linear, Lasso and Ridge) have similar performances in this study, implying that the overfitting methods added in Lasso and Ridge are not necessarily relevant here. The maximum R 2 difference between Lasso and Linear or Ridge is about 0.64 while this distance is about 0.005 between Ridge and Linear.

Figure 15b, showing the MAE distribution, has a different ranking with respect to the linear regressors. With this performance measure, Lasso yields better results than Linear and Ridge. Otherwise the conclusions for this performance measure are the same as for the R 2 score. This performance measure has however a more tangible interpretation: Of all the models tested, only the Gradient Boosting regressors predict the average waiting time with an average error of 10 minutes or less for every terminal of entry. The Gradient Boosting regressors have a MAE of 5 minutes or less for 28 terminals.

3) Training size analysis: In order to visualize the impact of the training set's size, the regressors were trained using training sets of different lengths (i.e. different starting dates) and their performance was tested on the data from the year 2018. As shown in Figure 16, decreasing the size of the training set does not have a major impact on the performance of the Random Forest regressor and the Gradient 17 shows a comparison per terminal of entry between four of the best performing models' mean absolute error and the standard deviation of the average wait time over the test year 2018. This plot shows that except for one exception (SNA Terminal C), the Gradient Boosting models mean absolute errors are significantly better than the standard deviation of the values to predict.

V. CONCLUSION

This paper is a first step in a systematic analysis of passenger wait times at customs across all airports of entry in the United States using publicly available data. This analysis makes it easier to uncover some long-term trends, i.e. an increase in the number of arriving flights coupled with a decrease in the number of open custom booths, while also enabling a per airport comparison. This analysis also laid the ground to implementing machine learning regression models in order to predict the average wait time per airport of entry. These models could be used to better anticipate the number of required booths once the number of arriving passengers is known.

APPENDIX

Figure 1 :

 1 Figure 1: Comparison of the total number of open booths (red) vs. the total number of arriving flights (blue) per day from January 2013 to January 2019

Figure 2 :

 2 Figure 2: Average wait time for all passengers per day across all airports from January 2013 to January 2019.

Figure 3 :

 3 Figure 3: Yearly comparison of the average wait time for all passengers per day across all airports

Figure 4 :Figure 5 :

 45 Figure 4: Boxplots per month of the hourly average wait time for all passengers across all terminals from 2013 to 2019

Figure 6 :

 6 Figure 6: Yearly comparison of the average wait time for all passengers per hour across all airports

Figure 7 :

 7 Figure 7: Average wait time distribution per day of the week from 2013 to 2019

  (a) On year 2017 (b) On year 2018

Figure 8 :

 8 Figure 8: Average wait time distribution for passengers from January to December for the years 2017 and 2018.

  (a) US passengers (b) Non US passengers

Figure 9 :

 9 Figure 9: Average wait time and standard deviation for passengers from 2013 to 2019 across all airports

Figure 10 :

 10 Figure 10: Airport comparison using boxplots over the year 2018

Figure 11 :

 11 Figure 11: Airport comparison using boxplots for the month of August from 2013 to 2019

Figure 12 :

 12 Figure 12: Airport comparison using boxplots for the month of February from 2013 to 2019

•

  Month of the year, • Day of the month, • Day of the week, • Hour of the day, • Number of passengers at this hour, • Number of open booths at this hour, • Number of flights at this hour. • Ratio of passengers per open booth Intuitively the waiting time at border security seems to depend on the state of the border area in the previous hour as well, e.g. if not all previously arrived passengers were processed, therefore the following features can be added to the base set: • Number of passengers at the previous hour, • Number of open booths at the previous hour, • Number of flights at the previous hour,

  (a) Comparison of the average wait time for US citizens per hour across airports (b) Comparison of the average wait time for non-US citizens per hour across airports. In red is indicated the median wait time for US citizens from Figure 13a. (c) Comparison of the ratio of average wait time for non-US citizens over US citizens. In red is indicated the interval [1, 2.5] for a better visualization.

Figure 13 :

 13 Figure 13: US vs non-US wait times comparison using boxplots over the year 2018

Figure 14 :

 14 Figure 14: Mean absolute error comparison between one-hot encoding for airports and different regressors per airports

  (a) R 2 score (b) Mean absolute error

Figure 15 :

 15 Figure 15: Box plot comparison of the performance measures for the five chosen benchmarks and for the six chosen regressors when predicting the average wait time for all passengers

  (a) Median mean squared error (b) Average mean squared error

Figure 16 :

 16 Figure 16: Evolution of the median and average regressors performance with the beginning of the training set

Figure 17 :

 17 Figure 17: Comparison per terminal of entry of the regressors' mean absolute error with the average wait time standard deviation over the year 2018

Table I :

 I Terminal of arrivals abbreviations

	Abbreviation	Airport (IATA)	Terminal name
	ATL -CE	ATL	Concourse E
	ATL -CF	ATL	Concourse F
	AUS -M	AUS	Main
	BOS -TE	BOS	Terminal E
	BWI -IA	BWI	International Arrivals
	CLT -M	CLT	Main
	DEN -I	DEN	International
	DFW -TD	DFW	Terminal D
	DTW -MT	DTW	McNamara Terminal
	DTW -NT	DTW	North Terminal
	EWR -TB	EWR	Terminal B
	EWR -TC	EWR	Terminal C
	FAT -M	FAT	Main
	FLL -T1	FLL	Terminal 1
	FLL -T4	FLL	Terminal 4
	GUM -MT	GUM	Main Terminal
	HNL -MOT	HNL	Main Overseas Terminal
	IAD -IA	IAD	International A
	IAH -I	IAH	IAB
	JFK -T1	JFK	Terminal 1
	JFK -T4	JFK	Terminal 4 (IAT)
	JFK -T5B	JFK	Terminal 5 (Jet Blue)
	JFK -T7	JFK	Terminal 7 (British)
	JFK -T8	JFK	Terminal 8 (American)
	LAS -T3	LAS	Terminal 3
	LAX -S2	LAX	Satellite 2
	LAX -S5	LAX	Satellite 5
	LAX -S7	LAX	Satellite 7
	LAX -T4	LAX	Terminal 4
	LAX -TBIT	LAX	Tom Bradley International Terminal
	MCO -A1	MCO	Airside 1
	MCO -A4	MCO	Airside 4
	MDW -MT	MDW	Main Terminal
	MIA -CT	MIA	Central Terminal
	MIA -NT	MIA	North Terminal
	MIA -ST	MIA	South Terminal
	MSP -T1L	MSP	Terminal 1 Lindbergh
	MSP -T2H	MSP	Terminal 2 Humphrey
	OAK -M	OAK	Main
	ORD -T5	ORD	Terminal 5
	PBI -M	PBI	Main
	PDX -M	PDX	Main
	PHL -TA	PHL	Terminal A
	PHX -M	PHX	Main
	RDU -T2	RDU	Terminal 2
	SAN -M	SAN	Main
	SAT -M	SAT	Main
	SEA -SS	SEA	South Satellite
	SFB -TA	SFB	Terminal A
	SFO -TA	SFO	Terminal A
	SFO -TG	SFO	Terminal G
	SJC -M	SJC	Main
	SJU -SJA	SJU	San Juan AA
	SLC -M	SLC	Main
	SMF -MT	SMF	Main Terminal
	SNA -TC	SNA	Terminal C
	SPN -IA	SPN	International Arrivals
	STL -M	STL	Main
	TPA -M	TPA	Main

ACKNOWLEDGMENTS

The authors would like to thank Nikunj Oza from NASA-Ames, the BDAI team from Verizon, Palo Alto as well as the French government for their financial support.


Barbara Gris

Chong Chen
email: chench@lsec.cc.ac.cn

Image reconstruction through metamorphosis

Keywords: Inverse problem, indirect registration, metamorphosis

The paper describes a method for reconstructing an image from noisy and indirect observations by registering, via metamorphosis, a template. The image registration part consists of two components, one is a geometric deformation that moves intensities without changing them and the other that changes intensity values. Unlike a registration with only geometrical deformation, this framework gives good results also when intensities of the template are poorly chosen. It also allows for appearance of a new structure. The approach is applicable to general inverse problems in imaging and we prove existence, stability and convergence, which implies that the method is a well-defined regularisation method. We also present several numerical examples from tomography.

Introduction

A key difficulty in shape based and/or spatiotemporal image reconstruction is to match an image against an indirectly observed target (indirect image registration). In the following, we give a brief overview of these notions along with a short survey of existing results.

Shape based reconstruction

The goal in shape based reconstruction is to recover shapes (morphologies) of interior sub-structures of an object whereas variations within these is of less importance. Such imaging studies arise in nano-characterisation of a specimen by means of electron microscopy [START_REF] Bladt | Electron tomography based on highly limited data using a neural network reconstruction technique[END_REF] or quantification of sub-resolution porosity in materials by means of x-ray phase contrast imaging.

In these imaging applications it makes sense to account for qualitative prior shape information during the reconstruction. Enforcing an exact spatial match between a template and the reconstruction is often unrealistic since shape information is almost always approximate, so the natural approach is to perform reconstruction assuming the structures are 'shape wise similar'.

Spatiotemporal imaging Here one images an object that changes over time, i.e., one needs to recover both the object structure and its time variation from noisy time series of indirect observations (measured data). An important case is when the object is the only entity that varies with time.

Spatiotemporal imaging arises when imaging moving organs in medicine. For more details on this problem, the reader is referred to [START_REF] Chen | A new variational model for joint image reconstruction and motion estimation in spatiotemporal imaging[END_REF] and the references therein. It is in particular relevant for techniques like positron emission tomography (PET) and single photon emission computed tomography (SPECT) that are used to image the distribution of injected radiopharmaceuticals (activity map). The latter is an inherently dynamic quantity since anatomical structures, like the heart and lungs, undergo motion during the data acquisition, see [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF] for a survey of organ motion models. Disregarding organ motion is known to degrade the spatial localisation of the radiotracer, leading to spatially blurred images. To exemplify the above mentioned issues, consider SPECT based cardiac perfusion studies and 18 F-FDG-PET imaging of lung nodules/tumours. The former needs to account for the beating heart and the latter for respiratory motion of the lungs and thoracic wall. Studies show a maximal displacement of 23 mm (average 15-20 mm) due to respiratory motion [START_REF] Schwarz | Implications of respiratory motion for the quantification of 2D MR spectroscopic imaging data in the abdomen[END_REF] and 42 mm (average 8-23 mm) due to cardiac motion in thoracic PET [START_REF] Wang | Cardiac motion of coronary arteries: Variability in the rest period and implications for coronary MR angiography[END_REF]. Finally, even when organ motion can be neglected, there are other dynamic processes, such as the uptake and wash-out of radiotracers from body organs. Visualising such kinetics of the radiotracers can actually be a goal in itself as is the case in dynamic PET, which is commonly used in pre-clinical imaging studies related to drug discovery/development [START_REF] Gullberg | Dynamic single photon emission computed tomographybasic principles and cardiac applications[END_REF].

Indirect image registration (matching)

The aim in image registration is to deform a template image so that it matches a target image. This becomes challenging when the template is allowed to undergo non-rigid deformations.

Diffeomorphic image registration is a well developed framework that meets this purpose. Here, the image registration problem is recast as the problem of finding a suitable diffeomorphism that deforms the template by a group action into the target image [START_REF] Younes | Shapes and Diffeomorphisms[END_REF][START_REF] Bauer | Overview of the geometries of shape spaces and diffeomorphism groups[END_REF]. The underlying assumption is that the target image is contained in the orbit of the template under the group action of diffeomorphisms. This can be stated in a very general setting where diffeomorphisms act on various deformable objects, like landmark points, curves, surfaces, scalar images, or even vector/tensor valued images.

The registration problem becomes even more challenging in indirect image registration, which is the case when the target is only known indirectly through measured data. Application of small diffeomorphic deformations in the context of indirect image registration is presented in [START_REF] Öktem | Shape based image reconstruction using linearised deformations[END_REF]. One can in a similar manner use large deformation diffeomorphic metric mapping (LDDMM) framework in indirect image registration as shown in [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF][START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF].

Overview of paper and specific contributions

The paper adapts the metamorphosis framework [START_REF] Trouvé | Metamorphoses through Lie group action[END_REF] to the indirect setting. Metamorphosis extends the LDDMM framework (diffeomorphometry) [START_REF] Younes | Shapes and Diffeomorphisms[END_REF][START_REF] Miller | Diffeomorphometry and geodesic positioning systems for human anatomy[END_REF] by allowing diffeomorphic changes to both the geometry of the template and its grey-scale values (intensities). As detailed below, an important feature of metamorphosis compared to diffeomorphometry is its capability to reconstruct images even when the template has erroneous grey-scale values that lack geometric features, e.g., in the form of a noisy background. We also show how indirect registration by metamorphosis formally defines a regularisation method for inverse problems in imaging, satisfying properties of existence, stability and convergence. Next, we study robustness with respect to choice of regularisation parameters. Finally, we provide an algorithm and its implementation1 that relies on a Fourier transforms for kernel convolutions and a Euler integration scheme for the temporal integration. Note that other approaches could be developed, like [START_REF] Neumayer | Regularization of inverse problems via time discrete geodesics in image spaces[END_REF] where the time discrete path method in [START_REF] Effland | Image extrapolation for the time discrete metamorphosis model: Existence and applications[END_REF] is extended to the indirect setting, or adaptation of the robust schemes for diffeomorphic registration developed in [START_REF] Mang | Constrained Hˆ1-regularization schemes for diffeomorphic image registration[END_REF][START_REF] Mang | A Lagrangian Gauss-Newton-Krylov solver for mass-and intensity-preserving diffeomorphic image registration[END_REF].

We start by recalling necessary theory from LDDMM-based indirect registration (section 3). Using these notions, we next adapt the metamorphosis framework to the indirect setting (section 4) and present our proof that this is a regularisation method (section 4.3). Its numerical implementation is outlined in section 4.4. Finally, we present several numerical examples from indirect registration in 2D tomography (section 5). In particular, we give a preliminary result for motion reconstruction when the acquisition is done at multiple time points. We also study the robustness of our methods w.r.t. the choice of regularisation parameters.

3 Indirect diffeomorphic registration

Large diffeomorphic deformations

Following [START_REF] Arguillere | Shape deformation analysis from the optimal control viewpoint[END_REF], (large) diffeomorphic deformations are defined by flows of time-varying vector fields. More precisely, let Ω ⊂ R d be a fixed bounded domain and let X := L 2 (Ω, R) represent grey scale images on Ω. Next, let V denote a fixed Hilbert space of vector fields on R d . We will assume V ⊂ C p 0 (Ω), i.e., the vector fields are supported on Ω and p times continuously differentiable. Finally, L 1 ([0, 1], V ) denotes the space of time-dependent V -vector fields that are integrable, i.e.,

ν(t, • ) ∈ V and t → ν(t, • ) C p is integrable on [0, 1].
The following (semi) norm on L 1 ([0, 1], V ) will be frequently used:

ν p := 1 0 ν(t, • ) p V dt 1/p .
Here, • V is the norm on the Hilbert space V of vector fields that is given by the inner product.

The following proposition shows that velocity fileds in L 1 ([0, 1], V ) can be used to generate flows in Diff p 0 (Ω) (set of p-diffeomorphisms that are supported in Ω ⊂ R d , and if Ω is unbounded, tend to zero towards infinity). This yields a method for generating diffeomorphisms.

Proposition 1. Let ν ∈ L 1 ([0, 1], V ) and consider the ordinary differential equation (flow equation):

     d dt φ(t, x) = ν t, φ(t, x) φ(0, x) = x for any x ∈ Ω and t ∈ [0, 1]. (1)
Then, (1) has a unique absolutely continuous solution φ(t,

• ) ∈ Diff p 0 (R d ).
The above result is proved in [START_REF] Arguillere | Shape deformation analysis from the optimal control viewpoint[END_REF] and the unique solution of ( 1) is henceforth called the flow of ν. We henceforth use the notation ϕ ν s,t : R d → R d , which refers to

ϕ ν s,t := φ(t, • ) • φ(s, • ) -1 for s, t ∈ [0, 1] (2)
where φ : Ω → R d denotes the unique solution to (1). As stated next, the set of diffeomorphisms that are given by solving the flow equation (1) forms a group that is a complete metric space [START_REF] Arguillere | Shape deformation analysis from the optimal control viewpoint[END_REF]. Its elements are called large diffeomorphic deformations. Proposition 2. Let V ⊂ C p 0 (Ω) (p ≥ 1) be an admissible reproducing kernel Hilbert space (RKHS) and define

G V := φ : R d → R d | φ = ϕ ν 0,1 for some ν ∈ L 2 ([0, 1], V ) . (3)
Then G V forms a sub-group of Diff p 0 (R d ) that is a complete metric space under the metric

d G (φ 1 , φ 2 ) := inf ν 1 : ν ∈ L 1 ([0, 1], V ) and φ 1 = φ 2 • ϕ ν 0,1 = inf ν 2 : ν ∈ L 1 ([0, 1], V ) and φ 1 = φ 2 • ϕ ν 0,1 .
The group G V may act on X through the geometric group action, which defines a deformation operator

W : G V × X → X where W(φ, I 0 ) := I 0 • φ -1 . (4)
We conclude by characterizing the regularity properties of flows of velocity fields as well as the group action in (4), the proof is given in [START_REF] Bruveris | Geometry of image registration: The diffeomorphism group and momentum maps[END_REF]. These results will play an important role in what is to follow.

Proposition 3. Assume V ⊂ C p 0 (Ω) (p ≥ 1
) is a fixed admissible Hilbert space of vector fields on Ω and {ν n } n ⊂ L 2 ([0, 1], V ) a sequence that converges weakly to ν ∈ L 2 ([0, 1], V ). Then, the following holds with ϕ n t := ϕ ν n 0,t :

1. (ϕ n t ) -1 converges to (ϕ ν 0,t ) -1 uniformly w.r.t. t ∈ [0, 1] and uniformly on compact subsets of Ω ⊂ R d .

lim

n→∞ W(ϕ n t , I 0 ) -W(ϕ ν 0,t , I 0 ) X = 0 for any I 0 ∈ X.

Indirect image registration

Image registration (matching) is the task of deforming a given template image I 0 ∈ X so that it matches a given target image I * ∈ X [START_REF] Brown | A survey of image registration techniques[END_REF][START_REF] Zitova | Image registration methods: a survey[END_REF].

The above task can also be stated in an indirect setting, which refers to the case when the template I 0 ∈ X is to be registered against a target I * ∈ X that is only indirectly known through data g ∈ Y where g = A(I * ) + e.

(

) 5
In the above, A : X → Y (forward operator) is known and assumed to be differentiable and e ∈ Y is a single sample of a Y -valued random element that denotes the measurement noise in the data. Access to a shape prior (template), as in indirect registration, can have profound effect in solving challenging inverse problem in imaging. As shown in [START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF], tomographic imaging problems that are otherwise intractable (highly noisy and sparsely sampled data) can be successfully addressed using indirect registration even when the template has a shape that is far from the true unknown target image that gave rise to measured data. The next step is to specify what is meant by deforming a template image and we consider diffeomorphic (non-rigid) deformations, i.e., diffeomorphisms that deform images by acting on them through a group action.

LDDMM-based registration Image registration by large diffeomorphic (non-rigid) deformations is typically formulated as minimising the functional

G V φ → γ 2 d G (Id, φ) 2 + W(φ, I 0 ) -I * 2 X for given γ > 0.
If V is admissible, then one can show that minimising the above functional on G V is equivalent to minimising the following functional on L 2 ([0, 1], V ): [40, Theorem 11.2 and Lemma 11.3]:

L 2 ([0, 1], V ) ν → γ 2 ν 2 2 + W(ϕ ν 0,1 , I 0 ) -I * 2 X given γ > 0.
Such a reformulation is advantageous since L 2 ([0, 1], V ) is a vector space, whereas G V is not. Hence, it is easier to devise algorithms for minimising a functional over L 2 ([0, 1], V ) rather than over G V . As shown in [START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF], one can extend the above to the indirect setting where the image target I * is not observed but only a data target g (LDDMMbased indirect registration). We simply need to replace the data fidelity term W(ϕ ν 0,1 , I 0 ) -I * 2 X with one that is adapted for indirect observations, which will involve the forward operator A and a data negative log-likelihood L : Y × Y → R. More precisely, the corresponding indirect registration problem can be adressed by minimising the functional

L 2 ([0, 1], V ) ν → γ 2 ν 2 2 + L (A • W)(ϕ ν 0,t , I 0 ), g .
The functional L : Y × Y → R is typically given by an appropriate affine transform of the data negative log-likelihood [START_REF] Bertero | Iterative image reconstruction: a point of view[END_REF], in which case minimising f → L A(f ), g amounts to seeking a maximum likelihood solution of ( 5). An interpretation of the above is that the template image I 0 , which is assumed to be given a priori, acts as a shape prior when solving the inverse problem in [START_REF] Bladt | Electron tomography based on highly limited data using a neural network reconstruction technique[END_REF] and γ > 0 is a regularisation parameter that governs the influence of this shape priori against the need to fit measured data. This interpretation becomes more clear when one re-formulates LDDMM-based indirect registration as The LDDMM-based indirect registration with geometric group action outlined in section 3.2 is remarkably stable when the template has correct topology and intensity levels [START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF]. Usage of the geometric group action, however, prevents creating or removing intensity, e.g., it is not possible to start out from a template with a single connected structure and deform it to a image with two separated structures. Figure 1 shows the effect that erroneous intensity has on LDDMMbased indirect registration. This severely limits the usefulness of LDDMMbased indirect registration, e.g., spatiotemporal images (movies) are likely to involve changes in both geometry (objects appear or disappear) and intensity.

             min ν∈L 2 ([0,1],V ) γ 2 ν 2 2 + L (A • W) φ(1, • ), I 0 , g d dt φ(t, x) = ν t, φ(t, x) (t, x) ∈ Ω × [0, 1], φ(0, x) = x x ∈ Ω. (6)
One approach to address this drawback of LDDMM-based indirect registration is to replace the geometric group action with one that alters intensities, e.g., a mass preserving group action [START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF]. This allows to correct motion artefacts in PET as shown in [START_REF] Gigengack | Motion correction of cardiac pet using mass-preserving registration[END_REF][START_REF] Gigengack | Mass-preserving motion correction of pet: Displacement field vs. spline transformation[END_REF][START_REF] Gigengack | Motion correction in dual gated cardiac pet using masspreserving image registration[END_REF] but not to model changes in intensities that do not come from a geometric change. Another is to keep the geometric group action, but replace LDDMM with a framework that model the transformation of an image by a combination of a diffeomorphic deformation and intensities changes. The latter leads to metamorphosis-based indirect registration, which we next describe.

The metamorphosis framework

Just like in LDDMM, diffeomorphisms are generated by flows [START_REF] Arguillere | Shape deformation analysis from the optimal control viewpoint[END_REF]. The difference in metamorphosis is that these diffeomorphisms are combined with intensities changes. Hence, image registration by LDDMM is restricted to targets that lie in the orbit of the template under the diffeomorphic group action whereas in metamorphosis, the ability to change intensities allows one to reach any target in X [START_REF] Trouvé | Metamorphoses through Lie group action[END_REF]. As such, metamorphosis extends LDDMM.

To proceed, we start by the abstract definition of metamorphosis.

Definition 1 (Metamorphosis [START_REF] Trouvé | Metamorphoses through Lie group action[END_REF]). Let V ⊂ C p 0 (Ω) be an admissible Hilbert space, G V is defined as in (3), and '.' denotes some group action of G

V on X. A metamorphosis is a curve t → (φ t , J t ) in G V × X. The curve t → f t := φ t .J t is called the image part, t → φ t is the deformation part, and t → f t is the template part.
The image part represents the temporal evolution that is not related to intensity changes (evolution of the underlying geometry), whereas the template part is the evolution of the intensity. Both evolutions, which are combined in metamorphosis, are driven by the same underlying flow of diffeomorphisms in G V .

A important case is when the metamorphosis t → (φ t , f t ) has a deformation part that solves the flow equation ( 1) and a template part that is C 1 in time. More precisely, let L 2 ([0, 1], X) denote the space of functions in X that are square integrable, i.e.,

ζ(t, • ) ∈ X and t → ζ(t, • ) X ∈ L 2 ([0, 1], R).
A natural norm on L 2 ([0, 1], X) is given by

ζ 2 := 1 0 ζ(t, • ) 2 X dt 1/2 .
We will henceforth use the notation

L 2 ([0, 1], V × X) := L 2 ([0, 1], V ) × L 2 ([0, 1], X) .
Bearing in mind the above notation, for given (ν,

ζ) ∈ L 2 ([0, 1], V × X) and I 0 ∈ X, define the curve t → I ν,ζ t that is absolutely continuous on [0, 1] and solves      d dt I ν,ζ t (x) = ζ t, ϕ ν 0,t (x) I ν,ζ 0 (x) = I 0 (x) with ϕ ν 0,t ∈ G V as in (2). (7)
The metamorphosis can now be parametrised as t → (ϕ ν 0,t , I ν,ζ t ).

Indirect registration

The indirect registration problem in section 3.2 can be approached by metamorphosis instead of LDDMM. Similar to LDDMMbased indirect image registration, we define metamorphosis-based indirect image registration as the minimization of the objective functional

J γ,τ ( • ; g) : L 2 ([0, 1], V × X) → R defined as J γ,τ (ν, ζ; g) := γ 2 ν 2 2 + τ 2 ζ 2 2 + L A W(ϕ ν 0,1 , I ν,ζ 1 ) , g (8)
Here, I ν,ζ t : Ω → R is given by ( 7) and we are given fixed regularisation parameters γ, τ > 0, measured data g ∈ Y , and initial template I 0 ∈ X that defines the initial condition I ν,ζ 0 (x) := I 0 (x). Hence, performing metamorphosis-based indirect registration of a template I 0 against a target indirectly observed through data g amounts to solving ( ν, ζ) ∈ arg min

(ν,ζ) J γ,τ (ν, ζ; g). (9)
As shown in proposition 4, the above optimisation problem has a solution assuming the data discrepancy and the forward operator fulfils some weak requirements. From a solution we obtain the following:

Initial template: I 0 ∈ X such that I ν,ζ 0 := I 0 . Reconstruction: Final registered template f ν, ζ 1 = W ϕ ν 0,1 , I ν, ζ 1 ∈ X.
Image trajectory: The evolution of both geometry and intensity of the template, given by t → W ϕ ν 0,t , I ν, ζ t .

Template trajectory: The evolution of intensities of the template, i.e., the part that does not include evolution of geometry:

t → I ν, ζ t .
Deformation trajectory: The geometric evolution of the template, i.e., the part that does not include evolution of intensity: t → W(ϕ ν 0,t , I 0 ).

Regularising properties

In what follows, X := L 2 (Ω, R) and Y a Hilbert space. We will here prove several properties (existence, stability and convergence) of metamorphosisbased indirect image registration. Taken together, these show that metamorphosis-based indirect image registration is a well defined regularisation method in the sense of [START_REF] Grasmair | Generalized Bregman distances and convergence rates for non-convex regularization methods[END_REF].

Proposition 4 (Existence). Assume A : X → Y is continuous and the data discrepancy L( • , g) : Y → R is weakly lower semi-continuous for any g ∈ Y . Then, J γ,τ ( • , g) : L 2 ([0, 1], V × X) → R defined through (8) has a minimiser in L 2 ([0, 1], V × X) for any I 0 ∈ L 2 (Ω, R).
Proof. We follow here the strategy to prove existence of minimal trajectories for metamorphosis (as in [START_REF] Charon | Metamorphoses of functional shapes in Sobolev spaces[END_REF] for instance). One considers a minimising sequence of J γ,τ ( • ; g), i.e., a sequence that converges to the infimum of J γ,τ ( • ; g) (such a sequence always exists). The idea is to prove that such a minimising sequence has a sub-sequence that converges to a point in L 2 ([0, 1], V × X), i.e., the infimum is contained in L 2 ([0, 1], V × X) which proves existence of a minima.

Bearing in mind the above, we start by considering a minimising sequence

(ν n , ζ n ) n ⊂ L 2 ([0, 1], V × X) to J γ,τ ( • ; g), i.e., lim n→∞ J γ,τ (ν n , ζ n ; g) = inf ν,ζ J γ,τ (ν, ζ; g). Since ν n n ⊂ L 2 ([0, 1], V ) is bounded, it has a sub-sequence that converges to an element ν ∞ ∈ L 2 ([0, 1], V ). Likewise, ζ n n ⊂ L 2 ([0, 1], X) has a sub- sequence that converges to an element ζ ∞ ∈ L 2 ([0, 1], X).
Hence, with a slight abuse of notation, we conclude that

ν n ν ∞ and ζ n ζ ∞ as n → ∞.
The aim is now to prove existence of minimisers by showing that

(ν ∞ , ζ ∞ ) is a minimiser to J γ,τ ( • ; g) : L 2 ([0, 1], V × X) → R.
Before proceeding, we introduce some notation in order to simplify the expressions. Define

I n t := I ν n ,ζ n t and ϕ n s,t := ϕ ν n s,t for n ∈ N ∪ {∞}. (10)
Hence, assuming geometric group action (4) and using (2), we can write

J γ,τ (ν n , ζ n ; g) = γ 2 ν n 2 2 + τ 2 ζ n 2 2 + L A I n 1 • ϕ n 1,0 , g
for n ∈ N ∪ {∞}. Assume next that the following holds:

I n 1 • ϕ n 1,0 I ∞ 1 • ϕ ∞ 1,0 as n → ∞. (11)
The data discrepancy term L( • , g) : Y → R is weakly lower semi continuous and the forward operator A : X → Y is continuous, so L( • , g) • A is also weakly lower semi continuous and then [START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF] implies

L A(I ∞ 1 • ϕ ∞ 1,0 ), g ≤ lim inf n→∞ L(A(I n 1 • ϕ n 1,0 ), g). (12)
Furthermore, from the weak convergences of ν n and ζ n , we get

γ 2 ν ∞ 2 2 + τ 2 ζ ∞ 2 2 ≤ lim inf n→∞ γ 2 ν n 2 2 + τ 2 ζ n 2 2 . ( 13
)
Hence, combining [START_REF] Chen | Gradient-flow-based semi-implicit finite-element method and its convergence analysis for image reconstruction[END_REF] and ( 13) we obtain

J γ,τ (ν ∞ , ζ ∞ ; g) ≤ lim n→∞ J γ,τ (ν n , ζ n ; g). Since (ν n , ζ n ) n ⊂ L 2 ([0, 1], V × X) is a minimising sequence, this yields J γ,τ (ν ∞ , ζ ∞ ; g) = inf (ν,ζ)∈L 2 ([0,1],V ×X) J γ,τ (ν, ζ; g), which proves (ν ∞ , ζ ∞ ) ∈ L 2 ([0, 1], V × X) is a minimiser to J γ,τ ( • ; g).
Hence, to finalize the proof we need to show that (11) holds. We start by observing that the solution of ( 7) can be written as

I n t := I n 0 (x) + t 0 ζ n s, ϕ n 0,s (x) ds for n ∈ N ∪ {∞}, (14)
and note that (t, x)

→ I n t (x) ∈ C([0, 1] × Ω, R). Next, we claim that I n 1 I ∞ 1 for some I ∞ 1 ∈ X, which is equivalent to lim n→∞ I n 1 -I ∞ 1 , J = 0 for any J ∈ L 2 (Ω, R). ( 15
)
To prove [START_REF] Effland | Image extrapolation for the time discrete metamorphosis model: Existence and applications[END_REF], note first that since continuous functions are dense in L 2 , it is enough to show (15) holds for J ∈ C 0 (Ω, R). Next,

I n 1 -I ∞ 1 , J = Ω t 0 ζ n s, ϕ n 0,s (x) -ζ ∞ s, ϕ ∞ 0,s (x) J(x)dsdx (16)
=

Ω t 0 ζ n s, ϕ n 0,s (x) -ζ n s, ϕ ∞ 0,s (x) J(x)dsdx (17)
+ Ω t 0 ζ n s, ϕ n 0,s (x) -ζ ∞ s, ϕ ∞ 0,s (x) J(x)dsdx. (18)
Let us now take a closer look at the term in [START_REF] Gigengack | Mass-preserving motion correction of pet: Displacement field vs. spline transformation[END_REF]:

Ω t 0 ζ n s, ϕ n 0,s (x) -ζ n s, ϕ ∞ 0,s (x) J(x)dsdx = Ω t 0 ζ n (s, x)J ϕ n 0,s (x) Dϕ n 0,s (x) dsdx - Ω t 0 ζ ∞ (s, x)J ϕ ∞ 0,s (x) Dϕ ∞ 0,s (x) dsdx = Ω t 0 ζ n (s, x) J ϕ n 0,s (x) Dϕ n 0,s (x) -J ϕ ∞ 0,s (x) Dϕ ∞ 0,s (x) dsdx - Ω t 0 ζ ∞ (s, x) -ζ n (s, x) J ϕ ∞ 0,s (x) Dϕ ∞ 0,s (x) dsdx = ζ n , J n -J ∞ -ζ ∞ -ζ n , J ∞ where J n ∈ L 2 ([0, 1], X) is defined as J n (s, x) := J ϕ n s,0 (x) Dϕ n s,0 (x) for n ∈ N ∪ {∞}. ( 19
)
By proposition 3 we know that ϕ n s,0 → ϕ ∞ s,0 and Dϕ n s,0 → Dϕ ∞ s,0 uniformly on Ω. Since J is continuous on Ω, we conclude that J n -J ∞ 2 → 0. When combined with the boundedness of ζ n , we get

ζ n , J n -J ∞ ≤ ζ n 2 • J n -J ∞ 2 → 0. Furthermore, since ζ n ζ ∞ , we also get ζ ∞ -ζ n , J ∞ → 0.
Hence, we have shown that [START_REF] Gigengack | Mass-preserving motion correction of pet: Displacement field vs. spline transformation[END_REF] tends to zero, i.e., lim

n→∞ Ω t 0 ζ n s, ϕ n 0,s (x) -ζ n s, ϕ ∞ 0,s (x) J(x)dsdx = 0.
Finally, we consider the term in [START_REF] Gigengack | Motion correction in dual gated cardiac pet using masspreserving image registration[END_REF]. Since ζ n ζ ∞ , we immediately obtain

Ω t 0 ζ n s, ϕ ∞ s (x) -ζ ∞ s, ϕ ∞ s (x) J(x)dsdx = ζ n -ζ ∞ , J ∞ → 0.
To summarise, we have just proved that both terms ( 17) and ( 18) tend to 0 as n → ∞, which implies that (15) holds, i.e., I n 1

I ∞ 1 . To prove (11), i.e., I n 1 • ϕ n 1,0 I ∞ 1 • ϕ ∞ 1,0 , we need to show that lim n→∞ I n 1 • ϕ n 1,0 -I ∞ 1 • ϕ ∞ 1,0 , J = 0 for any J ∈ L 2 (Ω, R), (20)
and as before, we may assume J ∈ C 0 (Ω, R). Using [START_REF] Grasmair | Generalized Bregman distances and convergence rates for non-convex regularization methods[END_REF] we can express the term in ( 20) whose limit we seek as

I n 1 • ϕ n 1,0 -I ∞ 1 • ϕ ∞ 1,0 , J ≤ I n 1 , J n (1, • ) -J ∞ (1, • ) + I n 1 -I ∞ 1 , J ∞ (1, • ) ≤ I n 1 • J n (1, • ) -J ∞ (1, • ) + I n 1 -I ∞ 1 , J ∞ (1, • ) .
Since I n 1 is bounded (because ζ n is bounded) and since I n 1 I ∞ 1 (which we shoed before), all terms above tend to 0 as n → ∞, i.e., [START_REF] Gullberg | Dynamic single photon emission computed tomographybasic principles and cardiac applications[END_REF] holds.

This concludes the proof of [START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF], which in turn implies the existence of a minimiser of J γ,τ ( • ; g).

Our next result shows stability in the sense that the solution to the indirect registration problem is (weakly) continuous w.r.t. variations in data.

Proposition 5 (Stability). Let {g k } k ⊂ Y and assume this sequence converges (in norm) to some g ∈ Y . Next, for each γ, τ > 0 and each k, define

(ν k , ζ k ) ∈ L 2 ([0, 1], V × X) as (ν k , ζ k ) = arg min (ν,ζ) J γ,τ (ν, ζ; g k ).
Then there exists a sub sequence of (ν k , ζ k ) that converges weakly to a minimiser of J γ,τ ( • ; g) in [START_REF] Bruveris | Mixture of kernels and iterated semidirect product of diffeomorphisms groups[END_REF].

Proof. J γ,τ ( • ; g k ) has a minimiser (ν k , ζ k ) ∈ L 2 ([0, 1], V × X) for any g k ∈ Y (proposition 4
). The idea is first to show that the sequences (ν k ) k and (ζ k ) k are bounded. Next, we show that there exists a weakly converging subsequence of (ν k , ζ k ) that converges to a minimiser (ν, ζ) of J γ,τ ( • ; g), which also exists due to proposition 4.

Since (ν k , ζ k ) minimises J γ,τ ( • ; g k ), we have by ( 8) that

ν k 2 2 ≤ 2 γ J γ,τ ( • ; g k )(ν k , ζ k ) ≤ 2 γ J γ,τ ( • ; g k )(0, 0) for each k. ( 21
)
Observe now that if ν = 0 and ζ = 0, then ϕ ν 0,1 = Id by (1) and I ν,ζ 1 = I 0 by [START_REF] Bruveris | Geometry of image registration: The diffeomorphism group and momentum maps[END_REF], so in particular W ϕ ν 0,1 , I ν,ζ 1 ) = I 0 whenever ν = 0 and ζ = 0. Hence, J γ,τ ( • ; g k )(0, 0) = L A(I 0 ), g k and, in addition, ν 2 = 0 and ζ 2 = 0, so (21) becomes

ν k 2 2 ≤ 2 γ L A(I 0 ), g k → L(A(I 0 ), g) as k → ∞. (22)
In conclusion, the sequence

(ν k ) k ⊂ L 2 ([0, 1], V ) is bounded. In a similar way, we can show that (ζ k ) k ⊂ L 2 ([0, 1], X) is bounded.
The boundedness of both sequences implies that there are sub sequences to these that converge weakly to some elements

ν ∞ ∈ L 2 ([0, 1], V ) and ζ ∞ ∈ L 2 ([0, 1], X), respectively. Thus, to complete the proof, we need to show that (ν ∞ , ζ ∞ ) ∈ L 2 ([0, 1], V × X) minimises J γ,τ ( • ; g), i.e., that J γ,τ (ν ∞ , ζ ∞ ; g) ≤ J γ,τ (ν, ζ; g) holds for any (ν, ζ) ∈ L 2 ([0, 1], V × X).
From the weak convergences, we obtain

γ 2 ν ∞ 2 2 + τ 2 ζ ∞ 2 2 ≤ γ 2 lim inf k ν k 2 2 + τ 2 lim inf k ζ k 2 2 ≤ 1 2 lim inf k γ ν k 2 2 + τ ζ k 2 2 . ( 23
)
The weak convergence also implies (see proof of proposition 4) that

W ϕ k 0,1 , I ∞ 1 W ϕ ∞ 0,1 , I ∞ 1 in X.
In the above, we have used the notational convention introduced in [START_REF] Chen | A new variational model for joint image reconstruction and motion estimation in spatiotemporal imaging[END_REF]. By the lower semi-continuity of L, we get

L A W(ϕ ∞ 0,1 , I ∞ 1 ) , g ≤ lim inf k L A W(ϕ k 0,1 , I k 1 ) , g k . (24)
Hence,

J γ,τ (ν ∞ , ζ ∞ ; g) = γ 2 ν ∞ 2 2 + τ 2 ζ ∞ 2 2 + L A W(ϕ ∞ 0,1 , I ∞ 1 ) , g . ≤ 1 2 lim inf k γ ν k 2 2 + τ ζ k 2 2 + lim inf k L A W(ϕ k 0,1 , I k 1 ) , g k ≤ lim inf k J γ,τ (ν k , ζ k ; g k ). (25) Next, since (ν k , ζ k ) ∈ L 2 ([0, 1], V × X) minimises J γ,τ ( • ; g k ), we get J γ,τ (ν ∞ , ζ ∞ ; g) ≤ lim inf k J γ,τ (ν, ζ; g k ) for any (ν, ζ) ∈ L 2 ([0, 1], V × X). Furthermore, J γ,τ (ν, ζ; g k ) → J γ,τ (ν, ζ; g), so J γ,τ (ν ∞ , ζ ∞ ; g) ≤ J γ,τ (ν, ζ; g) for all (ν, ζ) ∈ L 2 ([0, 1], V × X).
In particular, we have shown that (ν ∞ , ζ ∞ ) minimises J γ,τ ( • ; g).

Our final results concerns convergence, which investigates the behaviour of the solution as data error tends to zero and regularisation parameters are adapted accordingly through a parameter choice rule against the data error.

Proposition 6 (Convergence). Let g ∈ Y and assume

A W(ϕ ν 0,1 , I ν,ζ 1 ) = g for some (ν, ζ) ∈ L 2 ([0, 1], V × X).
Next, for parameter choice rules δ → γ(δ) and δ → τ (δ) with δ > 0, define

(ν δ , ζ δ ) ∈ arg min (ν,ζ) J γ(δ),τ (δ) (ν, ζ; g + e δ )
where e δ ∈ Y (data error) has magnitude e δ ≤ δ. Finally, assume that δ → γ(δ)/τ (δ) and δ → τ (δ)/γ(δ) are bounded, and

lim δ→0 γ(δ) = lim δ→0 τ (δ) = lim δ→0 δ 2 γ(δ) = lim δ→0 δ 2 τ (δ) = 0.
Then, for any sequence δ k → 0 there exists a sub-sequence δ k such that

(ν δ k , ζ δ k ) converges weakly to a (ν * , ζ * ) satisfying A W(ϕ ν * 0,1 , I ν * ,ζ * 1 ) = g.
Proof. Let (δ k ) be a sequence converging to 0 and, for each k, let us denote

g k := g + e δ k , ν k := ν δ k , and ζ k := ζ δ k .
Similarly to previous proofs, we will show that the sequences (ν k ) and (ζ k ) are bounded, and then that the weakly converging subsequence that can be extracted from (ν k , ζ k ) converges to a suitable solution. Define γ k := γ(δ k ) and τ k := γ(δ k ). Then, for each k we have

ν k 2 2 ≤ 1 γ k J γ k ,τ k ,g k (ν k , ζ k ) ≤ 1 γ k J γ k ,τ k ,g k ( ν, ζ) = 1 γ k γ k ν 2 2 + τ k ζ 2 2 + L(g, g k ) ≤ ν 2 2 + τ k γ k ζ 2 2 + δ k γ k .
From the assumptions on the parameter choice rules, we conclude that

(ν k ) ⊂ L 2 ([0, 1], V ) is bounded. Similarly, one can show that (ζ k ) ⊂ L 2 ([0, 1], X) is bounded.
From the above, we conclude that there is a subsequence of (ν k , ζ k ) that converges weakly to ( ν,

ζ) in L 2 ([0, 1], V ) × L 2 ([0, 1], V ). Then (see proof of proposition 4) L A W(ϕ ν 0,1 , I ν, ζ 1 ) , g ≤ lim inf k L A(W ϕ ν k 0,1 , I ν, ζ 1 ) , g k .
Furthermore, the above quantity converges to 0 since

L A W(ϕ ν k 0,1 , I ν k ,ζ k 1 ) , g k ≤ J γ k ,τ k ,g k (ν k , ζ k ) ≤ J γ k ,τ k ,g k ( ν, ζ) = γ k ν 2 2 + τ k ζ 2 2 + L(g, g k ) → 0 and k → ∞.
Hence, A W(ϕ ν 0,1 , I ν, ζ 1 ) = g.

Numerical implementation

In order to solve (9), we use a gradient descent scheme on the control variable (ν, ζ) ∈ L2 ([0, 1], V × X) with a uniform discretisation of the interval [0, 1] into N parts, i.e., t i = 1/N for i = 0, . . . , N and the gradient descent is performed on ν(t i , • ), ζ(t i , • ) , for i = 0, 1, . . . , N . we use a Euler scheme on this discretisation for computing the numerical integrations. The flow equation ( 1) is computed using the following approximation with small deformations:

ϕ ν t i ,0 ≈ ϕ ν t i-1 ,0 • Id - 1 N ν(t i-1 , • ) .
Algorithm 1 presents the implementation 2 for computing the gradient of J based on expressions from appendix A. Note that the number 'N ' is not a fixed value, which is just a value for numerical discretisation. As long as such value is sufficiently large, we can obtain a small linearised diffeomorphic deformation. Then the large diffeomorphic deformation can be obtained by compositions of these small linearised deformations [START_REF] Younes | Shapes and Diffeomorphisms[END_REF]. The computation of the Jacobian determinant Dϕ ν t i ,1 (x) at each time point is based on the following approximation similar to the one used in [START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF]:

Dϕ ν t i ,1 ( • ) ≈ 1 + 1 N div ν(t i , • ) Dϕ ν t i+1 ,1 • Id + 1 N ν(t i , • ) .
By simple derivations for the 2D case, we have

Dϕ ν t i ,1 ( • ) = 1 + 1 N div ν(t i , • ) + 1 N 2 | Dν(t i , • )| Dϕ ν t i+1 ,1 • Id + 1 N ν(t i , • ) .
It is easy to see that the error is about O(N -2 ) for the above approximation. Since all of the numerical implementations are based on the discretisations in temporal and spatial domains, the exact computation to the Jacobian determinant is hardly required. Even if we use the exact one, it might not increase the overall precision of the algorithm. Moreover, if the N is large enough, this error can be neglected.

Application to 2D tomography 5.1 The setting

Let X = L 2 (Ω, R) and its elements represent 2D images on a fixed bounded domain Ω ⊂ R 2 . In the application shown here, diffeomorphisms act on X through the geometric group action in (4). The goal is to register a given differentiable template image I 0 ∈ X against a target that observed indirectly as in [START_REF] Bladt | Electron tomography based on highly limited data using a neural network reconstruction technique[END_REF].

The forward operator The forward operator A : X → Y in 2D tomographic imaging is the 2D ray/Radon transform, i.e.,

A(f )(ω, x) = R f (x + sω)ds for ω ∈ S 1 and x ∈ ω ⊥ .
Here, S 1 is the unit circle, so (ω, x) encodes the line s → x + sω in R 2 with direction ω through x. The data manifold M is the set of such lines that are included in the measurements, i.e., M is given by the experimental set-up, so Y = L 2 (M, R). We will consider parallel lines in R 2 (parallel beam data), i.e., tomographic data g ∈ Y are noisy digitised values of an L 2 -function on this manifold. The forward operator is linear, so it is Gateaux differentiable and the adjoint of its derivative is given by the backprojection [START_REF] Natterer | Mathematical Methods in Image Reconstruction[END_REF][START_REF] Markoe | Analytic Tomography, volume 106 of Encyclopedia of mathematics and its applications[END_REF].

The squared 2-norm corresponds to the data likelihood when data is corrupted by additive Gaussian noise:

L : Y × Y → R with L(g, h) = g -h 2 2 .
The noise level in data is specified by the peak signal-to-noise ratio (PSNR), which is defined as PSNR(g) = 10 log 10 g 0 -g 0 2

e -e 2 for g = g 0 + e.

In the above, g 0 is the noise-free part and e is the noise component of data with g 0 and e denoting the mean of g 0 and e, respectively. The PSNR is expressed in terms of dB.

Algorithm 1 Computation of ∇ J (ν, ζ).

Require: ν(t i , • ) and ζ(t i , • ) with t i ← i/N for i = 0, 1, . . . , N .

1: for i = 1, . . . , N do Compute ζ(t i , • ) • ϕ ν 0,ti 2: temp ← ζ(t i , • ) 3:
for j = i -1, . . . , 0 do 4:

temp ← temp • Id + 1 N ν(t j , • ) 5:
end for 6:

ζ(t i , • ) • ϕ ν 0,ti ← temp 7: end for 8: for i = 1, . . . , N do Compute f ν,ζ (t i , • ) := I ν,ζ (t i , • ) • ϕ ν ti,0 9
:

I ν,ζ (t i , • ) ← I 0 + i-1 j=0 I ν,ζ (t j , • ) + 1 N ζ(t j , • ) • ϕ ν 0,tj 10
:

I ν,ζ (t i , • ) • ϕ ν 0,0 ← I ν,ζ (t i , • ) 11:
for j = 1, . . . , i do 12:

I ν,ζ (t i , • ) • ϕ ν tj ,0 ← I ν,ζ (t i , • ) • ϕ ν 0,tj-1 • Id -1 N ν(t j-1 , • ) 13:
end for 14: end for 15: for i = 1, . . . , N do Compute I 0 • ϕ ν ti,0 16:

I 0 • ϕ ν 0,0 ← I 0 • ϕ ν 0,0 = I 0 17: I 0 • ϕ ν ti,0 ← I 0 • ϕ ν ti-1,0 • Id -1 N ν(t i-1 , • ) 18: end for 19: for i = 1, . . . , N do 20: G(t i , • ) ← ∇(I 0 • ϕ ν ti,0 ) + ti-1 j=0 1 N ∇(ζ(t j , • ) • ϕ ν ti,tj ) 21: end for 22: Dϕ ν t N ,1 = Dϕ ν 1,1 = 1 Compute Dϕ ν ti,1
23: for i = N -1, . . . , 0 do 24:

Dϕ ν ti,1 ← 1 + 1 N div ν(t i , • ) Dϕ ν ti+1,1 • Id + 1 N ν(t i , • ) 25: end for 26: ∇ L f ν,ζ (1, • ), g ϕ ν t N ,1 ← ∇ L f ν,ζ (1, • ), g 27: for i = N -1, . . . , 0 do Compute ∇ L f ν,ζ (1, • ), g ϕ ν ti,1 28: ∇ L f ν,ζ (1, • ), g ϕ ν ti,1 ← ∇ L f ν,ζ (1, • ), g ϕ ν ti+1,1 • Id + 1 N ν(t i , • ) 29: end for 30: for i = 1, . . . , N do Compute ∇ J (ν, ζ) 31: ∇ ν J γ,τ (ν, ζ, g)(t i , , • ) ← 2γν(t i , , • ) - Ω K(x, • ) Det(dϕ ν ti,1 (x)) ∇ L f ν,ζ (1, • ), g ϕ ν ti,1 (x) G(t i , x)dx 32: ∇ ζ J γ,τ (ν, ζ)(t i , , • ) ← 2τ ζ(t i , , • ) + | Det(dϕ ν ti,1 ) ∇ L f ν,ζ (1, • ), g ϕ ν ti,1 (x) G(t i , • ) 33: end for 34: return ∇ J (ν)(t i , • ), ∇ J (ζ)(t i , • ) for i = 1, . . . , N .
Joint tomographic reconstruction and registration Under the geometric group action (4), metamorphosis based-indirect registration reads as 8), i.e., given regularisation parameters γ, τ ≥ 0 and initial template I 0 ∈ X we solve min

f ν, ζ 1 = W ϕ ν 0,1 , I ν, ζ 1 = I ν, ζ 1 • ϕ ν 1,0 where ( ν, ζ) ∈ L 2 ([0, 1], V × X) minimises (
(ν,ζ) γ 2 ν 2 2 + τ 2 ζ 2 2 + A f 1, φ(1, • ) -1 -g 2 2                      d dt f (t, x) = ζ t, φ(t, x) f (0, x) = I 0 (x) d dt φ(t, x) = ν t, φ(t, x) φ(0, x) = x. (26)
We will consider a set V of vector fields that is an RKHS with a reproducing kernel represented by symmetric and positive definite Gaussian. Then V is admissible and is continuously embedded in L 2 (Ω, R 2 ). The kernel we pick is

K σ : Ω × Ω → R 2×2 + K σ (x, y) := exp - 1 2σ 2 x -y 2 1 0 0 1 for x, y ∈ R 2 and σ > 0. ( 27
)
The kernel-size σ also acts as a regularisation parameter.

Overview of experiments

In the following we test various aspects of using metamorphoses based indirect registration for joint tomographic reconstruction and registration. The tomographic inverse problem along with characteristics of the data are outlined in section 5.1. The results are obtained by solving (26) via a gradient descent, see appendix A for the computation of the gradient of the objective. For each reconstruction, we list the number of angles of the parallel beam ray transform, the kernel-size σ in [START_REF] Natterer | Mathematical Methods in Image Reconstruction[END_REF], and the two regularisation parameters γ, τ > 0 appearing in the objective functional in [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF]. The first test (section 5.3) aims to show how metamorphoses based indirect registration handles a template with intensities that differ from those of the target. Section 5.4 considers the ability to handle an initial template with a topology that does not match the target. This is essential when one has simultaneous geometric and topological changes. As an example, in spatiotemporal imaging it may very well be the case that geometric deformation takes place simultaneously as new masses appear or disappear. Next, section 5.5 studies the robustness of the solutions with respect to variations in the regularisation parameters. Finally, section 5.6 shows how indirect registration through metamorphoses can be used to recover a temporal evolution of a given template registered against time series of data. This is an essential part of spatio-temporal tomographic reconstruction.

Sections 5.3 to 5.5 have a common setting in that grey scale images in the reconstruction space are discretised using 256 × 256 pixels supported in a rectangular region Ω. The tomographic data is noisy samples of the 2D parallel beam ray transform of the target sampled at 100 angles uniformly distributed angles in [0, π] with 362 lines/angle. Data is corrupted with additive Gaussian noise with differing noise levels. In the last section 5.6 we add a Poisson noise to the data in order to avoid inverse crime [START_REF] Wirgin | The inverse crime[END_REF] and also test robustness against misspecification of the noise model.

Consistent topology and inconsistent intensities

Here, topology of the template is consistent with that of the target, but intensities differ. The template, which is shown in fig. 2(a), is registered against tomographic data shown in fig. 2(c). The (unknown) target used to generate data is shown in fig. 2(b). Also, data has a noise level corresponding to a PSNR of 15.6 dB and the kernel size σ is set to 1/16:th of the side length of the image domain Ω. The final reconstruction is shown in fig. 2(h), which is to be compared against the target in fig. 2(b). Figure 2 also shows image, deformation and template trajectories.

We clearly see that metamorphosis based indirect registration can handle a template with wrong intensities. As a comparison, see fig. 1(c) for the corresponding LDDMM based indirect registration using the same template and data. Furthermore, the different trajectories also provides easy visual interpretation of the influence of geometric and intensity deformations.

Inconsistent topology and intensities

Here, both topology and intensities of the template differ from those in the target. The template, which is shown in fig. 3(a), is registered against tomographic data shown in fig. 3(c). The (unknown) target used for generating the data is shown in fig. 3(b). Also, data has a noise level corresponding to a PSNR of 10.6 dB and the kernel size σ is set to 1/16:th of the side length of the rectangular image domain Ω. The final reconstruction is shown in fig. 3(h), which is to be compared against the target in fig. 3(b). Figure 3 also shows image, deformation and template trajectories.

We clearly see that metamorphosis based indirect registration can handle a template where both intensities and the topology are wrong. In particular, we can see follow both the deformation of the template and the appearance of the white disc.

Robustness

Metamorphosis based indirect registration, which amounts to solving [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF], requires selecting three parameters: the kernel-size σ and the two regularisation parameters γ and τ . Here we study the influence of these parameters on the final registered image (reconstruction). We conclude that the reconstruction is very robust with respect to the choice of these three parameters. This is illustrated in figs. [START_REF] Bertero | Iterative image reconstruction: a point of view[END_REF] and 5 where the reconstructed image is in the left-most column. The metamorphosis framework provides, in addition to this reconstruction, a decomposition into a template part and a deformation part. However, unlike reconstruction, the task of decomposition is sensitive to the choice of these parameters.

We first study the influence of the choice of the two regularisation parameters γ and τ by considering the setup in section 5.3 for varying choices of these parameters. Two extremal solutions are valid for the change of the shape of the triangles: it can be generated by a change in the intensity value, or by a geometric change. The first solutions is the optimal solution for γ = 10 -1 and τ = 10 -5 (i.e., when the deformation part is more penalised than the change in the intensity value) and the second one for λ = 10 -5 and τ = 10 -1 . For intermediate choices of parameters, like γ = 10 -3 and τ = 10 -3 , the optimal solution is an intermediate one when these two effects combine. We present these results in fig. 4. We emphasis again that for all choices of parameters the reconstructed images are of good quality, only its decomposition into deformation and template parts varies.

Let us next study the influence of the choice of the kernel size σ. This parameter defines the characteristic scale for the allowed deformations: if σ is small, deformations of small areas of the template are favored, while if σ is large, bigger areas will be easily displaced. Apart from the visual perception, the reconstruction is quantitatively compared using structural similarity (SSIM) [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] and PSNR. Table 1 lists SSIM and PSNR values for reconstructions obtained using different choices of kernel size. We see that the reconstruction is of good quality for all choices of kernel size σ.

Interestingly, even if the reconstruction looks good, its decomposition into template and deformation is really different so the template and deformation parts seem to balance out each other in an non-intuitive way when they are combined to form the reconstruction. As an example, for extremal choices of the kernel size (like σ ≈ 1/3 of the domain size), false structures appear in the template part of the decomposition. This is clearly an issue if the decomposition is used for interpretation. However these artefacts appear only for extremal choices of kernel size σ, it does not appear for values of the kernel size that lie between 1/32:nd and 1/6:th of the image size. We currently lack an understanding of this phenomenon. σ 0.3 0.6 1 2 3 5 10 SSIM 0.660 0.703 0.737 0.769 0.766 0.764 0.682 PSNR -7.75 -7.03 -6.57 -6.36 -6.49 -6.66 -8.98 Table 1: SSIM and PSNR values for metamorphosis based indirect registration with varying kernel size σ and fixed regularisation parameters γ = τ = 10 -5 .

Spatio-temporal reconstruction

The goal here is to recover the unknown temporal evolution of a template matched against (gated) parallel beam 2D ray transform data acquired at 10 different time points (from t = 0.1 to t = 1), so the target undergoes a temporal evolution. At each of the 10 time points, we only have limited tomographic data in the sense that i:th acquisition corresponds to sampling the parallel beam ray transform of the target at time t i using 10 angles randomly distributed in [(i -1)π/10, iπ/10] using 362 lines/angles. We added a Poisson noise to these data. Similarly to previous experiments, the reconstruction space discretised as 256 × 256 pixel grey scale images supported in Ω.

The registration of the template I 0 against the temporal series of data g i , 1 ≤ i ≤ 10 at the 10 time points t i is performed by minimising the following functional with respect to one trajectory (ν, ζ) ∈ L 2 ([0, 1], V × X):

J γ,τ (ν, ζ; g 1 , . . . , g 10 ) := γ 2 ν 2 2 + τ 2 ζ 2 2 + 10 i=1 L A W(ϕ ν 0,t i , I ν,ζ t i ) , g i where t → I ν,ζ t , is the absolutely continuous solution to      d dt I ν,ζ t (x) = ζ t, ϕ ν 0,t (x) I ν,ζ 0 (x) = I 0 (x) with ϕ ν 0,t ∈ G V as in (2).
The target, the gated tomographic data, and the three trajectories (image, deformation and template) resulting from the metamorphosis based indirect registration are shown in fig. 6. We see that metamorphosis based indirect registration can be used for spatio-temporal reconstruction even when (gated) data is highly under sampled. In particular, we can recover the evolution (both the geometric deformation and the appearance of the white disc) of the target. As a comparison, fig. 6(f) presents reconstructions obtained from filtered back projection (FBP) and total variation (TV) [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF][START_REF] Sidky | Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[END_REF][START_REF] Chen | Gradient-flow-based semi-implicit finite-element method and its convergence analysis for image reconstruction[END_REF][START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF][START_REF] Chen | A new linearized split Bregman iterative algorithm for image reconstruction in sparse-view X-ray computed tomography[END_REF]. Here, data is a concatenation of the 10 gated data sets, thereby corresponding then sampling the ray transform using 100 angles in [0, π]. Note that the temporal evolution of the target is not accounted for in these reconstructions.

Conclusions and discussion

We have introduced a metamorphosis-based framework for indirect registration and showed that this corresponds to a well-defined regularisation method. We also present several numerical examples from tomography supporting this theoretical result.

In particular, section 5.6 illustrates that this framework can be used to recover the temporal evolution of a template from time series data, even when data at each time point is very incomplete. This approach does however assume access to an initial template. In spatio-temporal reconstruction, such an initial template is unknown and it needs to be recovered as well. One approach for doing this is by an intertwined scheme that alternates between to steps (similarly to [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF]): (i) given a template, estimate its evolution that is consistent with times series of data using the metamorphosis framework for indirect registration, and (ii) estimate the initial template from times series of data given its evolution. The method in section 5.6 solves the first of the above steps, which is the more difficult one.

Another topic is the choice of hyper-parameters. Our metamorphosisbased framework for indirect registration relies on three parameters with the kernel-size σ being the most important as shown in section 5.5. This parameter has a strong influence on the way the reconstructed image trajectory decomposes into a deformation and a template part. Clearly it acts as a regularisation parameter and a natural problem is to devise a scheme for choosing it bearing in mind the size (scale) of the features undergoing deformation. Unfortunately, similarly to direct registration by LDDMM, the choice of this parameter (and more generally choice of kernel for the RKHS V ) is still an open problem [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Chen | Indirect image registration with large diffeomorphic deformations[END_REF][START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. One way is to use a multiscale approach [START_REF] Bruveris | Mixture of kernels and iterated semidirect product of diffeomorphisms groups[END_REF][START_REF] Risser | Simultaneous multi-scale registration using large deformation diffeomorphic metric mapping[END_REF][START_REF] Sommer | A multi-scale kernel bundle for LDDMM: towards sparse deformation description across space and scales[END_REF] but a general method for selecting an appropriate kernel-size remains to be determined.

A Gradient computation

This section presents the computation of the gradient of J γ,τ ( • ; g), which is useful for any first order optimisation metod for minimising the functional J γ,τ ( • ; g) in [START_REF] Bruveris | Mixture of kernels and iterated semidirect product of diffeomorphisms groups[END_REF]. The computations assume X = L 2 (Ω, R) with

I 0 ∈ X ∩ C 1 (Ω,
Then .

∂ ζ L f ν,ζ t , g (ζ)(η) = ∇ L(f ν,ζ t , g), ∂ ζ f ν,ζ t (ζ)(η) = Ω t 0 ∇ L f ν,ζ 1 , g η(τ, ϕ ν t,τ (x)
In order to compute the differential of the discrepancy term with respect to ν, we start by computing the differential of f

• ,t ) ∇ L f ν,ζ t , g • ϕ ν • ,t ∇(I 0 •ϕ ν • ,0 )( • )+ • 0 ∇(ζ(τ, • )•ϕ ν • ,τ )( • )dτ , µ L 2 ([0,1],L 2 (Ω,R d )) = - Ω K(x, • )1 • ≤t Det(dϕ ν • ,t (x)) ∇ L f ν,ζ t , g ϕ ν • ,t (x) ∇(I 0 •ϕ ν • ,0 )(x)+ • 0 ∇(ζ(τ, • )•ϕ ν • ,τ )(x)dτ , µ L 2 ([0,1],V )
.

Figure 1 :

 1 Figure 1: Reconstruction (c) by LDDMM-based indirect registration using a template (a) with a geometry that matches the target (b), but with incorrect background intensity values. Target is observed indirectly through tomographic data (d), which is 2D parallel beam Radon transform with 100 evenly distributed directions (see section 5.1 for details). The artefacts in the reconstruction are due to incorrect background intensity in template.

  (a) Template. (b) Target. (c) Data (sinogram). Image (d) t = 0. (e) t = 0.2. (f) t = 0.5. (g) t = 0.7. (h) t = 1. Deformation (i) t = 0. (j) t = 0.2. (k) t = 0.5. (l) t = 0.7. (m) t = 1. Intensity (n) t = 0. (o) t = 0.2. (p) t = 0.5. (q) t = 0.7. (r) t = 1.

Figure 2 :

 2 Figure 2: Metamorphosis based indirect-matching of template in (a) against data in (c), which represents 2D ray transform of target in (b) (100 uniformly distributed angles in [0, π]). The second row (d)-(h) shows the image trajectory t → W(ϕ ν 0,t , f t (ν, ζ)), so the final registered template is in (h). The third row (i)-(m) shows the deformation trajectory t → W(ϕ ν 0,t , I 0 ), likewise the fourth row (n)-(r) shows the intensity trajectory t → f t (ν, ζ).

  (a) Template I0. (b) Target. (c) Data (sinogram). Image (d) t = 0. (e) t = 0.2. (f) t = 0.5. (g) t = 0.7. (h) t = 1. Deformation (i) t = 0. (j) t = 0.2. (k) t = 0.5. (l) t = 0.7. (m) t = 1. Intensity (n) t = 0. (o) t = 0.2. (p) t = 0.5. (q) t = 0.7. (r) t = 1.

Figure 3 :

 3 Figure 3: Metamorphosis based indirect-matching of template in (a) against data in (c), which represents 2D ray transform of target in (b) (100 uniformly distributed angles in [0, π]). The second row (d)-(h) shows the image trajectory t → W(ϕ ν 0,t , f t (ν, ζ)), so the final registered template is in (h). The third row (i)-(m) shows the deformation trajectory t → W(ϕ ν 0,t , I 0 ), likewise the fourth row (n)-(r) shows the intensity trajectory t → f t (ν, ζ).

Figure 4 : 3 σ

 43 Figure 4: Reconstruction results and their decomposition into template part and deformation part for various regularisation parameters γ and τ .

Figure 5 :

 5 Figure 5: Reconstruction results and their decomposition into template part and deformation part for various kernel size σ. As comparison, the image domain has size 32 in this example.

( a )

 a The temporal evolution of the target. (b) The (gated) tomographic data. Each data set is highly incomplete (limited angle). (c) Image trajectory (reconstruction), combines deformation and template trajectories.(d) Deformation trajectory, models mainly geometric changes. (e) Template trajectory, models mainly intensity changes.(f) FBP (left) and TV (middle) reconstructions from concatenated data (right).

Figure 6 :

 6 Figure 6: Reconstructing the temporal evolution of a template using metamorphosis. Target (a), data (b), and results (c)-(e), are shown at selected time points t = 0.2, 0.4, 0.6, 0.8, and 1.0. As a comparison we show reconstructions assuming static target obtained from concatenating the gated tomographic data (f).27

)dτ dx = Ω 1 0 1 τ

 11 ≤t | Det(dϕ ν τ,t (x))|∇ L(f ν,ζ t , g)(ϕ ν τ,t (x))η(τ, x)dτ dx = 1 • ≤t | Det(dϕ ν • ,t )|∇ L(f ν,ζ t , g))(ϕ ν • ,t ), η L 2 ([0,1],L 2 (Ω,R))

+ s 0 0 0 ∇- Ω 1 0 1 s≤t 0 ∇

 00110 ∇ζ(τ, • )•ϕ ν s,τ (x), dϕ ν s,τ (x) µ(s, x) dτ ds dx = -Ω t Det dϕ ν s,t (x) ∇ L f ν,ζ t , g ϕ ν s,t (x) • ∇(I 0 •ϕ ν s,0 )(x), µ(s, x) + s (ζ(τ, • ) • ϕ ν s,τ )(x), µ(s, x) dτ dsdx = Det dϕ ν s,t (x) ∇ L f ν,ζ t , g ϕ ν s,t(x)∇(I 0 • ϕ ν s,0 )(x) + s (ζ(τ, • ) • ϕ ν s,τ )(x)dτ , µ(s, x) dsdx = -1 • ≤t Det(dϕ ν

  R) and (ν, ζ) ∈ L 2 ([0, 1], V × X) . Furthermore, for each t ∈ [0, 1] we also assume t → ζ(t, • ) ∈ C 1 (Ω, R).In numerical implementations, we consider digitised images and considerations of the above type are not that restrictive.Let us first compute the differential of the data discrepancy term with respect to ζ using the notation f ν,ζ

		t	:= W(ϕ ν t,0 , I ν,ζ t ) = I ν,ζ t	• ϕ ν t,0 . As noted
	in (14), we have
	f ν,ζ t (x) = (I ν,ζ t	• ϕ ν t,0 )(x) = I 0 ϕ ν t,0 (x) +	t	ζ τ, ϕ ν t,τ (x) dτ.
			0

  ν,ζ 1 with respect to ν. Hence, let µ ∈ L 2 ([0, 1], V ) and x ∈ Ω. Then

	Using (29), we can compute the derivative of → L W(ϕ ν+ µ t	, I ν+ µ,ζ t	) at
	= 0:
	d d	L W ϕ ν+ µ t	, I ν+ µ,ζ t	=0 = ∇ L f ν,ζ t , g ,	d d	f ν+ µ,ζ t	=0
		= -					t	∇ L f ν,ζ t , g (x) • ∇I 0 ϕ ν t,0 (x) , dϕ ν s,0 ϕ ν t,s (x) (µ(s, ϕ ν t,s (x)))
			Ω			0
				s
		+			∇ζ(τ, • )•ϕ ν t,τ (x), dϕ ν s,τ ϕ ν t,s (x) µ s, ϕ ν t,s (x)	dτ ds dx
			0
		= -					t	Det dϕ ν s,t (x) ∇ L f ν,ζ t , g ϕ ν s,t (x) •
			Ω			0
								∇I 0 ϕ ν s,0 (x) , dϕ ν s,0 (x) µ(s, x)
	d d	f ν+ µ,ζ t	(x)	=0	= ∇I 0 ϕ ν t,0 (x) ,	d d	ϕ ν+ µ t,0	(x) =0
		+	0	t	∇ζ(τ, ϕ ν t,τ (x)),	d d	ϕ ν+ µ t,τ	(x) =0 dτ
			t
		= -			∇I 0 (ϕ ν t,0 (x)), dϕ ν s,0 ϕ ν t,s (x) µ(s, ϕ ν t,s (x)	ds
			0
			t						τ
					∇ζ τ, ϕ ν t,τ (x) ,	dϕ ν s,τ ϕ ν t,s (x) µ s, ϕ ν t,s (x) ds dτ
			0						t
			t
		= -			∇I 0 ϕ ν t,0 (x) , dϕ ν s,0 ϕ ν t,s (x) (µ (s, ϕ ν t,s (x)	ds
			0
				t		s
		-					∇ζ(τ, • )•ϕ ν t,τ (x), dϕ ν s,τ ϕ ν t,s (x) µ s, ϕ ν t,s (x)	dτ ds.
			0			0
									(29)

https://github.com/bgris/IndirectMatchingMetamorphosis

https://github.com/bgris/IndirectMatchingMetamorphosis

Acknowledgements

The work by Ozan 'Oktem, Barbara Gris and Chong Chen was supported by the Swedish Foundation for Strategic Research grant AM13-0049.