Keywords: Travel-time budgets, duration model, in-depth interview, perceptions of travel time, societal influences

This paper combines both quantitative and qualitative methods to provide supplemental information and results that clarify the reasons that individuals devote ample time to travel.

Through a quantitative analysis of travel-time budgets (TTB) in seven cities (Zurich, Bern, Geneva, Rennes, Lyon, Strasbourg, and Grenoble), the research initially establishes the relationship between the amount of daily travel time allocated and socio-demographic characteristics of individuals and households. The application of a duration model method additionally highlights the influence of diverse types of activities at the destination on TTB.

Our analysis confirms that the overall TTB of individuals has increased over roughly the past ten years, with Swiss cities possessing higher travel times than French ones. In many cities, around 20% of our samples also devote excessive amounts of time to travel. The duration model therefore indicates an atypical individual behavior not explained by the covariates of the model. Conversely, other dimensions could certainly explain the choices or obligations that influence these high TTB.

Accordingly, in-depth interviews are given to individuals that travel during the week for longer than 2 hours a day. Their profiles provide several explanations to the reasons behind their choice (or obligation) to spend major amounts of time in daily travel. Consequently, the allocation of travel time of a person is not solely determined by the activities at a destination, but also by individual aspects. Additionally, travel time is perceived as time when other activities can be performed simultaneously (ex. working on the train, listening to music, etc.).

Introduction

Travel time is perceived as a mobility indicator, yet certain attempt to quantify while others research its qualitative elements. Nonetheless, both undertake questions on the measurement, the universality, and the perception of time. Consequently, the analysis of travel time appears to provide an opportunity to reunite multiple disciplines. This is notably one of the aspects discussed at the conference "Mobilities and Temporalities" 1 in 2004 from the workgroup "Spatial Mobility and Societal Fluidity" 2 of the AISLF. [START_REF] Montulet | Introduction[END_REF] illustrates how issues of congestion and individual/social time, including its acceleration, can be articulated through the use of travel time, showing that one issue strongly correlates with the other. He also evocates the large complexity of mobility and temporality relationships that "proposes, illustrating a primary approach, to resort to the conjugation between social time, the organization of transport modes, and temporal attitudes" 3 . Measuring travel time is indeed necessary, for example, for the articulation of monetary costs, in order to better understand mobility behavior or improve our comprehension of urban dynamics. At the same time, certain challenges and limits often appear within quantitative operations, which could be complimented by qualitative dimensions.

From an economic perspective of transport, the study of travel time constitutes a central tool for the comprehension of mobility behavior with its relationship to a trip's finality, or the activities at the destination. This provides our initial point that demonstrates a greater issue regarding the creation of daily agendas, which appears all more complex in the current context of "deregulating" social and individual time. Our mobility analysis first concentrates on travel time and individual characteristics, then progresses to related dimensions, such as activities performed at a destination, residential location, and possible reallocation of travel time.

1 Mobilités et temporalités 2 Mobilité spatiale et fluidité sociale 3 p. 13 from [START_REF] Montulet | Introduction[END_REF].

Alternatively, the sociological perspective of travel time attempts to understand the individual reasons, both intentional and inadvertent, that influence a person to assimilate a grand portion of time to mobility. This section particularly targets persons possessing high travel-time budgets (TTB) to identify their opinions on their respective travel time and the influential factors that condition their perceptions and choices. In this perspective, individuals not only organize their travel time according to the activities at a destination, but also to certain supplemental aspects. These are not related to the trip itself, but rather the individual and their past mobility history, culture, convictions, etc. Research has previously evoked how an individual formulates their unique mobility through experiences and cognitive allocation of transport [START_REF] Flamm | Operationalising the Concept of Motility: A Qualitative Study[END_REF]. Choices are therefore not made purely on the satisfaction of a specific need, but take into consideration attitudes and personal attributes (Redmond and Mokhtarian, 2001). This double approach attempts to go beyond the usual conceptions of travel time, especially from certain economists that strictly relate mobility to monetary costs. The management of the temporal resource implies complex behavior, since travel time is the center of numerous interactions between individual and social choices. This unique role of travel time integrates the individual within a space filled with opportunities in multiple urban contexts. Consequently, managing mobility progressively appears to not be simply reduced to monetary or temporal costs, independent from other dimensions.

The objective of this article is to explore mobility that implies an atypical organization of the temporal resource. Sociological and econometric analyses are therefore proposed to comprehend the influential factors that cause certain individuals to accept their particularly high TTB. The quantitative approach is based on mobility studies of seven Swiss and French agglomerations, while the qualitative section uses interviews of individuals working or residing the Leman basin that possess high TTB.

After a general presentation of mobility and urban contexts of the chosen agglomerations (section 1), the article will present the implemented sociological and econometric methodologies (section 2). Section 3 will then present the results that aim to consolidate the findings of the two methodologies. Finally, the last section (section 4) will restate the complementarities in terms of results and methods, in order to highlight potential contributions identified between the two disciplinary fields.

Urban contexts and mobility disparities

Urban contexts and differentiated mobility

