Olivier Dézerald
email: olivier.dezerald@gmail.com

Céline Leroy

Bruno Corbara

Jean-François Carrias

Laurent Pélozuelo

Alain Dejean

Régis Céréghino

Food-Web Structure in Relation to Environmental Gradients and Predator-Prey Ratios in Tank-Bromeliad Eco systems

Little is known of how linkage patterns between species change along environ mental gradients. The small, spatially discrete food webs inhabiting tank-bromeliads provide an excellent opportunity to analyse patterns of community diversity and food-web topology (connectance, linkage density, nestedness) in relation to key environmental variables (habitat size, detrital resource, incident radiation) and predators:prey ratios. We sampled 365 bromeliads in a wide range of understorey environments in French Guiana and used gut contents of invertebrates to draw the corresponding 365 connectance webs. At the bromeliad scale, habitat size (water volume) determined the number of species that constitute food-web nodes, the proportion of predators, and food-web topology. The number of species as weil as the proportion of predators within bromeliads declined from open to forested habitats, where the volume of water collected by bromeliads was generally lower because of rainfall interception by the canopy. A core group of microorganisms and genera list detritivores remained relatively constant across environments. This suggests that (i) a highly-connected core ensures food-web stability and key ecosystem functions across environments, and (ii) larger deviations in food-web structures can be expected following disturbance if detritivores share traits that determine responses to environmental changes. While linkage density and nestedness were lower in bromeliads in the forest than in open areas, experiments are needed to confirm a trend for lower food-web stability in the understorey of primary forests.

Introduction

Biodiversity is not only the sum of coexisting species, it is also the diversity of interactions that connect these species. It is widely acknowledged that linkage patterns among species, rather than species richness alone, are a key component of the diversitystability relationship [START_REF] Gilbert | Connectance indicates the robustness of food webs when subjected to species Joss[END_REF]. For instance, for a fiXed number of species, food-web persistence and robustness (two concepts linked to network stability) are believed to increase with connectance, the proportion of all possible interactions that are realized [START_REF] Rall | Food-web connectance and predator interference dampen the paradox of enrichment[END_REF]. A rich body ofliterature has focused on the relationship between network structure and dynamics [START_REF] Kondoh | Foraging adaptation and the relationship between food-web complexity and stability[END_REF], and on how different types of interactions (either antagonistic or mutualistic) affect network properties [START_REF] Fontaine | The ecological and evolutionary implications of merging different types of networks[END_REF]. Most studies however focused on network structure regardless of variation in the abiotic environment. Consequendy, our understanding of the environmental determinants of network structure lags behind the increasingly vast knowledge of species richnessand interactions-mediated patterns [START_REF] Strong | Insects on Plants: Community Patterns and Mechanisms[END_REF][START_REF] Dunne | Food-web structure and network them-y: the role of connectance and size[END_REF]. Identifying linkages between ecological networks and physical environments could be relevant to predict the extent to which entire communities will respond to increasing environmental disturbances (e.g., climate change, habitat fragmentation, overharvesting, introduction of invasive species, pollution).

On a local scale, both the spatial/temporal changes in habitat conditions and species traits regulate community structure [START_REF] Townsend | The patch dynantics concept of stream community ecology[END_REF], thereby affecting ecosystem functions. For instance, the nature and extent of the riparian vegetation determines food quality and availability in still [START_REF] Declerck | Ecological characteristics of small farnùand ponds: Associations with land use practices at multiple spatial scales[END_REF] and running waters [START_REF] Minsball | The River Continuum Concept[END_REF], th us influencing the functional feeding group composition of macroinvertebrate communities [START_REF] Compin | Spatial patterns ofmacroinvertebrate functional feeding groups in streams in relation to physical variables and land-cover in Southwestern France[END_REF]. However, we know little about how linkage patterns among species (food-web topology) change against a background of environmental gradients over broad, regional scales. Most of what we know about environmental determinants of food-web structure cornes from experirnents in single locations [11], on a target fraction of the food web (e.g., vertebrates, macro-and micro-invertebrates excluding components such as microorganisms, bacteria, flagellates, rotifers and viroses; [START_REF] Paradise | Relationship of water and leaf litter variability to insects inhabiting treeholes[END_REF]).

It has proven challenging to study the effects ofboth abiotic and biotic factors on food webs. High species diversity and population densities in large ecosystems preclude the accurate characterization of trophic links (notably in continuous habitats), and because it is often difficult to manipulate entire ecosystems. To tackle these issues, we focussed on small, spatially discrete food webs that naturally span a broad range of environmental gradients. Bromeliaceae are flowering plants represented by 59 genera and sorne 3140 species native mainly to the Neotropics [START_REF] Givnish | Phylogeny, adaptive radiation, and historical biogeography in Bromeliaceae: Insights from an eight-locus plastid phylogeny[END_REF]. The interlocking leaves of tank-forming bromeliads form wells that collect water, leaf litter and other organic detritus. The detritus that enter the tank constitute the main source of nutrients for the aquatic food web [START_REF] Benzing | Bromeliaceae: profile of an adaptive radiation[END_REF]. The aquatic communities inhabiting tankbromeliads provide excellent opportunities to study food-web structure because they contain several trophic levels (from bacteria to predatory macroinvertebrates; [START_REF] Brouard 0 | Understorey environments influence functional diversity in tank-bromeliad ecosystems[END_REF]) and can be exhaustively sampled [START_REF] Romero | Food-web composition affects crossecosystem interactions and subsidies[END_REF].

The aim of this study was to determine whether spatial patterns in food-web structure can be predicted from a small set of environmental factors and/ or the richness of predators relative to their prey. To address this question, we sampled 365 tankbromeliads in a wide range of environments in French Guiana (plantations, pioneers growths, rock savannah, primary forest) and drew the 365 corresponding food webs using the gut contents of invertebrates as weil as observations of predator-prey encounters. Previous studies on tank-bromeliads concluded that detrital resources at the base of the food web, understorey light environments (energy available for algal production), habitat size and predation play key roles in shaping aquatic community composition [ 1 7, 18]. Detritus constitute the main source of energy for aquatic bromeliad communities, however, the high algal biomass found in sun-exposed bromeliads [START_REF] Brouard 0, Lejeune | Are algae relevant to the detritus-based food web in tank-bromeliads?[END_REF] may provide a complementary non-detrital resource to the upper trophic levels. If (i) species richness and abundance increase with bromeliad (habitat) size [START_REF] Armbruster | Factors influencing community structure in a South American tank bromeliad fauna[END_REF][START_REF] Jabiolj | Structure of aquatic insect communities in tank-bromeliads in a East-Amazonian rainforest in French Guiana[END_REF][START_REF] Cérégbino | Antplant mutualisms promote functional diversity in phytotelm communities[END_REF][START_REF] Leroy | Ants mediate foliar structure and nitrogen acquisition in a tank-bromeliad[END_REF][START_REF] Céréghino | Ants mediate the structure ofphytotelm communities in an ant-garden bromeliad[END_REF], and (ii) bromeliads from open area benefit from bath detrital and non-detrital resources and understorey bromeliad foodwebs are solely supported by detrital matter [START_REF] Brouard 0 | Understorey environments influence functional diversity in tank-bromeliad ecosystems[END_REF], then, for a given habitat size, the diversity of invertebrate functional feeding groups should increase from forest to open areas. One may expect shifts in food-web connectance, linkage density, and/or nestedness as species with particular traits are replaced or complemented by species with other traits when shifting from forest understorey to open areas. Conversely, if the resource that supports food webs does not differ from forest to open areas, then one should not observe significant shifts in foodweb structure and functions. To test these hypotheses, we used Linear Mixed Effect modelling to analyse patterns of community diversity and food-web topology in relation to key environmental variables and predator:prey richness ratios.

Materials and Methods

Ethics Statement

This study was conducted according to relevant national and international guidelines. Sample collections necessary to scientific research were authorized by the French Office National des Forêts (ONF) provided that their impact upon the environment is considered negligible.

Study Area, Bromeliads and Environmental Variables

The study was conducted in French Guiana, from March 2006 to October 2011. The climate is tropical moist with 3,000 -3,400 mm of yearly preCipitation mainly distributed over 280 days. There is a major reduction in rainfall between September and November and another shorter and more irregular dry period in March. The maximum monthly temperature averages 33.5°C (32.1-35.8°C), and the monthly minimum averages 20.3°C (19.7-2l 0 C).

We selected five sampling localities distributed across a southeast to north-west range (Fig. 1) and sampled 365 bromeliads (i.e. 365 food webs) in the understorey of primary and transitional forests, in pioneer growths, in a rock savannah, and in plantations. The main habitat characteristics of tank-bromeliad species found in five vegetation types at five localities (hereafter "vegetation types", within localities) are summarized in Table l. Further descriptions of the Nouragues (sampling period in April 2006), Petit-Saut (March 2009) and Kaw (Oct. 2008) localities and their bromeliads can be found in [START_REF] Brouard 0 | Understorey environments influence functional diversity in tank-bromeliad ecosystems[END_REF][START_REF] Cérégbino | Antplant mutualisms promote functional diversity in phytotelm communities[END_REF][START_REF] Bongers | Dynantics and plant-animal interactions in a neotropical rainforest[END_REF]. Saint-Elie (Oct. 2007) and Angoulème (April 20 1 0) are Citrus plantations. Bath epiphytie bromeliads and those that had taken root on the ground were included in the study. ln arder to prevent seasonality effects (i.e., dry vs rainy season) on water volumes and species abundance, the sampling took place at the transition between rainy and dry seasons, where invertebrate abundances reach a peak. We also sampled tank-bromeliads that were full of water. Species occurrence was not an issue as bromeliad invertebrates have fast life cycles (<2 months) with overlapping cohorts. To characterize habitat size, we emptied the wells in each plant by sucking the water out (see invertebrate sampling) and recorded the corresponding volume of water (WV, mL). The amount affine particulate organic matter (FPOM; 1000 to 0.45 j.UI1 in size) was expressed as preserved volume (mL) after decantation in graduated test-tubes [START_REF] Paradise | Relationship of water and leaf litter variability to insects inhabiting treeholes[END_REF]. Finally, percentages of total incident radiation (IR) above the bromeliads were calculated using hemispherical photographs and an image processing software (Gap Light Analyzer 2.0; [START_REF] Frazer | Gap Light Analyzer (GLA) 2.0: Imaging software to extract canopy structure and gap light transmission indices from true-colour fisheye photographs: users manual and program documentation[END_REF]), as described in [START_REF] Leroy | Ants mediate foliar structure and nitrogen acquisition in a tank-bromeliad[END_REF]. ln this study, we consider that ail bromeliads with an IR below 50% were located in partially shaded areas or forested environments while higher percentages defined sun-exposed areas.

