Resilient greedy routing on GPS-free surface sensor networks

We design a greedy routing scheme specifically for GPS-free large-scale wireless sensor networks deployed on surfaces of complex-connected three-dimensional settings. Compared with other greedy embedding–based surface network routing scheme, the proposed one is cut free such that no pair of nodes suffers a long detour to reach each other. The routing scheme is designed to be resilient to node or link failures especially under random node or link failure model where each node in a network has an equal and independent probability of failure during some time interval. The proposed algorithm is fully distributed and scalable to both the size and the topological complexity of a network. Each sensor node requires only limited and constant storage. Simulation results show the proposed routing scheme with a higher successful delivery ratio, a lower average stretch factor, and a lower normalized communication cost compared with other resilient routing methods.


Introduction
A wireless sensor network may be deployed on a twodimensional (2D) plane 1 (e.g. for crop sensing in fields or wildlife tracking on plains), or in a three-dimensional (3D) volume 2-4 (for underwater or space reconnaissance), or on a 3D surface 5,6 (e.g. with sensors mounted on ocean floor, mountain surface, or the surface of various man-made structures). Routing has always been a fundamental challenge in large-scale wireless sensor networks. A diversity of sensor network routing algorithms have been developed, but they primarily target at 2D planar or 3D volumetric settings. [7][8][9][10][11][12][13][14][15][16][17][18] This article focuses on designing a routing algorithm resilient to random node or link failures in a large-scale wireless sensor network deployed on the surface of a complexconnected 3D setting. Such 3D setting generally has multiple handles and forms a high genus surface. 19 Applications for such wireless sensor networks on high genus surfaces have been demonstrated for structural health monitoring (SHM), for example, at the Golden Gate Bridge 20 and National Stadium of China, 21 where sensors are deployed on megastructures to gauge changes in materials or geometric properties that could hinder the system's performance. Other applications can be found along the corridors of buildings for fire detection, and in water, sewer, or gas system for monitoring underground pipelines as introduced in Yu et al. 22 We consider densely deployed sensors that operate on high radio frequency and extremely low transmission power with short communication range. Only nearby sensors along a 3D surface can communicate with each other, whereas the wireless links connecting remote sensors across space are negligible and can be removed via a simple preprocessing. As a result, the dense sensors form a 3D surface network.