The following is based on mobility studies conducted in France and Switzerland4 , with our respective study concentrating on seven agglomerations: Bern (1994 and2000), Geneva (1994 and2000), Grenoble (1992and 2001), Lyon (1985[START_REF] Allison | Survival analysis using SAS -A practical guide[END_REF], Rennes (1991 and2000), Strasbourg (1988and 1997), and Zurich (1994and 2000). Large-scale studies are regularly performed, allowing the study of two consecutive dates and possible transformations in mobility practices for each city. The data available for each city offers several advantages, such as the contemporaneity of the studied periods and the convergence of objectives and study methodologies. Consequently, the databases constitute an adapted tool to study mobility at both aggregated and disaggregated levels (see table A-1 in the appendix).

As expected, certain adjustments are necessary for our respective study on urban mobility. A methodology comparing different databases induces and comprises of several limitations and inaccuracies. The treatment and homogenization of information necessitates the choice of certain constraints, yet the information used here principally remains at a comparable level of traditional mobility analyses.

The observed mobility analyses of various urban contexts permit an identification of certain conditions that differentiate mobility forms. The choice of the cities was performed while taking into account availability constraints and the following criteria:

• The form and dispersion of the suburbanization process

• The level of development for the network of heavy public transport modes

• The urban and transport policies that target the organization of city center accessibility

The transport systems and policies of the Swiss agglomerations are principally oriented towards public transport and walking. Contrarily, in the French cities, transport policies, such as automobile restriction and parking limitations in city centers, have only recently appeared compared to Swiss city policies.

Table 1 illustrates the means of TTB and the number of trips per individual in each agglomeration. The study of the means reveals that the TTB per individual in French cities, around 70 minutes, is less than the estimated 80 minutes in Swiss ones.

However, the analyses of the TTB distribution for each city (table A-2 in the appendix) illustrate the asymmetric distribution of TTB, resulting in an overestimation of the TTB mean for all cities. This non-normal distribution will be questioned and illustrated in the duration model section.

In cases of asymmetric distributions, the median is a relevant indicator, regardless of its responsiveness to rounding. The median TTB is 60 minutes for French cities and 65 minutes for Swiss ones. Since the median is less sensitive to extreme values of the distribution, the separation of TTB between the two countries is reduced. The distribution quartiles indicate that the large separation between the average Swiss and French TTB is partially due to the overrepresentation of the population characterized by particularly high TTB in Swiss cities. Moreover, an observable effect of agglomeration size is apparent by comparing the TTB levels of Rennes with Zurich or Lyon. The relative proximity of TTB medians strongly reduces the separation of TTB between these agglomerations. With the exception of Rennes and its lower median at 50-55 minutes, the majority of individuals in all the agglomerations do not have a TTB that exceeds 60-70 minutes. Finally, the median TTB difference between the two dates for each city is around 5 minutes. Accordingly, the TTB is characterized between the two dates by an increase in both Swiss and French cities. For each city, the increase in mean TTB is notably the result of an increase in persons with particularly high travel times of over 100 minutes. Figure 1 shows the increase of the number of people with over 100 minutes of daily travel time between the two observable dates.

Figure 1 Source: Enquêtes-Ménages-Déplacements and the Microrecensement

As the results of the duration models show, these individuals with extremely high daily travel times, or "intense commuters", do not necessarily distinguish themselves from the rest of the population by classic socio-economic characteristics or their mobility attributes.

Therefore, the objective of implementing two methodologies is to first clarify this phenomenon with the use of statistical data, then the use of sociological means, in order to illustrate the attributes of these intense commuters. The application of an econometric model on all the chosen cities allows the proposal of conclusions relative to intense commuters and their atypical organization of the temporal resource allocated to transport. Conversely, indepth interviews propose unique explanations regarding the behavior of intense commuters. They not only validate sociologically certain econometric results, but also exemplify certain limits of the quantitative approach.

Methodologies 5

Two different methodologies were implemented to understand the behavior of intense commuters. The disaggregated analysis of displacement durations integrates the econometric tool of duration models, while the sociological aspect utilizes in-depth interviews.

Duration model applied to travel-time budgets

The duration model aims to produce an analysis of daily transport durations. Models such as these are frequently used to analyze activity durations, but their application strictly refers to non-transport activities. Two types of results are expected from its implementation:

1.

The demonstration of influential variables through the use of personal and household socio-economic characteristics, residential location, and certain attributes on activity and mobility programs. This should clearly illustrate the relationships between TTB and certain variables.

2.

The identification and description of the temporal dynamic that partially dictates the organizational processes of daily travel time. This temporal dynamic is analyzed through the variation of sample size as the TTB gradually increases. As in the analysis of activity durations, a temporal dynamic specific to daily travel time can be identified. Its description will indicate the rhythm and possible cycles that characterize daily travel time through the course of the day.

In-depth interviews of intense commuters

The qualitative research targets individuals with high TTB, in order to understand the perceptions and experiences of travel time. A general TTB of 100 minutes was set to establish persons that were considered intense commuters. Certain additional criteria were imposed in the selection process:

The individual works or lives in the Leman region.

The individual works or goes to school at least three times during the week.

The sample respects a heterogeneity that reflects the regional population. This does not only consider socio-demographic aspects, but also transport modes used. The names have been changed to keep the interviewees anonymous. 7 We used the Michelin web site www.michelin.fr to calculate the distance by car between the two referred communes. The numbers themselves are not necessarily reliable for statistical calculations (ex. speed), considering that the exact address was never requested. However, the kilometers provide a way to compare relatively distance between interviewees.

Comparative results

Two types of results are shown through the use of the two different methodologies:

An atypical organization of travel times

The first result from the duration models of intense commuters is made by the nonparametric estimation of TTB distributions. The form of the estimated hazard curve is nonmonotonic and similar for each city8 . The conditional probability of an interruption in travel time allocation given the time passed first increases until 120 minutes, then decreases after 120 minutes.

The interpretation of this non-monotonic property of the estimated hazard curve becomes relevant when considering the estimated median residual survival lifetime(figure 2), or the travel time remaining or residual according to the time already devoted to mobility. The non-monotonic, estimated hazard derives from a residual travel time that first decreases, then increases. Therefore, according to the estimated hazard for Lyon 1995, during the process of travel time allocation, an individual that has already passed 30 minutes of travel time has "on average" a remaining travel time of around 45 minutes. Accordingly, a person with 60 minutes of travel time generally still has 30 minutes remaining.

In an intuitive sense, an individual has a remaining travel time that decreases from the beginning of the temporal process, which appears coherent with the hypothesis of the behavioral reduction of travel time (due to the lassitude or the saturation of preferences). However, once greater than 100 minutes, the hazard trend is opposite and the remaining time increases. In other words, individuals tend to devote greater supplemental time to transport after 100 minutes. The resistance to allocate greater time for travel appears to reduce after 100 minutes. Accordingly, the estimated hazard and residual survival median level suggests a transition in travel time allocation behavior for a population group with a TTB of around 100 minutes. The non-monotonic nature could be the result of the non-integration of an important explicative dimension. Conversely, the parametric analysis of the model, presented in the following section, was constructed to introduce individual and household characteristics while testing the different hazard forms. The performed tests that integrated individual attributes validated the non-monotonic nature of the hazard (table 3). Therefore, this only could derive from an unobserved heterogeneity or non-described elements in the data. This hazard form contradicts the simplistic economic vision that the reduction of travel time, with monetary costs, would be an optimization condition of mobility. In contrast, the econometric result illustrates a less direct mechanism. Other dimensions must be integrated, which does not represent time allocation behavior by means of travel time reduction. The change of the estimated hazard form indicates the existence of two phases in the allocation process of time to travel. Intuitively, the increasing conditional probability to interrupt the process refers to the behavioral reduction of the temporal cost of mobility. The change in the tendency would indicate an irregularity in the mechanism, almost as though certain individuals cannot, or choose not to reduce their travel times after a certain level. Those with TTB above 100 minutes may be faced with situations and constraints that do not allow them to easily reduce their travel time, or quite simply their behavior cannot be represented by a simplification of temporal use.

The estimated non-monotonic hazard refers to the economic notion of generalized cost. This illustrates that the economic objective of an individual is to not reduce each cost, but rather reduce their integral sum. This concept integrates the possibility to reduce one cost to the detriment of another.

Personal factors could partially explain this phenomenon in the mobility behavior of intense commuters. Individuals construct their unique mobility plan with the factors related to the trip (ex. financial costs, speed, duration) and complementary specific elements, such as environmental protection (***), personal space, past history, etc. (Flamm, 2003). Accordingly, in the analysis of intense commuters, subjective elements that are not associated with the trip must be taken into account.

The analysis of the in-depth interviews suggests that other dimensions, in addition to monetary or temporal constraints, could influence the allocation of time to travel. Generally, each individual trades of between travel time and a combination of additional elements (residential location, modal choice, etc.). Intense commuters distinguish themselves in their time allocation behavior, due to the lack of optimizing the temporal resource.

However, the difficulty is to test explicative dimensions of these atypical time allocation behaviors. The second series of results compares quantitative and qualitative results, in order to deduce similarities or divergences of interpretation.

The influence of socio-economic characteristics, mobility indicators, and activity program attributes

The second result of the duration models is an estimation of the effects of explicative variables from the three following dimensions: individual and household socio-economic attributes, mobility attributes, and activity programs. Due to space constraints, only the estimated models with the national samples are presented. The customary distribution forms have been tested9 . The distributions that produce the best adjustment are the loglogistic and generalized-gamma distributions (table 3). The loglogistic distribution is presented here for an easier interpretation. 4 indicates the effects of explicative variables on TTB. The estimated coefficients are interpreted after an exponential transformation, such as the travel time quotient of individuals characterized by the variable and the travel time of other individuals. For example, the associated coefficient of the professional activity variable is 0.114. The variable is significant at a 99% confidence level, and has a positive effect on TTB. More precisely, a worker has a mean TTB that is estimated at 12% greater than a person unemployed (e 0,114 = 1,12). From the econometric model, the different influential elements in travel time allocation are discussed simultaneously with the results from the sociological analysis. These different elements in time allocation to transport, from the econometric model, are discussed at once with the results from the sociological analysis.

Weekly regularity

By retaining the hypothesis that TTB are partially influenced by the activities at a destination, the multiplicity of these activities, according to the day during the week, would increase or decrease the travel time. Essentially, the two methodologies show the differences of TTB in relation to the day of the week, which forms certain regularity in the long term and creates a weekly cycle.

During the week, many of the interviewees have an increase in their TTB once or twice a week, due to activities outside the workplace. Frequently, these activities are premeditated and performed on precise days. For example, John increases his TTB by 20 minutes two times per week for soccer practice. After 50 minutes of travel by train, Thomas drives 40-50 minutes once a week for band rehearsal. Anne also adds 10 minutes of travel to drop off and pick up her laundry two times a week.

In addition, the notion of opportunity travel cost is perceived within the interviews. The intense commuters often see their "free time" for other activities clearly reduced. The constraint of "increased scarcity" with time consequently modifies certain mobility aspects. The first effect relates to the difficulty of returning back home when supplemental and occasional activities are planned. These supplemental trips are therefore frequently linked with regular ones to gain time for activities. This also influences the choice of transport modes. For example, instead of habitually taking the train to go to work, Rachel changes transport modes and travels by car when she has other activities planned at Lausanne after work, often for plays and musicals.

Another observation is shown through the influence of planning one's mobility for other days of the week. John and Thomas, for example, must leave work early to attend their soccer and band meetings. They then must recuperate the hours lost by working more on the other days of the week. In this case, secondary activities can even influence TTB for the entire week. These daily variations are observed in numerous quantitative studies10 . The daily cycle of a TTB model commonly shows an increase from Monday to Friday. In general, the analyses of TTB in the seven agglomerations also produce this same daily cycle.

Urban and individual contexts

Two types of contextual dimensions have been identified that can have an influence on mobility practices and TTB. First, the urban context constitutes a framework where daily activities are performed. Individual contexts then play a role in mobility behavior, via the combination of personal, household, and professional spheres.

Generally, the effects of urban spaces produce ambiguous results among the studied cities. Size and density factors are sometimes contradictory. Joly ( 2005) presents a rapid overview of literature of these effects of the urban context. Most often, an inconsistency between two urban systems appears between spread-out and dispersed cities, depending on transport systems that favor automobile use, and concentrated ones that possess a greater equilibrium among transport modes. These correspond to two principal profiles, with French cities commonly characterized by higher TTB in the city center, while lower TTB levels in Swiss cities.

The qualitative interviews clearly confirm this behavioral change, consistent with urban contexts and lifestyle, with both dimensions being influential. Compared with the three interviewees from the city center, individuals residing in the country are noticeably more stable and oriented towards family life. They also often live with several other people. According to their testimony, they sacrifice a part of their daily time for the advantages of the countryside, not only for themselves, but also for the other people in the household (ex. green and secure areas for children). Trade of exists between several dimensions, such as mobility costs (temporal and monetary), residential location and land price, access to property ownership [START_REF] Orfeuil | Les dépenses pour le logement et pour les transports des ménages franciliens[END_REF], lifestyle changes (Cao and Mokhtarian, 2005), and even personal situations of each member of the household [START_REF] Flamm | Comprendre le choix modal. Les déterminants des pratiques modales et des représentations individuelles des moyens de transport[END_REF].

Contrarily, the three interviewees of the city center prefer the lifestyle of those in the countryside, yet subsidiary factors restrain them to downtown areas. For example, George desires to live in the countryside, but he does not have enough money to buy a car. John and Cindy believe that their lives are too unstable for the moment, since their professional position is temporary and precarious. In the future, these individuals could presumably look for a certain lifestyle (ex. life in the countryside, car, family, etc.) after stabilization in their professional careers, without necessarily attempting to reduce their travel time.

The results presented in table 4 partially illustrate these lifestyle effects, but they are generally very difficult to quantify. Age, number of children, household structure, or zone of residential location are indicators used to assess lifestyle cycles and modes. Frequently, these attributes of travel time are different according to the level of responsibility of the members of the household, professional status, or the type of residence.

However, cultural influences are difficult to measure quantitatively, yet emerge in the interviews. Two interviewees noted that their decision to stay at their current residential location was based on family reasons, resulting in a weak accessibility to public transport and a high TTB compared to the distance covered. Both are 30-year old women, dependant of their parents, and are strongly engaged in daily family affairs. When asked about the reasons for their acceptance of a high TTB, they believe that their choice is greatly influenced by their African culture that promotes daily interaction with family members. Consequently, traditional indicators to indicate lifestyle, such as age or professional position, have a difficult time to grasp the mobility behavior of individuals strongly influenced by cultural or familial aspects at greater scales.

The allocation of travel time

Travel time is not solely used to physically move from an origin to a destination. Activities can also be performed during the trip. Several questions were asked in the interviews to first identify the activities during these trips, then understand the perceptions and experiences of the trip through the use of these activities.

Generally, public transport users have a greater variety of activities than automobilists. All interviewees using public transport read, work, sleep, socialize with others, or listen to music. Contrarily, automobilists are limited to listening to music or talking on the phone. In most cases, the interviewees believe that they partially take advantage of the travel time by simultaneously performing other activities.

However, according to all those interviewed, travel time is nevertheless considered to be time "lost". To a degree, activities help to compensate for the lack of free time, but all interviewees recognize that they could do more activities and be more efficient if they were not traveling at the same time. For example, certain interviewees have problems working in the train, due to crowds or the short travel time. Others are too tired for their return trip and are only productive in the morning. For automobilists interviewed, the concentration needed for highway travel prohibits them from completely taking advantage of the travel time.

Frequently, this incomplete compensation for travel time is nonetheless valorized in a different manner when an individual uses the travel time to escape from pressure or stressful situations. For instance, Caroline listens to the news or music on the radio when she drives and the activities only partially aid in compensating for lost time. However, instead of activities, her passion for driving makes the commuter trips pleasurable. Robert, on the contrary, takes advantage of his travel time differently. He listens to CDs disliked by his family in his car. The CDs improve his musical knowledge, which creates a personal space outside social obligations. Lastly, Thomas takes advantage of his time in the train to be away from his children. According to the interviewee, his children are young and require a lot of attention. The travel time on the train allows him to relax without the worries of familial responsibilities. In all these examples, the trip not only represents physical travel, but also a change in mental state.

For the interviewees, the sentiment of lost time in travel is correlated to the execution of other activities during the trip and the quality of these activities. On the other hand, in the case of intense commuters, these activities do not appear to completely compensate for the lack of free time. This mechanism of "reallocating" travel time provides a supplemental dimension to the trip beyond simple physical movement. To understand the choices behind transport, notably the choice in transport modes, these behavioral notions must be taken into consideration. In this perspective, the trip cannot be perceived as a simple function of monetary and temporal cost. In effect, an individual can take advantage of travel time to do other activities, which is clearly an issue for intense commuters.

To restate the hypothesis of [START_REF] Mokhtarian | How derived is the demand for travel? Some conceptual and measurement considerations[END_REF], the utility of transport can derive from three sources: the utility at the destination, the intrinsic utility of the trip, and the utility of the activities performed during the trip. The role of transport is not simply attributed to an activity. The conception of travel time is not perceived merely as a cost (or a component of a generalized cost), but as interstitial time considered as an activity in itself, a period of time devoted to execute another activity, and naturally the required time to reach a destination. Transport is therefore an activity with an intermediary role between an obligated and discretionary activity.

Personal convictions

Individual convictions emerge as an important role in transport behavior, notably in the choice of transport mode. This is clearly illustrated in two individuals that are multimodal. Rachel and Thomas first drive 10 to 15 minutes to a train station, then continue the rest of the trip by train.

According to the interviewees, the bad public transport accessibility in their villages forces them to do part of the trip by car. Both enjoy the tranquility of the countryside and live in houses outside the city. However, their past transport experiences were exclusively with public transport to commute to work. Neither owned a car before buying their house, because their prior apartments were always located in areas with dense public transport networks. They both consider the car "dangerous" and "bad for the environment", while public transport represents a "calm" and "tranquil" atmosphere.

The part of the trip by car is perceived as an obligation for the two interviewees. They voluntarily sacrifice more of their spare time in their trip with the use of public transport. Their reasons, however, are not uniquely related to a better quality of life in the countryside, but also other personal convictions and preferences. Therefore, in these instances, certain individual inclinations influence a person to use public transport, regardless of a possible increase in temporal or monetary costs compared to the automobile.

Conclusion

Beyond the results or trends proposed by the interviews, the mobilization of the two methods illustrates the complementarities of both approaches and perceptions. Frequently, the weakness of one method is compensated by the dynamics of the other. The will to modulate and synthesize, with the supplemental contribution to explain in detail, are perfectly harmonizing and necessary for one another.

Our study attempts to demonstrate a mixed method, by mobilizing two tools that address the same subject. This work shows certain complementarities for mixed methods [START_REF] Tucci | Mixed Methods and simulation research designs[END_REF]:

1. The approach by using qualitative interviews aids to identify unobserved heterogeneities within data and proposes other variables to explain them.

2. The results of a mixed approach can contribute by explaining previously misunderstood statistical results.

3. Qualitative research can help discover problems in the definitions and measures of quantitative instruments.

4. The quantitative approach can be used to examine the extent of results from a qualitative study and provide support for their generalizations.

The disaggregated model identifies a certain number of influential factors in the durations of daily travel, such as individual and household socio-demographic characteristics, mobility attributes, and residential locations. Conjointly, duration models estimation illustrates the particularity found at an individual level with regards to travel time management. Indeed, a significant part of our samples (around 20%) are characterized by high TTB, which paradoxically counters typical tendencies of increasing TTB. Despite the integration of certain explicative variables, the qualitative results compliment the quantitative study. The interviews provide a reminder that multiple dimensions interact and/or determine mobility behavior. Together, the two methods show that this specific population with high TTB may not wish, or have a harder capacity to reduce their travel time, proving that other dimensions must be integrated in interpreting TTB besides temporal and monetary costs. Therefore, by using a mixed method, this research on TTB demonstrates that one methodology is not necessarily independent from another. The mutual implementation of two methods permitted to surpass their respective limits by questioning and developing them together.

Particularities of duration data

The character of duration data12 imposes a certain number of restrictions on the usable model types. The major problem of applying classic models to duration data is the normality hypothesis of duration distribution, which is rarely observable. The linear regression is therefore not vigorous when contradicting the normality hypothesis. In addition, the statistical inference properties of the estimators are not valid without the normality hypothesis of the residuals.

11 For supplemental notions of other estimation methods, refer to Kalbfleich and Prentice (1980), [START_REF] Allison | Survival analysis using SAS -A practical guide[END_REF], [START_REF] Hosmer | Applied survival analysis[END_REF]Lemeshow (1999), or Lawless (2003).

Modeling the conditional probability of interruption

Classic methods state beforehand the probability of an event. However, they treat equally events occurring before or after an observation period. Duration models concentrate on modeling the conditional probability of interruption of a respective event. They therefore consider more complete information than classic models and provide unique information on the studied duration process. For example, linear or logistic regressions model the nonconditional probability of an event (the probability of an event, independently of the presented date of the event), while the duration model concentrates on its conditional probability.

Consequently, duration models introduce the temporal dynamic notion of a studied process. The probability of interruption can vary throughout the duration of the process. The analysis of the conditional probability of interruption is also explicit in providing information on the individual process of time management.

Characterizing the distribution of duration and the temporal dynamic

The objective is to model and study daily duration dedicated to transport and the probability of interrupting this process during its execution. Classic models are founded on the probability of execution, an estimated event, obtaining a higher or lower TTB, and its survival attributes:

[ ] ∫ ∞ = - = > = t su u f t F t T t S ) ( ) ( 1 Pr ) (
This corresponds to the probability that a process continues until date t.

The duration models estimate the probabilities of interruption at date t, knowing that the process continued until the respective date t. This conditional probability is noted h(t) and qualifies as a hazard function, defined as:

( ) Δ > Δ + < ≤ = + → Δ t T t T t P t h 0 lim ) (

Interpretation of the hazard and survival

The slope of the hazard function indicates variations of the slope and survival function. The variations of the hazard therefore correspond to modifications of the rhythm of the survival, the acceleration, or the deceleration of the process. The slope of the hazard curve h directly illustrates the variations of the descending rhythm of the survival function. The obtained log-logistic hazard shows a reversal of the slope. Consequently, this conditional probability increases over time, then decreases. The interruption is increasingly probable at the beginning of the process, then decreasingly probable for higher durations. The form of the hazard function therefore has important implications on the representation of the temporal dynamic of the process [START_REF] Hensher | Hazard-based duration models and their application to transport analysis[END_REF].

The estimated median residual survival lifetime

A unique indicator of duration analysis is the median, preferred over the mean because of the reduced bias formed through possible dissymmetric distribution. The non-parametric estimation produces for the population at risk for each date, the supplemental survival time reached by 50% of individuals. This is the median residual survival lifetime at date t,

showing the median remaining time for the population at risk for the date.

Qualitative interviews

We selected eleven people without using announcements to the general public (ex. local advertising). By asking various people known by the team, we successfully identified individuals that had high TTB. Listing and comparing certain socio-demographic criteria for each interviewee confirmed the diversity of the sample. Certain individuals also had unique attributes that enriched the results: Paul works in Geneva, but lives in France.

Sarah and Wanda spend a lot of time using public transport, without necessarily traveling long distances.

Cindy regularly uses a bicycle for her trips during the week.

Rachel and Thomas use both the car and public transport for the same trip to get to work.

Caroline carpools with a colleague 29

The interviews were either performed in person or via telephone. The interviews generally lasted between 45 and 60 minutes, with the questions focused on four phases: childhood, adolescence, after moving out of the parent's home, and the present. Several questions were asked on the perceptions and experiences of these trips, residential location, aspirations, other members of the household and their influences on the interviewee's mobility, the activities performed during the trip, and the regularity of these trips on a weekly or monthly basis.

The criteria of the city or countryside in residential location, positive or negative effects, and other subjective data all reflect the views of the interviewee.

Figure 2 :

 2 Figure 2: median residual survival lifetime in minutes

Table 1 :

 1 Mean and median of travel time budgets and number of daily trips

	French cities	TTB		N
	City	Mean	Median	Mean	Median
	Grenoble 1992	67.20	55	4.26	4
	Grenoble 2001	78.92	65	4.72	4
	Lyon 1985	67.79	60	4.04	4
	Lyon 1995	78.44	65	4.26	4
	Rennes 1991	57.44	50	4.04	4
	Rennes 2000	70.91	55	4.25	4
	Strasbourg 1988	69.94	60	4.67	4
	Strasbourg 1997	78.90	65	4.96	4
	Swiss cities	TTB		N
	City	Mean	Median	Mean	Median
	Bern 1994	74.13	60	3.78	4
	Bern 2000	85.41	69	4.05	4
	Geneva 1994	83.41	60	3.96	4
	Geneva 2000	84.06	65	4.22	4
	Zurich 1994	82.11	65	3.75	4
	Zurich 2000	87.35	70	3.96	4
	Source: Enquêtes-Ménages-Déplacements and the Microrecensement

  : Percentage of persons with TTB over 100 minutes

	40%	Date 1						38%
		Date 2				36%	35%
	35%
	30%						29%	30%
						27%
	25%	23%	22%		23%
	20%
			17%	17%	17%
		16%
	15%
				11%
	10%
	5%
	0%
		Grenoble	Lyon	Rennes	Strasbourg	Bern	Geneva	Zurich

Table 3 :

 3 Presentation of interviewees

	Name 6	Age Profession	Principal transport modes	Zone and type of residence	Household situation	TTB	Distance km) 7 home-work (in	Specific attributes
	John	30	Doctoral assistant	Train, tram	Downtown-Apartment	In couple	150-170	75
	Paul	38	Employee in a rental car agency	Car	Countryside-Apartment	In couple	180-195	55	French resident working in Switzerland
	Sarah	26	Hair dresser	Bus	Countryside-Residence	With brother parents and	120-130	13	African origin, very slow speed transport
	Cindy	26	Administrative assistant	Train, bus, bicycle	Downtown-Apartment	In couple	150-160	65	Bicycle use
	Robert	46	Doctor	Car	Countryside-House	Married with three kids	135-195	95	Long distance
	Anne	34 Ex pharmacist	Car	Countryside-Villa	In couple	100-125	65	Ancient commuter
	Rachel	50	Secretary	Train, metro, car	Countryside-Villa	Married	100-160	40	Multi-modal
									African
	Wanda	29	Doctoral assistant	Train, bus	Countryside-Apartment	With mother	210-270	65	origin, very slow transport
									speed
	George	29	Designer	Train, bus	Downtown-Apartment	Single	220-265	105
	Thomas	36	Editor	Train, car	Countryside-House	Married with two kids	160-215	75	Multi-modal
	Caroline	35	Employee in a cafeteria	Car	Countryside-House	In couple	85-115	65	Carpooling
	6
					12

Table 3 :

 3 Likelihood and likelihood ratio tests of a priori distributions

	Distribution	Log likelihood
		Swiss Sample	French Sample
	Weibull	-9587,08	-54824,51
	Log Normal	-9310,61	-51999,94
	Generalized Gamma	-9269,36	-51858,53
	LogLogistic	-9164,85	-51648,07
	Model 1 vs Model 2	Test Statistic
		Swiss Sample	French Sample
	Weibull vs GG	635,44 ***	5931,96 ***
	Log N vs GG	82,5 ***	282,82 ***
	Table

Table 4 :

 4 Parametric log-logistic models for the two national samples

	French sample	Swiss sample

Table A -

 A 3: Indicators of TTB distributions

				BTT
	Ville	Mean	CV	Q1	Median	Q3
	Grenoble 1992	67.20	78.32 35.00	55.00	85.00
	Grenoble 2001	78.92	81.43 40.00	65.00 100.00
	Lyon 1985	67.79	73.41 35.00	60.00	90.00
	Lyon 1995	78.44	76.44 40.00	65.00 100.00
	Rennes 1991	57.44	72.27 30.00	50.00	75.00
	Rennes 2000	70.91	97.00 35.00	55.00	85.00
	Strasbourg 1988	69.94	77.80 35.00	60.00	90.00
	Strasbourg 1997	78.90	79.27 40.00	65.00 100.00
	Bern 1994	74.13	90.27 31.00	60.00	94.00
	Bern 2000	85.41	88.21 39.00	69.00 110.00
	Geneva 1994	83.41 102.07 39.00	60.00 100.00
	Geneva 2000	84.06	90.82 38.00	65.00 108.00
	Zurich 1994	82.11	99.72 34.00	65.00 102.00
	Zurich 2000	87.35	84.80 40.00	70.00 114.00
	7.1 Duration models 11

Specifically, the studies are the Enquêtes-Ménages-Déplacements in France and the Microrecensement in Switzerland.

Refer to the appendix for more precise information on the two methodologies.

[START_REF] Littlejohn | La croissance des budget-temps de transport en question : Nouvelles approches[END_REF] show this result of the hazard curve, confirmed by the

95% confidence interval.

The Exponential, Weibull, Lognormal, Loglogistic, and Generalized-Gamma have been tested. Only Exponential, Weibull, Lognormal, and Generalized-Gamma were formally testable, because of their nested property (propriété d'emboitement).

Van der Hoorn (1979), Zahavi and Talvitie (1980),Prendregast and Williams (1981),Kumar and Levinson (1995)

The first particularity of the duration data shows that they are positive variables. In addition, despite not being the focus of this study, they could be reproached and characterized by explicative variables that change over time.

Acknowledgements

This article is based on a study conducted through a PREDIT study on "Questioning the Increase of Travel Time Budgets: New Approaches".

Appendix

Who?

All persons older than 5 years old in the residence Referred individuals (6 years +) from households (2 pers. for households > 4 members over 6 years hold)

Which trips?

All trips performed the day before the study (a weekday) All trips performed the day before or 2 days before the study (weekdays and weekends)

When?

One reference day over several months out of the year (October to May) One reference day over the entire year

Where?

The perimeter of the study best represents the agglomeration (defined by principal client)

The entire country can be considered an origin or a destination of a trip


M Guerquin-Kern

L Lejeune

K P Pruessmann

M Unser

Realistic Analytical Phantoms for Parallel Magnetic Resonance Imaging

Keywords: MRI, Fourier analytical simulation, Shepp-Logan, inverse crime

The quantitative validation of reconstruction algorithms requires reliable data. Rasterized simulations are popular but they are tainted by an aliasing component that impacts the assessment of the performance of reconstruction. We introduce analytical simulation tools that are suited to parallel magnetic resonance imaging and allow one to build realistic phantoms. The proposed phantoms are composed of ellipses and regions with piecewise-polynomial boundaries, including spline contours, Bézier contours, and polygons. In addition, they take the channel sensitivity into account, for which we investigate two possible models. Our analytical formulations provide well-defined data in both the spatial and k-space domains. Our main contribution is the closed-form determination of the Fourier transforms that are involved. Experiments validate the proposed implementation. In a typical parallel MRI reconstruction experiment, we quantify the bias in the overly optimistic results obtained with rasterized simulations-the inverse-crime situation. We provide a package that implements the different simulations and provide tools to guide the design of realistic phantoms.

I. INTRODUCTION

An active area of research in magnetic resonance imaging (MRI) is the development of reconstruction algorithms. In particular, the inverse-problem approach is getting popular [START_REF] Ribés | Linear inverse problems in imaging[END_REF], where one relies on an accurate model of the measurement process and possibly on additional information about the object being imaged.

In general, the development of any reconstruction approach requires that it be evaluated and compared to others. There are several reasons to rely on simulations in a first step

• saving the costs inherent to getting real scanner data • testing the suitability of the implemented discrete forward model • quantitatively evaluating the performance of the reconstruction software • providing reliable ground-truth data to compare with. However, for the results to be meaningful, simulations must be accomplished carefully. For instance, the inverse-crime situation, where exactly the same discrete model is used for simulation and reconstruction, leads to artificially good results. In the context of MRI, many developers of algorithms base their simulations on rasterized images. One should just be aware This work was supported by the Swiss National Competence Center in Biomedical Imaging (NCCBI). A preliminary version of this work was presented in [START_REF] Guerquin-Kern | Analytical form of Shepp-Logan phantom for parallel MRI[END_REF]. that such testing does not account for the full continuousdomain reality, because it neglects the aliasing that is inherent to spatial discretization. More realistic simulations are required to remove this bias and to ensure that the methods will perform adequately in practice.

A method to obtain resolution-independent simulations is to formulate the simulation analytically in the continuous domain. This approach goes back to Shepp and Logan [START_REF] Shepp | The Fourier reconstruction of a head section[END_REF], who introduced an ellipse-based phantom (SL) for X-ray tomography. For MRI, several analytical phantoms have been proposed. The first works, based on the SL phantom, are by Smith et al. [START_REF] Smith | Alternatives to the use of the DFT in MRI and spectroscopic reconstructions[END_REF], followed by Van de Walle et al. [START_REF] Van De Walle | Reconstruction of MR images from data acquired on a general nonregular grid by pseudoinverse calculation[END_REF]. More recently, Koay et al. [START_REF] Koay | Three-dimensional analytical magnetic resonance imaging phantom in the Fourier domain[END_REF] worked out the MR contribution of an ellipsoid for the 3-D extension of the SL phantom. Gach et al. [START_REF] Gach | 2D & 3D Shepp-Logan phantom standards for MRI[END_REF] adapted these elliptical phantoms specifically for MRI, introducing realistic physical parameters as well as T 1 and T 2 relaxation times. The family of analytical phantoms is extended by two recent works by Greengard and Stucchio [START_REF] Greengard | Spectral edge detection in two dimensions using wavefronts[END_REF] that use Gaussian functions, and Ngo et al. [START_REF] Ngo | Three dimensional digital polyhedral phantom framework with analytical fourier transform and application in cardiac imaging[END_REF] that introduce 3-D polyhedra.

The attractiveness of currently known analytical phantoms is limited for two reasons. First, the vast majority of currently available phantoms (except [START_REF] Ngo | Three dimensional digital polyhedral phantom framework with analytical fourier transform and application in cardiac imaging[END_REF]) use ellipses as basic elements. While such simple shapes have the advantage of mathematical tractability, they do not lend themselves well to the generation of images with realistic anatomical features. Secondly, to the best of our knowledge, no analytical phantom has been proposed that would take into account MRI receiving-coil sensitivities in the context of the simulation of parallel MRI experiments [START_REF] Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF].

In this work, we extend the class of available analytical phantoms by introducing regions parameterized by spline contours which are general enough to reproduce polygons and Bézier contours. Our shapes are well suited for the description of realistic anatomical regions [START_REF] Segars | A realistic spline-based dynamic heart phantom[END_REF]. To accurately simulate image formation in parallel MRI, we also make use of analytical models for the coil sensitivity maps. Specifically, we investigate the use of two classes of basis functions-polynomials [START_REF] Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF] and complex sinusoids-which both have the ability to generate maps that are physically realistic. These parametric forms are used to derive closed-form solutions for the MRI coil data. We have implemented and tested both models. Our conclusion is that the new sinusoidal one is preferable because it is better conditioned and robust to roundoff errors, while offering of the modeling flexibility that is required. The polynomial model works well with ellipses and polygons, but tends to display

