C Gicquel

L Hege

M Minoux

W Van Canneyt

A discrete time exact solution approach for a complex hybrid flow-shop scheduling problem with limited-wait constraints

Keywords: Hybrid flow-shop scheduling, batch scheduling, limited-wait constraints, discrete-time model, mixed-integer linear programming, valid inequalities

We study a real-world complex hybrid flow-shop scheduling problem arising from a bio-process industry. There are a variety of constraints to be taken into account, in particular zero intermediate capacity and limited waiting time between processing stages. We propose an exact solution approach for this optimization problem, based on a discrete time representation and a mixed-integer linear programming formulation. The proposed solution algorithm makes use of a new family of valid inequalities exploiting the fact that a limited waiting time is imposed on jobs between two successive production stages. The results of our computational experiments confirm that the proposed method produces good feasible schedules for industrial instances.

Introduction

Hybrid flow-shop scheduling (HFS) problems arise in manufacturing environments in which a set of jobs are to be processed following the same production flow in a series of stages, each stage comprising several parallel machines. Hence, the main decisions in operating such a system consist in assigning jobs to machines and in scheduling jobs on machines at each stage of the production process while optimizing a given objective function.

As mentioned in [START_REF] Ruiz | The hybrid flow shop scheduling problem[END_REF], in the standard form of the HFS problem, machines at a given stage are identical, a machine can process only one operation at a time, each job has to be processed on exactly one machine at each stage, setup and removal times are negligible, buffer capacities between stages are unlimited and there is no restriction on the waiting time of a job between two stages.

In the present paper, we investigate a complex hybrid flow-shop scheduling problem arising from a bio-process industry. This problem significantly deviates from the standard form recalled above due to the presence of several additional constraints and assumptions to be taken into account, in particular:

• a job may be processed simultaneously by more than one machine at some stages of the production process, i.e. there exist multiprocessor tasks.

• there is no buffer capacity between stages so that if no machine is available at the next stage, a completed job is "blocked" in the previous stage and must wait on the machine which processed it.

• due to unstable intermediate products, there is a limited waiting time between operations of a given job ("no-wait" or "limited-wait" constraints).

• time-consuming setup and removal operations are required before and after the processing of a job, leading to positive setup and removal times.

The main contributions of this paper can be summarized as follows. First, we report a real-world complex hybrid flow-shop scheduling problem found in a bio-process industry, thus making a step toward reducing the gap between the theory and practice of scheduling mentioned in [START_REF] Ruiz | Modeling realistic hybrid flowshop scheduling problems[END_REF] and [START_REF] Ribas | Review and classification of hybrid flow shop scheduling problems from a production system and a solutions procedure perspective[END_REF]. Second, we present an exact solution approach based on discrete time representation and mixed-integer linear programming. The proposed solution algorithm makes use of a new family of strong valid inequalities exploiting the fact that a limited waiting time is imposed on jobs between two successive production stages. As far as we know, this is the first time the issue of valid inequalities is addressed in the context of limited waiting time constraints. Computational experiments confirm that the proposed method consistently produces good feasible schedules for industrial-size instances involving up to 35 jobs. The rest of this paper is organized as follows. We provide in section 2 a detailed review of the related literature. In section 3, we describe the industrial scheduling problem in more detail. We then develop in section 4 the proposed exact solution approach and report in section 5 numerical results based on industrial data. Finally, we provide some concluding remarks in section 6.

Literature review

Hybrid flow-shop scheduling

Recent literature reviews on the HFS problem can be found in [START_REF] Ruiz | The hybrid flow shop scheduling problem[END_REF] and [START_REF] Ribas | Review and classification of hybrid flow shop scheduling problems from a production system and a solutions procedure perspective[END_REF]. Both papers use the standard α|β|γ structure to describe the scheduling problem with respect to shop configuration (α), job constrains (β) and objective functions (γ). References are then classified according to the characteristics of the production system and/or the type of approach used to solve the optimization problem. The reader is referred to these papers for comprehensive surveys on the HFS literature. In the present literature review, we focus on those papers in which one of the complicating features mentioned in the introduction is considered.

The situation where an operation at a certain stage requires more than one machine ("multiprocessor task") is addressed among others in [START_REF] Oguz | Heuristic algorithms for multiprocessor task scheduling in a two-stage hybrid flow-shop[END_REF] and [START_REF] Ying | Scheduling multistage hybrid flowshops with multiprocessor tasks by an effective heuristic[END_REF]. The authors of [START_REF] Sawik | An exact approach for batch scheduling in flexible flow lines with limited intermediate buffers[END_REF], [START_REF] Wardono | A tabu search algorithm for the multi-stage parallel machine problem with limited buffer capacities[END_REF] and [START_REF] Wang | A tabu search heuristic for the hybrid flowshop scheduling with finite intermediate buffers[END_REF] study the case of limited buffer capacities between production stages which leads to a "blocking scheduling problem" where a completed job remains on a machine and blocks it until a downstream machine becomes available. Contributions on hybrid flow-shop scheduling under no-wait constraints can be found for instance in [START_REF] Grabowski | Sequencing of jobs in some production system[END_REF] and [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF]. We point out here that the case of limited-wait constraints has only received limited attention: noticeable exceptions can be found in [START_REF] Yang | A two-machine flowshop scheduling problem with limited waiting time constraints[END_REF], [START_REF] Su | A hybrid two-stage flowshop with limited waiting time constraints[END_REF] and [START_REF] Akkerman | Influence of capacity-and timeconstraint intermediate storage in two-stage food production systems[END_REF]. Finally, HFS with positive setup and/or removal times is studied e.g. in [START_REF] Low | Simulated annealing heuristic for flow shop scheduling problems with unrelated parallel machines[END_REF]. However, as mentioned in [START_REF] Ruiz | Modeling realistic hybrid flowshop scheduling problems[END_REF], even if there are several papers addressing realistic extensions of the hybrid flow-shop scheduling, very few papers consider several complicating features jointly. Moreover, as can be seen from the detailed survey presented in [START_REF] Ruiz | The hybrid flow shop scheduling problem[END_REF], there seems to be no previous attempt in the Operations Research literature to solve the variant of the HFS problem studied in the present paper, i.e. the HFS problem with multiprocessor tasks, zero buffer capacity, limited waiting time between consecutive operations of a job and positive setup and removal times.

There is a wide variety of solution approaches for the HFS problem. The authors of [START_REF] Ruiz | The hybrid flow shop scheduling problem[END_REF] propose a classification into three broad classes: exact methods, heuristics and metaheuristics. Exact solution approaches for the HFS problem mostly rely on problem-specific Branch & Bound algorithms where nodes correspond to partial schedules and lower bounds are computed by exploiting specific properties of the HFS problem. But, as can be seen from the surveys presented in [START_REF] Ruiz | The hybrid flow shop scheduling problem[END_REF] and [START_REF] Kis | A review of exact solution methods for the nonpreemptive multiprocessor flowshop problem[END_REF], research in this area has focused on simplified versions of the problem by considering either problems with a restricted number of processing stages (typically 2) or problems close to the standard form of the HFS. Moreover, existing Branch & Bound algorithms appear to be limited to situations where the objective is to minimize makespan or mean flow time. In the present paper, we handle a HFS problem involving a general number of processing steps, several complicating job constraints and a more complex criterion, namely total weighted tardiness. This is why we did not consider the development of an exact solution approach based on a specific Branch & Bound algorithm but rather chose to use a standard Branch & Bound algorithm with a tight mixed-integer programming linear formulation.

Batch scheduling in chemical processes

A variant of the HFS problem has been widely studied in the context of chemical process industries where the problem is often referred to as the "scheduling problem for multi-product multi-stage batch plants". Such plants usually produce similar chemical compounds by performing on each batch of product an identical sequence of processing operations defined by the product recipe. Each batch to be processed corresponds to a job in the HFS problem and has to be assigned to a production unit and scheduled on this unit at each stage of the production process.