An overview of resilient routing
Node and link failures are unavoidable in a large-scale distributed sensor system. For example, nodes may die out of battery and communication links can be temporarily or permanently disabled due to attacks or interference. We can model network failures with two failure models: the independent node failure model and the geographically correlated failure model. In the independent node failure model, each node in a network has an equal and independent probability of failure during some time interval. 23 Independent node failure can represent failure due to energy dissipation or localized environmental effects. In the geographically correlated failure model, all nodes within a certain range or hop count fail. 23 Geographically correlated failure may happen due to natural environmental effects (such as rain fades) or some human activity within a geographic region.
Resilient routing has been discussed in conventional 2D sensor networks. For example, Directed Diffusion 24 is a well-known routing algorithm that establishes gradients and uses gradual reinforcement of better paths to allow path recovery, enabling a system to be robust to a local network and sink dynamics. Multipath routing represents an efficient routing strategy to increase network resilience to node and link failures. Algorithms have been proposed for Internet multipath routing [25][26][27][28][29] and sensor network multipath routing. 23,30 In Ganesan et al., 23 the authors consider two different approaches to construct multipaths between two nodes. One is node-disjoint multipath, where the alternate paths do not intersect with each other. The other approach builds many braided paths-partially disjoint alternate paths. In Lou, 30 a branch-aware flooding scheme is utilized to construct a spanning tree and discover a set of node-disjoint paths for each sensor node back to the base station. Later, a tree structure is applied for multipath routing. 31,32 One recent advance in multipath routing is homotopy-based routing. 33 The authors apply the concept of homotopy to effectively evaluate the diversity of alternative paths of a network deployed on 2D plane with uncovered holes. The simplest case to illustrate the concept of paths with different homotopy types is a network deployed on a planar region with a lake inside. Two paths connecting a pair of source and destination are homotopic to each other if they are on the same side of the lake; otherwise, they belong to different homotopy types because one path cannot continuously deform to the other without crossing the lake. The most recent work in multipath routing is introduced in Liu et al. 34 with a security disjoint routing-based verified message (SDRVM) scheme proposed to improve the network performance with comprehensive consideration of universal resistance to attacks, energy efficiency, and transmission delay. The main strategy is to divide nodes into several sets based on their remaining energy and assign different duties to each set.
Multipath routing schemes are effective for geographically correlated failures. Alternative paths, especially those disjoint with the failed one, has a high probability to avoid a network failure region. Homotopy-based multipath routing schemes choose alternative paths topologically different from the failed one. These homotopy different paths provide a high probability of routing success even under massive geographically correlated failures. However, under the independent node failure model, the performance of multipath routing schemes decrease dramatically with a slight increase in node failure probability. The reason is that an alternate path is always longer than the original one for any multipath method, so the failure probability of the alternate path is always higher than the original one considering each node has an equal probability of failure. It is obvious that multipath routing schemes, especially those choosing alternative paths based on disjoint or different homotopy types, have intrinsically a very low resilience to independent node failures. If the probability of node failure is not low, switching to an alternate path does not help.
Instead of choosing an alternative path, a local detour is more efficient under the independent node or link failure model. We can consider dead nodes or regions with broken links as holes. To avoid torrents of triggered global updates of network topology, a practical strategy is to apply Face routing and its alternatives 7-14 to locally detour packages around the boundary of each hole. Specifically, a node always forwards a packet to one of its neighbors closest to the destination of the packet. When greedy forwarding of a packet fails due to holes, the packet can be forwarded along the boundary of holes in either clockwise or counterclockwise direction until greedy forwarding is achievable.
However, Face routing and its alternatives require all sensor nodes to know their own positions. For networks deployed on indoor, underground or underwater surfaces, GPS signal is forbidden. They also require networks to be planar. For networks deployed on surfaces of complex-connected 3D settings, when greedy forwarding of a packet is stuck at a local minimum including a hole, as shown in Figure 1(a), there is no deterministic path to successfully forward the packet to the destination based on local information only. We propose a fully distributed algorithm to apply such local detour strategy to achieve efficient resilient routing for GPS-free sensor networks deployed on surfaces of complex-connected 3D structures. Routing is based on local information only. Network dynamic information only needs to be updated locally.

Our approach
A few greedy routing methods 19,22 have been developed recently to address the GPS-free surface sensor network routing problem. In Yu et al., 19,22 the authors propose to cut a closed high genus surface network open and then map the cut-open one to a planar rectangle with two pairs of boundaries. Furthermore, the planar rectangle can be mapped to a planar annulus with one pair of boundaries aligned and glued together. Each sensor node stores a set of planar coordinates and uses it to enable greedy routing. However, these methods suffer a large stretch factor in average. The averaged stretch factor is defined as an average of the ratios of the hop counts of greedy forwarding paths to the shortest ones. When the original surface network is cut open and mapped to either a planar disk, rectangle, or annulus, nodes located at the two sides of the cut-open boundaries need to take a long detour to reach each other instead of directly crossing the boundary. Figure 1(b) shows an example of a pair of such nodes.
The problem of a large stretch factor is intrinsic to the methods proposed in Yu et al. 19,22 They need to cut the originally closed surface network open and then map the network to plane. We propose a cut-free embedding of the original network in the plane. Locally, any region of the network is mapped to plane one-to-one and continuously. The key of the seamless embedding is to extract a triangular mesh structure from the connectivity graph of the original network and compute a set of special edge lengths for the mesh. Based on the computed edge lengths, the original surface network can be embedded in the plane seamlessly. A packet is greedily forwarded based on the planar coordinates stored at each node. Later, the network dynamically maintains the mesh structure. Each node periodically updates its neighboring nodes and local mesh structure due to a potential node or link failure. A face of the mesh is always an n-polygon with n ø 3. We consider all faces with n . 3 holes. When a packet is routed to hit a hole, the packet is forwarded along the boundary of the hole in either counterclockwise or clockwise direction until greedy forwarding is achievable.
The simulation shows the proposed routing scheme with a higher successful delivery ratio, a lower average stretch factor, and a lower normalized communication cost compared with other resilient routing methods. Such results indicate the proposed routing scheme has a strong resilience to network failures, a low cost to detour, and long durability of the network system.
Our contributions. Our contributions can be summarized as follows: We design a greedy routing scheme specifically for GPS-free wireless sensor networks deployed on surfaces of complex-connected 3D settings. Compared with other greedy embedding-based surface network routing scheme, the proposed approach is cut free such that no pair of nodes suffers a long detour to reach each other. The routing scheme is designed to be resilient to node or link failures especially under random node or link failure model. Simulation results show a high successful delivery ratio, a low average stretch factor, and a low normalized communication cost. The proposed algorithm is fully distributed and scalable to both the size and the topological complexity of a network.
The rest of this article is organized as follows: section ''Uniformization metric'' gives a brief introduction of uniformization metric-the special set of edge lengths we compute for a surface network. Section ''Algorithms'' elaborates the proposed resilient routing algorithms. Section ''Simulations'' presents simulation and comparison results. Section ''Conclusion'' concludes the article.