∈ R 2 spatial coordinates (XY plane) k ∈ R 2 k-space coordinates (XY plane) ω ∈ R 2 Fourier angular frequency ρ(r) ∈ R + object (proton density) in space m S (k) ∈ C k-space observation from channel S S(r) ∈ C spatial sensitivity of the channel f (ω) ∈ C function f in the Fourier domain 1 R (r) ∈ {0, 1} characteristic function of a region R ∂R contour of a region R Jn ∈ R R n-
i = x i δ k,l ∈ {0, 1} Kronecker's delta (1 if k = l and 0 otherwise) j ∈ C imaginary unit such that j 2 = -1 Multi-Index Notations for α ∈ N d z α = z α i i ∈ R |α| = α i ∈ N p! = p i ! ∈ N p q = p i q i = p! (p-q)!q! ∈ N b p=a = b 1 p 1 =a 1 b 2 p 2 =a 2 . . .
numerical instabilities with Bézier contours when the order becomes too high. This document is organized as follows: in Section II, we present the different models considered for the parallel MRI measurement process, the analytical phantom, and the coil sensitivities. We motivate and compare the polynomial and the proposed sinusoidal models. In Section III, we propose the main theoretical elements that make the analytical MRI simulation possible, deferring the more technical considerations until Appendices A, B, and C. Finally, we present in Section IV the experiments that validate our implementation of the theoretical tools and an application that quantifies the bias of rasterized simulations on linear and nonlinear reconstructions, in a typical parallel MRI setup.