Recent reviews on the chemical engineering literature on short-term scheduling of batch processes are provided in [START_REF] Floudas | Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review[END_REF] and [START_REF] Mendez | State-of-theart review of optimization methods for short-term scheduling of batch processes[END_REF]. One of the main features used in these papers to distinguish among scheduling problems is the so-called inventory storage policy. Namely, unlike discrete manufacturing industries, chemical industries are commonly characterized by the use of buffer tanks to store intermediate liquid or powder products between production stages. Moreover, it is usually impossible to mix in the same storage tank different batches of the same product. The authors of [START_REF] Sundaramoorthy | Modeling of storage in batching and scheduling of multistage processes[END_REF] propose a classification of storage policies in multi-stage batch plants based on the available storage capacity and the amount of time a batch can wait after its processing at a given stage before undergoing the next production operation. They identify three types of storage policies: (1) unlimited intermediate storage, (2) limited intermediate storage (i.e. there is a limited number of available storage tanks which can be either dedicated to a given type of product or shared among all products), (3) no intermediate storage (i.e. there is no intermediate storage tank but a batch, once processed by a production unit, can wait in this unit until a unit in the next stage becomes free to receive it). In terms of the time a batch can wait between two production operations, the authors of [START_REF] Sundaramoorthy | Modeling of storage in batching and scheduling of multistage processes[END_REF] distinguish between unlimited and limited waiting time policies and classify the "zero-wait" constraint as a special case of the limited waiting time policy. Contributions on multi-stage batch plants scheduling with limited intermediate storage capacity and/or zero waiting time can also be found among others in [START_REF] Moon | Mixed-integer linear programming model for short-term scheduling of a special class of multipurpose batch plants[END_REF], [START_REF] Lee | List-based treshold-accepting algorithm for zero-wait scheduling of multiproduct batch plants[END_REF], [START_REF] Liu | Scheduling multistage batch plants with parallel units and no interstage storage[END_REF], [START_REF] Liu | An effective hybrid particle swarm optimization for batch scheduling of polypropylene processes[END_REF]. In what follows, we study a multi-product multi-stage batch scheduling problem with no intermediate storage capacity and non-zero limited waiting time between operations. To the best of our knowledge, the case of non-zero limited waiting time has only been considered in [START_REF] Sundaramoorthy | Modeling of storage in batching and scheduling of multistage processes[END_REF] where, in contrast with the model investigated in the present paper, an approach based on a continuous-time formulation is proposed.

Optimization models for sequential batch scheduling can be broadly classified into two main groups based on the time representation. Discrete time models rely on a division of the scheduling horizon into a number of time intervals of identical duration and allow events such as the beginning or ending of the processing of a batch on a unit to happen only at the boundaries of these time intervals. Examples of work based on such an approach can be found in [START_REF] Kondili | A general algorithm for shortterm scheduling of batch operations-I. MILP forumation[END_REF], [START_REF] Shah | A general algorithm for short-term scheduling of batch operations-II. Computational issues[END_REF], [START_REF] Blomer | LP-based heuristics for scheduling chemical batch processes[END_REF], [START_REF] Maravelias | Minimization of the makespan with a discrete-time state-task network formulation[END_REF] and [START_REF] Burkard | Review, entensions and computational comparison of MILP formulations for scheduling of batch processes[END_REF]. One of the main advantages of discrete time models is that scheduling constraints have to be monitored only at specific and known time points. This reduces the problem complexity and makes the model structure simpler and easier to solve, particularly when limitations on renewable resources availability or on storage capacity must be taken into account. However, to achieve a suitable accuracy of the solution, it is usually needed to use a sufficiently small time interval, resulting in large-size combinatorial optimization problems. In contrast, continuous time models avoid an external discretization of the scheduling horizon and allow events to take place at any time instant. A first type of approach following this direction is based on the concept of time-slots which stands for a set of predefined time intervals with unknown durations. The production schedule is obtained by assigning a batch to be processed to each time slot on each production unit. The authors of [START_REF] Liu | Scheduling multistage batch plants with parallel units and no interstage storage[END_REF], [START_REF] Moon | New MILP models for scheduling of multiproduct batch plants under zero-wait policy[END_REF] and [START_REF] Lamba | Scheduling parallel production lines with resource constraints. 1. Model formulation[END_REF] among others use slot-based models to solve short-term batch scheduling problems. A second type of approach is based on sequence-based models where the batch sequence on each unit is given by a set of binary variables defining precedence relations between batches. Contributions on sequence-based models can be found for instance in [START_REF] Sundaramoorthy | Modeling of storage in batching and scheduling of multistage processes[END_REF], [START_REF] Cera | A mixed-integer linear programming model for short-term scheduling of single-stage multiproduct batch plants with parallel lines[END_REF] and [START_REF] Mendez | Optimal scheduling of batch plants satisfying multiple product orders with different due-dates[END_REF]. Continuous-time models usually lead to smaller combinatorial optimization problems as compared to discrete time models. However resource constraints (such as the limitations on manpower availability or on storage capacity arising in our problem) are difficult to handle with continuous time formulations. Namely, as mentioned in [START_REF] Mendez | State-of-theart review of optimization methods for short-term scheduling of batch processes[END_REF], they require "the definition of more complicated constraints involving many big-M terms, which tends to increase the model complexity and the integrality gap and may negatively impact on the capabilities of the method". This is why we propose a discrete time representation to model and solve the batch scheduling problem under study.

Problem description

We now provide a detailed description of the industrial case study which provided the initial motivation for this work. The optimization problem consists in scheduling production for a bio-process in which fermentation techniques are used.

Bio-processes are processes that use living cells or microorganisms (bacteria, yeasts, fungi...) to obtain products such as antibiotics, antibodies or enzymes. Bio-processes often involve a processing step called fermentation. Fermentation consists in placing the living organisms together with nutrients in an appropriate environment where temperature, pressure and oxygen content are controlled so that their metabolism produces the expected material. These biochemical reactions result in the formation of many undesired byproducts. Therefore, after the fermentation is completed, it is necessary to carry out recovery operations to separate the desired end product from the other residues and isolate it in its pure form. Optimization of production scheduling in such bio-processes have been studied e.g. in [START_REF] Yuan | On-line application oriented optimal scheduling for penicillin fed-batch fermentation[END_REF] for penicillin production and in [START_REF] Samsatli | An optimization based design procedure for biochemical processes. Part II: detailed scheduling[END_REF] for intra-cellular enzyme production.

In the considered case study, the production process comprises four production stages: (1) fermentation step 1, (2) fermentation step 2, (3) broth preparation, (4) recovery. The first step of fermentation consists in placing the microorganisms in a small tank (called a fermenter) where they start growing and reproducing. The cells are then transferred in a larger fermenter, containing previously sterilized nutrients and water. During this second step of the fermentation, the cells continue to grow and multiply, while producing the expected end product. When the fermentation is complete, the mixture of cells, nutrients, residues and end product, called the broth, is transferred in a preparation tank where it undergoes some transformations prior to recovery operations. The final step is the recovery where the desired end product is separated from the broth, using mainly filtration and purification techniques.

There is a variety of aspects that need to be considered when developing scheduling models for the above specific batch process. The main features are as follows:

1. This is a multi-product multi-stage sequential process. Production is represented by a predefined number of fixed-size batches which have to be processed under the same production flow according to similar recipes. There are four production stages with several identical production units at each stage. Recipes differ only with respect to the value of processing times and production unit requirements.

2. There is no intermediate storage capacity so that after being processed at a given stage, a batch is either transferred immediately to a production unit in the next stage or must wait in the current production unit until one becomes available at the next production stage. Moreover, waiting times between processing stages are limited. More precisely, in the application considered here, there are no-wait constraints between stages 1-2 and stages 2-3 and limited-wait constraints between stages 3-4.

3. There are some multiprocessor tasks at stage 3. Namely, depending on its size, a batch may require simultaneously one or two preparation tanks during the broth preparation.

4. There are sequence-independent setup and removal times at stages 1 and 2. Indeed, fermenters need to be sterilized and filled up with the appropriate cell culture medium prior to processing of a new batch and some cleaning operations are required in fermenters after the batch has been transferred to the next production stage. Moreover, there also are sequence-independent removal times at stage 3 due to the fact that during recovery operations, the volume of a batch is progressively transferred from the preparation tanks to the recovery unit. More precisely, the preparation tanks in which the batch is processed remain occupied during its recovery and are freed only when the filtration process for the entire batch has been completed.

5. The setup of a fermenter at stage 2 of the production process requires a number of difficult manual operations to be carried out by the plant staff, especially during the first period of setup. Due to manpower restrictions, it is possible to start the setup of at most one fermenter in stage 2 during any specified time period.

