A now well established and growing community in the surroundings of the OASIS climate model coupling library is perpetually recreating its coupled systems. The coupling library maintenance and upgrading activity, constant but modest, cannot alone explain this popularity. We assume that the constant growth of the OASIS-based model livestock is the consequence of a quasi organic development of climate models in a favourable direction. A three year long period of on-site dedicated supports to the OASIS coupling helped to understand how and why our coupler was used. We rely on this experience to better define the role that modularity is playing in the building of our community. We propose to call OASIS a community coupler, and examine how it partakes of an ethical politics of scientific research, which argues against its dismantling

1-Rationale

Used by more than seventy climate modelling groups in the five continents, the OASIS coupling library [Valcke et al. 2021] plays a crucial role for the development and efficient handling of climate models. The library allows synchronized exchanges of coupling information between components of systems which aims at representing earth climate. Recently, a set of on-site dedicated supports was organised within the IS-ENES-3 and ESIWACE-2 EU infrastructure projects. This gave us the opportunity to help designing, upgrading or enhancing the implementation of the OASIS library and setting up a tailored and computationally efficient coupled systems. But at the same time, this work was an opportunity to try to understand how and why our coupler was used. In this document, we propose to rely on this experience to try to better define the role that our coupler is playing in our community and to argue against its dismantling.

A dedicated support program was already led with success through a networking activity during the first fold of the IS-ENES project [Maisonnave et al. 2013]. The first twelve IS-ENES-1 OASIS dedicated user supports (ODUS) had facilitated the implementation and the use of coupled models in Europe, identifying or resolving practical issues during coding, or simply by discussion. Focusing and isolating the work of both model and coupler specialists on a given time period (one month), the ODUS gave enough time to precisely identify issues, sometimes involving more than one modelling group at the same time. Sometimes, a feedback effect of the work led to the library enhancement, with an added benefit for the whole community. This first ODUS program gave us a better idea about present and future model community requirements and contributed to tailor the OASIS following versions in that direction, diverting from trends driven by exogenous desires and futile fashions. The 2009-2013 ODUS program also contributed to disseminate OASIS coupling best practices through laboratories and created an onset of OASIS experts network. An extension of this activity, though lighter, was lonely sustained by CERFACS from 2013 onwards [Oouchi et al. 2015;Will et al. 2017].

Taking benefit of the same chain of infrastructure projects, the original 90's version of the coupler was enhanced, by adding interpolation (SCRIP) [Jones 1999] and parallel communication (MCT) [Jacob at al 2005] libraries. A comprehensive documentation was maintained. Access to beginners was facilitated by training sessions. The management of these activities at a continental level advertised the tool on a set of workshop and conferences.

But in our sense, this activity, constant but modest, cannot alone explain the success of our library. From our point of view, the constant growth of the OASIS-based model livestock is less the outcome of an imperialism whose tanks remain to be detected by the spy satellites, than the consequence of a quasi organic development of climate models in a favourable direction.

Of course, such a bold postulate needs to be corroborated but the full demonstration, if reachable, would require, in addition to our practical experience, a larger theoretical base in fields such as science epistemology, or even sociology. At our climate and computing scientist level, just let us stress a selected number of practical features of the OASIS software that should, in our sense, facilitate the climate modelling process and explain its success.

Climate modelling is fundamentally a process in motion. Mainly because climate research is a dynamical process per se, which modifies the object that it aims to represent, but also because of the nature of the representation tool: a computing code, sensible to the evolution of the underlying hardware. Climate is a complex topic, which has been studied for decades. The complexity of a climate model grows during time, since it is integrating more and more complex processes but also because it is costly to simplify it by rewriting the original code (legacy) too often1 .

One of the main benefit of the OASIS library lies in the initial assumption that has governed its conception: a climate system is made of more than a single atmosphere model, and more components need to be added at its interface rather than prescribing fixed boundary conditions. The modularity induced by this definition has its counterpart: the arbitrary concept of internal limits within the system, which makes more difficult to think and represent a continuum. The representation of any physical interface is made more complex because it is split between at least two models. These two drawbacks tend to dissuade from more detailed studies of the physical interface such as [ Pelletier et al. 2021]. In our sense, this is the main limit of an OASIS induced representation of the climate system, which present document does not deny. But we assume that it is still possible to make progress in our science within this limit, which will not affect the main conclusions of this document.

The benefit of modularity is high regarding code development, since components developed on another laboratory are easily reusable, even more if these components are already endowed with an OASIS interface. In addition, updates in one module do not affect the whole system. These benefits reduce the development cost, limit duplicated developments in the community, thus save time for more precise and complementary studies using the same model (e.g. coupled systems using the same NEMO ocean propose a comprehensive understanding of this model behaviour in atmosphere coupled mode). Since the OASIS tool is long lasting2 , it is reliable and has incited the community to bet on the long term availability of the model they develop and of the paradigm on which it is based. This has strengthened the possibility of long term studies with the same model.

Another benefit of the OASIS interface is its distributed writing and maintenance. Users of one or several modules are defining the needs, and can participate to the coding of the interface but also of parts of the OASIS library itself, all written in an accessible and standard language: FORTRAN. A guarantee of the library continuous developments is ensured by the conceiver's laboratory (CERFACS), while the community involvement ensures adequacy to the needs. A "weak" copyleft, induced by the GNU Lesser General Public License that governs the OASIS use, secures the possibility for the community to freely adapt the library to their needs while preventing the (benign) risk of seizures.

Last but not least, OASIS gives you the possibility to build your own coupled system. As stated by Michel de Montaigne in his "Essais", "j'y veux pouvoir quelque chose du mien": I would have a power of introducing something of my own 3 . This freedom that our tool offers is granted back by a willing participation to the tool development and the community permanence. With a better knowledge of the tool, freedom of development increases. In interaction with a community, possibilities widen. One would see a positive feedback to this mechanism, and another possible element of this supposed quasi organic process that led the development of the OASIS community.

Apparently, a now well established community in the surroundings of the OASIS library is perpetually recreating its coupled systems, spreading them on a larger extend, and also contributes to the coupling tool update, thus its relevance and its ability to convince more users. Can we say more about the mechanisms at work behind this development ?

2-A community coupler

The 2019-2021 OASIS dedicated support program was the ideal view point to try to identify our community needs and practices. The reader will find the comprehensive description of the eight ODUS provided in the three yearly reports made available during the IS-ENES-3 project The description is intentionally reduced in this document, where we focus on summarising the practices that help to understand how OASIS is used, and then try to deduce from that how it can be adopted by a larger community.