II. MODELING

In this section, we present the MRI measurement model and building blocks that are used to define our phantom. The main notations adopted are summarized in Table I.

A. Parallel MRI

We use the well-established linear model for parallel MRI that relates the object ρ to the k-space signal m Sn observed by each receiving coil, via the Fourier integral m Sn (k) = S n (r)ρ(r)e j2πk.r dr,

where S n accounts for the sensitivity map of the n-th receiving channel. The sensitivity map is defined as S n (r) = B x (r) -jB y (r), where, by the principle of reciprocity, B(r) = (B x (r), B y (r), B z (r)) corresponds to the magnetic field generated at point r by a unit-value steady current in the coil. The Biot-Savart law relates this field to the coil geometry through the relation

B(r) ∝ coil du × (u -r) u -r 3 . (2)

B. Analytical Phantom

We mathematically define the phantom ρ as a simple function, involving R regions R i of constant intensity ρ i

ρ(r) = R i=1 ρ i 1 Ri (r). (3)
The term region refers to a connected and bounded set. The symbol 1 R denotes the characteristic function of a region R.

Such a phantom has a limited spatial support ( R i=1 R i ) that we call a region of interest (ROI).

This model allows us to render realistic phantoms of two kinds

• piecewise-constant phantoms that mimic segmented data with sharp contours (e.g., the SL brain phantom) • textured phantom via a triangular-mesh approach. We investigate the first approach in this paper. The contours that are considered are ellipses, polygons, and quadratic-spline curves. We show in Figure 1 three such phantoms that we use in our experiments.