Uniformization metric
In this section, we introduce the uniformization metric, which is essentially a special set of edge lengths we compute to support the proposed resilient routing algorithm. Denote S a closed surface embedded in R 3 and L a loop on S. L is surface separating if it can be expressed as the symmetric difference of boundaries of topological disks embedded in a surface as shown in Figure 2 with L 1 and L 2 ; otherwise, it is non-separating as shown in Figure 2 with L 3 .
The genus of S, equivalent to the number of handles of S, is the maximum number of disjoint non-separating loops L 1 , L 2 , . . . , L g in M; that is, any L i and L j have no topological intersection if i 6 ¼ j, and Mn(L 1 [ Á Á Á L g ) is connected. The genus number is the most basic topology information of a surface. For example, a torus is a genus one surface, and a double torus as shown in Figure 2 is a genus two surface.
A special set of edge lengths to embed S to plane without cutting is called uniformization metric. Uniformization metric exists for general closed surfaces according to Uniformization Theorem. 35 Informally speaking, let S with Riemannian metric induced from the Euclidean metric of S. We denote the Riemannian metric of S as g. There exists a unique metric represented as g = e 2u g, where u : S ! R is a scalar function defined on S. We call g a uniformization metric of S. It can be verified that g is also a Riemannian metric on S and g induces constant Gaussian curvature on S.
Specifically, based on surface uniformization metric, S with genus one can be periodically embedded into Euclidean plane and S with genus larger than one can be periodically embedded into the hyperbolic plane. If S is a discrete surface, for example, a triangular mesh, we have the following discrete definitions of Riemannian metric and Gaussian curvature.
Given a triangular mesh denoted as M, we denote a vertex set as V , an edge set as E, and a face set as F. e ij represents the edge connecting vertices v i and v j , and f ijk denotes the face formed by v i , v j , and v k . The edge lengths of a triangular mesh are sufficient to define a discrete Riemannian metric on the discrete surface l : E ! R + as long as, for each face f ijk , the edge lengths satisfy the triangle inequality: l ij + l jk . l ki .
The discrete Gaussian curvature denoted as K i on a vertex v i 2 V can be computed from the angle deficit where u jk i represents the corner angle attached to vertex v i in the face f ijk , and ∂M represents the boundary of the mesh. Since u jk i can be computed directly from edge lengths of f ijk , it is obvious that discrete Gaussian curvatures are determined by the discrete metrics.
Uniformization metric also exists on a discrete surface. It is a special set of edge lengths such that the discrete surface can be embedded in the plane based on the edge lengths with vertex Gaussian curvatures being zero everywhere. Discrete surface Ricci flow is an efficient tool to compute discrete surface uniformization metric. The algorithm of discrete surface Ricci flow is fully distributed, so it is well suited for wireless sensor network applications. We refer readers to Jin et al. 36 for details of discrete surface Ricci flow.

Algorithms
This section describes the proposed routing scheme in three steps. Section ''Preprocessing'' gives the preprocessing step. Section ''Computing virtual planar coordinates'' introduces the algorithm to compute virtual planar coordinates for each sensor node of a giving network. Based on the computed virtual planar coordinates, resilient greedy routing is introduced in section ''Resilient greedy routing.'' Later, we analyze the time complexity and communication cost in section ''Time complexity and communication cost.'' Preprocessing Given a wireless sensor network deployed on a high genus surface, we apply the algorithm proposed in Zhou et al. 37 to extract a triangular mesh denoted as M from the connectivity graph of the network based on locally measured distances between nodes within onehop communication range. Vertices of the triangular mesh are the set of sensor nodes. An edge between two neighboring vertices indicates the communication link between the two sensors. The process is local and has no constraint on communication models.
For densely deployed sensor nodes, we assume each face of the extracted mesh an n-polygon with n = 3. However, we allow the mesh with non-triangular faces, n-polygons with n.3, due to the randomness of sensor deployment. We consider these non-triangular faces as holes. We will show later that these holes won't affect the computation of uniformization metric and resilient greedy routing.
The topology of the surface of a 3D setting, that is, the genus number, could be available before the deployment of sensor nodes. Otherwise, we can still have the following simple algorithm to automatically detect the surface network topology.
The algorithm starts from one randomly chosen vertex v i of M. v i can be the one with the smallest node ID. v i initiates a package with three counters v n , e n , and f n and sets them as 1, 0, 0, respectively. v i marks itself and sends the packet to one of its direct neighbors, denoted as v j . v j marks itself and adds v n by one. v j checks its neighboring edges. For one with both ending vertices marked (e.g. e ij ), v j adds e n by one. v j also checks its neighboring triangles. For one with three ending vertices marked, v j adds f n by one. v j then sends the package to one of its unmarked neighbors. Once every vertex has been marked, v n , e n , and f n in the packet have counted the total numbers of vertices, edges, and faces of M, respectively. The number of holes of M, denoted as b, can be counted in a similar way. A non-triangular face initiates a package with a counter b set to one. The package is forwarded to neighboring faces until all faces in M have been visited.
Once v n , e n , f n , and b are flooded through the network, each node can easily compute the genus number of M based on Euler-Poincare´Theorem.
Computing virtual planar coordinates With the extracted triangular mesh and known topology of the surface, we apply discrete surface Ricci flow 36 to compute the special set of edge lengths, the uniformization metric of M. As we discussed in section ''Preprocessing,'' we allow the mesh with non-triangular faces, n-polygons with n.3, due to the randomness of sensor deployment. For a face with n.3, the algorithm simply adds one virtual vertex and virtual edges to connect the virtual vertex with all other vertices of the face. The face with n.3 is then split into n triangle faces. We refer readers to Jin et al. 36 for theoretical background introduction of discrete surface Ricci flow.
Here, we only provide the algorithm of computing uniformization metric using discrete surface Ricci flow. The algorithm is fully distributed, and each node only needs to exchange information with its one-range neighbors.
Specifically, each edge e ij has its initial edge length set to unit one. Each vertex v i is associated with a unit circle with radius denoted as g i = 1. The intersection angle of the two unit circles centered at v i and v j is set as the edge weight f ij at e ij . In each iteration, each v i exchanges its current g i with its direct neighbors and updates its adjacent edge lengths fl ij je ij 2 Eg according to the following Euclidean cosine law for genus one surface l 2 ij = g 2 i + g 2 j + 2g i g j cos f ij or the hyperbolic cosine law for genus larger than one surface cosh l ij = cosh g i cosh g j + sinh g i sinh g j cos f ij With the updated edge lengths, v i updates its corner angles fu jk i jf ijk 2 Fg according to the following inverse Euclidean cos law for genus one surface u jk i = cos À1 l 2 ki + l 2 ij À l 2 jk 2l ki l ij or the inverse hyperbolic cos law for genus larger than one surface u jk i = cos À1 cosh l ki cosh l ij À cosh l jk sinh l ki sinh l ij Then, v i computes its current discrete Gaussian curvature K i (equation (1)). If for every v i , its current Gaussian curvature K i is less than a threshold, discrete surface Ricci flow converges. Otherwise, each v i updates its u i : u i = u i + d( K i À K i ), where u i = log g i and d is the step length.
With the computed edge lengths when discrete surface Ricci flow converges, any node can initiate an embedding process. We assume the node with the smallest node ID, denoted as v i . Its planar coordinates denoted as uv(v i ) are set to (0, 0). Then, it chooses one of its direct neighbors, denoted as v j . For genus one surface M, the planar coordinates of v j is set to uv(v j ) = (0, l ij ). Otherwise, the planar coordinates of v j is set to uv(v j ) = ( tanh (l ij =2), 0). For vertex v k , adjacent to both v i and v j , it calculates the intersection points of two circles centered at v i and v j with radii l ik and l jk , respectively. For genus one surface M, the two circles are Euclidean circles; otherwise, they are hyperbolic circles. A hyperbolic circle can be easily converted to a Euclidean circle according to equation (3) given in Appendix 1, so the computation of the intersection points of two hyperbolic circles is equivalent to the computation of the intersection points of two Euclidean circles. v k then chooses one intersection point as its planar coordinates such that v i , v j , and v k satisfy right-hand rule in plane. With f ijk embedded in plane, we continue to embed adjacent triangles with f ijk into plane. Such breadth-first embedding continues until each vertex of M has at least two sets of planar virtual coordinates.

Resilient greedy routing
With the stored virtual planar coordinates at each node, the whole network is embedded in plane seamlessly. Locally, any region of the network is mapped to plane one-to-one and continuously. The algorithm then removes all the virtual vertices added in section ''Computing virtual planar coordinates'' and considers faces with n.3 holes. For edges shared by only one triangle face, they are labeled as boundary edges, otherwise, non-boundary edges. Vertices connecting with boundary edges are labeled as boundary vertices.
Given a pair of source and target nodes denoted as S and T , respectively, considering both S and T having more than one set of planar coordinates, S calculates the Euclidean/Hyperbolic distances of all combinations and chooses the closest pair as the source and target coordinates. S attaches the target coordinates to a packet and then applies greedy perimeter stateless routing (GPSR) 8 to forward it. Specifically, in greedy forwarding mode, a packet is forwarded to next hop with shortest Euclidean/Hyperbolic distance to T . When greedy forwarding encounters holes or local minimums, the packet is changed to perimeter mode. In perimeter model, if the current stuck node is a boundary node, the packet is forwarded along the boundary in counterclockwise direction until greedy forwarding is achievable. Otherwise, the packet is stuck at a local minimum. It is forwarded to the next node of the current triangle face following counterclockwise direction.
At the same time, each node periodically updates its neighboring nodes and local mesh structure due to potential node or link failures. However, the resulting mesh may exist some dangling vertices and edges as shown in Figure 3 that violate a triangular mesh structure. A dangling vertex is a vertex belonging to more than one boundaries. The existence of a dangling vertex confuses the orientation of boundaries. We remove dangling vertices by vertex split operation. As shown in Figure 3(a), v 1 is connected with boundary vertices v 2 , v 4 , v 5 , and v 7 that belong to two disconnected boundaries with boundary edges marked with red color. We virtually split v 1 and denote the newly added node as v 0 1 . v 0 1 shares the same planar coordinates as v 1 . v 1 is then connected with v 7 and v 2 , and v 0 1 is connected with v 4 and v 5 . Similarly, a dangling edge is an edge belonging to more than one boundary and not shared by any triangle face. We remove dangling edges by edge removal operation. As shown in Figure 3(b), e 01 belongs to two disconnected boundaries. We remove this edge such that the previously disconnected boundaries are merged into one. The graph now is a triangular mesh with each face a triangle. A boundary vertex or edge belongs to exactly one boundary. Therefore, boundary edges and vertices along newly generated holes can be easily detected and labeled.

Time complexity and communication cost
Given a network with n nodes, the time complexity to extract a triangular mesh M from the connectivity graph of the network is O(n). The time complexity to detect the network topology is O(n) too. The time complexity to compute uniformization metric of M using discrete surface Ricci flow is measured by the number of iterations. It is given by ÀC( log e=l), where C is a constant, e is the threshold of curvature error, and l is the step length of each iteration. 36 We set e and l to 1e À 5 and 0:1, respectively. The total number of iterations of each network model is no more than hundreds in our simulations. The time complexity to seamlessly embed M in the plane is also linear to n.
The majority of communication cost, measured by the number of messages, comes from computing uniformization metric using discrete surface Ricci flow. It is O(ÀC( log e=l))nd, where d is the average number of vertex degree since each vertex needs to exchange messages with its direct neighbors in each iteration.
In section ''Simulations,'' we apply the normalized communication cost to measure the consumed energy of a network to conduct a routing. We calculate the average number of messages exchanged in the network during one communication, whether it succeeds or not. The results show that the communication cost of our proposed resilient greedy routing algorithm is O(n) for any pair of source and target nodes based on the precomputed virtual planar coordinates. It is well scalable to the size of a network.

Simulations
The proposed resilient routing algorithm can be applied for a wide range of large-scale sensor networks densely deployed on surfaces of coal mine tunnels for disaster prevention and rescue, or corridors of buildings for fire detection, or sewer or gas systems for monitoring underground pipelines. These surfaces generally have complex shapes and multiple handles, that is, genus numbers. To this end, we create several representative network models as shown in Figure 4 to carry out simulations. These network models have various genus numbers, ranging from one to four. Their sizes also vary from 500 to 1100 nodes.
In our simulations, we check two different failure models of a network: independent node failure, where each node in network has an equal and independent probability of failure during some time interval; and geographically correlated failure, where all nodes within a fixed radius or hop count fail simultaneously and the center of failure region is along the primary path. We assume the source and target nodes can locate anywhere in a network. We use three metrics to evaluate the performance of a resilient routing algorithm. The first one is the successful delivery ratio. It measures the probability of an algorithm successfully finding a path when the nodes or links of a network are not reliable. The successful delivery ratio indicates the resilience of a routing algorithm. The second metric is the average stretch factor. We divide the end-to-end hop count from source to target including the cost of forwarding messages and local detours, by the hop count of a shortest path without failure triggered. The third metric is the normalized communication cost. We calculate the average number of messages exchanged in a network during one communication regardless of the communication succeeds or not. This metric reflects the consumed energy of the network to conduct one (successful/failed) communication and the durability of the system. For those broadcasting-based methods, broadcasting dominates the communication cost of a newly established pair of source-target nodes. However, different source-target pairs sharing the same source node need only one-time broadcasting to build various paths from one source to different target nodes. We define source-to-pairs ratio as the reciprocal of the number of source-target pairs sharing the same source. To be fair, we set two different source-to-pairs ratios: 1 : 100 and 1 : 400 in our simulations.
Since the primary target of the proposed method is network resilience under independent node failure model, we compare our method with other well-known resilient routing methods that choose local detour instead of disjoint alternative paths, including Directed Diffusion, 24 Braided Multipath Routing, 23 and Robust Cooperative Routing. 38 Besides, SINUS method 22 is the most recent greedy routing algorithm specifically designed for networks on high genus surfaces. We implement it with GPSR module to enhance its resilience against node failures and consider it as another comparison method. We randomly choose 20, 000 source-target pairs on each network model. Simulation results given in sections ''Independent node failure'' and ''Geographically correlated failure'' show that our method achieves the best performance under both independent node and geographically correlated failure models.

Independent node failure
The first row of Figure 5 compares the successful delivery ratio of different methods as a function of the probability of independent node failure for networks shown in Figure 4. As the probability increases, the successful delivery ratios of Directed Diffusion, Braided, and Robust Cooperative Routing methods decrease dramatically. When the probability of independent node failure is 15%, more than half of the sampled sourcetarget pairs of these methods fail. However, SINUS with GPSR and our method perform much better. As the probability of independent node failure continuously increases, our method shows a bigger advantage. The reason is that SINUS generates a planar embedding with uneven density. Failure nodes can easily cause isolated regions around the source or target nodes such that it is impossible to find a successful detour with GPSR in this situation.
For the same set of network models, the second row of Figure 5 shows the average stretch factors of successful deliveries under different probabilities of individual node failure. As we can see, our method has a very low stretch factor, while SINUS with GPSR generates much longer paths in the network. We have discussed the reason in section ''Introduction.'' SINUS cut the original surface network open along computed virtual boundaries to an annulus and then embed the annulus to plane. A pair of nodes with the shortest path crossing a virtual boundary in the original network needs to take a much longer detour to reach each other based on the computed virtual planar coordinates. We set the maximum probability of individual node failure for comparison of the average stretch factor as 10% À 20% only. If the failure probability is higher, the successfully delivered paths of comparison methods are quite short since longer paths are more likely to fail. To avoid such biased influence on the comparison of the average stretch factor, we compare different methods when there is still a significant number of successfully delivered paths.
We set the probability of independent node failure as 10% when measuring the communication cost. and 1 : 100, respectively. As we can see, even with a smaller source-to-pairs ratio, Directed Diffusion, Braided, and Robust Cooperative Routing spend much more energy than our method. When the ratio is larger, the energy consumption of these methods is more dramatic. Compared with SINUS with GPSR, our method consumes around 1=3 less energy for successful delivery in Network I and II. The difference is bigger when a network has more nodes and higher genus number. Besides, our method wastes much less for unsuccessful communications. Therefore, our method benefits the system a lot, especially for those time-sensitive applications that involve frequent communications in a network.

Geographically correlated failure
The first row of Figure 7 compares the successful delivery ratios of different methods as a function of the radius of a failure region. As the radius of a failure region increases, the delivery ratios of Directed Diffusion, Braided, and Robust Cooperative Routing decrease dramatically. However, our method consistently performs the best with the highest successful delivery ratio, even when the failure region is huge. The only exception is the Network II model. When the radius of a failure region is larger than 3 in Network II model, both SINUS with GPSR and our method drop dramatically. The reason is that Network II model has narrow handles with fewer nodes covered compared with other network models. When the primary path of a pair of source and target nodes passes along a handle, a large failure region can easily break the handle such that a local detour cannot find a successful alternative. It is obvious that when the radius of a failure region becomes huge in geographically correlated failure model, a routing method with high resilience needs to take the global topology information of a network into consideration and choose a path with different homotopy type, that is, a path along different handles. The second row of Figure 7 shows the average stretch factors of successful deliveries under different radii of a failure region. Since the successful delivery ratio of Robust Cooperative Routing method drops to zero when the radius of a failure region is 1, we remove it from the comparison. As the radius of a failure region becomes larger, the successful delivery ratio of Directed Diffusion and Braided methods drops dramatically, and the successfully delivered paths are short ones since the longer paths all fail. To avoid a biased influence on the computation of average stretch factor, we set the maximum radius of a failure region as 2 only. With the increase in the radius of a failure region, the stretch factor of a successful routing path gets  longer for all methods. However, our method still has a consistently lower stretch factor.
We set the radius of a failure region as 1 when measuring the communication cost. Figure 8 shows the normalized total/successful/failed communication cost under geographically correlated failure model when the source-to-pairs ratio is 1 : 400 and 1 : 100, respectively. Since the successful delivery ratio of Robust Cooperative Routing method drops to zero when the radius of a failure region is 1, we remove it from the comparison. The results are similar to those under the independent node failure model. Our method consistently performs the best among all.

Conclusion
We have proposed a cut-free resilient greedy routing scheme for GPS-free large-scale wireless sensor networks deployed on surfaces of complex-connected 3D settings. It is designed to be resilient to node or link failures especially under random node or link failure model. Simulation results have shown that the proposed scheme achieves a higher successful delivery ratio, a lower average stretch factor, and a lower normalized communication cost compared with other resilient routing methods.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors are partially supported by NSF CNS-1320931.