Aquatic lnvertebrate Communities

For bath ethical (extensive sampling could destroy local populations) and legal reasons (the Nouragues Research Station is located in a nature reserve), we decided to use a non-destructive sampling technique for ail studied bromeliads. To sample the water retained in the tanks, we used 5-mL and 1 0-mL automatic micropipettes with the end trimmed to widen the orifice. We homogenized the water within leaf wells by sucking in and out with a pipet, before sucking out the content. Although less efficient than plant dissection, we and other researchers have already successfully used this technique [START_REF] Céréghino | Ants mediate the structure ofphytotelm communities in an ant-garden bromeliad[END_REF][START_REF] Jocqué | How effective are non-destructive sampling methods to assess aquatic invertebrate diversity in bromeliads?[END_REF]. It was consistendy used for ail of the samples and most of the water (>95%) was collected. Bath early and late instars of prey and predator invertebrates were captured, so we were confident that our technique was efficiendy implemented. The samples were preserved in the field in 4% (fmal concentration) formalin. lnvertebrates were sorted in the laboratory and preserved in 70% ethanol. They were identified to genus, species or morphospecies and enumerated (species lists in [START_REF] Jabiolj | Structure of aquatic insect communities in tank-bromeliads in a East-Amazonian rainforest in French Guiana[END_REF][START_REF] Céréghino | Ants mediate the structure ofphytotelm communities in an ant-garden bromeliad[END_REF]). Species abundance data (individuals per plant) were used to calculate evenness (Simpson index) and entropy (Shannon-Weaver index) for each invertebrate community. The Simpson index was calculated as D = 'L pt The Shannon-Weaver index (hereafter Shannon index) was defined as H= -'L Pi log(b) p;, where Pi is the proportional abundance of species i and b is the base of the logarithm (natural logarithm in this case). Species richness and abundance being components of these two indices, they were not so • so • 40 . u so-30 . 1 ).

doi:l 0.1371/journal.pone.0071735.g001 considered alone in subsequent modeling. Macroinvertebrate taxa were partitioned into predators (carnivorous species which attack and consume living macroinvertebrates) and prey (species which sift fme particulates and microorganisms from the water column and/ or gather FPOM and associated microorganisms from the accumulated debris), and these categories were used to calculate Predator:Prey Richness ratios (number ofpredatory taxa/number of prey taxa, hereafter PPR).

Food Web Characterization

The diet of the various invertebrate species that make up food webs was determined by dissecting the entire guts. Twenty ta 50 individual guts from each taxon were dissected across bath localities and vegetation types in arder ta encompass the variability in diets of omnivores and predators. The gut contents of predators were placed into a drop of water on a glass slide, spread out, and analysed using a binocular microscope (Leica® MZ 12.5) and an Optiphot-2 Nikon® microscope whenever necessary. Most of the prey items could be identified from the guts of predators by comparing the chitinous parts (head capsules or legs, setae ... ) with specimens of bromeliad invertebrates archived in our collection (Univ. of Toulouse III). Only the gut contents of piercers (Heteroptera Velüdae, Diptera Tabanidae) could not be identified visualiy; in this case, we relied on observations of arranged encounters in test tubes.

The microorganisms found in the guts of detritivores were identified ta a coarse taxonomie leve! using an Optiphot-2 Nikon microscope at x600 magnification. Determinations were based on earlier descriptions ( [START_REF] Carrias | A preliminary study of freshwater protozoa in tank bromeliads[END_REF][START_REF] Carrias | An ant-plant mutualism induces shifts in the protist community structure of a tankbromeliad[END_REF]), and on recent microbiological samples (Carrias J-F, unpublished data). The presence of particulate organic matter was also recorded in the gut contents, and we used an ocular micrometer ta distinguish fme particulate organic matter (FPOM, 1000 ta 0.45 11m in size) from coarse particulate organic matter (CPOM, > 1000 llffi). Gut contents and other observations were used ta build 365 interaction matrices (one per bromeliad) that described 365 "connectance webs" sensu [START_REF] Post | The long and short food chain length[END_REF] (see Fig. 1 ).