C. Sensitivity Models

For computations, we need to parameterize the complex sensitivity maps. It is commonly admitted that they are smooth and slowly-varying spatially. It is therefore possible to generate physically-realistic sensitivity maps using a reasonably small number of lowpass basis functions. Here, we discuss two models that are well-suited for this task. They both relate linearly the parameters to the complex sensitivity values. Moreover, their corresponding MRI models involve the Fourier integrals of monomials over the regions of the phantom.

Definition II.1. For α ∈ N d and R a region of R d f α R (ω) = R r α e -jω.r dr. (4)
Here, we adopted the multi-index notation r α defined in Table I.

1) Polynomial Sensitivity: This model, first proposed in [START_REF] Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF] to represent the local behavior of the sensitivity, assumes that the coil sensitivity S is represented by a polynomial of degree D inside the ROI as

S(r) = D d=0 |α|=d s d,α r α , ∀r ∈ ROI. ( 5
)
As the degree D increases, the model will reproduce sharper transitions. The number of polynomial coefficients is N p = (D + 1)(D + 2)/2. The corresponding MR response is given by

m S (k) = R i=1 ρ i D d=0 |α|=d s d,α f α Ri (-2πk). (6)
2) Sinusoidal Sensitivity: Alternatively, the coil sensitivity is defined by the linear combination of complex exponentials

S(r) = v s v e jr.v , ∀r ∈ ROI. ( 7
)
We propose to constrain the problem to the angular frequencies v on a Cartesian grid with spacings that correspond to twice the considered field of view (FOV). The low-frequency properties are ensured by only considering the L × L angular frequencies around the origin (see Figure 2). Similarly to the effect of the polynomial degree D, an increase in the parameter L allows one to reproduce sharper transitions. The number of coefficients is given by N s = L 2 . The corresponding MR response is given by

m S (k) = R i=1 ρ i v s v f 0 Ri (-2πk -v). (8)
3) Comparison: In order to evaluate and compare the ability of the two models to describe realistic sensitivity maps, we considered a 256 × 256 rasterization of the SL phantom and the 27,648 pixels of its ROI. Using Biot-Savart's law (2), we simulated the complex sensitivity maps of a 24channel circular head coil array (FOV=28cm, distance to center=17cm, radius=5cm) distributed around the phantom. Then, the parameters of the two models where selected to fit the maps. In Figure 3, the average fitting properties of the two models are presented as a function of the number of parameters.

We observe that the fitting accuracy of both models rapidly increases with the number of parameters, with a sensible advantage for the sinusoidal model. The downside is an increased condition number for the fitting operations. With respect to that criterion, the sinusoidal model behaves also better. The maximal spatial errors are comparable for both models.

III. ANALYTICAL MRI MEASUREMENTS

A. Overview of Analytical Fourier Computations

In this section, we present the theoretical tools that are necessary to derive the analytical expression of the MRI measurements. Proofs and additional calculation details are provided in Appendices A, B, and C.

The models presented in the previous section allow us to decompose the analytical MRI measurements into Fourier integrals of the sensitivity over the regions that compose the phantom. Depending on the type of region or sensitivity model, we propose tailored methods to decompose the analytical response as a sum of special functions that can be computed accurately and rapidly. In Figure 4, we present the roadmap of these decompositions that are defined and worked out in the sequel.