The batch scheduling problem can thus be stated as follows: given a set of batches with specified due-dates, we seek to determine an assignment of batches to units and a sequence of the processing of the batches in units at each stage so that all operating constraints are satisfied and total weighted tardiness is minimized. According to the α|β|γ notation used in [START_REF] Ruiz | The hybrid flow shop scheduling problem[END_REF] and [START_REF] Ribas | Review and classification of hybrid flow shop scheduling problems from a production system and a solutions procedure perspective[END_REF], the scheduling problem considered in the present paper can be noted:

F F 4 , ((P M (i) ) 4 i=1 )|block, no -wait, size j3 , ST SI , R SI |T w .

An exact solution approach

We now discuss a mathematical programming model to describe and solve the batch scheduling problem stated in section 3. We assume here that time is uniformly discretized, i.e. that the horizon is divided into planning periods of identical duration and that events such as the beginning or the end of the processing of a batch occur only at the beginning of a planning period. Moreover, we make use of the fact that identical parallel production units are available at each stage to reduce the model complexity. Namely, we do not treat the detailed assignment of batches to units at each stage s but rather use an aggregate representation of production capacity and consumption at stage s. This yields the detailed MILP formulation described in the next subsection.

MILP formulation

We wish to schedule production of a number of batches indexed b = 1...B over a scheduling horizon divided into t = 1...T periods of identical duration. Each batch is to be processed according to a production flow involving s = 1...S processing stages, with U s identical production units available at stage s. We denote M the amount of available manpower. Parameters relative to batch b are denoted as follows:

-DD b : due-date of batch b, w b : weight of batch b, k bs : number of units needed to process batch b at stage s (k bs > 1 for multiprocessor tasks), α τ bs : amount of manpower needed by batch b at stage s τ periods after the beginning of setup operations for this batch, σ bs : setup time required on the production units of stage s before batch b can be transferred to stage s, π bs : production time of batch b at stage s, ρ bs : removal time (due to cleaning or transfer operations) needed on stage s after batch b has been transferred to stage s + 1, δ bs : maximum amount of time a batch can wait after completion of its processing at stage s before being transferred to stage s + 1. We use the following binary decision variables:

y t bs =
1 if processing of batch b at stage s starts in period t 0 otherwise.

Using these notation, the batch scheduling problem can be formulated as follows:

min B b=1 w b T t=DD b -π bs +1 (t + π bS -DD b )y t bS (1)
T t=1 α τ bs y t-τ +1+σ bs b,s

y t bs = 1 ∀b, ∀s (2)
B b=1 k bs t+σ bs τ =1 y τ bs - t-ρ bs τ =1 y τ b,s+1 ≤ U s ∀s = 1...S -1, ∀t (3)
B b=1 k bS t+σ bS τ =t-π bS -ρ bS +1 y τ bS ≤ U S ∀t (4)
≤ M ∀t (5)
y t bs + t+π bs τ =1 y τ b,s+1 ≤ 1 ∀b, ∀s, ∀t (6)
y t bs = y t+π bs b,s+1 ∀(b, s) s.t. δ bs = 0, ∀t (7)
y t bs + T τ =t+π bs +δ bs y τ b,s+1 ≤ 1 ∀(b, s) s.t. δ bs > 0, ∀t (8)
y t bs ∈ {0, 1} ∀b, ∀s, ∀t (9)
The objective, minimizing the total weighted tardiness, is expressed by [START_REF] Ruiz | The hybrid flow shop scheduling problem[END_REF]. Constraints [START_REF] Ruiz | Modeling realistic hybrid flowshop scheduling problems[END_REF] ensure that processing of batch b at stage s starts only once during the scheduling horizon. Constraints (3)-( 4) are production capacity constraints: they state that the number of units required at stage s in period t cannot be greater than the total number of units available U s . Note that a batch b requires k s production units at stage s over the set of periods between the beginning of setup operations (i.e. σ bs periods before its processing at stage s starts) and the end of removal operations (i.e. ρ bs periods after it has been transferred to stage s + 1 or its processing at stage S is finished). Constraints [START_REF] Ying | Scheduling multistage hybrid flowshops with multiprocessor tasks by an effective heuristic[END_REF] ensure that, in each period, the total demand on manpower is bounded by the available amount of manpower. Constraints [START_REF] Sawik | An exact approach for batch scheduling in flexible flow lines with limited intermediate buffers[END_REF] guarantee that every batch b undergoes operations following the prescribed recipe, i.e. that processing of batch b at stage s + 1 cannot begin before processing at stage s is completed. Constraints [START_REF] Wardono | A tabu search algorithm for the multi-stage parallel machine problem with limited buffer capacities[END_REF] guarantee that no-wait restrictions are respected. Similarly, [START_REF] Wang | A tabu search heuristic for the hybrid flowshop scheduling with finite intermediate buffers[END_REF] guarantee that the maximum waiting time between operations at stages s and s + 1 is respected. The binary character of variables y t bs is enforced by (9).

Valid inequalities

The formulation proposed in subsection 4.1 enables us to solve exactly only small instances of the scheduling problems. A possible explanation for this lies in the observation that the linear relaxation of (1)-( 9) only provides a poor approximation to the exact optimal integer solution values. In order to circumvent this difficulty, we investigate a way of strengthening the formulation based on the use of a new family of valid inequalities. More precisely, we exploit a feature of the considered scheduling problem which has received limited attention in the literature, namely the presence of limited-wait constraints, to derive the following valid inequalities.

Proposition 1

All solutions of (1)-( 9) satisfy:

∀(b, s) s.t. δ bs > 0, ∀(t 1 , t 2 ) s.t. t 2 ≥ t 1 t 2 τ =t 1 y τ bs + t 1 +π bs τ =1 y τ b,s+1 + T τ =t 2 +π bs +δ bs y τ b,s+1 ≤ 1 (10)
Moreover, (10) can be interpreted as clique constraints associated with maximal cliques in an underlying incompatibility graph.

Before proceeding to the proof, we briefly explain the idea underlying [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF]. We consider a batch b at stage s where the waiting time before transfer to stage s + 1 is assumed to be bounded by δ bs . Inequalities [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF] state that starting the processing of batch b at stage s within the time interval [t 1 , t 2 ] is incompatible with starting the processing of batch b at stage s + 1 before period t 1 + π bs or after period t 2 + π bs + δ bs . In other words, if the starting date of processing of batch b at stage s is fixed within time interval [t 1 , t 2 ], there is a rather tight time window within which processing of batch b at stage s + 1 may start, namely [t 1 + π bs + 1; t 2 + π bs + δ bs -1].

Proof

We now provide the proof for proposition 1. We arbitrarily choose a batch b and a stage s such that δ bs > 0 as well as two time periods (t 1 , t 2 ) s.t. t 2 ≥ t 1 and show that the corresponding inequality of type ( 10) is valid. The proof relies on the fact that there are a number of explicit or implicit binary exclusion constraints in formulation (1)-( 9) which corresponds to the edges of an underlying incompatibility graph. We first derive from formulation (1)-( 9) the following binary exclusion constraints:

y t bs + y τ bs ≤ 1 ∀t, ∀τ (11)
y t b,s+1 + y τ b,s+1 ≤ 1 ∀t, ∀τ (12)
Indeed, ( 11) and ( 12) are induced by ( 2), ( 13) are induced by ( 6) and ( 14) by ( 8). Now we consider the incompatibility graph G = (V, A) defined as follows:

-there is a node v = (s, t) ∈ V for each binary decision variable y t bs , t = 1..T and a node v = (s+1, t) ∈ V for each binary decision variable y t b,s+1 , t = 1..T . -there is an edge a ∈ A between any two nodes if the corresponding decisions are incompatible, i.e. if in system ( 11)-( 14), there is a binary exclusion constraint linking the corresponding decision variables.