2.1-Usage

During the first year of the project, we could upgrade the OASIS3-MCT library owned by the users to the current version 4 and make available new coupler functionalities in their coupled systems, such as the parallel computing of interpolation weights. Interfaces were modified to A rough classification of these tasks we achieved during the ODUS could give a rather comprehensive view of the coupling library potentialities. The classification would discriminate between (i) the set up of new coupled system, (ii) the upgrade of the porting of an existing framework and (iii) the use of new OASIS functions.

(i) New coupled system set up

The coupled systems we had to study were mainly built with well established components. The ODUS work may consist in adding a new interface, if the component is a newcomer in the OASIS zoo, e.g. adding the MFWAM wave model to NEMO ocean for a future Météo-France weather forecast system. It may also consist in an upgrade of the existing interface, if the model is used in a new configuration, e.g. a two way coupled AGRIF zoom included in the original OpenIFS-NEMO ocean-atmosphere model at GEOMAR. More exceptionally, an existing model can be split into independent pieces, that are coupled in a second step with an OASIS interface. Computing performance is the key motivation to such operation. This was done at MetOffice, for both SI3 sea-ice and TOP-MEDUSA bio-geo-chemistry (BGC) modules of the NEMO framework.

The extreme case of the neXtSIM sea-ice coupling at NERSC is worth to be noticed. The discretisation of the existing NEMO ocean model was simplified (removal of the tri-polar folding zone) and proposed as a new sea-ice model grid, the coupling of the two models being ensured, without error prone interpolation, by the existing NEMO ocean interface.

These examples, experienced within the ODUS program, make clear that modularity is the primary need that favours the use of the library OASIS. The addition of a new component is done quickly, thanks to the low intrusiveness of the OASIS API. When it is done, the prime core model boundary conditions can be provided from static data or coupled fields. They are ready to be filled by coupling fields coming from a future component, if this should have been interfaced with OASIS API routines.

(ii) Upgrading and porting of an existing coupled system

As already noticed, the non intrusiveness of the OASIS API facilitates the standard exchange of one component, particularly when a new version has to be included in the user system, as it was experienced at NERSC or ETHZ. Porting of the whole coupled system can be much more challenging, particularly when made on non CPU architectures. Not because of the compatibility of the library itself, not used on GPU and already compatible with vector machines, but because of the heterogeneity of the computing units that the coupled system were gathering on the supercomputers, e.g. IO servers on CPU and computing processes on vector engines at DWD. Assuming that the trend of hardware manufacturing will keep favouring, in the near future, the assembling of specialised chips [ Thompson & Spanuth 2018 ], the ODUS has helped validating the hypothesis that a modular system (i) can be ported on such architectures and (ii) makes the most of the heterogeneous computing power by distributing each element of the modular system on the kind of hardware where it is the most efficient.

(iii) Use or implementation of new OASIS functions

The third class of ODUS tasks we can identify gathers (i) the use of new OASIS functionalities, related to the core features of the library, such as its interpolations (locally conservative interpolation for runoff, GEOMAR [ Voldoire 2019, Maisonnave 2020 ]) or additional tools (new load balancing tool, MetOffice) and (ii) the use of existing functions, but in innovative ways : parameter tuning of the Gaussian distance weighted interpolation (GEOMAR), joint use of heterogeneously (scalar/vector machines) compiled OASIS libraries (DWD), validation of the pyOASIS API in a simplified coupled system, used during a post-processing step to calculated in parallel (OpenMP) weights and addresses (SMHI). Doing that, we contribute to test functions of the library not often used by the community, and make sure that the last OASIS version can handle them (bug reports, commits into the repository).

More scarcely, bug can be corrected and new functions can be added by users to the OASIS library itself. And this is particularly facilitated by the readiness of the FORTRAN language and the profusion of the inlined commentaries. Our implementation of an efficient 3D coupling (without vertical interpolation) and the replacement of the average weighting operation by a maximum between source grid points, shows at least that adapting the library to user needs is possible. Of course, this was done by experienced users and we have to admit that, in history, only a few number of users are at the origin of OASIS modifications. And no new experimental observation can corroborate the assumption that the opportunity to modify the OASIS library is intensively taken in the community.

2.2-Impact on community

The modular design of OASIS made coupled systems is what greatly improves the collaboration of laboratories willing to share the same modules. Rather than developing an new model, it is easier to work on a component developed elsewhere by better specialists, even if this was done on different coupled architectures, e.g. OpenIFS-NEMO in global or regional mode, with or without AGRIF zoom, with or without a BGC component, etc. During the IS-ENES-3 ODUS program, we observed that the OASIS interfacing of a set of new components, or the upgrade of this interface in existing components (see Table 2) has facilitated their redeployment to another laboratory of the community (see Table 3). The last axis of what could be considered as a community effect is the enhancement of the coupling library by its users. The existence of this axis remains unclear, since the library upgrade mechanism is centralised and marginally affected by user contributions. However, user questions are compiled and addressed in the library website forum, and the users and steering committees are consulted by the OASIS development head. During our ODUS experience, library bugs that slows down, sometimes prevents, the implementation of a coupled system, were reported. But their number is modest (less than a dozen), even though they could have there importance regarding to the overall reliability of the coupler and the trust that the community can award to it. Most of them are ready to be integrated to the next release. However, the inclusion of library deeper modifications (3D coupling, source maximum, single precision computations) is less certain.

Despite this restriction, that will have an importance for the future of the OASIS library, the ODUS experience has confirmed that both an ensemble of intrinsic properties and an ensemble of usages (which proceed from these intrinsic properties) allows to call our coupler a community coupler. In that multiple sense that (i) it is a widely spread tool for building (coupled systems), (ii) that this technical building must be done in a community, (iii) that the built product can be adapted to other needs for a larger community, (iv) that this sharing also propagated building techniques that also contributes to build an even larger community, (v) that finally forms a meta-community which participates to the tool enhancement.

3-Recommendations

The maintenance of the OASIS framework itself may be subject to a weakening in the years to come. The main CERFACS' contribution is languishing, and its planed partition from EU infrastructure projects will necessarily reduce the development and maintenance effort. We hope to have convinced the reader of the importance of maintaining a community coupler. In our opinion, this requires to focus our attention and some part of our working efforts in prospective actions.

Despite its evergreen community of users and its growing popularity, OASIS is materially relying on a code which is at risk. As previously emphasised, the endogenous code development by the community itself is the weak link in the chain of subtle actions and feedbacks that sustains the library success. Ways of empowering this feedback must be found to have a chance to adapt the code to our science evolution. But in case of failure, can the solution come from another coupling software developed in the community ?