B. Elliptical Regions

Let us consider an elliptical region E parameterized by its center r c , the angle θ formed between its semimajor axis A and the abscissa, and its semiminor axis B. The linear transformation

r → u = D -1 R T (r -r c ) , (9)
with D = diag(A, B) and R the rotation matrix of angle θ, maps E into a unit disk, that is to say,

E = {u | u ≤ 1}.
The Fourier transform of the unit disk involves the functions

G n (x) = J n ( x )/ x n , (10)
where J n denotes the n-th order Bessel function of the first kind [START_REF] Vembu | Fourier transformation of the n-dimensional radial delta function[END_REF].

Using the sinusoidal sensitivity model, the integral f 0 E can be worked out [START_REF] Van De Walle | Reconstruction of MR images from data acquired on a general nonregular grid by pseudoinverse calculation[END_REF] as

f 0 E (ω) = 2π|D|e -jω.rc G 1 DR T ω , (11)

MRI model

(1)

Sensitivity • Polynomial • Sinusoidal k-space trajectory Phantom (3) • Ellipse • Bézier
Surface Integrals ( 6) and ( 8) Change of Coordinates ( 9) and ( 14)

f α R in (4)
Ellipse Partial Derivatives (12)
Polynomial Recursion [START_REF] Sorets | Fast fourier transforms of piecewise constant functions[END_REF] Bessel functions [START_REF] Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF] Sinus. where |D| represents the absolute value of the determinant of matrix D.

When considering the polynomial sensitivity model, we suggest to first consider the change of variables [START_REF] Ngo | Three dimensional digital polyhedral phantom framework with analytical fourier transform and application in cardiac imaging[END_REF], rather than computing f α E directly. We write that

E u α e -jω.r dr = 2π|D|j |α| e -jω.rc ∂ |α| G 1 ∂x α DR T ω . ( 12
)
The interesting point is that the partial derivatives ∂ |α| G 1 /∂x α can be decomposed recursively as a sum of G n thanks to the property

∇G n (x) = -xG n+1 (x). ( 13
)
The coefficients of the polynomial in terms of the new coordinates (9) are required to satisfy

S(r) = D d=0 |α|=d s d,α r α = D d=0 |α|=d t d,α u α . ( 14
)
They can be computed by inverting the matrix that relates the N p coefficients to the sensitivity values at N ≥ N p randomly chosen points in terms of the new coordinates.

The MR contribution of such an elliptical contour is presented in Table II.

C. Piecewise-Quadratic Contours

In this section, we first provide relations for the computation of the d-dimensional Fourier transform of a monomial delimited by a connected subset B of R d . With methods that are similar to the ones used in [START_REF] Sorets | Fast fourier transforms of piecewise constant functions[END_REF], we show how to decompose the d-dimensional Fourier integral into a sum of integrals over the contour ∂B. These summed integrals are of reduced dimensionality. In a second step, we show how quadratic-spline curves involve a family of 1-D integrals.

1) Fourier Transform of Monomials over a Connected Set: We show that the surface integral f α B in (4) can be decomposed into a sum of contour integrals.

Definition III.1. We define

g α B (0) = ∂B r α+e k 1 + α k e k .ndσ, for any k, (15)
g α B (ω) = ∂B r α e -jω.r ω 2 ω.ndσ, ∀ω ∈ R d \ {0} ,( 16
)
where n stands for the outward-pointing unit normal of boundary element dσ. Note that g α B is not continuous at the origin ω = 0.

Theorem III.2. For ω ∈ R d \ {0} and α ∈ N d , f α B (ω) = j α m=0 -jω ω 2 α-m |α -m|! α m g m B (ω), (17)
and

f α B (0) = g α B (0). ( 18
)
The consequence of Theorem III.2 is that the d-dimensional integral f α B can be decomposed into a sum of (d -1)dimensional integrals. By recursion, the f α B can be computed via 1D integrals. The proof is provided in Appendix A.

Note that the case ω = 0, which corresponds to the calculation of the moments of the region, has been worked out first by Jacob et al. in [START_REF] Jacob | An exact method for computing the area moments of wavelet and spline curves[END_REF] for parametric 2-D spline contours.

Sensitivity model

Ellipse Quadratic Bézier curve Sinusoidal [START_REF] Gach | 2D & 3D Shepp-Logan phantom standards for MRI[END_REF] 2π|D|e j2πk.rc

v sve jv.rc G 1 (DR(-θ)(2πk + v)) v svf 0 B (-2πk -v) Polynomial (14) 2π|D|e 2πjk.rc D d=0 |α|=d j |α| t d,α ∂ |α| G 1 ∂x α (2πDR(-θ)k) D d=0 |α|=d s d,α f α B (-2πk) c 0 c 1 c 2 c 3 c 4 c 5 r 0 r 1 r 2 r 3 r 4 r 5 r (t) r(t) r ⊥ (t) n(t)
Fig. 5. Example of a quadratic-spline-defined region with N = 6 control points. The boundary is described counter-clockwise and the normal vector is outward-pointing. The contour is represented piecewise by quadratic Bézier curves.

2) Parameterization of a Contour in 2-D: The region B is defined by its boundary, the contour ∂B. In 2-D, a convenient way to parameterize the contour is by the use of a B-spline generating function ϕ such that

∀r ∈ ∂B, ∃t ∈ R, r(t) = p∈Z c p ϕ(t -p). (19)
The considered contour is closed. Consequently, the vectorvalued function r must be periodic. In addition, the number N of coefficients c p that characterize the curve must be finite. The simplest way to satisfy these constraints is to impose that the sequence of coefficients c p be N -periodic. This enforces the N -periodicity of r.

If we note ϕ p the N -periodized version of ϕ, the contour is parameterized either globally as

∀t ∈ [0, N [ , r(t) = N -1 q=0 c q ϕ p (t -q) (20)
or piecewise, with

0 ≤ t = n + λ < M , n ∈ 0 . . . N -1 and λ ∈ [0, 1[, as r(λ + n) = N -1 q=0 c n-q ϕ p (λ + q). (21)
3) Decomposition of the Contour Integrals: We introduce the notation z ⊥ for the vector perpendicular to z with same norm and pointing outwards the region B at the considered point (see Figure 5). We write r (t) = ∂r ∂t (t). The piecewise representation of the contour (21) can be exploited to decompose the contour integral of interest, for instance [START_REF] Gautschi | Computational aspects of three-term recurrence relations[END_REF] or [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], which leads to ∂B F(r).ndσ = N -1 q=0 1 0 F (r(q + λ)).r ⊥ (q + λ)dλ. [START_REF] Wilm | Higher order reconstruction for mri in the presence of spatiotemporal field perturbations[END_REF] 4) Quadratic Bézier Curves: In the sequel, we focus on contours represented by linear and quadratic B-splines. The former describe polygons while the latter give a piecewise description of quadratic Bézier curves. Three equivalent piecewise representations can be useful and are given in Table III with their relationships.

Definition III.3.

h (m) (a, b) = 1 0
λ m e -jλ(a+λb) dλ.

(

) 23
Proposition III.4. For ω ∈ R d \ {0} and a contour ∂B parameterized piecewise by r(λ + n) = r n + λβ n + λ 2 γ n , with n ∈ N and λ ∈ [0, 1], we have that

g α B (ω) = 1 ω 2 N -1 n=0 e -jω.rn 2|α|+1 i=0 d n,i h (i) (ω.β n , ω.γ n ) (24) while g α B (0) = N -1 n=0 2|α|+3 i=0 d n,i h (i) (0, 0), (25)
where, with the notation

σ p = |p 2 | + 2|p 3 |, the symbol d n,i stands for p1+p2+p3=α r p1 n β p2 n γ p3 n α! p 1 !p 2 !p 3 ! ω. β ⊥ n δ σp,i + 2γ ⊥ n δ σp,i-1 ,
(26) and d n,i stands for

p1+p2+p3=α+e k r p1 n β p2 n γ p3 n α! p 1 !p 2 !p 3 ! e k . β ⊥ n δ σp,i + 2γ ⊥ n δ σp,i-1 . (27)
The values h (m) (a, b) follow a three-term recurrence relation [START_REF] Gautschi | Computational aspects of three-term recurrence relations[END_REF]. More details on their numerical computation are given in Appendix B.

Note that the piecewise parameterization of the contour of a polygon corresponds to the particular case of a quadratic parameterization with β n = r n+1 -r n and γ n = 0. Such simpler polygonal models with homogeneous sensitivities have been considered in prior work [8, Prop. 3.2] using a similar formulation.

IV. EXPERIMENTS A. Implementation Details

Our implementation uses Matlab 7.12 (Mathworks, Natick). The experiments run on a 64-bit 8-core computer, 8GB RAM, Mac OS X 10.6.7.

We implemented the analytical computations as described by the scheme in Figure 4, with double float precision. For

λ 2 2 cn + 1 2 + λ -λ 2 c n-1 + 1 2 (1 -λ) 2 c n-2 rn = 1 2 (c n-1 + c n-2 ) Polynomial (piecewise) rn + λβn + λ 2 γn βn = 2 (c n-1 -rn) Bézier curve (design) (1 -λ) 2 rn + 2λ(1 -λ)c n-1 + λ 2 r n+1 γn = r n+1 + rn -2c n-
(m) (a, b)| ≤ 1/(m + 1
) is sometimes violated for orders m ≥ 2 and large values of the first argument. This prevented us to present valid simulations of piecewise quadratic contours combined with a polynomial sensitivity. Given the comparison of the two models in Section II-C3, we considered the sinusoidal model with parameter L = 7, that is N s = 49 in Figure 3, which lead to accurate representations of the physical sensitivities and numerically tractable inversions.

As an alternative to our analytical method, we consider the traditional simulation procedure that consists in i) sampling the phantom with a grid of a given size and ii) resampling the DFT of this discrete image according to the desired k-space trajectory. We call this procedure a rasterized simulation. It is expected to be consistent with our analytical method only when considering an infinitely dense sampling.