A set C ⊂ V is called a clique if each pair of nodes in C is connected by an edge. A maximal clique is a clique which is not properly contained in another clique. Each maximal clique in G gives rise to a valid inequality called a maximal clique constraint stating that the sum of corresponding binary variables should be less than or equal to 1. We consider the node set C s = {(s, t 1 ), ..., (s, t 2 )}: thanks to constraints [START_REF] Yang | A two-machine flowshop scheduling problem with limited waiting time constraints[END_REF], we know that there is an edge between any two nodes of C s , i.e. C s is a clique. Similarly, the node set C s+1 = {(s + 1, 1), ..., (s + 1, t 1 + π bs ), (s + 1, t 2 + π bs + δ bs ), ..., (s + 1, T )} is a clique (see constraints [START_REF] Su | A hybrid two-stage flowshop with limited waiting time constraints[END_REF]). We first show that the node set C s ∪ C s+1 is a clique. Namely, let v = (s, τ ), τ ∈ [t 1 ; t 2 ] be an arbitrarily chosen node in C s and v ′ = (s + 1, θ), θ ∈ [1; t 1 + π bs ] ∪ [t 2 + π bs + δ bs ; T ] be an arbitrarily chosen node in C s+1 : -either θ ∈ [1; t 1 + π bs ] and there is an edge between v and v ′ thanks to constraints [START_REF] Akkerman | Influence of capacity-and timeconstraint intermediate storage in two-stage food production systems[END_REF], -or θ ∈ [t 2 + π bs + δ bs ; T ] and there is an edge between v and v ′ thanks to constraints [START_REF] Low | Simulated annealing heuristic for flow shop scheduling problems with unrelated parallel machines[END_REF]. Thus, there is an edge between any pair of nodes (v, v ′ ) such that v ∈ C s and v ′ ∈ C s+1 : C s ∪ C s+1 is a clique. We now show that the node set C s ∪ C s+1 is a maximal clique. Namely, let us consider a node w not included in C s ∪ C s+1 : -either w = (s, τ ) with τ < t 1 : C s ∪ C s+1 ∪ {w} is not a clique as w is not linked to node (s + 1, τ + π bs + 1) which belongs to C s ∪ C s+1 , -or w = (s, τ ) with τ > t 2 : C s ∪ C s+1 ∪ {w} is not a clique as w is not linked to node (s + 1, τ + π bs + δ bs -1) which belongs to C s ∪ C s+1 , -or w = (s + 1, τ ) with τ ∈ [t 1 + π bs + 1; t 2 + π bs + 1]: C s ∪ C s+1 ∪ {w} is not a clique as w is not linked to node (s, τπ bs -1) which belongs to C s ∪ C s+1 , -or w = (s+1, τ ) with τ ∈: [t 1 +π bs +δ bs -1; t 2 +π bs +δ bs -1]: C s ∪C s+1 ∪{w} is not a clique as w is not linked to node (s, τπ bsδ bs + 1) which belongs to C s ∪ C s+1 . As a consequence, C s ∪ C s+1 is not contained in another clique, i.e. it is a maximal clique providing the multiple exclusion constraint [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF].

We note that the maximality property is important since this guarantees that inequalities [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF] are the strongest possible clique constraints based on ( 11)- [START_REF] Low | Simulated annealing heuristic for flow shop scheduling problems with unrelated parallel machines[END_REF]. The computational results displayed in section 5 will confirm that significant strengthening of the linear relaxation is obtained when adding these constraints to the model.

Cut & Branch algorithm

The number of valid inequalities of type [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF] grows very fast with problem size. Namely, for each batch b and each stage s such that δ bs > 0, there are 2 T valid inequalities of type [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF]. Therefore, it is not possible to include all of them in the formulation. They can however be generated as needed according to a cutting-plane strategy. This is what we did in our implementation of the solution procedure. More precisely, we devised the following cutting-plane algorithm (CPA) which we use at the root node of the Branch & Bound search tree to generate violated valid inequalities of type [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF] and add them to the initial formulation ( 1)-( 9).

Algorithm (CPA)

Step 1. Solve the linear relaxation of the problem and denote ỹ the optimal solution of the relaxed problem.

Step 2. For each batch b and each stage s such that δ bs > 0, -for t 1 = 1..T , for t 2 = 1...T , compute:

LHS(t 1 , t 2 ) = t 2 τ =t 1 ỹτ bs + t 1 +π bs τ =1 ỹτ b,s+1 + T τ =t 2 +π
) = max (t 1 ,t 2 ) LHS(t 1 , t 2 ). -if V max > 1,
1) add the valid inequality corresponding to (t max

1 , t max 2 ) to formulation (1)-(9). 2) let LHS(t max 1 , τ ) = LHS(τ, t max 2 ) = 0 for τ = 1...T . 3) let nc = nc + 1.
-else, let test = 0.

Step 3. If no violated inequalities are found in step 2, stop. Otherwise go back to step 1.

Algorithm (CPA) is intended to generate at each iteration the K most violated inequalities of type [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF] for each batch b and each stage s such that δ bs > 0. However, preliminary computational experiments showed that, if nothing is done to prevent it, there is a rather high probability that a large number of similar cutting planes will be generated. This is explained by the fact that a valid inequality of type [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF] is often violated by the current continuous solution due to the values of a small subset of the T + δ bs -1 variables involved in its expression.

For instance, we often have the case where, for the valid inequality corresponding to (t max 1 , t max ) > 1 due to the fact that

t max 2 τ =t max 1 y τ bs + t max 1 +π bs τ =1 y τ b,s+1 > 1 and T τ =t max 2
+π bs +δ bs y τ b,s+1 = 0. In this situation, any valid inequality (t max 1 , θ 2 ) such that θ 2 > t max 2 will also be violated with a value LHS(t max ) and might belong to the K cuttingplanes generated during the current iteration of algorithm (CPA). This will lead to the unnecessary generation of a large number of similar cutting planes, each one having a rather low contribution to closing the integrality gap.