From at least the end of the EU PRISM infrastructure project onwards [ Valcke et al. 2006 ], a centripetal force contrives to multiply the number of couplers or coupling frameworks. This objectively contributes to divide and weaken a common effort, and to forever defer the possibility of a worldwide models interoperability. From now on, these softwares are devoted to the coupling of a dedicated number of models, following national or laboratory coopted criteria of choice. This latter promotes integrated coupled models, making practically impossible any change in the laboratory pre-defined module list. Some of these softwares, making use of object oriented languages, jeopardise their source openness. And to crown it all, a lingering trend towards high performance, which now verges on gigantism, darkens the climate model computing codes (computing oriented languages, opaque libraries, etc), thus their tied proprietary coupled libraries.

Though justified at short term, these deliberate choices denies some long term tendencies and recent scientific clarifications that should make us ruling such dangerous policies out :

1-It is now obvious that the CPU multi purpose era is ending and most of the computers and supercomputers available include specialised processors such as GPU [ Thompson & Spanuth 2018]. This heterogeneity, and the short lifetime of this market guided hardware, necessarily suppose the adaptability and portability of our software and should prevent to orient the community corpus towards esotericism. At the opposite, a legible and simple coding syntax must facilitate the involvement of the climate scientists in their software library development and upgrading, including for couplers.

2-Despite the unquestionable media coverage of the Global Warming Issue, external societal trends strongly suggests that the expansion phase of our community is also behind us. The conjunction of the rarefaction of public investment to poorly profitable sciences4 , from one hand, but also from the now multi-decadal summons to energy savings and more recent stresses on raw material production, thus semiconductor industry, on the other hand, should stimulate more imaginative approaches of climate model development, less blindly grounded on cheap labour force and unlimited computing power.

3-Despite our thirst of explicit phenomena representation, tradeoff with parametrisation in climate models is not questionable, because of the absence of any target resolution where one could argue that everything of interest to climate is resolved [ Balaji et al. 2022 ]. The race to more resolution, thus uncontrolled growth of computing resources, encouraged by GPU migration, is unsustainable, converges to make our codes and couplers opaque and reserved to a selected number of trustee laboratories. This finally sends most of our scientists away from their models understanding and evolution.

At the opposite, modularity that our community coupler brings can be the pivotal concept for sustainability. 1986 ] during times of more severe restrictions in the German Democratic Republic (GDR), and for a totally different purpose. Motorbikes conceived by this Bauhaus follower were built in such a way that the maintenance, the repairing and the component replacement (particularly in case of significative technical innovation) were facilitated. Focused on function, not fashion, the conception was combining standardisation (of the components) and variety of choice (of the overall system). Tunable, the vehicle was fitting the users needs at its best possible. Despite the unsteadiness of the GDR economy, it remained in production for decades and some specimen are still in use today5 . This is this kind of modularity that a coupling library should reflect to make it more durable. For that, the involvement of the user in the model maintenance and evolution is of paramount importance. To face limitations that will totally redefine our working habits, an Openness. Rather than a Weltanschauung that self consistent systems implicitly carries. To avoid as much as possible preconceived ways of thinking, to facilitate implementation by and for scientists,

this is what has to be constantly stimulated. To favour a practice, a pedagogy, instead of the passive use of an object, which maintenance is untrusted to its authors and which comprehension is more and more out of reach of its users. A practice is made possible by the modularity at the basis of the concept of community coupler, which facilitates module switch and comparison, freedom of choice, in the path of fulfilling true needs by the comprehension of means. A continuous back and forth between practice and theory founds, not a topic, but a subject: manners of defining climate complex systems and their module interfaces.

This idea makes possible to think a seamless poetics6 -ethics-politics: a community coupler is the place where points of view are put together, and selected by the community itself in a visible process, through an informed self selection, in short, an ethical politics of scientific research. At the opposite of a Research of Prestige. A modernity, against all obscurantisms, and despite all difficulties. Though this latter was the only reason to do it, we would recommend to maintain, by any mean, our community coupler.

  [ Maisonnave 2019 ; Maisonnave & Kjellsson 2021 ; Maisonnave 2022 ], jointly with the two reports of the ESIWACE-2 funded supports [Maisonnave & Bourdallé-Badie 2022; Maisonnave & Berthet 2022].

  allow single precision computations (ETHZ), concurrent coupling of ocean and ice (MetOffice) or full ocean zoom coupling (GEOMAR), with a significant performance improvement. During the 2020-2021 period, we modified the call of OASIS API routines in models (ocean-ice and runoff mapper-ocean), to enhance the physical interface (NERSC) or make the most of a new OASIS functionalities, e.g. the locally conservative interpolation (GEOMAR). Computing performance was checked, with more accuracy since the new OASIS event timeline[ Maisonnave et al. 2020 ] was made available. Finally, during the last year, a complex and computationally efficient coupled system, involving the XIOS I/O server, was ported on a non scalar processor (DWD) but the ODUS also helped to validate interpolation choices (at DWD and SMHI), to check the functionning of the newly developed pyOASIS API (SMHI) and to design an interface in a component (wave model) recently made available in the community (Météo-France). More advanced coupling were prepared for a Eulerian version of the neXtSIM sea-ice model (NERSC) and the coarsening of the TOP-MEDUSA biogeo-chemistry model (MetOffice). In the first case, the ORCA global grid was split into two independent hemispheric components, as a simplified basis of the neXtSIM future discretisation, having an interpolation-less coupling with the NEMO ocean in mind. In the second case, an exact coarsening of the NEMO ocean tracers was implemented in the OASIS library to allow performing the future BGC component computations at a lower resolution.

Table 1 :

 1 Granted laboratories, task, coupled system implemented or modified and interaction with coupling library (Enhanced, bug Fixed, Version upgraded on existing model)

	Laboratory, Country	Main tasks			Coupled system components	Coupler
					(Regional or Global)
	ETHZ, CH	Model upgrade, porting		Atmosphere -Land (R)	V
	MetOffice, UK	Interface design, performance	Ocean -Sea Ice (G)	F
	GEOMAR, DE	Interface design, performance	Atmosphere -Ocean (G/R)	V,F
	NERSC, NO	Interface design		Ocean -Sea Ice (R)	none
	GEOMAR, DE	Interpolation setup		Atmosphere-Ocean-Runoff	V
					(G/R)
	NERSC, NO	Model upgrade, interface design	Ocean-Sea Ice (G)	none
	DWD, DE	Porting			Atmosphere -Ocean (R)	V, F
	Météo-France, FR	Interface design, porting	Ocean -Wave (G)	none
	SMHI, SE	Workflow upgrade, interpolation	Atmosphere-Ocean-Runoff (G) E
		setup
	MetOffice, UK	Interface	design,	interpolation	Ocean-Biogeochemistry (G)	E
		development