For reconstructions, we consider an optimization problem of the form

x = arg min x m -Ex 2 2 + λP(x), (28)
where x represents an image, x is the reconstructed one, m is the concatenated scanner data vector, E is the encoding matrix, and P is a regularization function. With N receiving channels and M k-space measurements at positions {k m } m=1..M , the MRI encoding matrix E is formed as with I n representing the n × n identity matrix, the symbol ⊗ standing for the Kronecker product, and s n being the nth coil sensitivity map vector in the same way as x. The encoding matrix E 0 corresponds to the same MRI scan with a single homogeneous receiving coil and is defined as

E = (I n ⊗ E 0 ) [diag (s 1 ) , . . . , diag (s N )] T , (29)
E 0 = [v 1 , . . . , v M ] T . ( 30
)
There, v m are vectors such that, for a pixel of coordinates p, (v m ) p = exp (-2jπk m .p).

We used two types of regularizations in our experiments

• P(x) = x 2 2
corresponds to a Tikhonov regularization, which leads to linear reconstructions that we implemented with the conjugate-gradient method.

• P(x) = x TV is the isotropic total-variation pseudonorm, which leads to a nonlinear reconstruction problem. This reconstruction scheme is often used in compressedsensing research and is particularly suited for dealing with piecewise-constant images such as our phantoms. We implemented it using the iteratively reweighted leastsquares algorithm, also known as the additive form of the half-quadratic minimization [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], [START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF].

B. Validation of the Implementation 1) Simple Example with Homogeneous Sensitivity:

As first validation, we consider the simple phantom composed of a rectangular region that is represented in Figure 1. Under a proper change of variables, it yields a square and its Fourier transform is given by a product of sinc functions. This phantom is composed of a polygon and consequently falls in the category of the spline-defined contours. We test the accuracy of our proposed simulation method and of the rasterized approach against the closed-form solution. To do so, we consider the MR response associated with a homogeneous receiving coil sensitivity and a 256 × 256 Cartesian k-space sampling. The simulation errors are reported in Tables IV andV.

As expected, the error of rasterized simulations decreases when the sampling density increases. Meanwhile, the accuracy of our analytical implementation is as good as the machine double float precision would allow. Thus, we conclude that we can indistinctly use the closed-form ground truth or our proposed analytical model in the conditions of Section IV-B1. 2) Validation with non-Homogeneous Sensitivity: We now use our analytical phantom as a gold standard to evaluate the accuracy the measurements obtained from rasterized simulations. We consider the SL and brain phantoms. The single sensitivity map is computed using Biot-Savart's law and is approximated on the support of each phantom with the sinusoidal model. The k-space is on a 128 × 128 Cartesian grid. Errors are reported in Table VI.

We observe that the errors decrease with the same trend as in the rectangle case, which strongly suggests that our gold standard is accurate. Meanwhile, for a given sampling density, the errors occurring with the SL phantom are consistently larger than the ones corresponding to the brain phantom. This is explained by the fact that the SL phantom presents edge transitions of larger intensity. When the benefits of an analytical model are forsaken, the MRI data are generated from a rasterized version of the phantom and the sensitivity, using the (non-necessarily uniform) discrete Fourier transform (DFT)

m M (k) = |M|F e -2πjMk , (31)
with Mk ∞ ≤ 1/2 and

F e 2jπν = p∈Z 2 f (p)e -2jπp.ν = q∈Z 2 f (ν + q) . ( 32
)
The right-hand side of (32) can be worked out using Poisson's summation formula. The terms with q = 0 represent the aliasing that occurs with rasterized simulations. Due to the intrinsically discontinuous nature of the phantom ρ, the Fourier transform f decreases slowly, leading to significant aliasing artifacts. However, as the sampling density increases (Tr M → 0), the impact of aliasing is reduced, as we saw in Section IV-B. Let us define an ideal anti-aliasing filter h in the Fourier domain as

ĥ(ν) = 1 if ν ∞ ≤ 1/2, 0 otherwise. ( 33
)
For normalized frequencies ν such that ν ∞ ≤ 1/2, the analytical simulation (unaliased) is characterized as the DFT of the samples of the lowpass-filtered continuous signal

f (ν) = p∈Z 2 (h * f ) (p)e -2jπp.ν , (34)
where (h * f ) represents the spatial continuous convolution of h and f . When using a full Cartesian k-space sampling, the classical approach to reconstruction is to perform an inverse DFT. In this case, the samples of the signal f will be perfectly recovered out of the rasterized simulation (32) which is not desired because it conceals the existence of the Gibbs phenomenon due to the antialiasing filter (see, for instance, [START_REF] Archibald | A method to reduce the Gibbs ringing artifact in MRI scans while keeping tissue boundary integrity[END_REF]). By contrast, the data provided by our analytical model lead to a fairer reconstruction where the Gibbs phenomenon appears. This effect is illustrated in Figure 6.

Counterintuitively, the reconstructions out of rasterized simulations lead to aliasing effects that have a positive impact on visual quality. This situation, which occurs when the same model is used for both simulation and reconstruction, is sometimes referred to as "inverse crime". It arises because of the artificially imposed consistency between the computational forward models used for simulation and reconstruction. In such an inverse-crime situation, the continuous nature of the underlying physical model is not taken into account.

2) Impact of Rasterized Simulations on Reconstruction: We consider a plausible pMRI setup. It involves an array of 8 receiver coils that are uniformly distributed around the phantom. The corresponding sensitivity maps are computed according to Biot-Savart's law. Spiral and EPI k-space trajectories are considered, both supporting a 256 × 256 image with reduction factor R = 4. The simulated raw data are generated using our analytical method as well as 256 × 256 and 512 × 512 rasterized simulations (see Section IV-C1). The same realization of complex Gaussian noise is added to the simulated data with different intensities, according to three scenarios: very low noise (40dB SNR), normal data (30dB SNR), and very noisy data (20dB SNR). Reconstructions are performed using quadratic (Tikhonov linear solution) and TV regularizations. The reconstruction algorithms use the same forward model, in the form of the same encoding matrix E. The experiments only differ in terms of the input data. The regularization parameter is tuned to optimize the SER with respect to the ground-truth phantom (256 × 256 rasterization of the phantom). We report our results in Table VII for the spiral trajectory and in Table VIII for the EPI experiments. Reconstructed images are shown in Figures 7 and8, together with their error maps, in order to illustrate the impact of the inverse-crime situation (the 256×256 rasterized simulation) in the different scenarios.

The reconstructions in the spiral experiment are penalized compared to the EPI ones, in the sense that the high-frequency corners of the k-space are not sampled which leads to slightly inferior resolution. This explains that, all other parameters remaining constant, the EPI reconstructions outperform the spiral ones qualitatively and quantitatively. We observe that the reconstructions from rasterized simulations consistently outperform the ones obtained from analytical measurements. While large differences can occur between the inverse-crime scenario (the 256 × 256 rasterized simulations) and the analytical simulation data, the 512 × 512 simulations yield much closer performance, with at most a 0.6 dB SER difference. This is explained by the reduced aliasing artifacts when doubling the sampling density (see Section IV-C1). As expected for this type of piecewise-constant phantom, the TV reconstructions consistently outperform the linear ones. Whatever the simulation method is, TV brings a significant improvement in the very noisy scenario. However, for the other scenarios (SNR 30dB and 40dB), the improvement over linear reconstruction is modest when using the analytic measurements, whereas it is artificially spectacular using the 256 × 256 rasterized simulations. We believe that our quality assessment, obtained analytically, offers fairer predictions of the practical worth of a reconstruction method than its overly optimistic rasterized version.

V. CONCLUSION

We proposed a method to develop realistic analytical phantoms for parallel MRI. Our analytical phantom approach offers strong advantages for the quantitative validation of MRI and pMRI reconstruction softwares: it is flexible enough to represent general imaging targets, it provides highly accurate representation of the physical continuous model and avoids overly optimistic reconstructions. This kind of framework is also applicable to the assessment of advanced MRI reconstruction methods such as autocalibrating parallel imaging, B0 correction [START_REF] Man | Multifrequency interpolation for fast off-resonance correction[END_REF], motion correction [START_REF] Zaitsev | Magnetic resonance imaging of freely moving objects: prospective real-time motion correction using an external optical motion tracking system[END_REF], [START_REF] Ooi | Prospective real-time correction for arbitrary head motion using active markers[END_REF], or higher order field imaging [START_REF] Wilm | Higher order reconstruction for mri in the presence of spatiotemporal field perturbations[END_REF]. Implementations of the phantom are made available to the community.

C. Proof of Proposition III.4