Three topological descriptors were used ta describe food webs: connectance, linkage density, and nestedness. ln the literature, two measures of the connectance are widely used: connectance and normalized connectance. Bath refer ta the proportion of ali possible interactions that are realized, but normalized connectance is rather used for food webs comprising Jess than 20 species [1 J.

However, in this study, the two measures were strongly correlated (Pearson's correlation coefficient= 0.966, p<O.OOO 1) and we therefore used the connectance as defmed by [START_REF] Morin | Community ecology[END_REF]: C= 2xLinks Species x (Species-1)

Linkage density was calculated as the ratio between the number of links and the number of species. Levels of nestedness were expressed as matrix temperatures, with values ranging from 0 (perfecdy nested) ta 100 (see [START_REF] Bascompte | The nested assembly of plant-animal mutualistic networks[END_REF]) were calculated from the 365 interaction matrices (binary data) using the BINMATNEST3 program [START_REF] Rodriguez-Gironés | A new algorithm to calculate the nestedness temperature of presence-absence matrices[END_REF]. Departure from random was tested against a number of 1000 random matrices using the most conservative nul! mode! in BINMATNEST3 (namely "mode! 3''). A matrix temperature is affected by bath its structure (the fil!) and its size (the number of rows and colunms). The extent ta which the observed temperature departs from random depends on how random matrices are built. Unlike the nul! models 1 and 2, the mode! 3 is a packing algorithm that generates random matrices that are Jess influenced by the structure of the input matrix, thus aliowing for comparisons with other studies [START_REF] Rodriguez-Gironés | A new algorithm to calculate the nestedness temperature of presence-absence matrices[END_REF]. More specifically, the nul! mode! 3 generates random matrices where each cel! has a probability of containing a "l" equal ta the averaged probability of occupancy of rows and colunms of the input matrix. Therefore, the fil! of each cel! depends on bath the fraction of "!s" in rows and columns in which the cel! is included.

Data Analysis

Ta analyze the relationships between food-web descriptors, environmental variables and PPRs, we used linear mixed effect modeling. Ali variables were log-transformed ta fit a normal distribution. Because a given vegetation type was found at only one locality, we qualified this structure as the "vegetation type nested within locality". ln contrast, a given bromeliad species could be sampled at many localities or many vegetation types, therefore the variable "bromeliad species" was not nested within the former or the latter. Since "locality", "vegetation type", "vegetation type nested within locality" or "bromeliad species" could be included as random factors, we performed a mode! selection based on the more conservative Bayesian Information Criterion (BIC) of the full models (models considering ali environmental variables) for each dependent variable. Then, for a selected random factor (the remaining random effects being dropped from the madel), the relationship between dependent and independent variables were explored using a stepwise backwards removal procedure and only the final models containing significant variables were presented. Departures from homoscedasticity and the normality of the residual errors were evaluated graphicaliy for each fmal mode!. Ali statistical analyses were evaluated under a 95% confidence leve! and were conducted using R software V. 2.14.1 [START_REF]R: A language and environment for statistical computing[END_REF] and the associated packages for Linear Mixed Effect Modeling (nlme, lmer4).

Results

Abiotic Environments and Food-web Composition

There were large differences in WV and FPOM among bromeliad species, even at a given vegetation type (Table 1). The two bromeliads growing in rock savannah (A. aquilega and C. berteroniana) showed the highest and lowest mean values for FPOM (10.84 mL vs. 0.74 mL, respectively), and the highest and one of the lower mean values for WV (949.23 mL vs. 40.62 mL, respectively) compared ta the others tank-bromeliad species. Th us, even though the understorey of primary forests consistendy received lower incident radiation (i.e. V. splendens and G. lingulata) and higher litter inputs than sun-exposed areas, the amount of FPOM inside the tank ofbromeliads from primary forests was not necessarily higher than those from rock savannah. However, bromeliads from primary forests had higher FPOM:WV ratios (0.05--0.09 mL FPOM/mL WV) than bromeliads from open areas (0.01--0.02 mL FPOM/mL WV).

Food webs comprised 8 ta 20 taxa (median= 11 taxa), including macroinvertebrates, rotifers, and miscellaneous microorganisms (bacteria, cyanobacteria, fungi, algae, heterotrophic flagellates, and ciliates) identified in the gut contents. The predators belonged ta the Diptera Culicidae (Toxorhynchites purpureus), Corethrellidae (Corethrella sp.) and Tabanidae, and ta the Odonata (one unidentified Coenagrionidae species). Detritivores mosdy consisted in Diptera Culicidae (Culex spp. and Jt!Yeomyia spp.), Limoniinae, Chironomidae (Ta'!)ltarsus sp.) and Psychodidae (Telmatoscopus spp.). The list of taxa and the corresponding functional groups is provided in Fig. 2.

Model Selection and Influence of Abiotic and Biotic variables on Food Webs

Ali random factors were selected in at !east one mode! (Table 2), suggesting that different spatial scales could provide complementary exp1anations ta the observed pattern. More specificaliy, "vegetation type" and "bromeliad species" were included as random factors in the Simpson index (BIC= -314. [START_REF] Romero | Food-web composition affects crossecosystem interactions and subsidies[END_REF]) and connectance models (BIC= -1108.648), respectively. "Locality"

D Figure 2. Examples of connectance webs under two contrasted environ mental conditions: open areas (A, B) and forest areas (C, D).

The upper trophic level (predators) is at the top of the graphs, and the lower (algae, detritus) at the bottom. Numbers and abbreviations are for: Wyeomyia spp. (1; filter-feeder); Culex spp. (2; filter-feeder); Forcypomiinae sp2 (6; filter-feeder); Bezzia sp. (7; predator); Tanypodinae (9; predator); Chironominae (10; detritivore); Tanytarsinii (11; detritivore); Corethrella sp. (12; predator); Telmatoscopus sp1 (14; detritivore); Limoniinae (17; detritivore-shredder); Cyphon sp. (22; shredder); Coenagrionidae (24; predator); Aulophorus superterrenus (25; detritivore); Hydracarina (28; detritivore); Elpidium sp. (34; detritivore-scraper); Bacteria (1); Cyanobacteria (Il); Fungi (Ill); Algae (IV); Flagellate (V); Ciliate (VI); Rotifera (VIl); Fine Particulate Organic Matter (FPOM). doi:1 0.1371/journal.pone.0071735.g002 was selected in the Shannon index, linkage density and nestedness models (BIC= 21.83, BIC= -586.99, and BIC= -546.25, respectively). The subsequent stepwise backward removal procedures for each food-web descriptor were thus performed using their selected random factors. The computation oflikelihood ratios demonstrated that we could not discriminate the madel with "locality" from the one with "vegetation type" as random factor for the Shannon index madel (L = 0.034, p = 0.92). Therefore, the variable selection for the Shannon index madel was also performed with "vegetation type" as random factor. In general, the inclusion of random factors with higher BIC values (e.g., "vegetation type", "vegetation type nested within locality" or "bromeliad species") only slighdy affected the outcomes of the models in terms of significant and non-significant interactions among variables. More specifically, interactions that were barely significant with the selected random factor became non-significant when other random factors were included. The Shannon and Simpson indices were both positively correlated with WV (p<O.OOOl) and with PPR (p =0.0009 and p = 0.0027, respectively; Table 3), showing that invertebrate community diversity increased with habitat size and the relative number of predators. Moreover, these two models showed a negative interaction between FPOM and PPR (p = 0.036 and p = 0.04 7, respectively). However, when "vegetation type" was included as random factor in the Shannon index mode!, both water volume and PPR remained significant but the FPOM:PPR interaction was not significant (p = 0.081 ). Last, the variable "IR" was not significandy correlated with the Shannon index, whatever the random factor included (p = 0.6--0.7 5). Finaliy, the two negative interactions IR:WV and PPR: FPOM become significant (p = 0.046 and p = 0.0 14, respectively) for the Simpson index when "Locality" was included as a random factor.

Connectance showed a negative and significant correlation with WV and PPR (p = 0.0079 and p<O.OOO 1, respectively), and there was a trend for decreasing connectance with increasing FPOM amount (p = 0.039; Table3). Incident radiation was not significandy correlated with connectance (p = 0.24). When "Locality" was included, the two negative interactions IR:WV and WV:FPOM turned to be significant (p = 0.0002 and p = 0.009, respectively) in the connectance model. Both the linkage density and nestedness models showed a significant negative interaction between WV and IR (p=0.008 and p=0.0071, respectively). Linkage density and nestedness were positively correlated with WV and IR, and negatively correlated with PPR (Table 3). Linkage density was positively correlated with the IR:FPOM interaction (p = 0.034), but did not correlate with FPOM alone (p = 0.064). Finaliy, levels of nestedness (Nj ranged from 21.54 to 49.79 (median N= 39.98). Overali, 44.4% of our interaction matrices were significandy different from the null matrices obtained at random (p=O.OOOl to 0.99, median p=0.08). Nevertheless, this percentage increased with the number of taxa; for instance, 66.1% and 88.8% of ali of the interaction matrices comprising :::::13 and :::::17 taxa, respectively, were significandy different from random.

Discussion

Based on an unprecedented number of replicates un der natural conditions, our models have proven informative in assessing whether environmental factors and biological interactions explain sorne patterns of food-web structure. Overali, variations in foodweb structure across environments and spatial scales were mosdy due to interactions between habitat size and the distribution of predators (Kendali's test, positive correlation between WV and PPR, p=O.l79,p<0.0001). There was however a trend for more feeding links in open areas compared to forest understorey. Finaliy, a core of generalists remained relatively constant across large regional scales, so that patterns of food-web connectance, linkage density and nestedness were mosdy related to the addition or loss of predators.

Water volume was one of the most significant variables in ali models, suggesting that habitat size was the primary factor controlling species composition and abundance patterns and, subsequendy, food-web structure. Larger habitats are more easily colonized by immigrants, resulting in positive species-area relationships [START_REF] Drakare | The imprint of the geographical, evolutionary and ecological context on species-area relationships[END_REF]. An increase in habitat area also fosters functional diversity [START_REF] Takimoto | Environmetal deterntinants of food-chain length: a meta-analysis[END_REF], while aliowing for a better partitioning of food resources by coexisting species [START_REF] Richardson | Nutrients, diversity, and community structure of two phytotelm systems in a lower montane forest, Puerto Rico[END_REF]. Hence, at any given locality, both the number of taxa and individuals per plant increased with habitat size (there was a positive correlation between WV and PPR, Kendall's test=O.l79, p<O.OOOl), which is consistent with previous fmdings on bromeliad ecosystems [START_REF] Jabiolj | Structure of aquatic insect communities in tank-bromeliads in a East-Amazonian rainforest in French Guiana[END_REF] and smali wedands in general [START_REF] Oertli | Does size matter? The relationship between pond area and biodiversity[END_REF]. However, there was a trend for bromeliad with higher incident radiation to accumulate more water (positive WV to IR relationship, r = 0.32, p<0.05, logtransformed data) certainly because in open areas there are fewer overhanging trees to keep most of the rain from reaching the bromeliads. lndeed, based on measures made in the understorey of the primary forest in Petit-Saut, we estimated that 30--38% of the rain is intercepted by the canopy. Hence, for a given bromeliad species and size, the containers hold more water at sun-exposed areas. Nevertheless, the complex interaction between habitat size and incident radiation (significant interaction between IR and WV in the linkage density and nestedness models, but not in the connectance mode!) precludes a precise distinction of their relative impacts on food webs. We cautiously suggest an indirect effect where IR mediates food-web features through bromeliads' carrying capacity for aquatic invertebrates. ln aquatic ecosystems, detrital inputs form a strong trophic link between plant production, decomposer microorganisms, and larger metazoans [START_REF] Vondracek | Land use, spatial scale, and stream systems: Lessons from an agricultural region[END_REF]. ln tank-bromeliads, detritus constitutes the main source of nutrients for the aquatic food web [START_REF] Benzing | Bromeliaceae: profile of an adaptive radiation[END_REF]. Debrischewing invertebrates process incoming litter. Smali particles, including faeces, are then washed into the plant pools where the FPOM is further processed in the gut of invertebrate collectors and fùterers. Dead organisms, litter, and faecal particles, which collect in the leaf bases, are utilized by bacteria and other rnicroorganisms. The nature and extent of the vegetation that surrounds these systems was therefore expected to have a strong influence on foodweb structure through food quality and availability. On one hand, algae were found to account for more than 30% of the total microbial diversity and biomass in sun-exposed bromeliads [START_REF] Marino | Habitat size deterntine algae biomass in tank-bromeliads[END_REF][START_REF] Brouard 0, Lejeune | Are algae relevant to the detritus-based food web in tank-bromeliads?[END_REF], which usually receive lower amounts ofleaflitter [START_REF] Brouard 0 | Understorey environments influence functional diversity in tank-bromeliad ecosystems[END_REF]. We suggest that such a complementary "green" food web could contribute ta reducing functional redundancy in FPOM-poor systems. Overall, there was a positive correlation between WV and FPOM (Kendall's test= 0.349, p<O.OOOl). However, higher FPOM concentrations (i.e. the amount of FPOM in relation ta WV) decreased the amount of open water inside of the tank, which had a negative effect on aquatic organism diversity. High FPOM concentrations may have the effect of clogging habitats in small freshwater ecosystems and, ta a certain extent, space availability [START_REF] Jabiolj | Structure of aquatic insect communities in tank-bromeliads in a East-Amazonian rainforest in French Guiana[END_REF]. Since FPOM was not consistendy significant in our models, we assume that resource availability is not a limiting factor in bromeliad systems. However, the influence of resource quality (detritus vs. algae) on individual food apportionment should deserve more attention in future studies.

The PPR was a highly significant variable in ail models. Connectance, linkage density and nestedness (matrix temperature) significandy decreased with increasing proportions of predators. A decrease in connectance and linkage density can be due ta a gain of specialists and/ or a loss of generalists [START_REF] Gilbert | Connectance indicates the robustness of food webs when subjected to species Joss[END_REF]. ln our study, predators were generalist species in that they fed on numerous prey species and multiple trophic levels, but they could be considered as "node specialists" since they interacted mainly with their prey while prey established links with bath basal species (algae, rotifers, ciliates ... ) and predators. The extent ta which predators are specialized on specifie prey plays a great role in generating patterns of community diversity from site-ta-site ta large regional scales [START_REF] Ryberg | Predators alter the scaling ofdiversity in prey metacommunities[END_REF][START_REF] Chase | Predators temper the relative importance of stochastic processes in the assembly of prey metacommunities[END_REF]. We found that the Shannon and Simpson indices (entropy and evenness, respectively) were positively correlated with the PPR at relatively small scales (from individual bromeliads ta local scales), but overall, bromeliad food webs comprising less than 13 species were primarily composed of microorganisms and primary consumers, and were characterized by high levels of diffuse interactions where ail species were more or less closely linked ta each other. lndeed, a core of highlyconnected species remained relatively constant at lower trophic levels (see also [START_REF] Valladares | Forest fragmentation leads to food web contraction[END_REF]), bath in terms of species identity and ecological function, between individual bromeliads and from open ta closed areas. For instance, the dominant generalist detritivores belonging ta the Culicidae (U)eomyia sp., Culex sp., and Arwpheles sp.), Chironomidae and Oligochaeta were found in 310, 215 and 133 plants out of 365 respectively. The main top-predators, namely Corethrellidae, Ceratopogonidae and Toxorl!Jnchites sp., occurred in 122, 105 and 83 plants respectively. It is likely that this highly-connected core ensures food-web stability and key ecosystem functions (e.g., decomposition, nutrient cycling) across environments. If detritivores share traits that determine responses ta environmental changes in any given area (e.g., resistance ta dessication, dispersal ability), then larger deviations in food-web structures can be expected following a major disturbance. Beyond 13 species, community diversity and food-web nestedness increased as a result of the addition of specialist predators and generalist detritivores. Our observations suggest that predators could enter the food webs at an average water volume of 49.9 mL, which corresponded ta 1 predator species for 3 prey species on average. Regardless of the locality, larger bromeliad hosted a higher proportion of predators. In addition ta their intrinsic value for biological diversity, predators (i.e., Toxorl!Jnchites sp., Corethrella sp. and Coenagrionidae) thus played a key role in generating foodweb patterns. This assumption is supported by the fact that, whilst predators are less frequent than detritivores, the predator:prey ratio was a significant variable in ail three food-web topology models (connectance, linkage density and nestedness). We suggest that the addition of species at higher trophic levels contributes ta departure from random interactions through the partition of prey resources between coexisting predators. Our results are also in line with studies which demonstrated that specialized species tend ta be dependent on a core of densely connected generalist species, and that natural communities display non-random interaction patterns [32,43--45]. Therefore, although we relied on unweighted trophic links to analyse spatial patterns in food-webs, this study provides strong empirical support on how food-web structure change along environmental and predator richness gradients. We conclude that at the bromeliad scale, habitat size determines the number of species that constitute food-web nades and the proportion of predators (positive relationship), and subsequendy, food-web topology. Hence, the number of species as weil as the proportion of predators within assemblages globally decline from open to closed (forest) habitats, where the volume of water collected by bromeliads is generally lower because of rainfall interception by the canopy. The extent to which stability is related to linkage patterns (notably nestedness) in networks such as food webs remains however unclear and obviously deserves further research [START_REF] Tylianakis | Conservation of species interaction networks[END_REF]. While linkage density and nestedness were lower in forested

Figure 1 .

 1 Figure 1. Distribution of the sampling locations in French Guiana (see Table1).

Table 1 .

 1 Main habitat characteristics (mean±standard error) of tank-bromeliad species found in five vegetation types.

	Locality	Vegetation types	Species	N	IR	wv	FPOM	#Taxa
	Nouragues	Rock savannah	CB	29	73.3:t2.4	40.62:t3.75	0.74:t0.08	10.55:t0.15
			AA	31	66.97:t2.73	949.23:t102.64	10.84:t1.26	15.87:t0.42
		Transitional forest	VP	30	25.12:t0.29	73.2:t9.11	1.7:t0.28	10.3:t0.35
			AB	26	25.69:t0.47	137.85:t21.1	1.08:t0.25	11.23 :t 0.32
		Primary forest	vs	26	18.75:t0.4	48.54:t5.03	4.29:t0.42	1 0.35 :t 0.25
			Gl	19	15.9:t0.6	17.46:t2.53	1.04:t0.16	9.21 :t0.29
	Petit-Saut	Pioneer growth	AM	63	39.33 :t 2.63	84.38:t 10.1	6.16:t1.26	12.49 :t 0.28
		Primary forest	vs	34	16.64 :t 0.38	26.18:t3.84	3.46:t0.71	1 0.68:t0.23
	Kaw	Pioneer growth	AM	45	35.1 :t2.98	92.46:t 11.82	4.3:t0.68	11.78 :t 0.28
	Angoulème	Citrus plantation	AM	35	33.03 :t 2.95	31.01 :t4.54	2.32:t0.31	11.8:t0.34
	Saint-Elle	Citrus plantation	AM	27	32.07:t1.3	56.41 :t 8.29	4.58:t0.86	14.22:t0.36

CB: C. berteroniana; AA: A. aqui/ega; VP: V. pleiosticha; AB: A. bromeliifolia; VS: V. splendens; Gl: G. lingulata; AM: A. mertensii. N = number of plants sam pied, IR= incident radiation (%), WV = water volume extracted (ml), FPOM =fine particulate organic matter (ml after decantation in test-tubes), #Taxa= number of taxa per food web. doi:1 0.1371/journal.pone.0071735.t001

Table 2 .

 2 Bayesian Information Criterion (BIC) of the full models with four different random effects for each food-web descriptor. LocalityNegetation type' means that the vegetation types are nested within the localities. Bold characters highlight the lowest BICs. doi:1 0.1371 /journal.pone.0071735.t002

	Random effeds

Table 3 .

 3 Models evaluating the patterns of community diversity (Shannon's entropy and Simpson's evenness) and food-web structure (connectance, linkage density, nestedness) in relation to environmental variables and their interactions.

	Random

doi:1 0.1371 /journal.pone.0071735.t003

Acknowledgments

We are grateful to the members of the Laboratoire EnoironTiiffTimt de Petit Saut (HYDRECO) and to the members of the Nouragues Research Station for technical support and the use of their facilities. The English text was proofread by Andrea Yockey-Dejean and Kurtis Trzcinski. Two anonymous reviewers provided insightful comments on earlier version of this paper.

Author Contributions

Conceived and designed the experiments: RC BC CL ADJ-FC. Performed the experiments: RC BC CL AD J-FC LP. Analyzed the data: OD. Contributed reagents/materials/analysis tools: RC BC J-FC AD CL LP. Wrote the paper: OD RC CL.


F Mora

L Aveneau
email: lilian.aveneau@xlim.fr

O Apostu
email: oana.apostu@xlim.fr

D Ghazanfarpour
email: djamchid.ghazanfarpour@xlim.fr

Lazy visibility evaluation for exact soft shadows

Keywords: Analytical visibility, exact soft shadows Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Visible surface algorithms, Shadowing

This paper presents a novel approach to compute high quality and noise-free soft shadows using exact visibility computations. This work relies on a theoretical framework allowing to group lines according to the geometry they intersect. From this study, we derive a new algorithm encoding lazily the visibility from a polygon. Contrary to previous works on from-polygon visibility, our approach is very robust and straightforward to implement. We apply this algorithm to solve exactly and efficiently the visibility of an area light source from any point in a scene. As a consequence, results are not sensitive to noise, contrary to soft shadows methods based on area light source sampling. We demonstrate the reliability of our approach on different scenes and configurations.

Introduction

Shadows and especially soft shadows play an important role in visual perception as they unveil spatial relationships between objects. This makes soft shadows essential to realistic rendering. However, they are very expensive to compute since they require solving the visibility of an area light source from any visible point in a scene. The accuracy of any solution to this problem determines the soft shadows quality and, as a consequence, the realism of the scene. Dedicated algorithms were developed to reduce the computational complexity inherent to soft shadows. Interactive solutions exist, but they make crude assumptions leading to approximate or plausible soft shadows. They are forced to sacrifice the correctness of the shadows in order to reach interactive frame rates. Geometry-based approaches can lead to fine results, using silhouette edges [LLA06, FBP08] for example. While the accuracy of the visibility computation is crucial for soft shadows, the wide majority of algorithms are based on sampled light sources. As a consequence, whatever their accuracy, they are still sensitive to noise according to the sampling strategy and density.

In this paper, we present a new high quality soft shadow algorithm based on exact visibility computed from the area light source. Our method relies on an equivalence relation allowing to group lines according to the geometry they intersect [Pel04]. We build on this theoretical framework to provide a novel and practical approach to lazily compute the exact visibility from a polygon. This visibility data allows to solve efficiently and exactly the visibility of an area light source for each point to be shaded, taking advantage of the visibility coherence between neighbour points. Since the direct illumination is analytically evaluated using exact visibility, the shadows are optimal and noise-free. While exact approaches suffer from implementation and robustness issues, our method is easy to implement, scalable and numerically robust. In addition, our results show that it outperforms an optimized ray tracer.

The paper is organized as follows: The first section outlines the works related to this paper. The second section presents the geometrical and mathematical framework we use. Moreover, we demonstrate an essential theorem to our algorithm. The third section details our exact from-polygon visibility algorithm which is the basis to solve the visibility of an area light source. At last, in section 4, we test our algorithm's behaviour in different configurations. Before con-cluding, we discuss these results, we analyse some limitations and we consider several perspectives.

Previous work

Many contributions on soft shadows can be found in the literature. A survey has been presented by Hasenfratz et al. [HLHS03]. Recently, Eisemann et al. [EASW09] provide another interesting lecture. In this section, we only focus on widely used or recent physically-based algorithms able to compute very high quality or exact soft shadows. In addition, since the core of our method relies on an exact from-polygon visibility algorithm, we also focus on this topic.

Soft shadows in ray and beam tracing

Ray tracing easily supports soft shadows using shadow rays and stochastic sampling of the light sources. As a counterpart, the results are sensitive to noise. This requires increasing the number of samples in order to avoid visual artifacts. But this drastically slows down the rendering. Packet [WBWS01] and frustrum [RSH05, WIK * 06] traversal techniques improve ray-tracing efficiency. They take advantage of the spatial coherence and have efficient implementations, using SIMD instructions. Moreover, they can be designed for a specific graphic architecture [BW09]. Despite these improvements, soft shadows remain a challenging task. If the sample distribution can be improved using importance or interleaved sampling, the results remain sensitive to noise.

Beam tracing [HH84] methods are suited for soft shadows computation. A single beam is cast to the light to solve exactly its visibility. As a consequence, results are not sensitive to noise. Beam tracing is usually considered as lacking robustness, which can be improved by adaptive methods [START_REF] Ghazanfarpour | A beam tracing with precise antialiasing for polyhedral scenes[END_REF]. Overbeck et al. [ORM07] propose a robust beam-triangle intersection to handle more complex models with soft shadows.

Packets, frustrum or beam techniques take advantage of the spatial coherence from each point to be shaded. Our algorithm uses the spatial coherence from the light source. This is an important difference since the shading of the points is not an independent operation anymore.

Silhouette based soft shadows

Akenine-Möller and Assarsson [START_REF] Akenine-Möller T | Approximate soft shadows on arbitrary surfaces using penumbra wedges[END_REF][START_REF] Assarsson | A geometrybased soft shadow volume algorithm using graphics hardware[END_REF] propose the soft shadow volume algorithm on graphics hardware. It extends the shadow volume algorithm [Cro77] using penumbra wedges. A wedge is the place where the light source view can change because of a silhouette edge. Given a point inside a wedge, the relevant silhouette edge is projected onto the light to determine the light samples visibility. However, overlapping silhouette edges lead to overestimating the shadow and silhouette edges are underestimated since they are computed only from the center of the light source. As a consequence, soft shadows are approximated.

The soft shadow volume algorithm is extended to offline ray-tracing by Laine et al. [LAA * 05]. All the silhouette edges are taken into account and the algorithm deals properly with the overlapping ones. Special rules determine the light samples affected by each silhouette edges. Then a single ray is shot to recover the samples visibility. Lehtinen et al. [LLA06] remark that this algorithm can be overconservative and less efficient than a classical ray-tracer. They change the data structure used to store the silhouette edges, improving both the algorithm efficiency and its memory consumption.

Forest et al. [START_REF] Forest | Accurate shadows by depth complexity sampling[END_REF] build on the works of [AAM03] and [LLA06]. They remove the shooting ray step and provide a dedicated implementation on graphics hardware allowing several frames per second on moderate scenes using an average of 16 light samples per light sources.

Research on silhouette based algorithms has lead to very efficient solutions. However, they require triangles connectivity with consistent orientation. Since these methods are object-based, visibility computations are accurate. But they are still subject to noise because of the light source sampling. The sampling strategy has a significant impact as shown in [START_REF] Forest | Accurate shadows by depth complexity sampling[END_REF].

Other methods

Laine et al.

[LA05] change the pre-processing order of ray tracing. Instead of finding a triangle hit by a shadow ray, they search all the shadow rays intersecting a given triangle. Notice that the same idea was used to provide interactive GPU algorithms based on triangles [ED07, SEA08] or extended to silhouette edges [JHH * 09]. However, approximations are made to reach such frame rates. In this work, we also focus on the sets of lines intersecting triangles. However, we use an analytical representation of those sets.

Exact from-polygon visibility

Exact solutions to visibility problems are important to better understand its underlying nature, and design new algorithms. However, exact from-polygon (and thus from-region) visibility is very complex, because it is a four-dimensional problem [DDP02]. There are few solutions in the literature. They all rely on the same theoretical framework described by Pellegrini [Pel04] in the Plücker space, a 5D space of lines. His analysis provides a worst case complexity of O(n 4 log n) in space and time, where n is the number of triangles. This underlines the complexity inherent to the visibility in space. To remain practicable, a very special care must be given to the algorithm design. As an example, a first approach is proposed by Mount and Pu [MP99]. But their experiments do not exceed 14 random triangles, because, as stated by the authors, they reach the theoretical complexity upper bound.

The two first practicable algorithms are provided by Nirenstein et al. [NBG02] and Bittner [Bit02] (see [START_REF] Aveneau | Coherent and exact polygon-to-polygon visibility[END_REF] for details about the differences between the two solutions), and further improved by Haumont et al.

[HMN05] and Mora et al. [START_REF] Aveneau | Fast and exact direct illumination[END_REF] respectively. It is worth underlying a major difference between the two approaches initiated by Nirenstein and Bittner respectively :

• Nirenstein and Haumont design their algorithms to solve submitted to COMPUTER GRAPHICS Forum (12/2013).

whether two polygons are mutually visible. Thus the algorithm output is a boolean value. Whenever possible, the computation of the whole visibility set is avoided since the process can be terminated as soon as a visibility is found.

• Mount, Bittner and Mora algorithms aim to compute and store all the visibility information. In this case, the output is a partition of the Plücker space encoded in a Binary Space Partitioning tree.

Despite this difference, they all rely on Computing Solid Geometry in the Plücker space which is computationally very expensive. As a consequence, they are only used as a preprocess step. In addition 5D CSG operations are very complex to implement and prone to numerical instability. Thus, the reliability and the scalability of these solutions is restricted.

The soft shadow algorithm presented in this paper relies on a novel approach to capture the visibility from a polygon. It also leads to the construction of a BSP tree in the Plücker space. But in contrast to the previous works, the tree is built lazily and at run time, not as a pre-process step. Visibility is computed on-demand when and where it is required, according to the image resolution. It does not rely on any expensive and complicated 5D CSG operation. This makes the algorithm very easy to implement, efficient and computationally robust, whereas previous methods are complex and suffer from numerical instability.

Geometrical basis

In this section we introduce the geometrical knowledge underlying the work presented in this paper. In particular, we demonstrate a result on the orientation of the lines stabbing two polygons. This result is important for our algorithm.

Plücker's coordinates

The Plücker space is a five dimensional projective space, denoted P 5 , well known in computer graphics as an efficient solution for dealing with real 3D lines [Sho98]. We only recall here the properties used in this paper. Let's consider an oriented 3D line l going through two distinct points p and q of coordinates (px, py, pz) and (qx, qy, qz) respectively. The line l maps to the Plücker point denoted π l ∈ P 5 . The six coordinates of π l , denoted (l 0 , l 1 , l 2 , l 3 , l 4 , l 5 ) are defined by: l 0 = qxpx l 3 = qz pyqy pz l 1 = qypy l 4 = qx pzqz px l 2 = qzpz l 5 = qy pxqx py Any Plücker point is in bijection with its dual hyperplane:

h l (x) = l 3 x 0 + l 4 x 1 + l 5 x 2 + l 0 x 3 + l 1 x 4 + l 2 x 5 = 0
with x ∈ P 5 with coordinates (x 0 , x 1 , x 2 , x 3 , x 4 , x 5 ). Thus, any line l in the 3D space can be mapped to the Plücker space as a point π l or its dual hyperplane h l . We now recall the definition of the so-called side operator: side(l, r) = l 3 r 0 + l 4 r 1 + l 5 r 2 + l 0 r 3 + l 1 r 4 + l 2 r 5 where l and r are two 3D lines, (l 0 , l 1 , l 2 , l 3 , l 4 , l 5 ) and (r 0 , r 1 , r 2 , r 3 , r 4 , r 5 ) are respectively the coordinates of their Plücker point, π l and πr. The sign of the side operator dertermines the relative orientation of the two lines. In particular, two lines are incident (or parallel) if their side operator equals zero. We can notice that: side(l, r) = h l (πr). This leads to the geometrical interpretation of the side operator as depicted in Figure 1.

On the lines intersecting the same triangles

a πa l 0 h l 0 l 1 h l 1 l 2 h l 2 R 3 P 5
Figure 2: Left: A triangle in 3D space. l 0 , l 1 , l 2 are the lines spanning the triangle's edges. a is a line stabbing the triangle. Right: h l 0 , h l 1 and h l2 are the Plücker hyperplanes mapped from l 0 , l 1 and l 2 . πa is the Plücker point mapped from the stabbing line a. πa has a consistent orientation with respect to l 0 , l 1 and l 2 . From a geometrical point of view, πa lies at the intersection of the halfspaces induced by h l 0 , h l 1 and h l2 . Thus, these 3 hyperplanes provide an analytical representation of all the lines stabbing the triangle.

A direct application of the side operator is an easy and robust line-triangle intersection test. A line intersects a triangle if its orientation is consistent with respect to the lines spanning the triangle edges. This is depicted in Figure 2. Beyond this intersection test, we can notice that the 3 hyperplanes related to the triangle edges divide the Plücker space into cells and one of them contains all the lines (i.e. their Plücker points) stabbing the triangle, while the other cells contain the lines missing the triangle.

Pellegrini [Pel04] develops a more general approach. He uses the Plücker space as a framework to provide theoretical bounds on various problems involving lines. Let S be a set of several triangles (or convex polygons) and L S be the lines spanning the triangles edges in S. Each line in L S can be mapped to Plücker space as a hyperplane. This builds an arrangement in Plücker space: A decomposition of the space into cells by a set of hyperplanes. All the points in a same cell satisfy the following property: They all have the same sign with respect to its bounding hyperplanes. This is illustrated in Figure 3 with 2 triangles. Thus, in Plücker space, all the lines (i.e. their Plücker points) belonging to the same cell intersect the same subset of triangles in S. Notice that this subset can be empty if the lines miss all the triangles in S.

The decomposition of Plücker space into cells allows to group lines together according to the subset of triangles they intersect. This defines an equivalence relation on lines. Each cell corresponds to an equivalence class (sometimes called isotopy class [Pel91] or orientation class [CEG * 96]). To sum up, the Pellegrini approach allows an exact and analytical representation of all the sets of lines generated by a set of triangles S, using the set of lines L S .

Orientation of lines intersecting two polygons

We focus on the set of lines intersecting two convex polygons. We prove the following theorem: Theorem 1 Let A and B be two convex polygons with n and m vertices respectively. We define vvset

= {v i j , i ∈ [1, n], j ∈ [1, m]}

the set of the lines v i j defined by one vertex of A and one vertex of B. Let l be any line and q a line intersecting A and B

: side(l, x) ≥ 0, ∀x ∈ vvset ⇒ side(l, q) ≥ 0, side(l, x) ≤ 0, ∀x ∈ vvset ⇒ side(l, q) ≤ 0.
In other words, if all the vertex-to-vertex lines of two polygons have a positive (resp. negative) orientation with respect to any line l, all the stabbing lines of the two polygons will have a positive (resp. negative) orientation with respect to l. This is also illustrated by Figure 5. We refer the reader to the Appendix at the end of this paper for a detailed proof of the theorem 1.

Using this theorem, we can determine if the set of lines stabbing two polygons has a consistent orientation with respect to a given line. The visibilty algorithm presented in the next section uses this essential result.

It can be proved that the Plücker points of the lines in vvset are the vertices of the smallest convex polyhedron in the Plücker space containing the stabbing lines of A and B. A demonstration can be found in [CAF06, ACFM11], but it requires advanced knowledge in geometric algebra. For the understanding and the correctness of this work, the theorem demonstrated in this paper is sufficient.

Algorithm

Since equivalence classes are continuous sets of lines that hit/miss the same subset of triangles, they represent coherent paths through the scene independently of any viewpoint. Therefore, two lines belonging to the same equivalence class are spatially coherent. We use this property to build a Plücker space partition representing the visibility of an area light source.

Overview

We consider a convex area light source L, a triangle T , and define their occluders O as the triangles intersecting their convex hull. To represent the light source visibility, we focus on the sets of lines intersecting L and which either miss all the occluders in O, or hit at least one occluder in O. Thus, we define:

• A visible class: Any equivalence class representing a set of lines that do not intersect any occluders. • An invisible class: Any equivalence class representing a set of lines that intersect at least one occluder. • An undefined class: An equivalence class that is not yet found as visible or invisible.

Our algorithm builds a BSP tree in Plücker space, providing a hierarchical representation of the equivalence classes generated by the occluders. Each leaf represents one of these three classes. The algorithm is lazy: The BSP tree is grown on-demand depending on when and where visibility information is needed. The construction only relies on two operations: Inserting an occluder into the tree and growing the submitted to COMPUTER GRAPHICS Forum (12/2013).

tree. The former allows to find the leaves affected by an occluder in the tree, while the latter grows the tree by replacing a leaf with an occluder's equivalence classes.

l 0 l 1 l 2 h 0 h 1 h 2 o
Figure 4: Let l 0 , l 1 , l 2 be the lines defined by the edges of an occluder o, and h 0 , h 1 , h 2 their dual hyperplanes in Plücker space. We can build a BSP tree whose leaves are the four equivalence classes generated by the triangle: Three visible classes (left leaves) and one invisible class (right leaf). We note bsp(o) such a representation.

Growing a BSP tree:

Let o be an occluder, we note bsp(o) the BSP representation of the equivalence classes generated by o (see Figure 4). The tree grows each time an occluder is merged: The visible/invisible classes generated by the occluder are added to the tree. Merging o into the tree consists in replacing a leaf by the root of bsp(o). As a consequence, each inner node contains a hyperplane corresponding to an occluder's edge. Those are the 3 cases illustrated by Figure 5.

The BSP tree construction is driven by the visibility queries in order to compute only the required equivalence classes.

The following section presents this process and gives the details of our visibility algorithm. Next, we describe how it is used in our soft shadow framework.

Visibility algorithm

We The first one with a positive orientation with respect to h, the second one with a negative orientation. Figure 6 gives an illustration. Notice that the orientation of p is coherent with the orientation of h. As a consequence, if L is split by p, we can compute the two relevant polygons L ∩ p + and L ∩ p -so that they represent respectively the positive and negative subset of lines with respect to h. classes related to the visibility of L from xyz. It subdivides L into several convex fragments so that each of them belongs to a single equivalence class. Because it is lazy, it combines the query (lines 25-35) and the growing (lines 36-44) of the data structure in a single step. Algorithm 1 starts with all the occluders associated to a single undefined leaf. For each inner node, we compute the plane defined by xyz and the line stored in the node (line 26), as illustrated by Figure 6. Then, L is tested against this plane to determine the orientation of the lines stabbing L and xyz (line 27). If L lies in the positive (resp. negative) halfspace of the plane, then all the lines have a positive (resp. negative) orientation, and the algorithm continues in the left (resp. right) subtree (line 28 or 30). Otherwise, L is split against the plane and the algorithm continues recursively in both subtrees with the relevant parts of L (line 33). When a fragment reaches a leaf, two alternatives occur:

• The leaf has no occluders, therefore it is either a visible or an invisible class. In the former case, the fragment is a convex part of L which is visible from xyz, thus it is returned (line 37). Otherwise, the fragment is invisible from xyz and it is discarded. • The leaf has some occluders, the class is undefined and we cannot answer the query without further developing the tree.

In the latter case, the algorithm chooses a random occluder (line 39) ro among the occluders associated with the current leaf. This occluder is used to grow the tree by replacing the leaf by bsp(ro), the BSP representation of the four classes generated by ro (line 40). Next, the remaining occluders are inserted into bsp(ro) which is now a new part of the tree. The procedure at line 9 achieve an occluder insertion as described in Figure 5. When all the occluders have been inserted in bsp(o), the algorithm continues from the root of bsp(o) until it finds the equivalence class for the current light fragment.

Key points

This section underlines significant points of the algorithm.

Random selection of the occluders:

The algorithm's efficiency is related to the balance of the tree. To develop the BSP tree, the algorithm chooses an occluder randomly. Obviously, some choices may lead to a more balanced tree than others. However this is not predictable. In fact, we have tested different heuristics, trying to make a "good choice". Unfortunately, all of them achieved poor improvements compared to the extra computational cost. Moreover, their behaviour can be very different according to the nature of the scene. We have reached the conclusion that a random choice gives better results, and more important, it has a consistent behaviour independent of the rendered scene. This is similar to the choice of the pivot in the well known quicksort algorithm: Although a random pivot is not the optimal choice, it leads to the best performance in practice.

Conservative insertion:

If all the lines stabbing an occluder o do not have the same submitted to COMPUTER GRAPHICS Forum (12/2013).

orientation with respect to a hyperplane, o is inserted in both subtrees of the relevant node. As a consequence the insertion is conservative. This step could be computed exactly as in [HMN05,MA05]. But, as explained in section 2, this is expensive and prone to numerical errors. Our algorithm avoids these problems to remain simple and robust.

Tree growth:

Notice that all the queries will not develop the BSP tree. We can expect some queries to develop new equivalence classes, at least at the begining since the BSP tree is empty. We can expect queries to take advantage of the previous computations thanks to the visibility coherence. This is a key point of the algorithm's efficiency.

Soft shadows framework

To illustrate the efficiency and the reliability of our exact visibility algorithm, we plug it into a ray-tracing rendering software for computing high quality soft shadows. Algorithm 2 describes the process.

Algorithm 2

The following pseudocode illustrates how our visibility algorithm is plugged in a ray tracer software to analytically compute soft shadows • Using the primary rays, all the image points are grouped together with respect to the triangle they belong to. This builds a list of visible triangles (line 1). • Multiple lights are handled successively and independently (line 2) • For each visible triangle, an empty BSP tree is created and associated with its set of occluders (line 8). • Next, for each image point, the algorithm 1 is used to compute the visible parts of the light (line 10). In our framework, we consider area light sources with a uniform emission, therefore we compute direct illumination analytically (line 11) by integrating over the visible parts of the light [NN85].

This framework is designed to allow an efficient implementation. The loop order is chosen to build one BSP tree per light and per visible triangle. Since visible triangles are successively handled, BSP trees are developed successively and independently. This improves the memory coherence and avoids switching many times between BSP trees. This also limits the memory consumption since each BSP tree is deleted as soon as all its related image points are shaded. Moreover, the implementation can be easily multithreaded: A thread gets a visible triangle from the list, shades its image points and starts over until the list becomes empty. In this case, the triangle list access (line 6) has to be protected.

Implementation details

In Section 4, occluders are defined as the geometry intersecting the convex hull of a triangle and a light source. In practice, the occluders selection relies on a shaft culling approach as described in [Hai00]. Let's consider the bounding box of a visible triangle and the bounding box of an area light source. The shaft is defined as the convex hull of the two bounding boxes plus the triangle's plane support and the light's plane support. Any triangle intersecting the shaft is considered as a possible occluder. While this definition can lead to a conservative occluder set, it can be computed efficiently. At last, compared to Algorithm 1, a stack is managed to avoid recursive system calls. Except for this detail, the implementation is straightforward, it uses single floating point precision and does not use SIMD instructions.

Results

All tests were run on a 2.67 GHz Intel Core i7 920 processor with 4GB of memory. For comparison purpose, all pictures were rendered at 1280 × 720 pixels with one primary ray per pixel. Four sets of results are presented, three sets testing the global performance of our soft shadow framework, and a last one giving an insight on Algorithm 1 behaviour.

Comparisons on time and quality

The first set of results compares our method to ray-traced soft shadows both at comparable time and comparable quality. The ray tracer implementation is similar to [WBWS01]: It relies on an optimized SAH kd-tree, uses SIMD instructions to trace four rays simultaneously, and supports multithreading to render several parts of the picture in parallel. Ray-traced soft shadows are computed using groups of 4 shadow rays, and an uncorrelated stratified sampling of the area light sources. Since both our method and the raytracer support multithreading, all computations are run using 4 threads.

We use four scenes to test our method in different configurations. Despite its moderate geometrical complexity (26 673 triangles), the T-Rex scene is challenging for our approach because it presents difficult and complex shadows due to a long rectangular light source. This means that the light source visibility is complex and this is precisely what is computed by our algorithm. with the statue of Neptune (115 737 triangles) and Conference (282 873 triangles) are significant and detailed models with different kinds of shadow complexities. At last, Soda Hall (2 147K triangles) is used to test the scalability of our approach on a massive model with heterogeneous geometry. Figure 7 details these results, and presents the time spent in soft shadows computation for each method. At comparable time, ray traced soft shadows are always noisy. At comparable quality (i.e. the noise is not noticeable anymore in the stochastic shadows), our method is always faster. In addition, this is very noticeable on a complex case such as the T-rex scene, because the stochastic approach requires a very high number of shadow rays to remove almost all the visible noise. In contrast, the exact visibility algorithm produces high quality results in a few seconds.

Figure 8 presents the computation measurements (memory and time consumption) for our exact visibility algorithm. Since we compute a BSP tree in Plücker space for each visible triangle from the camera, the memory consumption varies during the process according to the building of the BSP tree. As a consequence, our results report the maximum memory load reached by our soft shadow framework. In any case, the memory footprint is low. The lazy evaluation of the visibility driven by the visibility queries allows to focus submitted to COMPUTER GRAPHICS Forum (12/2013). Figure 8: The time and memory consumption using our algorithm. The Max Size column is the maximum memory load reached during the process. The Shaft Culling column gives the time percentage spent to select the occluders. The Init column gives the time percentage spent to initialize each BSP tree. The Queries column gives the time percentage spent to query the BSP trees for shading each image point. The Total column gives the time in seconds for the whole process. MV Time column presents the result obtained using a Modified Version of our framework, where all queries are prevented from taking advantage of the others.

Memory

The last column gives the falling off factor between our framework (Time) and its modified version (MV Time).

on the equivalence classes related to the shadows and avoid the computation of useless visibility data. In addition, our framework is designed to avoid building too many BSP trees simultaneously (at most one per thread) allowing to keep the memory consumption low. The total time can be subdived in three steps: The occluder selection, the BSP tree initialization (including the computation of the hyperplanes and the vvset for each occluder), and the visibility queries used for shading the image points. We can notice that the computation time is clearly dominated by the visibility queries i.e. calls to the BSP_query function (see Algorithm 1) which is the core of our method.

Aside the comparison with the stochastic approach, these results show that we can make the most of from-polygon visibility coherence to design an efficient and robust algorithm, in contrast with previous works on this topic. As an example, we were unable to process the scenes presented in this paper with a method as described in [START_REF] Aveneau | Fast and exact direct illumination[END_REF]. The scenes are too complex for such an approach, which relies on 5D CSG operations. This becomes numerically unstable and leads to degenerate results. As expected, the sensitivity of our algorithm to the visual complexity of the light source is confirmed. For instance, the T-Rex scene required roughly the same amout of time as the Sponza and Neptune scene or the Conference scene.

About visibility coherence:

A crucial property of our approach is its ability to take advantage of the visibility coherence between image points. To test this ability, the same pictures were rendered again using a modified version of our visibility algorithm: Between each visibility query (i.e. between each BSP_query call), the related BSP tree is reset to its root node associated with its former occluder set (this additional operation is excluded from the timings). In such a case, each query is "the first one" and we prevent all queries from taking advantage of the previous ones. Figure 8 presents (two last columns of the table) the computation times compared to the timings obtained with the true version: The loss of efficiency is considerable. This demonstrates the capacity of the visibility algorithm to ben-efit from the visibility coherence, which is a key point to its efficiency.

Increasing the area of the light

The soft shadows complexity also depends on the size of the area light source. This is intrinsic to the soft shadows problem and it will inevitably affect our algorithm since it relies on the visibility coherence of the light source. If its area increases, the visibility coherence may decrease and lead to a loss of efficiency. As a consequence the second part of our tests investigates our framework's behavior when the area light source increases. These tests are run on the Conference model, selected for its significant number of triangles as well as the wide range of shadows cast in the scene. The tests start using a small area light source whose size is progressively increased until it becomes 100 times larger. For each size of light, the time and maximum memory used by our method are measured. Figure 9 sums up the results. As expected, it shows a loss of efficiency, both in time and memory, as the light source size grows (left and right graph of Figure 9). By extent, there is inevitably a critical light size where the time and the memory consumption would become a problem, in particular the memory since it is a limited resource. However, our tests show that, even with the largest area light source, we are far from such a point. In addition, further comparisons using our ray-tracer show that our approach remains fast. At comparable quality, 1024 samples are required for the largest light and the computation takes 91 seconds against 39 seconds using our algorithm (2.33 times slower). Independently of the increase in the noise, the break-even occurs for a light source 184 times larger than the smallest one (92 versus 93 seconds using 1024 samples).

Increasing the number of lights

Generally, a scene has several light sources. As detailed by Algorithm 2, our implementation supports multiple lights which are handled successively. Thus, it is interesting to test our framework behaviour in such a case. The Conference model is used again and rendered with 2 to 36 area light sources. All the lights have the same size and cast roughly the same "amount" of shadow. Figure 10 presents the results. The left graph shows that the time consumption is linear with respect to the number of area light sources. This is the expected behaviour since our framework evaluates the contribution of each light source one after another. In addition this makes the memory load independent from the number of lights, as shown on the right graph. However we can notice a significant growth of the memory load when the number of lights is more than 12. Indeed, the 12 first lights are above the conference table while the 24 other lights are mainly above the chairs, casting more complex shadows. The memory growth is independent from the number of lights in the scene, but it is coherent with the increase in the visibility complexity.

The previous set of results shows that oversized area light sources can become a limitation for our approach, mainly because of the memory consumption. The present tests give a solution: Any huge light source, even with a critical size, could always be treated as the union of several smaller lights. And the result is always guaranteed noise free, which is specialy interesting with huge light sources since they require a high number of shadow rays using a stochastic method.

Focus on the BSP_query behaviour

Previous results demonstrate the global efficiency and robustness of this work, but they do not highlight the visibility algorithm (Algorithm 1) behaviour, i.e. how the visibility data and the computational cost evolve with respect to the visibility queries. Figure 11 gives such an insight. It focuses on the construction of a single BSP tree from the Conference model. This tree was selected because it is representative of the BSP_query behaviour. It has 2 896 occluders and a significant number of queries (17 878) are performed. In particular, most of the image points are located in the umbra or penumbra, which represent the most complex cases for any visibility algorithm.

At first, the tree grows quickly because there is no visibility data and the algorithm has to develop it to answer the visibility queries. The timings show the extra computational cost required for this construction. As a second step, the tree growth slows down drastically because the previously computed equivalence classes can be re-used and only need to be completed from time to time. As a consequence the computational cost falls down. This is the global behaviour of the visibility algorithm. In addition, notice that the image points are shaded in the scanline order. This allows handling consecutive points which are likely to share the same visibility data. This is noticeable locally: An "expensive" query is al- The abscissa corresponds to the number of queries. The left ordinate is the time in milliseconds and is related to the continuous curve, while the right one is the memory consumption in MB and is related to the dashed curve. Above the graphic, a half-tone illustration of the light visibility for each query. Black means invisible while white is fully visible. ways followed by "cheaper" queries, taking advantage of the previous computational effort.

Discussions and future work

In this section we discuss our algorithm and point out some limitations and some issues we would like to address in the future. The lazy construction of the visibility data is an important feature of our algorithm. This requires to query the visibility from some points to drive the computations. On one hand, this allows our approach to fit the image resolution and to focus the computational effort where and when it is needed. On the other hand, it may be inappropriate to solve some visibility problems. For example it is not suitable to prove that the light source is invisible from a polygon or to compute exact Potentially Visible Sets. Indeed, it is not possible to query all the points on a polygon or on the PVS boundaries. However, using an adaptative sampling of the surfaces for example, our approach can still be used to query the visibility from each sample and produce an accurate solution to these problems.

In this paper, a BSP tree is built for each triangle with image points to shade and makes the most of the visibility coherence between those points. This is an easy solution to group image points. However, if the geometric resolution is very high with respect to the image resolution, it can lead to very few image points per triangle and thus, a loss of efficiency, because the visibility data may be dropped before being re-used. To overcome this problem we plane to develop an approach independent from the geometry. We are thinking about mapping image points into a Bounding Volume Hierarchy built to balance the number of image points per leaf. Next, the light source visibility could be computed per leaf, using its bounding box faces to apply our visibility algorithm. This would probably require to solve the self occlusions that may occur inside a bounding box. As a future work this is our main idea since it would solve any problems related to the image versus object resolution. For example, it could even handled micro-polygons.

We also think our implementation could be improved using a parallel computing architecture such as CUDA or OpenCL to benefit from a higher number of threads on the graphic hardware. This will require fixing some technical issues, in particular about memory management since the lazy construction of data structures requires dynamic allocations from the graphic device. When available, they are made from the global memory, which could lead to a loss of efficiency.

In this work, the visibility data is dropped as soon as it is not needed anymore. Instead, we are thinking it could be saved to be used again for other renderings. In addition, we plane to simplify and compress the visibility data. This may further improve the memory footprint and the visibility queries efficiency.

At last, we also plane to introduce depth information in the visibility data, which means grouping lines according to the first triangle they intersect. This would allow extending our approach to other problems and applications such as ambient occlusion.

Conclusion

In this paper, we have presented a new approach to compute analytically from-polygon visibility. The core of this method is an original algorithm based on equivalence classes of lines. It lazily encodes the visibility from a polygon, and allows taking advantage of the visibility coherence between successive visibility queries. In addition, the algorithm implementation is easy, and uses only basic geometric operations. As a consequence, it is fast and robust in contrast to previous works on exact from-polygon visibility.

Using this core algorithm, we have proposed a noise-free solution to analytically compute high quality soft shadows. We have tested our approach in various configurations, including different models, different sizes and different numbers of area light sources. These results demonstrated the reliability of our method.

To conclude, we hope this work illustrates that exact visibility, especially exact from-polygon visibility, is not necessarily incompatible with robustness and performance. Moreover, since visibility is a recurrent question in computer graphics, soft shadows is not the only problem that could benefit from the core of our method. As discussed in the previous section, we are interested in several issues and we intend to continue this work.

Acknowledgments: Sponza Atrium by Marko Dabrovic. Neptune model is provided courtesy of Laurent Saboret (IMATI, INRIA) from the AIM@SHAPE Shape Repository.

Appendix

To prove the theorem 1, we use a vector notation for Plücker's coordinates: Given two points p and q, the Plücker coordinates of (pq) is the vector couple ( u, v) such as:

(pq)    u = -→ pq = (l 0 , l 1 , l 2 ) v = -→ op × -→ oq = (l 3 , l 4 , l 5 )
where o is the origin in 3D space. For simplicity, we will omit the symbol , simply denoting the vector couple (u, v).

Using this notation, we now formulate the side operator as the sum of two dot products. For example, given two lines l and l ′ with Plücker coordinates (u, v) and (u ′ , v ′ ) respectively, then : side(l, l ′ ) = u • v ′ + v • u ′ . a and b define any line q stabbing A and B, and the Plücker point of q is (ab, oa × ob). Thus the relative orientation of a line q and any line l (u, v) is:

Proof

side(q, l) = u • (oa × ob) + v • ab.

Since the orientation of two lines is relative, we can assume without loss of generality that l goes through the origin of the 3D space (otherwise it is always possible to translate l).

Thus v (the last 3 Plücker coordinates of l) becomes zero and:

side(q, l) = u • (oa × ob) = u • (∑ n i α i oa i × ∑ m j β j ob j ) Since the cross product is distributive over addition: side(q, l) = u • (∑ m j ∑ n i (α i oa i × β j ob j )) Next, by compatibility with scalar multiplication: side(q, l) = u • (∑ m j β j ∑ n i α i (oa i × ob j )) At last, since the dot product is distributive over addition: side(q, l) = ∑ m j β j ∑ n i α i u • (oa i × ob j ) = ∑ m j β j ∑ n i α i (side((a i b j ), l)) Since α i ≥ 0 and β j ≥ 0, the sign of side(q, l) depends only on the sign of side(a i b j , l). In particular, if they are all positive (resp negative) then side(q, l) is positive (resp negative).

Figure 1 :

 1 Figure 1: Above: The relative orientation of two lines is given by the sign of the side operator. Under: Since side(l, r) = h l (πr), applying the side operator comes down to testing the position of the Plücker point of one line against the dual hyperplane of the other line

Figure 3 :

 3 Figure 3: Left: 2 triangles and 3 lines a, b, c in various configurations. Right: The arrangement of hyperplanes (illustrated by 6 2D lines) mapped from the 6 triangles edges. They divide the Plücker space into cells. Filled cells are set of lines intersecting at least one triangle. πa, π b and πc are the Plücker points mapped from a, b and c. They are located in the cells they belong to, according to the triangle(s) they stab. For example, π b has a consistent orientation with respect to the 6 hyperplanes since b intersects the two triangles. Its relevant cell holds all the lines intersecting the two triangles.

Figure 5 :

 5 Figure 5: 2D illustration of an insertion using Theorem 1. An occluder o intersects a subset of all the lines originating from the light L. This subset may have a negative (left example), a positive (right example) or a mixed (center example) orientation with respect to any hyperplane h from an inner node. Theorem 1 allows to define this orientation using the vvset lines, i.e. the lines defined by one vertex of o and one vertex of L. As a consequence the occluder o is inserted into the right or left or both subtrees of the node."

Figure 6 :

 6 Figure 6: The point xyz and the line h define a plane that sets apart the lines stabbing L and having a positive or negative orientation with respect to h

Figure 7 :

 7 Figure 7: From top to bottom: T-Rex, Sponza with Neptune, Conference and Soda Hall. The middle column shows the results computed with our algorithm. The left column presents the same pictures computed at comparable time and the right one at comparable quality, both using shadow rays. As an indication of the performance of our comparison method, we have rendered the same pictures under the same circumstances (same computer, 4 threads), using Mental Ray c and obtained the following results: Trex (512 samples -4"43) Sponza-Neptune (256 samples -3"52), Conference (256 samples -2"20), Soda (256 samples -3"04). These timings concern only the shadow rays.

Figure 9 :

 9 Figure 9: Increasing the light size. Left: The time consumption. Right: The maximum memory load reached during the computations. Center: A picture from our tests with the largest area light source (100 times larger than the smallest one. This corresponds to a square whose side is exactly the width of the conference table).

Figure 10 :

 10 Figure 10: Increasing the number of lights. Left: The time consumption. Right: The maximum memory load reached during the computations. Center: A picture from our tests with 36 area light sources.

Figure 11 :

 11 Figure11: Lazy construction of a representative BSP tree. The abscissa corresponds to the number of queries. The left ordinate is the time in milliseconds and is related to the continuous curve, while the right one is the memory consumption in MB and is related to the dashed curve. Above the graphic, a half-tone illustration of the light visibility for each query. Black means invisible while white is fully visible.

:

  Let A and B be two convex polygons with n and m vertices in the 3D space. Using barycentric coordinates, any point a on A (resp. b on B) can be defined by a linear combination of its vertices (a 1 , • • • , an) (resp. b 1 , • • • , bm)1, with ∀i β i ≥ 0.

  The visibility algorithm (Algorithm 1) finds the equivalence Algorithm 1 Visibility algorithm: Given a point xyz on T , the algorithm answers the query "which parts of L are visible from xyz ? L is used to drive the BSP tree construction and may be split into several fragments representing homogeneous sets of lines from xyz. The fragments reaching visible classes are the visible parts of L from xyz

	1: Node {
	2:	// if it is an inner node
	3:	Plucker hyperplane //a line from an occluder's edge
	4:	Node left, right // left and right children
	5:	// if it is a leaf
	6:	Enum class // visible, invisible or undefined
	7:	Occluder[] occluders //occluders in the leaf if it is undefined
	8: }
	9: insertOccluder(Node n, Occluder o)
	10: if n is a leaf then
	11:	if n is not an invisible class then
	12:	n.occluders ← n.occluders ∪ o
	13:	end if
	14: else if orientation( o.vvset, n.hyperplane) > 0 then
	15:	insertOccluder( n.le f t, o )
	16: else if orientation( o.vvset, n.hyperplane ) < 0 then
	17:	insertOccluder( n.right, o )
	18: else
	19:	insertOccluder( n.le f t, o )
	20:	insertOccluder( n.right, o )
	21: end if
	22:
	23: BSP_query(Node n, Polygon L, Point xyz ) return Polygons
	24: loop
	25:	while n is not a leaf do
	26:	Plane p ← makePlane( xyz, n.hyperplane )
	27:	if position( p, L ) > 0 then
	28:	n ← n.le f t
	29:	else if position( p, L ) < 0 then
	30:	n ← n.right
	31:	else
	32:	// split the light and work recursively
	33:
	34:	end if
	35:	end while
	36:	if n.occluders is empty then
	37:	return (n.class is visible) ? L : ∅
	38:	else
	39:	ro ← random occluder in n.occluders
	40:	n ← root of bsp(ro)
	41:	for o in n.occluders, o = ro do
	42:	insertOccluder( n, o )
	43:	end for
	44:	end if
	45: end loop

return BSP_query(n.le f t, L ∩ p + , xyz) ∪ BSP_query(n.right, L ∩ p -, xyz)

  The modified Sponza Atrium

	Comparable Time	Our Method	Comparable Quality
	7 seconds / 32 shadow rays	6.5 seconds	1 min 22 seconds / 512 shadow rays
	7 seconds / 32 shadow rays	7 seconds	23 seconds / 256 shadow rays
	6 seconds / 32 shadow rays	6 seconds	24 seconds / 256 shadow rays
	8 seconds / 32 shadow rays	5 seconds	38 seconds / 256 shadow rays

submitted to COMPUTER GRAPHICS Forum (12/2013).