Table 2 :

 2 List of OASIS compatible models which interface was upgraded or created during the ODUS Models are not the only things which exchanges are facilitated by the modularity. The ODUS program may have speed up the spreading of coupling solutions or best practices that are not currently well known (locally conservative runoff, ocean-ice coupling in NEMO, load balancing measurement …) but we guess that this diffusion is implicit in the communities we visited for the reason that coupling techniques are exchanged at the same time that new components or component versions.

	Ocean	Atmosphere	Others
	NEMO v3.6	COSMO-CLM v5	CLM v5 (land)
	NEMO v4	OpenIFS cy40	SI3 v4 (sea ice)
	NEMO v4.2	OpenIFS cy43r3	neXtSIM (sea ice)
	HyCOM v2.2	ICON-NWP v2.6.4	MF-WAM (waves)
			TOP-PISCES (bio-geo-
			chemistry)

Table 3 :

 3 Communities benefitting of the interfaced components developed during the ODUS

  First because a standardly interfaced code is more easily shared or recycled, which should contribute to save labour force, if correctly organised. But modularity also makes the component choice more accurate, thus less expensive, by tuning (i) the number of component in the system, e.g. from GCM to ESM in [ Séférian et al. 2019 ], (ii) the component complexity, e.g. full atmosphere model or simple boundary layer [ Lemarié et al. 2021 ], or (iii) the component resolution, e.g. by coarsening the resolution of the BGC of the NEMO ocean [ Maisonnave et al. 2021 ]. Metaphors can sometimes help to understand more quickly. So let forgive us the a priori strange comparison of the OASIS community coupler to the solution designed by [ Dietel

One could also mention here, as a possible explanation of the code growing intricacy, our proclivity to add new parametrisations to the model, for the only reason that they are supposed to enhance its realism. In that representation, the removal of any part of the code would take away the model from some kind of truth

The reusability is increased by the backward compatibility of the OASIS interface

See e.g. in our community the growing importance of Climate Change impact and mitigation studies compared to model infrastructure support

Although one of our reviewers questioned the possibility to offer us a ride on one of these Bauhaus inspired motorbike which is rusting in his garage. "Arcs triomphaux, pointes du ciel voisines / qui de vous voir le ciel même étonnez / Las, peu à peu cendre vous devenez / fable du peuple et publiques rapines !"

"Poetics" to be understood in its classical sense of production of works

Acknowledgments

The OASIS dedicated support has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824084 (IS-ENES-3) and No 823988 (ESIWACE-2).

• Balaji, V., Couvreux,F., Deshayes, J., Gautrais, J., Hourdin, F. & Rio., C, 2022: Are general circulation models obsolete? Proceedings of the National Academy of Sciences of the United States of America, 119 47, e2202075119 • Will, A., Akhtar, N., Brauch, J., Breil, M., Davin, E., Ho-Hagemann, H. T. M., Maisonnave, E., Thürkow, M., & Weiher, S., 2017: The COSMO-CLM 4.8 regional climate model coupled to regional ocean, land surface and global earth system models using OASIS3-MCT: description and performance , Geosci. Model Dev., doi:10.5194/gmd-10-1549-2017


Fang Yang

Li Chai
email: chaili@wust.edu.cn

Da Chen
email: chenda@ceremade.dauphine.fr

Laurent Cohen
email: cohen@ceremade.dauphine.fr

Geodesic via Asymmetric Heat Diffusion Based on Finsler Metric

Keywords: Geodesic, Minimal Path, Heat Diffusion, Finsler Metric

Current image segmentation involves strongly non-uniform, anisotropic and asymmetric measures of path length, which challenges available algorithms. In order to meet these challenges, this paper applies the Finsler metric to the geodesic method based on heat diffusion. This metric is non-Riemannian, anisotropic and asymmetric, which helps the heat to flow more on the features of interest. Experiments demonstrate the feasibility of the proposed method. The experimental results show that our algorithm is of strong robustness and effectiveness. The proposed method can be applied to contour detection and tubular structure segmentation in images, such as vessel segmentation in medical images and road extraction in satellite images and so on.

Introduction

Geodesic refers to the shortest path connecting two points in metric space and plays an important role in image processing and computer vision. It can be applied to tasks such as contour detection, tubular structure segmentation, surface remeshing and so on [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF]. In theory, the shortest path problem can be described by a static Hamilton-Jacobi equation or an anisotropic Eikonal partial differential equation (PDE). Numerical methods to compute the discrete geodesic distance on smooth surfaces and images can be classified into two classes: exact methods based on geometry and approximation method based on PDE. Generally, the approximation methods via solving the Eikonal PDE are widely used.

In terms of computation, it is difficult to get the analytical solution of the Eikonal equation directly. Therefore, an iterative relaxation scheme such as Gauss-Seidel is useful to get the numerical solution. In [START_REF] Hysing | The eikonal equation: numerical efficiency vs. algorithmic complexity on quadrilateral grids[END_REF], the authors summarize popular numerical solutions to the Eikonal equation on Cartesian grids. The most popular algorithms are the Fast Marching [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF] and Fast Sweeping Method [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF] etc. By comparing the errors, speed, accuracy and robustness of different Supported by the National Science Foundation of China (grant 61625305).

algorithms, the authors conclude in [START_REF] Hysing | The eikonal equation: numerical efficiency vs. algorithmic complexity on quadrilateral grids[END_REF] that the Fast Marching Method outperforms the other methods. However, the Fast Marching Method do not reuse its information, i.e., for the same data, once the initial point is changed, the geodesic distance should be recomputed from scratch.

Recently, the geodesic method based on heat flow has become popular because of its efficiency, robustness and insensitivity to noise. It is originated from the Varadhan's formula [START_REF] Varadhan | On the behavior of the fundamental solution of the heat equation with variable coefficients[END_REF]:

φ(p 0 , p x ) = lim t→0 -4t log u p0 (p x , t), (1)
which shows that within a small time t, the geodesic distance φ(p 0 , p x ) between two points p 0 and p x can be approximated via the heat kernel u p0 (p x , t).

Inspired by the classical result of Varadhan, Crane et al. proposed a different method for computing the geodesic distance on a Riemannian manifold, which is called the heat method [START_REF] Crane | Geodesics in heat: A new approach to computing distance based on heat flow[END_REF]. Instead of using the Varadhan's formula directly, the heat method is based on solving a pair of standard linear elliptic problems.

The key observation in [START_REF] Crane | Geodesics in heat: A new approach to computing distance based on heat flow[END_REF] is to divide the process of distance computation into two steps: the first step determines the direction along which the distance increases; the second step recovers the distance via solving a Poisson equation by using the normalized gradient of the heat flow. Each step relates to a standard problem in numerical linear algebra. Moreover, the sparse systems from elliptic PDEs can be solved in a very short time that is close to linear time. Therefore, the heat method can be facilitated by using existing algorithms and software, thus improves the efficiency and robustness of distance computation. In addition, for every single data, the Laplacian can be precomputed, i.e. recomputation is not necessary even the initial point is changed.

Though fast and efficient, the heat method described in [START_REF] Crane | Geodesics in heat: A new approach to computing distance based on heat flow[END_REF] only involves the simplest case -the diffusion coefficient is a constant. To expand the scope of use of the heat method, Yang and Cohen [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF] proposed to use the isotropic/anisotropic heat diffusion to obtain geodesic in images and on surfaces. In their method, the diffusion coefficient is no longer a constant, but an isotropic scalar or an anisotropic tensor computed from the image. Based on [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF], the authors propose a 2D+Radius heat method for segmenting the centerline and boundary of the tubular structures simultaneously [START_REF] Yang | Tubular structure segmentation based on heat diffusion[END_REF]. Besides, the diffusive nature of the heat equation causes instant smoothing, which is the reason why the heat method is robust to noise.

Because of the advantages of the heat method, based on Crane's work, a variety of tasks in computational science and data analysis are developed. For example, Belyaev and Fayolle [START_REF] Belyaev | On variational and pde-based distance function approximations[END_REF] reinterpret the heat method and achieve more accurate results by either iterating the heat method or by applying more advanced descent strategies. Zou et al [START_REF] Zou | Iso-level tool path planning for free-form surfaces[END_REF] use the heat method for efficient tool path planning; Solomon et al [START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF] use the heat method to facilitate solving optimal transport problems on geometric domains; Lin et al [START_REF] Lin | Geodesic distance function learning via heat flow on vector fields[END_REF] apply this approach to vector-valued data in the context of manifold learning.

The heat method is popular, however, the diffusion coefficient and tensor are constructed on Riemannian manifold, which can not provide asymmetric met-ric. While current image processing applications involve strongly non-uniform, anisotropic and asymmetric measures of path length, which can not be efficiently solved on Riemannian manifolds. Hence, in this paper, we are motivated to apply asymmetric Finsler metric to the heat method.

Finsler geometry can be described as Riemannian geometry without the quadratic restriction [START_REF] Chern | Finsler geometry is just riemannian geometry without the quadratic equation[END_REF]. Riemannian distances are defined by a position dependent inner product, while Finslerian distances are defined by a direction dependent inner product, computed from a position dependent norm. Since the level sets determined by norms can be more complex than those (ellipsoidal) given by an inner product, Finsler metric is more suitable for modeling complex diffusion profiles [START_REF] Astola | Finsler geometry on higher order tensor fields and applications to high angular resolution diffusion imaging[END_REF].

In [START_REF] Melonakos | Finsler active contours[END_REF], Melonakos et al.proposed a Finsler active contour method by replacing the original potential in [START_REF] Caselles | Geodesic active contours[END_REF] into a potential including the normal of the curves. Then Zach et al [START_REF] Zach | Globally optimal finsler active contours[END_REF] proposed a global optimal method for the Finsler active contours. Mirebeau [START_REF] Mirebeau | Efficient fast marching with finsler metrics[END_REF] proposes an efficient discretization scheme for the Fast Marching Method on Finsler manifold. Chen et al [START_REF] Chen | Global minimum for a finsler elastica minimal path approach[END_REF] introduces a minimal path model with Finsler metric and uses the scheme in [START_REF] Mirebeau | Efficient fast marching with finsler metrics[END_REF] to compute the distance. In [START_REF] Chen | Fast asymmetric fronts propagation for image segmentation[END_REF], the authors extend the framework of front propagation from Riemannian manifold to Finsler manifold. The Rander's metric used in their work [START_REF] Chen | Fast asymmetric fronts propagation for image segmentation[END_REF] prevents the fronts leaking problem during the fronts propagation.

The asymmetric metric is successfully applied in the Eikonal PDE, but as mentioned before, the numerical scheme Fast Marching Method is more sensitive to noise than the heat method. In addition, the Fast Marching Method do not reuse information. Therefore, we would like to develop the asymmetric heat method.

The contributions of this paper lie in that:

1. We extend the heat method based on Riemannian metric to Finsler metric; 2. We build two kinds of Rander's metric to detect contours and centerlines; 3. The asymmetric heat diffusion improves the performance comparing with the traditional heat method on image segmentation.

2 Background and Related Work

Finsler metric and minimal path model

Similar to the geodesic on Riemannian manifold, the minimal path on Finsler manifold is also obtained by minimizing the length between two fixed points.

Let Ω ⊂ R N , N ∈ {2, 3} denote the image domain, which is equipped with a Finsler metric F(x, v) > 0, where x ∈ Ω denotes the position and v ∈ R N denotes the orientation. Generally, the Finsler metric

F(x, v) is asymmetric, i.e., F(x, v) = F(x, -v).
The length of a Lipschitz continuous curve γ(t), 0 ≤ t ≤ 1 on a Finsler manifold is defined as:

L F (γ) = 1 0 F(γ(t), γ (t))dt. ( 2
)
Note that because of the asymmetry of F(x, v), the curve on a Finsler manifold and its reverse curve are usually not the same.

The geodesic distance defined in (2) can be generalized to the distance from any point p ∈ Ω to a set of points S ⊂ Ω, thus defines the geodesic distance map:

φ S (p) := inf{L F (γ); γ(0) = p, γ(1) ∈ S}, (3)
which is the unique viscosity solution to the following Eikonal PDE:

F * (∇φ(p)) = 1, ∀p ∈ Ω \ S, φ(p) = 0, ∀p ∈ S, (4)
where F * is the dual metric of F:

F * (x, u) = sup v =0 u, v F(x, v) . ( 5
)
Once the geodesic distance map φ is obtained, the shortest path from p to S can be tracked by solving an ordinary differential equation (ODE) [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using voronoi's first reduction of quadratic forms[END_REF]:

dγ(t) dt := ∇F * (∇φ(γ(t))). (6)
In this paper, we will utilize the Rander's Finsler metric, which is defined by the combination of a symmetric Riemannian metric and an asymmetric linear part:

F(x, v) = v M (x) + ω(x), v . (7)
According to [START_REF] Mirebeau | Efficient fast marching with finsler metrics[END_REF], the dual to F is also a Rander's metric F * with parameters M * and ω * :

F * (x, v) := v M * (x) + ω * (x), v , (8)
where M * and ω * are some algebraic expressions of M and ω [START_REF] Mirebeau | Efficient fast marching with finsler metrics[END_REF]:

η := 1 -ω, M -1 ω , M * := (M -1 ω)(M -1 ω) T + ηM -1 η 2 , ω * := - M -1 * M -1 ω η .
Hence, the backtracking ODE on the Rander's metric becomes [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using voronoi's first reduction of quadratic forms[END_REF]:

dγ(t) dt = M * (γ(t))∇φ(γ(t)) ∇φ(γ(t)) M * (γ(t)) + ω * (γ(t)). (9)
Compared with the Riemannian minimal path model [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF], the Rander's Finsler minimal path is with an additional linear part, which is able to deal with the asymmetric information of the images.

Finsler Heat Equation

In this paper, we are dedicated to developping the heat method on the Rander's Finsler manifold. Let us consider the following energy functional:

E(u) := Ω ∇u + g 2 D , ( 10
)
where D is the conductivity matrix and is positive definite by construction:

D = M * -ω * ⊗ ω * , (11)
and we set g = D -1 ω * .

The gradient descent of ( 10) is interpreted as:

∂ t u = ∆ F u, (12)
where ∆ F is the Finsler Laplace operator [START_REF] Akagi | The cauchy problem for the finsler heat equation[END_REF][START_REF] Ohta | Heat flow on finsler manifolds[END_REF], and ( 12) is called the Finsler heat equation.

Let us recall the heat method proposed in [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF], where the initial condition for the heat equation is u p0 (p x , t = 0) = δ p0 . Therefore, the initial condition for the Finsler heat method should be the same as in [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF]. The Varadhan's formula (1) presents that within a small time t → 0, the heat kernel can be used to approximate geodesic distance. It is proved in [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF] that the Varadhan's formula is applicable in Riemannian cases because that the geodesic distance is only related to the coefficients of the second-order derivative of the heat kernel.

In [START_REF] Ohta | Heat flow on finsler manifolds[END_REF], the Finsler Gauss kernel G F * is presented by the following exponential form:

G F * (x, t) := (4πt) -N 2 exp (- F * (x) 2 4t ). ( 13
)
Let t → 0, put the Finsler heat kernel (13) on the r.h.s of the Varadhan's formula [START_REF] Akagi | The cauchy problem for the finsler heat equation[END_REF], according to the L'Hospital's rule, we have:

lim t→0 -4t log(G F * ) = F * (x). ( 14
)
Equation ( 14) is already well understood in the case of homogeneous heat diffusion on the whole domain R 2 , where F * (x) = x and the corresponding distance is Euclidean.

Construction of Rander's Metric and Algorithm of the Asymmetric Heat Method

As mentioned in section 2.1, the Rander's metric is composed by a Riemannian symmetric quadratic form plus a linear part. In this section, the construction of the Rander's metric is presented. In addition, we should pay attention that in order to ensure that the metric F is always positive, the following constraint should be satisfied [START_REF] Chen | Global minimum for a finsler elastica minimal path approach[END_REF][START_REF] Chen | Fast asymmetric fronts propagation for image segmentation[END_REF]:

ω(x), M -1 (x)ω(x) < 1, ∀x ∈ Ω. ( 15
)
Different specific metrics are designed for detecting the boundaries of shapes and centerlines of tubular structures in images.

Rander's Metric for Boundary detection

In [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF] and [START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF], the authors use either the magnitude of the gradient as the diffusion coefficient or the structure tensor to enhance the flow-like structures or edge features. Inspired by their work, we adapt these metrics to our method for detecting boundaries.

Denote by g(x) the gradient of the image: g(x) = ∇f (x), where f (•) is a scalar function computed from the image I. The definite positive symmetric tensor field M (x) can be built by the magnitude of g(x):

M I (x) := ( g(x) + )I d , (16)
where M I (•) represents the tensor field in isotropic case, is a small constant that prevents the diffusion coefficient from being 0, and I d is a 2 × 2 identity matrix. M (x) can be also built by the tensor product:

M A (x) = λ 1 (x)e 1 (x)e T 1 (x) + λ 2 (x)e 2 (x)e T 2 (x), (17)
where M A (•) denotes the tensor field in anisotropic case, e 1 (x) = g(x) and e 2 (x) = g(x) ⊥ are the eigenvectors of M A (x), λ 1 (x) ≥ λ 2 (x) ≥ 0 are the corresponding eigenvalues. The anisotropy A is defined as:

A(x) = λ 1 (x) -λ 2 (x) λ 1 (x) + λ 2 (x) . (18)
The linear part can be obtained by using the orthogonal of the normalized gradient of the image:

ω(x) = τ ∇I(x) ∇I ∞ ⊥ , ( 19
)
where τ is an adjustable coefficient which is used to satisfy [START_REF] Melonakos | Finsler active contours[END_REF]. Note that in [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF], M A behaves better than M I when the feature of interest is very curved. However, in the proposed method, the difference between M A and M I is weakened because of ω(x).

Figure .1 shows the experimental results on a synthetic image under the same initial condition but different metrics. In (a), only ( 16) is used as the diffusion coefficient, without adding any additional part, which leads to a shortcut during backtracking from the endpoint to the source point. In (b) and (c), since ω(x) and -ω(x) defined in [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF] are used as the linear part in the Rander's metric respectively, the heat flows more on the edge of the central black curve. Hence, during the backtracking process, the extracted curves are along the edge of the shape. Note that the closed boundary of the black central curve can be obtained by simply merging the results of (b) and (c).

Figure .1 illustrates clearly that the additional linear part ω(x) forces the heat to flow more along the direction of ω(x), i.e., ω(x) enhance the heat diffusion on the edges. 16), (b) by using Rander's metric whose linear part ω is defined in ( 19) and (c) by using the Rander's metric, the linear part -ω. In the bottom row, the red and blue points denote the source point and endpoint respectively, the green curves are the geodesic curves.

Figure.

2 presents an example of using different metrics on a 2D Gaussian image. First we show the Euclidean distance with the anisotropy A = 0 and its corresponding geodesic curve, which is in fact a straight line (a). In (b), we use the diffusion tensor defined in [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using voronoi's first reduction of quadratic forms[END_REF] as the symmetric part M A (x) with the anisotropy A ≡ 0.9 everywhere. We can see that the heat flows along the direction of e 2 (x) and the geodesic line is curved because of the anisotropic diffusion. In (c), we add the linear part ω(x). Though the geodesic curves in (b) and (c) make no difference visually, it is not hard to distinguish the distance maps from (b) and (c), which are symmetric and asymmetric respectively. In (b) the distance map is symmetric because the diffusion tensor is with quadratic form, while in (c), because of ω(x), the heat diffusion is becomes asymmetric hence flows more on the direction of ω(x). In (d), we use -ω(x) as the linear part, which is opposite to ω(x), so the heat flows in the inverse direction of (c).

The results in Fig. 1 and 2 not only demonstrate the validity of the metrics designed for boundary detection, but also testify the effectiveness of the Rander's metric. In addition, both examples illustrate that F(x, v) = F(x, -v).

Rander's Metric for Centerline extraction

For the purpose of extracting the centerline, the vessel enhancement method can be considered, e.g. Hessian-based vesselness measures [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] Here we plan to construct the symmetric part of the Rander's metric based on the eigenvalues 17), (c) by using the Rander's metric with the linear part ω defined in ( 19) and (d) by using the Rander's metric with linear part -ω. In the bottom row, the red and blue points denote the source point and endpoint respectively, the blue curves are the geodesic curves.

and eigenvectors of the Hessian matrix because of its simplicity. In [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], the authors introduced a multi-scale vessel enhancement method by interpreting geometrically the eigenvalues of the Hessian matrix. The local orientation of the image can be estimated by using the Hessian eigenvectors, which allows us to find out the position of the tubular structures. First we convolve the image I with a Gaussian kernel G(x, σ):

G(x, σ) = 1 (2π)σ exp(- x 2 2σ 2 ) ( 20
)
Since σ can be considered as an estimator of the width of the tubular structures, we use different σ in different images. Let I G be the product after convolution, we obtain the symmetric Hessian matrix H by computing the second derivative of I G :

H = ∂ xx ∂ xy ∂ yx ∂ yy I G (21)
Then we decompose H to get the corresponding eigenvectors v 1 (x), v 2 (x) and eigenvalues ϑ 1 (x), ϑ 2 (x) in order to form M (x):

M = ϑ 1 (x)v 1 (x)v T 1 (x) + βϑ 2 (x)v 2 (x)v T 2 (x) ( 22
)
where β is a constant controlled by the users to adjust the anisotropy A.

In order to prevent the heat from leaking out the region of tubular structure, we can use the eigenvector v 1 (x) as the linear part of the Rander's metric:

ω = τ v 1 ( 23
)
where τ is an adjustable constant which plays the same role as in [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF] to ensure that ( 15) is satisfied.

Figure .3 shows an example on a synthetic U-tube image. In this example, we set σ in (20) to be 4, β in ( 22) to be 0.1 and τ = 1. In (a), we use the isotropic heat diffusion. The diffusion coefficient is simply obtained by using the graylevel of the image: α = (1 -(I -I x0 )) n , n = 3 in this case, and M I = αI d . Clearly, there is a shortcut connecting the source point and the endpoint. In (b), we use the anisotropic heat diffusion where M is the diffusion tensor. The geodesic curve goes along the U tube structure this time. However, our aim is to extract the centerline of the structure. In (c), we use the Rander's metric in heat diffusion by adding ω [START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF]. From the result, we can see that the asymmetric heat diffusion forces the heat to flow more on the main direction v 1 , In the bottom row of Fig. 3, although that both curves in (b) and (c) travel along the U tube structure, we can still find that the result by using asymmetric diffusion is closer to the centerline.

Algorithm of the Asymmetric Heat Method

The basic steps of the asymmetric heat method can be described as follows:

In our method, the numerical solution to the Finsler heat equation and the time of diffusion (Step.2) are determined according to a backward discretization scheme designed for anisotropic heat diffusion, see [START_REF] Fehrenbach | Sparse non-negative stencils for anisotropic diffusion[END_REF] for details.

Algorithm 1 The Asymmetric Heat Method

1. Construct the asymmetric metric F(x, u), Sect.3.1 and Sect.3.2 ; 2. Integrate the heat flow ∂tu = ∆F u from some source points for some fixed time t; 3. Approximate the distance φ using Varadhan's formula; 4. Backtrack the geodesic from a given endpoint to the source point using an ODE [START_REF] Chern | Finsler geometry is just riemannian geometry without the quadratic equation[END_REF].

Experiments and Discussions

Experiment data and settings

We testify the effectiveness of the proposed method on both synthetic and real images. Firstly, we use a 411×411 synthetic image to testify the noise-insensitive performance of the heat method, see Fig. 4. A 12% of pepper and salt noise is added to the original image. The metric is constructed by using ( 16) and ( 19), with τ = 1.

For detecting the boundary of shapes in images, we choose some natural images from the BSDS500 dataset [START_REF] Arbeláez | Contour detection and hierarchical image segmentation[END_REF], see Fig. 5. The gray images are used to construct metrics and initiate the heat flow.

For extracting the centerline of tubular structure, we use real road and vessel images, see Fig. 6. We compare the result obtained by asymmetric diffusion with isotropic and anisotropic diffusion in [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF]. For the isotropic diffusion, we preprocess the images with a sigmoid function:

f (x) = 1 - 1 1 + exp(k(I -α)) , (24)
where k = 10 and α = 0.4 in our tests, and f (x) is used to compute the diffusion coefficient M I in [START_REF] Mirebeau | Efficient fast marching with finsler metrics[END_REF]. For anisotropic and asymmetric diffusion, the quadratic form M A is constructed by using [START_REF] Varadhan | On the behavior of the fundamental solution of the heat equation with variable coefficients[END_REF] with a fixed σ = 8 and β = 0.1.

To display the results clearly, we use either green or blue curves to represent the extracted geodesic curve. The source points and endpoints for each image are given by the users. The red point denotes the source point and the blue ones are the endpoints.

Finally, we evaluate the quality of our algorithm on the images of BSDS500 dataset. Table .4.2 lists the recall of the detection results by symmetric diffusion method [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF] and the proposed method. Note that instead of all the edges, only the geodesic curves between the two chosen points are used to calculate the recall.

Results and analysis (1) Noisy case

In this test, the original synthetic image is corrupted by some pepper and salt noise. As we know that these kind of noise is impulse noise, it affects the image gradient a lot. However, as we can see in Fig. 4(a), most part of the extracted geodesic still goes along the edge of the center curve. Moreover, by adding the additional linear part ω and -ω [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF], both geodesics travel exactly on the boundary of the curve and complete the contour of the center shape. For the result by the fast marching method (d), we use the same metric as in (b), but the result is not as good as in (b) or (c). This test demonstrates that the heat method is robust to noise, since the exponential kernel (13) causes instant smoothing. (

2) Contour Detection

The boundary and contour detection results are displayed in Fig. 5. The comparison experiments are performed by using symmetric heat method in [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF] with isotropic metric M I or anisotropic metric M A , as shown in Column.1 and Column.3 in Fig. 5. The results of the proposed method are shown in Column.2 and Column.4. By using symmetric diffusion [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF], the detected edges do not include some main details, e.g., some elephants and the tree of elephant, the long pole, haystack and head of lakeman. As mentioned before, because of the diffusive nature of the heat kernel, the heat method can smooth the sharp changes in images, i.e., noise, spurious parts. Hence in such cases, the spurious part of the edges can not be fully detected by using the symmetric heat diffusion. However, after adding the linear part ω, all the missing details are included. Apparently, the shapes of all images are irregular and asymmetric. Because of the symmetry of the Riemannian metric, the geodesic connecting two points of using Riemmanian metric will surely be the shorter part on the contour of the shape, such as mushroom, redbird and so on. However, by adding the asymmetric linear part ω, the detected curves can go in the opposite way and thus complete the whole contour of the shape, such as hawk, swan, mushroom and birds etc. In addition, for river, waterfall and bear, the background is complex, therefore, the symmetric method fails to extract the edges of these shapes directly and causes some shortcuts. When adding the linear part ω, the curves are forced to go along the edges of the shapes, even the details can be detected.

The reason why asymmetric diffusion using Rander's metric works better is not only because that during the process of diffusion, the heat are more concentrated on the edges, but also on the way of backtracking, the vector field ω in ( 9) is taken into account. However, it is worth mentioning that the proposed method to detect the boundaries are somehow sensitive to the gradient of the image, e.g., the beak of the swan in Row.2, the contour of the bird in Row.5.

The quantitative results of contour detection in Fig. 5 (from left to right, top to bottom) is shown in Table .4.2. Obviously, in most cases, the recall of detection results using asymmetric diffusion are higher than the symmetric ones in [START_REF] Yang | Geodesic distance and curves through isotropic and anisotropic heat equations on images and surfaces[END_REF]. (3) Centerline Detection

The results of detecting centerlines in tubular structures are shown in Fig. 6. For the road image in Row.1, the curves extracted using different metrics are all smooth, despite that there is a lot of noise around the road. This example demonstrates again that the heat method is not sensitive to noise. The Column.1 of Fig. 6 displays the results by using isotropic heat diffusion. As we can see, the extracted lines are not on the centerline, especially for the curved part, there is much deviation from the center. This case can be improved by using the anisotropic diffusion, which helps the heat to concentrate inside the tubular structures. However, the centerlines of the very curved parts and junctions are still not very well detected (Column.2). Finally, we combine the quadratic form with the linear part ω (23) to construct the Rander's metric. The results by using the asymmetric diffusion are shown in Column.3. Thanks to ω the extracted geodesic curves are located much more on the centerlines visually.

Conclusion

In this paper, we extend the traditional symmetric heat method to asymmetric cases in order to overcome the difficulty on measuring asymmetric path length. Two kinds of Rander's metric are designed for contour detection and centerline extraction. Experimental results show that asymmetric heat diffusion is effective and robust. Thanks to the asymmetric metric, the proposed method performs better than the traditional heat method. However we only consider the Rander's metric in this paper. In future, we will design different Finsler metrics to apply to the heat method according to different image processing tasks.

Fig. 1 .

 1 Fig. 1. From top to bottom are the distance maps and their corresponding geodesic respectively, from left to right are the results: (a) by using Riemannian metric M defined in (16), (b) by using Rander's metric whose linear part ω is defined in (19) and (c) by using the Rander's metric, the linear part -ω. In the bottom row, the red and blue points denote the source point and endpoint respectively, the green curves are the geodesic curves.

Fig. 2 .

 2 Fig. 2. From top to bottom are the distance maps and their corresponding geodesic respectively, from left to right are the results: (a) Euclidean, (b) using anisotropic Riemannian metric (17), (c) by using the Rander's metric with the linear part ω defined in (19) and (d) by using the Rander's metric with linear part -ω. In the bottom row, the red and blue points denote the source point and endpoint respectively, the blue curves are the geodesic curves.

Fig. 3 .

 3 Fig.3. From top to bottom are the distance maps and their corresponding detected geodesic curves respectively, from left to right are the results by using: (a) isotropic heat diffusion, (b) anisotropic heat diffusion with metric (17), (c) asymmetric heat diffusion with the linear part ω defined in[START_REF] Weickert | Coherence-Enhancing Diffusion Filtering[END_REF].

Fig. 4 .

 4 Fig. 4. Experiment on a synthetic noisy image, from left to right are the results by using: (a) isotropic heat diffusion, (b) asymmetric heat diffusion with linear part ω defined in (19), (c) asymmetric heat diffusion with the linear part -ω, (d) asymmetric fast marching method (the same metric as used in (b)) .

Fig. 5 .

 5 Fig. 5. Experiment on several real natural images (from left to right and top to bottom): elephant, river, hawk, swan, lakeman, mushroom, redbird, waterfall, bear and bird. The red points are the source points and the blue ones are the endpoints, the blue or green curves are the detected geodesic curves.

Fig. 6 .

 6 Fig. 6. Experiment on real satellite road image and vessel images, from left to right are the results by using: (a) isotropic heat diffusion, (b) anisotropic heat diffusion, (c) asymmetric heat diffusion.

Table 1 .

 1 Quantitative results of symmetric and asymmetric heat methods

	Image	elephant river	hawk swan lakeman
	Symmetric recall	0.365	0.214	0.863 0.628 0.653
	Asymmetric recall 0.892	0.878	0.827 0.852 0.936
	Image	mushroom redbird waterfall bear bird
	Symmetric recall	0.637	0.774	0.033 0.797 0.913
	Asymmetric recall 0.921	0.836	0.866 0.453 0.824

Acknowledgements

The authors would like to thank Dr. Jean-Marie Mirebeau for his insightful suggestions on the numerical solutions to asymmetric heat diffusion. The authors would also like to thank Dr. Xin Su for his useful comments that allowed us to improve this paper.