Proof: We rewrite g α B using Characterization [START_REF] Wilm | Higher order reconstruction for mri in the presence of spatiotemporal field perturbations[END_REF] with F (r) = r α e -jω.r ω 2 ω for ω = 0 and F (r) = r α+e k 1+α k e k for ω = 0. The piecewise parameterization of the contour (Table III) is then used, and by virtue of the multinomial theorem, we expand the terms r α and r α+e k .

Fig. 1 .Fig. 2 .

 12 Fig.1. Phantoms parameterized by elliptical and Bézier-defined regions. From left to right: rectangle phantom used for validation, the Shepp-Logan phantom, and a proposed brain phantom. The PDF and SVG versions of these phantoms are available online at http://bigwww.epfl.ch/algorithms/mriphantom/.

Fig. 3 .

 3 Fig.3. Fitting properties of the two sensitivity models as a function of the number of parameters. From left to right: approximation signal to error ratio (SER) in dB, maximal absolute error, and condition number of the fitting matrix.

Fig. 4 .

 4 Fig. 4. Roadmap of analytical Fourier computations.

C. Applications 1 )

 1 Investigation of Aliasing Artifacts: Let us consider the function f (u) = Sρ(Mu) which depends on the spatial sampling step matrix M. According to (1), the analytical MR data are given by m S (k) = |M| f (-Mk).

Fig. 6 .

 6 Fig. 6. 64 × 64 SL full Cartesian sampling reconstructions. From left to right: analytical simulation, rasterized "inverse crime" simulation, and corresponding line profiles.

Fig. 7 .

 7 Fig. 7. Reconstructed Brain phantoms and error maps for the spiral SENSE experiment.

TABLE II CLOSED

 II -FORM MR CONTRIBUTION OF ELLIPTICAL AND QUADRATIC BÉZIER REGIONS FOR THE PROPOSED SENSITIVITY MODELS.

  Our package also includes graphical tools to design the analytical phantoms. For purposes of adequate visualization, export to the popular vector-graphics formats SVG 1.1 and PDF (via the PGF/Tikz L A T E X package) is supported. The package is distributed4 in order to provide sensitivity fitting, phantom-design interface, analytical simulation tools, and to allow replication of the experiments of this section.Unlike the sinusoidal model which is very robust to numerical errors, our current implementation of the three-term recurrence relation (see Appendix B) leads to instabilities when using the polynomial model. The theoretical relation |h

	1
	efficient computations of the error function of a complex
	variable, we coded the critical parts of erfz in C++/MEX,
	with POSIX multithreading, following Marcel Leutenegger's
	recommendations 1 . The code implementing Theorem III.2
	utilizes Matt Fig's npermutek 2 . The rasterization of spline-
	defined regions, which is performed without approximation,
	partly relies on Bruno Luong's fast MEX implementation of
	insidepoly 3 .

TABLE V ERRORS

 V OF THE RASTERIZED SIMULATIONS FOR THE RECTANGLE PHANTOM.

	Sampling density	256	352	400	512	704	800	1024	1408	1600	2048
	NRMSE	5.58e-02	2.51e-02 2.01e-02	1.25e-02	7.45e-03	6.04e-03 3.85e-03	1.59e-03	1.27e-03 8.32e-04
	max. error in k-space	5.5e-03	2.0e-03	1.5e-03	1.1e-03	5.5e-04	5.3e-04	3.6e-04	1.2e-04	1.0e-04	5.8e-05
	max. error inverse DFT	5.5e-01	1.9e-01	1.6e-01	1.0e-01	6.2e-02	5.4e-02	3.9e-02	1.5e-02	1.2e-02	6.7e-03

TABLE VI ERRORS

 VI OF THE RASTERIZED SIMULATIONS FOR THE BRAIN AND SL PHANTOMS VERSUS OUR ANALYTICAL SIMULATIONS.

		Sampling density	128	176	256	352	512	704	1024	1408	2048	2816
		NRMSE	1.45e-01 9.26e-02	5.45e-02	3.48e-02 2.13e-02	1.02e-02	6.70e-03	4.06e-03 2.03e-03	1.49e-03
	Brain	max. err. k-space	1.1e-02	5.7e-03	3.3e-03	2.6e-03	1.4e-03	6.6e-04	4.1e-04	2.4e-04	1.5e-04	9.5e-05
		max. err. inverse DFT	2.3e-01	1.4e-01	1.1e-01	7.9e-02	5.1e-02	2.0e-02	2.1e-02	1.4e-02	4.8e-03	5.6e-03
		NRMSE	2.76e-01 1.79e-01	9.74e-02	5.38e-02 2.85e-02	2.01e-02	1.28e-02	6.23e-03 3.34e-03	2.03e-03
	SL	max. err. k-space	2.9e-02	1.6e-02	8.8e-03	4.9e-03	2.6e-03	1.7e-03	1.0e-03	6.1e-04	3.0e-04	1.7e-04
		max. err. inverse DFT	4.7e-01	3.0e-01	1.6e-01	1.1e-01	6.5e-02	3.9e-02	3.3e-02	1.3e-02	7.0e-03	5.0e-03

TABLE VII RECONSTRUCTION

 VII BIAS OF RASTERIZED SIMULATIONS ON QUADRATIC AND TV REGULARIZATION RECONSTRUCTIONS OF THE BRAIN PHANTOM FOR THE SPIRAL SENSE EXPERIMENT. OPTIMIZED SER AND CORRESPONDING BIAS ARE SHOWN IN DB.

	Channel data SNR	40dB	30dB	20dB
	Sampling density	256	512	256	512	256	512
	Linear	SER Bias	24.61 19.92 5.07 0.37	20.31 17.99 2.56 0.24	14.09 13.45 0.75 0.11
	TV	SER Bias	33.75 20.80 13.45 0.49	27.60 20.26 7.75 0.42	19.61 17.72 2.43 0.54

TABLE VIII RECONSTRUCTION

 VIII BIAS OF RASTERIZED SIMULATIONS ON QUADRATIC AND TV REGULARIZATION RECONSTRUCTIONS OF THE BRAIN PHANTOM FOR THE EPI SENSE EXPERIMENT. OPTIMIZED SER AND CORRESPONDING BIAS ARE SHOWN IN DB.

	Channel data SNR	40dB	30dB	20dB
	Sampling density	256	512	256	512	256	512
	Linear	SER Bias	36.25 20.77 16.02 0.54	26.30 19.79 6.95 0.44	16.73 15.31 1.61 0.19
	TV	SER Bias	42.25 20.98 21.85 0.58	32.75 20.73 12.57 0.55	23.92 19.29 5.02 0.39

ACKNOWLEDGMENTS We thank Matt Fig and Bruno

Luong to make their efficient Matlab implementations of npermutek and insidepoly available. We are grateful to Jeff Fessler for his code of NUFFT and interesting comments during ISBI 2008 and 2010. We also acknowledge Dimitri Van De Ville and Isik Karahanoglu who took part in preliminary works.

APPENDIX

A. Proof of Theorem III. [START_REF] Ribés | Linear inverse problems in imaging[END_REF] First, we consider the case ω = 0. Proof: Take ψ(r) = r α and ϕ(r) = (e k .r) 2 /2 = r 2 k /2. Then, ∇ (r α ).e k = α k r α-e k , ∇ϕ = r k e k , and ∆ϕ = 1. Using Green's first identity yields (1 + α k )f α B (0) = ∂B r α+e k e k .ndσ.

For the case ω = 0, we need an intermediate result.

Proof: Use Green's first identity with ψ(r) = r α , ϕ(r) = -e -jω.r , and the fact that x.∇r α = i x i α i r α-ei .

Let us continue the proof of Theorem III.2 by mathematical induction on n = |α|. For n = 0, α = 0 and the result holds true according to Lemma A.1. When considering n = 1, α = e i , and Lemma A.1, we obtain f ei B (ω) = jg ei B (ω) + ωi ω 2 g 0 B (ω). This is true for all i, hence the result holds true for n = 1. Now, we suppose the result holds true at order n and we consider α such that |α| = n + 1. From Lemma A.1, we have that

Since |α -e i | = n, we substitute f α-ei B using the induction hypothesis and, after simplifications, we obtain

with d i (x) = x i /|x| for x = 0 and d i (0) = 0. By permutation of the sums and noting that i d i (x) = 1 for x = 0 and i d i (0) = 0, we get This is valid for all α such that |α| = n + 1. Hence, the induction hypothesis was proved at order n + 1 assuming it holds true at order n.

B. Characterization and Computations of a Family of 1-D Integrals

Proposition A.2. For m ∈ N, h (m) follows the recursion rule a+λb) dλ by parts and identify h (m+1) , h (m) and h (m-1) if m > 0.

Corollary A.3. For small values of a and b, one can rely on the backward iteration starting from a higher order M > m to get accurate results

Proposition A.4. For b nonzero and m ≥ 1, the forward iteration is used

Proof: From Proposition A.2 with m = 0 and b = 0, we get h (0) (a, 0) = e -ja/2 sinc(a/(2π)). In the case b = 0, we define t = λ + a 2b such that λ(a + bλ) =

where the lower incomplete gamma function is defined as γ(s, z) = z s 1 0 λ s-1 e -λz dλ. Proof: Note that e -jλ(a+λb) = e -jλa ∞ n=0 (-jλ 2 b) n /n!. By virtue of Fubini's theorem, we get

(-jb) n h (m+2n) (a, 0)/n!.

Identify γ(m + 2n + 1, ja) to (ja) m+2n+1 h (m+2n) (a, 0).