1 , θ 2 ) = LHS(t max 1 , t max
To avoid this, during a given iteration of (CPA), if a violated valid inequality is generated for the pair (t 1 , t 2 ), we prevent the algorithm from generating any other cutting plane corresponding to (t 1 , τ ) or (τ, t 2 ). This is achieved by setting to 0 the corresponding entries in table LHS.

Computational results

We discuss the results of some computational experiments carried out to evaluate the impact of the formulation enhancements proposed in subsection 4.2 and to estimate the size of the largest instances of the problem that could be solved to optimality within reasonable computation time limit.

Test instances

We use the data from our industrial case study to create instances of various size.

Time is discretized in periods of identical duration (typically 4 hours) and the scheduling horizon comprises 2 to 4 weeks. The plant configuration is the same for all instances: there are S = 4 processing stages with U 1 = 2, U 2 = 5, U 3 = 5 and U 4 = 1.

The plant produces P = 6 different kinds of end products, i.e. there are P = 6 different recipes. Each recipe p corresponds to a given set of values for setup, processing and removal times (σ s , π s , ρ s ) as well as for production unit requirement (κ s ). Table 1 provides a detailed description of the recipes where all time values are expressed as integer multiples of time periods.

There are no-wait constraints after stages 1 and 2 (i.e. ∀b, δ b1 = δ b2 = 0) and limited-wait constraints after stage 3 (i.e. ∀b, δ b3 = 9). Manpower availability is restricted to M = 1 team and manpower consumption arises

p σ 1 π 1 ρ 1 κ 1 σ 2 π 2 ρ 2 κ 2 σ 3 π 3 ρ 3 κ 3 σ 4 π 4 ρ 4 κ 4 1 2 5 1 1 3 24 2 1 0 3 4 2 0 4 0 1 2 1 1 1 1 2 4 2 1 0 2 2 1 0 2 0 1 3 2 5 1 1 2 19 2 1 0 2 2 2 0 2 0 1 4 2 2 1 1 2 15 3 1 0 2 2 1 0 2 0 1 5 2 2 1 1 3 13 2 1 0 2 4 1 0 4 0 1 6 2 2 1 1 3 15 2 1 0 2 2 1 0 2 0 1
only during the first period of setup operations at stage 2 of the process. We thus have: α τ bs = 1 if s = 2 and τ = 1, α τ bs = 0 otherwise. The instances tested differ with respect to the following characteristics:

• the length of the scheduling horizon: T ∈ {100, 120, 140, 160}.

• the number of batches involved in the scheduling problem: B ∈ {20, 25,

• the batch parameters. For each batch b, a recipe p is randomly chosen from a discrete uniform DU(1, P ) distribution. Due-dates occur every 20 periods, i.e. ∀b, DD b ∈ {20, 40, 60, 80, 100, 120, 140, 160}. For each batch b, a due-date DD b is randomly chosen from a discrete uniform distribution over the set of possible due-dates and a weight w b is randomly chosen from a discrete uniform DU(1, 2) distribution.

We generated two sets of instances:

• a set (set A) of 20 instances of various size to get a first evaluation of the proposed solution approach,

• a larger set (set B) of 50 instances involving 30 batches and 140 periods to confirm the conclusions of the first part of the computational study.

Table 2 provides for each set of instances the number of planning periods in the scheduling horizon T , the number of batches B as well as the resulting number of binary variables V ar and constraints Const in the initial formulation (1)-( 9).

Results for set A instances

We solved each instance of set A with a standard MILP software (CPLEX 11.1) using either the initial formulation (1)-( 9) or the strengthened formulation obtained after running algorithm CPA. All tests were run on a Pentium 4 (2.8 GHz) with 504 Mo of RAM, running under Windows XP. We used the default settings of CPLEX MILP solver. This means that some cutting planes, among which are clique cuts, cover cuts and Gomory fractional cuts, are added automatically to the model.

The computational results obtained with the two formulations are displayed in table 3. We provide:

• Cut CP A : the number of violated valid inequalities of type [START_REF] Wang | No-wait flexible flowshop scheduling with no-idle machines[END_REF] generated by algorithm CPA and added to the strengthened formulation.

• Cut ST D : the number of cutting-planes generated by CPLEX standard cutting-plane algorithm at the root node of the Branch and Bound tree.

• Gap 0 : the initial gap, i.e. the relative difference between the lower bound provided by the linear relaxation of the problem and the best integer solution found after 4 hours of computation. We use the value of the lower bound obtained after the cutting-plane generation carried out by CPLEX algorithm at the root node of the Branch & Bound tree has stopped.

• Nodes: the number of nodes of the search tree explored before a guaranteed optimal solution is found or the time limit of 4 hours is reached.

• CP U IP : the time in seconds needed to find a guaranteed optimal solution when one has been found. CP U IP comprises the time needed to strengthen the formulation using algorithm CPA.

• Gap F : the gap between the best lower bound found and the best integer feasible solution obtained after 4 hours of computation.

Results from table 3 show that using the strengthened formulation enables us to significantly enhance the efficiency of the Branch and Bound procedure embedded in CPLEX. Namely:

• for small size instances (I1-I10), the average computation time to obtain a guaranteed optimal solution is significantly reduced from 2368s with the initial formulation to 850s with the strengthened formulation.

• for industrial size instances (I11-I20), the number of instances for which a guaranteed optimal solution could be found within 4 hours of computation is increased from 4 with the initial formulation to 7 with the strengthened formulation. Moreover, the average remaining gap after 4 hours of computation is decreased from 6.5% with the initial formulation to 2.7% with the strengthened formulation.

The main explanatory factor for this lies in the observation that the lower bounds provided by the strengthened formulation appears to be stronger that the ones provided by the initial formulation. Indeed, the initial gap Gap 0 is on average 8.5% with the initial formulation as compared to 6.1% with the strengthened formulation. As a consequence, the number of nodes explored by the Branch and Bound procedure before a guaranteed optimal solution is obtained or the computation time limits are reached is divided on average by a factor of 4.6 (from 3010 to 655). However, in 5 out of the 20 instances of set A (A4, A5, A6, A14 and A20), the computation time needed to obtain a guaranteed optimal solution is greater with the strengthened formulation than with the initial formulation. This might be explained by the fact that the number of cutting planes generated by algorithm CPA is much larger that the number of cutting planes generated by the standard CPLEX algorithm. As a consequence, the linear program to be solved at each node of the Branch & Bound tree is larger and its resolution requires more time. In most cases, this increase of computation time at each node of the Branch & Bound tree is compensated by a significant decrease in the number of nodes to be explored before optimality is reached. But for a small subset of instances (especially the smaller ones), this compensation does not occur and the computation time with the strengthened formulation is larger than with the initial formulation. To make sure that the formulation strengthening obtained thanks to algorithm CPA has a statistically significant positive impact on the computation time and solution quality, we carried out additional computational experiments. We considered a larger set of industrial size instances and used a statistical test to check wether the observed improvements are significant.

Results for set B instances

The second set of instances comprises 50 industrial size instances involving 30 batches and 140 periods. The results obtained with the initial and strengthened formulations are displayed in Table 4 and5. We first provide:

• #Feas: the number of instances for which a feasible solution could be obtained within 4 hours of computation,

• #Opt: the number of instances for which a guaranteed optimal solution could be obtained within 4 hours of computation.

Moreover, for each instance, we computed:

• ∆Gap 0 = Gap 0 (IF )-Gap 0 (SF ) the difference between the initial gaps obtained with the initial and the strengthened formulations,

• ∆Nodes = Nodes(IF ) -Nodes(SF ) the difference between the numbers of nodes explored with the initial and the strengthened formulations,

• ∆CP U IP = CP U IP (IF )-CP U IP (SF ) the difference between the computation times with the initial and the strengthened formulations,

• ∆Gap F = Gap F (IF ) -Gap F (SF ) the difference between the final gaps obtained with the initial and the strengthened formulations.

For each of these performance measures, we provide in table 5 the mean, minimum and maximum values.

Results from table 4 confirm that using the strengthened formulation enables us to significantly enhance the efficiency of the Branch and Bound procedure embedded in CPLEX. Namely:

• the number of instances that could be solved to proven optimality is doubled while using the strengthened formulation, • the mean decrease in computation time (CP U IP ) obtained while using the strengthened formulation is 3043s,

• the mean decrease in the final gap (Gap F ) obtained while using the strengthened formulation is 12.45%.

Moreover, we conducted a t-test (see e.g. [START_REF] Crawley | Statistics: an introduction using R[END_REF]) to check whether the observed mean values of ∆CP U IP and ∆Gap F are statistically significantly different from 0.

The results of these tests indicate that:

• the mean decrease in the computation time obtained while using the strengthened formulation belongs to the interval [1834; 4252] with a 95% probability.

• the mean decrease in the final gap obtained while using the strengthened formulation belongs to the interval [4.9%; 19.9%] with a 95% probability.

This confirms that the use of algorithm CPA to strengthen the initial formulation leads to a statistically significant improvement of the computational efficiency of the Branch and Bound procedure embedded in CPLEX solver.

Conclusion

We investigated a scheduling problem arising from a bio-process industry. The problem is a complex variant of the hybrid flow-shop problem in which there are multiprocessor tasks, zero buffer capacity, limited waiting time between consecutive operations of a job and positive setup and removal times. We proposed a mixed-integer linear programming formulation based on a discrete time representation. We then derived a new family of strong valid inequalities exploiting the fact that a limited waiting time is imposed on jobs between two successive production stages and used these valid inequalities to devise a Cut & Branch algorithm. The results of our computational experiments confirm that the proposed method consistently provides good feasible schedules for industrial size instances with a reasonable computation effort.

Among the possible research directions suggested by the present work, it might be worth investigating the development of a heuristic solution approach so as to be able to provide good feasible schedules for industrial size instances within a shorter computation time.

σ

  bs +π bs +ρ bs τ =1

  bs +δ bs ỹτ b,s+1 . -let nc = 0 and test = 1. -while nc < K and test = 1, -look for V max = LHS(t max 1 , t max

2

 2

2 )

 2

2

 2

2

 2

Table 1 :

 1 Data relative to recipes

Table 2 :

 2 Data relative to instances

	Instances T	B	Var Const
	A1-A5	100 20 4000	4580
	A6-A10	120 25 6000	6640
	A11-A15 140 30 8400	9340
	A16-A20 160 35 11200 12320
	B1-B50	140 30 8400	9340

Table 3 :

 3 Results for set A instances Gap 0 Nodes CP U IP Gap F Cut CP A Cut ST D Gap 0 Nodes CP U IP Gap F

			Initial formulation			Strengthened formulation
	Cut ST D A1 172	8.0	1796	577	0	1508	8	6.5	434	493	0
	A2	175	16.7	2460	1148	0	1258	9	14.1	817	987	0
	A3	42	0	0	290	0	1292	33	0	0	289	0
	A4	109	0	0	119	0	1349	16	0	0	409	0
	A5	247	1.0	3	120	0	1621	5	0.5	0	166	0
	A6	95	0	0	928	0	1284	17	0	0	1391	0
	A7	244	7.3	3494	3864	0	2000	12	5.9	684	2162	0
	A8	199	3.1	469	790	0	1852	10	2.4	6	403	0
	A9	202	2.0	472	3780	0	1573	13	1.4	24	512	0
	A10	374	4.4	27554 12063	0	2416	11	2.9	836	1692	0
	A11	214	4.5	1059	8342	0	2253	33	2.9	114	1741	0
	A12	598	6.3	7421	14400	5.2	3766	4	1.8	1879	5833	0
	A13	249	3.5	207	1600	0	2293	14	2.2	1	559	0
	A14	138	48.5	1096	3323	0	1208	16	47.8	862	5472	0
	A15	401	12.4	9103	14400	11.6	2507	7	11.9	2950	14400	9.4
	A16	451	8.0	2848	14400	7.4	4418	18	2.6	3703	8143	0
	A17	679	14.6	1554	14400	13.5	4720	9	12.1	373	14400	11.9
	A18	411	17.6	618	14400	16.8	3115	15	5.8	420	14400	5.2
	A19	334	11.1	60	14400	10.8	2323	27	1.5	5	7314	0
	A20	327	0	0	12356	0	2608	34	0	10	12765	0

Table 4 :

 4 Results for set B instances

		Initial formulation	Strengthened formulation
	#Feas	45	49
	#Opt	10	20

Table 5 :

 5 Results for set B instances

		Mean value Min/Max values
	∆Gap 0	11.9	[-0.05; 98.6]
	∆Nodes	2952	[0; 11345]
	∆CP U IP	3043	[-2150;13181]
	∆Gap F	12.45	[-1.61;100]

Corresponding author. E-mail: celine.gicquel@ecp.fr . Tel: (+33) 1 41 13 16 30. Fax: (+33) 1 4113 12 72.


D Murray

J Baril

X Granier

Shape Depiction for Transparent Objects with Bucketed k-Buffer

Keywords: (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing, and texture

Figure 1: Enhancing surface with Mean Curvature Transparency and Bucketed k-Buffer. On the left of the image is classical rendering (OIT)

(a)). Observe how object details on each surface is highlighted to provide a better understanding of the composition of the object. The model runs on our system at 64 FPS with OIT and 24 FPS with our solution, for a resolution of 1000 by 1000 pixels and 16 layers of transparency. Model courtesy of dessindus.blogspot.fr.

Introduction

The simulation of a lighting model for visualizing threedimensional scenes is an essential tool for understanding the content of the scene. While current real-time rendering techniques are convincing for everyday objects, some models suffer from lack of real landmark and thus, from grasp of the shape. Actually, transparency is a common tool for analysing complex object subtleties, e.g. Computer-Aided Design. However, with raster-based applications, ambiguities from using classical alpha-blending are two-fold: (i) a correct composition involves to render the scene back-to-front; (ii) composition of multiple non-related fragments introduces a bias in shape perception. The back-to-front composition issue (i) is solved by Order-Independant Transparency (OIT) methods. However, the second one (ii) has no obvious solution.

In this paper, we provide customizable tools that help to work around shape perception ambiguities during alpha-blending. Our predicate is that shape perception relates on expressive rendering methods. Among the large set of techniques available in this domain, we are particulary interested in shape depiction algorithms that aim at stylizing the rendering of a shape according to its features. Basically, these features can be classified into three categories: surface orientation, surface curvature and inflection points. While the first one can be considered as supplied during rendering process, other ones should be computed during the rendering in submitted to Eurographics Symposium on Rendering -Experimental Ideas & Implementations (2016) order to benefit from screen-space formulation. Screen-space approaches for feature extraction have proven useful for exploring dense meshes obtained by photogrammetry using Radiance Scaling [GVP * 12], with applications in Cultural Heritage. We consider this as we focus on the exploration of complex models (e.g. CAD), which require efficiently computing curvature at a scale adapted to the current viewpoint. The computation of surface properties is linked to derivative considerations: curvature is a derivative of the normals and we need a derivative of the curvature to find inflection points. Using discrete differential geometry operators, the derivative computation in screen-space involves to have access to the neighborhood of a point, i.e. a fragment located at a pixel position. For opaque object, depth-test operation is enough to eliminate unnecessary fragments. Thus, with one fragment per pixel, finding neighbor information is straightforward. For transparent objects, a pixel may contain multiple fragments. Since a differential operator must be applied per fragment to perform shape depiction locally, the complexity of using such an operator is overweighted by the complexity of finding neighbor information. Thus, our main contribution is two-fold. First, we introduce the Bucketed k-Buffer, a data structure suited to find fragment's neighborhood efficiently. Second, we present a customizable stylisation that modulates shading and transparency per fragment, according to its curvature and its depth, which can be adapted for various kinds of application. * 13] proposed an altered A-Buffer, in which fragments are connected to their neighboring fragments. According to the connexity and the logical depth of a fragment, transparency is modulated such that creases are emphasized to nearly opaque, smoothly decreasing to nearly transparent otherwise. Günther et al. [GSE * 14] propose an enhancement with surface patches to make part of a surface transparent if an object is behind. However, to our knowledge, there is no method that computes screenspace geometric features per fragment to modulate shading and transparency in real-time. Based on a tailored data structure that allows such computation, our method provides a way to address existing shape depiction techniques on transparency shapes.

Related Work

Solution overview

The goal of the presented method is to enhance shape depiction for transparent object. For such intention, our process is composed of three main steps as shown in Figure 2: 1. Scene discretization into a screen-space representation. and data. B is used to group fragments belonging to the same shape per bucket. Details of the BKB's data structures are presented in Section 4. Note that the BKB contains geometric and color data per fragment. The second step is a stylization process which extracts geometric features from the BKB (i.e. based on the normals) to enhance shape depiction (i.e. by modulating shaded color). The stylization is customizable by providing user-defined transfer functions indexed by the fragment's depth and curvature (see Section 5). The last step consists in reordering fragments back-to-front in order to realize a correct alpha-blending.

Every stage of our system is performed in real-time on modern graphics hardware. Results are presented in Section 6.

Bucketed K-Buffer for efficient neighbor query

Leaving aside commutative OIT which excludes per fragment operations, we focus on per pixel fragment list. We seek a way to efficiently access a fragment's neighborhood. First, we introduce the concept behind our solution, that is, the elaboration of the Bucketed k-Buffer. Then, we provide several practical details about the method.

Definition of the structure

Our structure relies on a k-Buffer basis. Indeed, with current GPU, Depth Peeling is obsolete. A-Buffer relies on a linked list which implies unbound memory. Moreover, all fragment information is stored in the same buffer. For a better flexibility on data structure, the dual pass k-Buffer [START_REF] Kubisch | Order independent transparency in opengl 4[END_REF] technique is well-suited. Indeed, by decorraleting depth and color, we can manipulate only depth or color information separately. Furthermore, storing only the k nearest fragments in k layers ensures a total control over memory occupancy. Those two points are key elements for the choice of our structure.

In such structure, a naive approach to find the best neighbor candidate is to iterate through each neighboring pixel's fragment list. However, this implies costly memory access. Vardis et al. [START_REF] Vardis | A multiview and multilayer approach for interactive ray tracing[END_REF] use buckets to decompose a fragment list such that fragments are grouped by depth range. In the same way, we propose to reorganize the fragment lists so that fragments belonging to a same object are regrouped in a bucket. That way, only fragments in the corresponding bucket needs to be traversed, thus avoiding many memory access. Using such partition, finding a one-ring neighborhood would go down, in theory, from a complexity of O(k 2 ) per pixel to O(kn), where k is the number of fragments per pixel and n is the number of fragments in a bucket, so we should have a gain as soon as n < k.

The proposed algorithm shares the same concepts as the dualpass k-Buffer in the sense that we use two geometry passes to store depth and data in separate buffers. However, in order to facilitate access to neighbors, our data structures differ in the data we store and we use additionnal computing passes to arrange the data according to our needs. The remainder of this section describes the pipeline to build the BKB data structure in a sequential order. In the following, W and H denote respectively the width and the height of the rendering window. Z insertion Such as the dual pass k-Buffer [Kub14], we store fragment's depth per pixel in the first geometry pass. Relying on the benefit of such approach, the n nearest fragments from the camera, n ≤ k, are stored and ordered by depth in the buffer Z. For each fragment, we keep not only its depth (z) but also its shape ID (b) (see Figure 3 (1)), where b is a integer representing the identifier of the shape where the fragment is from. Thus, the Z data structure is of size

z 0 , b 2 z 1 , b 1 z 2 , b 2 z 3 , b 0 b 0 , z 3 b 1 , z 1 b 2 , z 0 b 2 , z 2 b 0 , 0, 0 b 1 , 1, 1 b 2 , 2,
W × H × k, containing pairs (z, b).
Z reordering During this pass, we reorder per pixel fragment pairs of the Z data structure along b instead of z (see Figure 3 (2)). Sorting fragment pairs in this manner allows grouping fragments that belong to the same shape consecutively in order to gather them into buckets.

Bucket creation Per pixel groups of fragments, which share the same shape ID, are implicitly defined by the reordering of Z in the previous pass. However, in order to have a direct access to the range of these groups, the B data structure is built in a computing pass. Such groups are defined as buckets: a bucket contains the shape ID (b, i.e. a bucket ID) and the indices (i, j) of the first and last fragments of a group. The bucket buffer B contains W × H × k triplets (b, i, j) (see Figure 3 (3)).

Since n fragments are stored in Z during the first pass, n ≤ k, no more than k buckets are stored per pixel. However, storing bucket information in a table indexed by b can lead to indices that are out of range if b > k. Thus, B is considered as a per pixel hash table in which each bucket is indexed by H(b), where H is a hash function defined as:

H(b) = b mod k
Using such hash function may lead to collision in the hash table. Thus, during insertion of a bucket b in buffer B, we check that the slot at index H(b) in the hash table is free. If the slot is unoccupied, the bucket information is stored at this index. If the slot is occupied, b is incremented until the index H(b) points to a free slot. Note that a free slot is always available since the number of buckets per pixel cannot be greater than k, the size of the hash table For a generated fragment, the shape ID is used to locate, with a linear search, the associated bucket (1) and extract its range [START_REF]Eurographics Symposium on Rendering -Experimental Ideas & Implementations[END_REF]. The current fragment is located within these range. To find its index, a depth comparison is used. Fragment's data, color and normal, are stored at this index (3). For neighbor query, the fragment (x, y, z) is used for a look-up in B(x + dx, y + dy) and Z(x + dx, y + dy).

incremental process is used to find a bucket, except that we check if the slot contains the right bucket. Organized in that way, B, in combination with Z, is used as an index lookup to perform efficient fragment query as shown in Figure 4.

Data storage

To store fragment's data, a second geometry pass is performed to fill the buffer D. The fragment's data are its color (c) and normal (n). Only fragments pre-selected during the first geometric pass must be stored. Thus, for each generated fragment, we perform a fragment query as shown in Figure 4. If the query is successful, data are stored to the given index, or rejected otherwise. Buffer D thus contains W × H × k pairs (c, n).

Implementation details

As presented in the previous section, we manipulate three 3D structures. Each of them will be able to contain W × H × k elements, where W and H correspond to the viewport size and k is the number of layers in the buffer.

During the first pass, the sorted-insertion by depth is done using atomic comparison. As atomic comparisons rely on a bitwise AND operator, depth must be used as Most Significant Bit, even if we store both depth and shape number. These two values are packed on 32 bits with 16 bits each, in the buffer Z. Using depth as MSB ensures that we keep the k-nearest fragments.

For the construction of the buckets, depth and shape number are swaped inside Z. This way, shape number is now used as MSB instead of depth and we can easily reorder our list. After reordering, buckets are constructed as described in the previous section. Using this structure requires 32 bits per bucket : 8 bits for the first fragment position, 8 bits for the last one and 16 bits for the shape number.

Finally, to insert the data in the buffer D, we rely on a binary search done directly inside the correct bucket. This bucket is located using the process in Figure 4 with a condition on shape number. Once we know the correct fragment's layer, we can store the fragment color after shading to use only 32 bits for color. The shading relies on the Bidirectionnal Reflectance Distribution Function (BRDF) proposed by Disney [BS12] for direct illumination and Image-Based Lighting to approximate indirect illumination. We store as well the fragment's normal, packed on 32 bits (15 bits for each x and y components and 2 for the sign of z). This results in using 64 bits per fragment for the data structure.

Feature-driven stylisation

In this section, we demonstrate how to use the BKB for an efficient feature-driven stylization. First, we present how to find fragment's neighborhood to compute curvature per fragment. Then, an application of stylization is shown adapted from existing shape depiction techniques. * 13] propose a method to find the neighbor of a fragment. A comparison is performed between a fragment and all the fragments contained in an adjacent pixel. The comparison function (ε) is based on normal and depth differences and the best neighbor candidate is the fragment which obtains the minimum value. We use a similar approach with the following criteria:

Feature extraction

Carnecky et al. [CFM

εn(i, j) = 1 -n i • n j εz(i, j) = |z i -z j | 2 ε(i, j) = εn(i, j) + εz(i, j)
To ensure that we correctly identify the neighborhood, each fragment labeled as the best candidate does a reverse check to determine if this candidate also views the current fragment as a neighbor. If one of the neighbors is not found, we use a special value indicating a missing neighbor. Figure 5 illustrates the interest of this reverse check. Comparing to the work of Carnecky et al. [CFM * 13], a crucial point in our method is the lower complexity of this algorithm. Actually, thanks to the BKB, fragments belonging to the same object are grouped together. Thus, as explained in the previous section, the number of pairs to test is restricted to a subset of candidates, i.e. the ones with the same shape number. We access to this subset using the same process as data storage, presented in Figure 4, except that the query is done on the neighboring pixel's lists (buckets and fragments).

The neighbor information is used to compute screen-space curvature per fragment, such as Vergne et al. [VPB * 09]. In this approach, principal curvatures, κ 1 and κ 2 , are extracted from surface normals. Instead of directly using of the mean curvature (κ 1 + κ 2 )/2, we prefer, for practical reason, the following mapping varying between -1 and 1:

κ = tanh(κ 1 + κ 2 )
For the rest of this paper, κ is considered as a curvature value. In our convention, -1 corresponds to a convetixity et 1 to a concavitiy. Contrary to Vergne et al. [VPB * 09], our method does not suffer from undefined behavior at the object's boundaries. By taking advantage of potential missing neighbors, we can handle those aera: the curvature of a fragment with a missing neighbor is set to κ = -1.

Stylization

For the stylization, our method modulates both shaded color and transparency. These modulations are controlled by transfer functions: an 1D texture for color and a 2D one for transparency.

The shaded color is modulated with a curvature-indexed transfer function, presented in Figure 6 we propose to modulate opacity with a transfer function indexed by absolute mean curvature and depth. A value of 1 corresponds to fully opaque and 0 to fully transparent. An example of transfer function can be seen in Figure 7. We choose to use the absolute mean curvature so that either convexities or concavities would have the same opacity. Indeed, based on our observations, both are crucial to fully grasp shape charactestics. Our system allows the use of any transfer function to modulate opacity. Figure 8 illustrates how different functions allow to emphazise specific parts of the scene. For instance, the function used in (a) and (c) results in nearly transparent fragments on flat surfaces close the viewport, while curved surface are nearly opaque. The opacity quickly grows as the distance to the camera increases. This behavior allows to place more focus on the front part of the scene. On the contrary, the function in (b) places less focus in front but depicts shape on a wider range of depth. Finally, (d) presents another behavior, where front parts are more opaque than far ones. Both use the function presented in Figure 6(a) to modulate shading so that convexities are enlightened while concavities are darkened.

Results and performances

To illustrate the benefits of extending shape depiction to transparent surfaces, Figure 9 with transparency, enhanced as well with our method. OIT allows only to guess the internal structure of the car while Mean Curvature Shading enhances surface features only on the outer parts. On the contrary, our method enhances both internal and external parts, with a focus on the front of the scene.

To measure the impact on performances of our method, we use a synthetic test configuration, in which rendering time is only dependent of the number of objects. The scene is composed of 32 stacked full-screen quads so that each of the 32 layers are always filled. Each quad is located at a random depth, and is associated with a random object. We then vary the maximum number of object, as described in Figure 10. We compare the result of our solution with the brute force approach. Note that the brute-force approach has a constant rendering time, as complexity is only dependent on the number of fragments. The latter roughly corresponds to the method of Carnecky et al. [CFM * 13]. Implementation and renderings were done using the Open Inventor R 9.6 SDK by FEI TM with an NVIDIA GTX 980. Results are shown on Figure 10 for 32 layers and a resolution of 512x512 (around 8,4 million of fragments). If our solution is obviously much slower that simple transparency (10 ms in this configuration), we can see that it improves performances for transparent shape depiction as soon as there is more than one bucket. With more than five buckets, we halve the time required, going down to a fifth of the performances of the brute force approach when there are many objects.

These observations fit rather well with theorical gain of Section 4, going from O(k 2 ) to O(nk), but only when n is at least 2. Indeed, when n ≤ 2, the cost of accessing all the fragments in a bucket becomes nearly equal to the cost of accessing the bucket, whereas for n > 2, the bucket cost becomes quickly insignificant with regards to the cost of accessing all fragments. However, this gain is valid only for neighbor query, as our structure still requires two sorting pass. In our algorithm, we use insertion sorting which has a complexity of O(k 2 ) in the worst case, O(k) when fragments are nearly-sorted, which is often the situation with our structure.

The last aspect concerns memory consumption. The proposed solution requires up to 128 bits per fragment (32 bits for depth and Figure 10: Cumulative rendering times of a fixed number of fragments (32) for different bucket configurations (colored quads) and 512 by 512 pixels. Full lines depict cumulative times of our method, the dot line represents the total time of the brute force approach. As a reference, a dual pass k-Buffer alone with constant transparency requires around 10 ms per frame.

bucket id, 64 bits for color and normal, and 32 for bucket informations) whereas a classic dual pass k-Buffer requires 64 bits per fragment. For a Full HD resolution, this over is around 32 MByte per layer versus 16 MByte for a simple k-Buffer.

Conclusions and limitations

We have presented a new approach to propose one solution to the bias in shape perception introduced by compositing multiple nonrelated fragments. We convey geometric properties of transparent surfaces with Mean Curvature Transparency based on the Bucketed k-Buffer. The BKB allows an easy access to the fragments of a shape. Moreover, we provide a customizable tool to enhance scene perception for transparent surfaces. Based on a user-defined transfer function, it is possible to choose which features should be submitted to Eurographics Symposium on Rendering -Experimental Ideas & Implementations (2016)

highlighted. The choice is highly dependent of which features need to be enhanced, as well as the complexity of the scene.

Our approach is dedicated to scenes with multiple objects: on a unique object, such an approach introduces too much complexity. Our structure performs better than existing ones, even if it implies a overhead in memory. However, our data structure may be optimized for more applications than neighbor query. In particular, we could use it with editing techniques like Fast-Edit, or extend the amount of optical phenomena that we currently simulate, like refraction, sub-surface scattering or translucency with a limited impact on performances.

On the other side, in a scene with many small objects, results can become less legible. The current stylization is based on the feedback provided by our team. A further user-study, with a representative panel, would be necessary to accurately validate the perceptual gain of our method. In future work, we also plan to investigate the impact of using line-based rendering ( [DFRS03, OBS04, JDA07, KST08], introduced in Section 2) to enhance third order features on transparent surfaces.

  Shape depiction techniques aim at conveying a message based on shape properties. Some previous methods propose to modify the reflection pattern according to constraints such as geometric features, material properties or stylisation. Gooch Shading [GGSC98] relies on a cool to warm color transition based on normal consideration giving information about surface orientation. With the same idea, Sloan et al. [SMGG01] propose to use precomputed Lit Sphere addressed with surface normals. Barla et al. [BTM06] extend this idea by adding depth information to mimic aerial perspective. Exaggerated Shading [RBD06] uses normals at multiple scales to enhance the relief of a shape. Using higher order consideration, Kindlmann et al. [KWTM03] modulate the shading between concavities and convexities based on object-space curvature. Vergne et al. [VPB * 09] propose to compute image-space curvature to shift the normals of a surface, thus exaggerating the reflection patterns. The screen-space approach has the additional benefit of exaggerating curvature magnitudes toward object contours, which increases surface legibility. Later, Vergne et al. [VPB * 10] propose a generalisation of this technique for any model of illumination and any type of material. Line-based rendering have also been used to convey geometric informations. Initiated by Saito and Takahashi [ST90], line-based techniques may be used to depict purely intrinsic geometric properties, for instance ridges and valleys (e.g. [OBS04]) or surface inflections (e.g. [KST08]). It can also enhance view-dependent features like occluding contours, Suggestive Contours [DFRS03] and later Suggestive Highlight [DR07], or Apparent Ridges [JDA07], an view-dependent proposition of [OBS04].

2 .Figure 2 :

 22 Figure 2: Pipeline used for feature-based stylization with Bucketed K-Buffer. The discretization pass fills the buffer Z, the bucket buffer B and the data buffer D with their respective information (see Section 4). The stylization extracts features and applies modulation on fragment's data (see Section 5). Finally, the blending pass consists in traditional OIT back-to-front composition.The first step consists in capturing the scene into the bucketed k-buffer (BKB). The BKB is a discretization of the scene such that fragments are organized pixelwise into three 3D data structures: Z, B and D. The roles of Z and D are nearly the same as data structures used in dual pass k-buffer[START_REF] Kubisch | Order independent transparency in opengl 4[END_REF] to capture, respectively, depth

Figure 3 :

 3 Figure3: Description of the discretization process. In the first pass (1), Z is filled with a depth-ordered fragment list, containing depth (z) and shape ID (b) for each fragment. In the second one (2) Z is reordered by shape Id. Then (3), B is filled with bucket's information : shape Id (b), first and last fragment index (i, j).

  .

Figure 4 :

 4 Figure4: Fragment query for data storage pass. For a generated fragment, the shape ID is used to locate, with a linear search, the associated bucket (1) and extract its range[START_REF]Eurographics Symposium on Rendering -Experimental Ideas & Implementations[END_REF]. The current fragment is located within these range. To find its index, a depth comparison is used. Fragment's data, color and normal, are stored at this index (3). For neighbor query, the fragment (x, y, z) is used for a look-up in B(x + dx, y + dy) and Z(x + dx, y + dy).

Figure 5 :

 5 Figure 5: Importance of doing a reverse check. (a) shows a simple check, which result in biased results at creases. (b) shows a reverse check, with creases identified as such.

  (a). This function is inspired by Mean Curvature Shading [KWTM03] and Radiance Scaling [VPB * 10], with bright convexities and dark concavities. It provided convincing results for our application, as we wanted to modulate a physically-based shading. However, for other applications, a different transfer function could be used in our pipeline for other stylisations as shown in Figure 6(b). The modulation of the transparency extends the work presented by Hummel et al. [HGH * 10] as we use curvature instead of normalvariation considerations. Based on the idea of X-Toon [BTM06],

Figure 6 :

 6 Figure 6: (a) The transfer function we use to modulate shading brighten convexities and darken concavities. (b) Another example : bright curved surfaces, dark flat ones.

Figure 7 :

 7 Figure 7: Example of different transfer functions. White corresponds to opaque (opacity of 1) and black to transparent (opacity of 0). Both follow the same idea: more transparency on flat and close surfaces, more opacity on curved surfaces as well as far ones.

Figure 8 :

 8 Figure 8: Examples of Mean Curvature Transparency. (a,b) The transfer functions used are respectively the one presented in Figure 7 (45 FPS versus 120 without modulation). (c,d) The transfer functions are shown in the top left corner (10 FPS, versus 20 without modulation). Times are obtained with a resolution of 1200 by 800 for the bike 1000 by 1000 pixels for the UNC Powerplant and 16 layers.

Figure 9 :

 9 Figure 9: (a) OIT using a k-Buffer with a constant opacity value of 0.6 (50 FPS). (b) Mean Curvature Shading only. (c) Mean curvature transparency (20 FPS). The transfer function used is on the top-right corner. Model courtesy of dessindus.blogspot.fr. Times are obtained with a resolution of 1000 by 1000 pixels and 16 layers.

Acknowledgements

We thank FEI TM for funding this work and providing a licence of Open Inventor R , and the members of the Open Inventor R team for their support and feedback, particularly Fedy Abi-Chahla for his useful advices on OIT. We are grateful to the GrabCAD community (Car prototype and Suzuki R1, credits go to their respective designer) and the AIM@SHAPE VISIONAIR repository for the Circular box model.