Jeanpierre Guédon

Nicolas Normand

THE MOJETTE TRANSFORM: the First Ten Years

In this paper the Mojette transforms class is described. After recalling the birth of the Mojette transform, the Dirac Mojette transform is recalled with its basic properties. Generalizations to spline transform and to nD Mojette transform are also recalled. Applications of the Mojette transform demonstrate the power of frame description instead of basis in order to match different goals ranging from image coding, watermarking, discrete tomography, transmission and distributed storage. Finally, new insights for the future trends of the Mojette transform are sketched.

Introduction

1.1 The Word "Mojette" "Mojette" is a well known word in Poitiers; in old french, it describes the class of white beans. Until recently, these white beans were the standard tool for a child to start computing additions and subtractions. They were also used for computing the number of victories for card games like "Aluette" still played with middle age spanish cards. Mojette is the name of our transform to remember first that when only adds are invoked, the computations can be easily made, second that by sharing the information pot each player will get a part of it.

The First Papers

The first paper on the Mojette transform has been published in 1995 after one year of hard work. At the moment, we were looking for something that we still did not find: a discrete mathematical tool that can split the Fourier plane into radial and angular sectors. The initial application was the psychovisual image encoding that mimics the human channelized vision. Hopefully, what we found was so much practical that many applications did appeared later on. The first talk on the Mojette transform (denoted MT in the following) was published in SPIE Visual Communications and Image Processing [START_REF] Guédon | Psychovisual image coding via an exact discrete radon transform[END_REF]. The fact that a novel transform was welcomed by the community was of prime importance to continue this kind of research. Even if all the theory was not yet there, the structure of the transform was already explained by its first inverse version using a recursive implementation. The first laboratory talk was given on June 22 d 1995.

Related Works

The first two precursors of the Mojette transform were Myron Katz and Gabor Herman. They were both looking at discrete geometric tools to properly inverse the Radon transform at the beginning of the X-ray scan. The concept of discrete angles that we used later on in the MT was taken from Myron Katz' work [START_REF] Katz | Questions of uniqueness and resolution in reconstruction from projections[END_REF]. Gabor Herman presented iterative algorithms to solve the equation f = Rp where f is the image, R is a discrete Radon projection operator and p is the set of projections obtained with a constant bin width [START_REF] Gabor | Image Reconstruction from Projections -Implementation and Applications[END_REF]. In this regard, mixing both works gives the Mojette transform.

The mood for a discrete Radon transform use only came back in the nineties; except for pioneer works as Attila Kuba who was starting to see the generality of the problem using discrete geometry in 1984 [START_REF] Kuba | The reconstruction of two-directionally connected binary patterns from their two orthogonal projections[END_REF]. Guédon [START_REF] Guédon | Bandlimited and haar filtered back-projection reconstuction[END_REF][START_REF] Guédon | Les problèmes d'échantillonnage dans la reconstruction d'images à partir de projections[END_REF] starts with a spline version of the Radon transform in order to fill the gap between discrete and continuous points of view. L. Dorst added a major point by linking the Radon (or slope) transform to mathematical morphology [START_REF] Dorst | Morphological signal processing and the slope transform[END_REF] as independently we did obtained our first reconstruction theorem by this link to the two pixel structuring element. Notice that a big step has been made in 1999 when Jean-Marc Chassery did organize a workshop on discrete tomography that was using discrete geometry for mixing the kind of solutions for a given polyominos problem [START_REF] Del Lungo | The number of convex polyominoes reconstructible from their orthogonal projections[END_REF] with the complexity of others very close problems [START_REF] Gardner | Geometric Tomography. Encyclopedia of Mathematics and its Applications[END_REF][START_REF]Discrete Tomography: Algorithms and Complexity, number 97042[END_REF]. However, it must be clear here that the Mojette transform does not exactly belongs to the "discrete tomography" community in the sense the word was defined by Larry Shepp in 1994, i.e. reconstructing an object from two to four projections.

A novel situation arises now, when many and very interesting papers are published to express different properties of transforms very close to the MT as [START_REF] Svalbe | Farey sequences and discrete radon transform projection angles[END_REF] and unread papers come back to the community as [START_REF] Matus | Image representation via a finite radon transform[END_REF].

Paper Organization

The paper is split into theoretical results and applications. The second section of the paper reviews the properties of the Mojette transform. The simple definition of the direct Mojette transform (that will be described by Dirac-Mojette afterwards) is first presented and its fast inversion (same order of complexity than the direct operator) follows. The strong relationships with mathematical morphology are then presented because they constitute both the core of the proofs of the reconstruction theorems we obtained and the geometric way to cope with the transform. The generalization to spline Mojette transform and other kernels are then presented to show the power of the relationships between continuous and discrete words. In other words, this link allows for translating various problems such as tomography into the discrete geometry world. Finally, extensions to higher dimensions demonstrate two notions : the order of complexity of the transform is still linear in the number of pixels, voxels, ixels (information elements) and linear into the number of hyperplanes. The third section present the variety of applications already found for the Mojette transform. The key point that is used here is the frame description notion (instead of the classical vectorial basis). Its first obvious success (because of its direct implementation) was to use the Mojette transform as a tool for multiple description which has evolved to a kind of new standard for communications [START_REF] Benoît Parrein | Multiple description coding using exact discrete radon transform[END_REF]. A step further in this area was to add the hierarchical (or multiresolution) data description with the concept of multi-layer buffering. Another completely different application was to master the noise properties into the Mojette domain in order to allows the Mojette transform for tomographic reconstruction. Even if this sounds obvious from the Radon inheritance, this works is the last one we were able to perform. Back to the image domain, the Mojette transform has been successfully implemented to solve for new techniques as watermarking as well as texture analysis for image and video as presented in this section. The final part of the paper gives some insights for the next ten years.

Basic Mojette Transform Properties

Direct and Inverse Dirac-Mojette Transform

The mojette transform is derived from the Radon transform [START_REF] Radon | Über die bestimmung von funktionen durch ihre integralwerte längs gewisser mannigfaltigkeiten[END_REF]. The 2-D transformed domain consists in projections (in the Radon sense) where each calculated element called a bin is the sum of ixel values. However, from an original block f (k, l) of information elements (denoted ixel for information element), the Mojette transform gives a linear set of projections proj θ (m) only for specific angles of the form θ = atan(q/p) where p ∈ Z, q ∈ Z + and are relatively prime (GCD(p, q) = 1). The Mojette transform set is defined by Mf(k, l) as a set of I projections M p,q f (k, l), where (k, l) belongs to the ixels block and ∆(m) = 1 if m = 0, 0 otherwise, as follows:

Mf = {M pi,qi f, i = 1, . . . , I} = {proj pi,qi , i = 1, . . . , I} (1)
M p,q f (k, l) = proj p,q (m) = proj θ (m) = f (k, l)∆(m + qk -pl), (2)
A bin value is then simply the summation of every ixel intersecting the line m = -qk + pl. As described in Fig. 1, the major difference with the Radon transform is that the bin spacing on a projection depends on the projection angle.

The well-posed nature of this linear transform and its explicit discrete nature were the reason of the new name not to confuse with the numerous regularized versions of the Radon transform (ill-posed inverse operator). A direct consequence of this sampling is that the number of bins of the projection indexed by (p, q) depends both of the projection angle and of the block dimensions. For instance, a P × Q rectangular block has its (p, q) projection composed of B bins with B = (Q -1)|p| + (P -1)|q| + 1. The order of complexity of the direct transform is obviously O(IN ) i.e. linear in the number of (p)ixels N and linear in the number of projections I.

The inverse transform can be implemented with the same order of complexity. This is an important characteristic of the Mojette transform which comes from the underlying mathematical morphology properties as described in the next paragraph. To get such a result, the reader has to pay attention to the projection of the "corners" of the shape: they exhibit some one to one ixel-bin relationships. Thus, these specific bins can be exactly backprojected in the image as soon as we keep for each bin the number of contributed ixels during the inverse transform. The tool for the proof of this algorithm (given in [START_REF] Normand | La transformée mojette : une représentation redondante pour l'image[END_REF]) is now presented.

Mathematical Morphology

Reconstruction conditions were demonstrated for any convex image [START_REF] Normand | La transformée mojette : une représentation redondante pour l'image[END_REF] by associating a projection angle (p,q) with a 2 pixels structuring element (2PSE) as a point couple {O, (p, q)} giving the same direction axis. The set of ixels that are projected in one bin m is the intersection of the image with the (p, q)line defined by: -qk + pl = m. If the image is convex, this set is connected relative to the (p, q)-neighborhood (all ixels are reachable from the others using (p, q) displacements). Thus, an opening (an erosion followed by a dilation) with the 2PSE {O, (p, q)} will empty the set if it is composed of exactly one ixel but will leave it unchanged in all other cases. The back projection in direction (p, q), by removing all ixels in a one to one correspondence with a bin, results in the opening of the image with the 2PSE:

I = I • {O, (p, q)} .
An ixel disappears from I only if it has no (p, q)-neighbor.

The final reconstruction result is not affected by the order bins are back projected. So the overall reconstruction can be viewed as an iterative series of openings with the 2PSEs:

I k+1 = I k • {O, (p 1 , q 1 )} • . . . • {O, (p I , q I )} .
The reconstruction stops if no ixels can be removed by any projection (I k+1 = I k or I k+1 = ∅). This occurs when every ixel in the image has at least one (p, q)-neighbor in every (p, q) direction. The smallest non-null image with that property is built by the series of dilations with the 2PSEs.

The reconstructibility criteria of Katz finds here another expression only using dilations. A convex image is not reconstructible if and only if the dilation result by 2PSE is not included in the image support (even if one pixel is concerned). This general result is expressed in the following theorem by the equivalence of two statements. This inversion criterion is an extension of Katz' criterion for rectangular images. In that particular case, R can be included in the image if its width (1 + i |p i |) and height (1 + i |q i |) are respectively smaller than the width P and height Q of the image. Thus, for a P × Q image, a complete reconstruction is achievable if [START_REF] Katz | Questions of uniqueness and resolution in reconstruction from projections[END_REF]:

P ≤ i |p i | or Q ≤ i |q i | (3)

Spline Mojette Transform

The previous Mojette transform is denoted Dirac-Mojette transform because the underlying 2D ixel representation is a Dirac field:

f (x, y) = k l f (k, l)δ(x -k)δ(y -l) . ( 4
)
However, very similar angles like (p, q) = (1, 1) and (p, q) = (30, 29) will do not have similar bin values: a direct inter-correlation computation between the two projections will be poor because the number of bins and the number of pixels contributing to each bin will be very different. However, Philippé and Guédon [START_REF] Philippé | Correlation properties of the mojette representation for non-exact image reconstruction[END_REF] have shown that the autocorrelation of the Mojette projection corresponds to the projection of the autocorrelation of the image. Summing ixel intensities weighted by the relative length of the ray inside each ixel, defines the spline 0 -Mojette transform. The corresponding continuous image is given by:

f (x, y) = k l f (k, l)β 0 (x -k)β 0 (y -l) , ( 5
)
where β 0 (x) = 1 when |x| < 1 2 , β 0 ( 1 2 ) = 1 2 and 0 elsewhere. It can be simply shown that the spline 0 -Mojette transform corresponds to the filtering of the Dirac-Mojette transform with a finite impulse response (FIR) filter. The FIR filter definition w p,q (m) corresponds to a trapezoidal shape, i.e. the projection of the shape of the ixel (a square in 2D) as described in [START_REF] Guédon | Spline mojette transform. application in tomography and communication[END_REF].

M 0 f = [M D f ] * w p,q , (6)
Where the kernel w can always be decomposed with small filters composed of unitary values. In other words only additions can be used for implementation of this convolution and the recursive filter that implements the inverse operator (going from spline 0 to Dirac Mojette data) will only use subtractions.

The generalization to higher spline order Mojette transform follows the same pattern. From a β n spline definition of the original function given as:

f (x, y) = k l f (k, l)β n (x -k)β n (y -l) , (7)
with β n (x) = β n-1 (x) * β 0 (x), the spline n -Mojette transform is defined as:

M n f (k, l) = [M n-1 f ] * w p,q (b) . ( 8
)
The main originating reason to use spline functions instead of Dirac fields is to complied with the generalization of the well-known Shannon-Wittaker sampling theorem obtained by Unser and AlDroubi [START_REF] Aldroubi | Families of wavelet transforms in connection with Shannon's sampling theory and the Gabor transform[END_REF][START_REF] Unser | Polynomial spline signal approximations: Filter design and asymptotic equivalence with shannon's sampling theorem[END_REF][START_REF] Unser | B-Spline signal processing: Part I -Theory[END_REF][START_REF] Unser | B-Spline signal processing: Part II -Efficient design and applications[END_REF]. This theorem allows to combine the projection decomposition with other image needs like wavelet or geometric (e.g. rotation) operators.

The main use of spline Mojette transform is to introduce a controlled redundancy inside each projection via the kernel w. Direct applications to error correcting codes have already been obtained [START_REF] Souhard | Simulation d'une chaîne de communication adaptée à la transmission d'images fixes sur canal réel[END_REF] in the presence of noise inside the projection data. This intra-redundancy can be managed at the same time that the inter projections redundancy to answer consistency of distributed information projections as shown in the applications.

nD Mojette Transform

Going from 2D results to 3D then nD is an indefectible temptation to generalize the obtained results. This was done by still considering the Mojette transform as giving (n -1)D hyperplanes from the nD discrete hypervolume. It can be easily shown that from N ixels distributed inside a convex set of dimension n and a choice of I hyperplanes, the direct and inverse Mojette transform still exhibit a complexity order of O(IN ). This is quite obvious for the direct transform. The way of having this very low cost algorithm for the inverse transform is to use the same mapping than previously used in 2D between ixels and bins. As a matter of fact, this mapping is dimension independent.

Another result that can be seen as independent of the dimension is the discrete central slice theorem (CST). The continuous version of the CST has been extensively used for Radon inversion as well as for image processing. It says that the 1D Fourier transform of a projection of angle θ equals the slice of the 2D Fourier transform of the image at angle θ + π/2. The discrete version of this theorem was established in 2D by Dudgeon and Mersereau using the Z transform [START_REF] Dudgeon | Multidimensional Digital Signal Processing[END_REF]. It has been checked with the 2D Mojette transform and extended to higher dimensions in a straightforward manner [START_REF] Verbert | An exact discrete backprojector operator[END_REF].

Mojette Transform and Redundancy

The Mojette transform matrix M is only filled by 1 and 0 values. Only additions and subtractions are required for M and M * . It is also the case for M -1 with the inverse Transform algorithm for exact values [START_REF] Normand | La transformée mojette : une représentation redondante pour l'image[END_REF]. The matrix M * M is Toeplitz-block-Toeplitz. Since only the additive group structure is needed for its definition, replacing the natural addition by modulus addition will not change any property of the transform. This can be done as soon as the initial values of the ixels are quantized onto an interval [0, 2b[ where b represents the number of bits of an ixel. In this case, a bin will belong to the same interval (by modulus addition), so do the reconstructed ixel. This represents an important matter for not having an overbinary representation in the Mojette domain as well as to use the same type of elements in both domains.

The computation of a redundancy factor as Red = #Bins/#Ixels -1, gives the true weight for managing a frame. However, this Red indicator does not explain at all the stability of the transform under inconsistency (e.g. noise) considerations. In such a case, the number of projections as well as the difference between the initial shape and the dilations of two ixel structuring elements (2ISE) will better explain the power of recovering inconsistency. For instance, for a 5 × 5 convex, the sets S 0 = {(5, 3)}, S 1 = {(3, 2), (2, 1)} and S 2 = {(1, 0), (-1, 1), (1, 1), (2, 1)} all fulfill the reconstruction conditions and give Red= 8/25. The corresponding reconstructed shapes are illustrated in Fig. 3.

The importance of splitting the initial information into (projection) pieces can be pictured as in the previous figure or computed as the conditioning factor of the matrix M * M . In this example, when the first degenerated information set gets a false bin value it will not been detected at the reconstruction step, the second set will detect the error but not correct it as the third set does.

The computation of a redundancy factor as Red= #Bins/ #Ixels -1, gives the true weight for managing a frame. However, this Red indicator does not explain at all the stability of the transform under inconsistency (e.g. noise) considerations. In such a case, the number of projections as well as the difference between the initial shape and the dilations of two ixel structuring elements (2ISE) will be better explained the power of recovering inconsistency. For instance, for a 5x5 convex, the sets S 0 ={(5,3)}, S 1 ={(3,2),(2,1)} and S 2 ={(1,0),(-1,1),(1,1),(2,1)} all fullfil the reconstruction conditions and give Red= 8/25. The corresponding reconstructed shapes are illustrated in figure 2.

Figure 2 : Redundancy of three sets of projections under a 5x5 f(k,l) lattice

The importance of splitting the initial information into (projection) pieces can be pictured as in the previous figure or computed as the conditionning factor of the matrix M * M. In this example, when the first degenerated information set gets a false bin value it will not been detected at the reconstruction step, the second set will detect the error but not correct it as the third set does.

MOJETTE TRANSFORM AND FAREY SERIES nn

Applications

The power of the Mojette transform lies in it abilities to manage inconsistency. The first way to have inconsistent data is to have not enough data whereas it exists some intrinsic correlation inside the data as shown in the last paragraph of section two. In this case, only some partial exact descriptions of the information are available: this corresponds to the multiple description problem. This information can also already be sorted in an hierarchical manner and treated with the same formilsm through the use of the Mojette transform. The second case is the shape description problem Fig. 3. Redundancy of three sets of projections under a 5 × 5 f (k, l) lattice

Applications

The power of the Mojette transform lies in it abilities to manage inconsistency. The first way to have inconsistent data is to have not enough data whereas it exists some intrinsic correlation inside the data as shown in the last paragraph of section two. In this case, only some partial exact descriptions of the information are available: this corresponds to the multiple description problem. This information can also already be sorted in an hierarchical manner and treated with the same formalism through the use of the Mojette transform. The second case is the shape description problem through some of its projections, where the Mojette transform becomes a geometric tool when allied with the mathematical morphology. The third case of inconsistent data arises when metadata is added inside the data as for the watermarking techniques. Then the two coexisting kinds of data must have different behaviors (one must be shown while the other should lie under the expressed data until a specific algorithm can reveal it. The fourth case is when enough projections are available but are corrupted by noise such that only an approximation can be computed as for tomography.

The Use of Redundancy: the Multiple Description Tool

The multiple description problem arises in the seventies and was studied by Ozarow and others as an information theory problem for telephony over two channels [START_REF] Ozarow | On a source coding problem with two channels and three receivers[END_REF]. It has been developed for a decade for Internet transmission and now for distributed storage. The paradigm is as follows: given an initial information set, computes I different descriptions such that for a lower number of descriptions there is still always an approximation that can be computed and used usually for real-time considerations.

Many kind of solutions have been employed ranging from Turbo-codes, Solomon codes, Tornado codes and of course Mojette projections. The major property of these codes is generally to be Minimum Description Separable (denoted MDS), that is without redundancy as for Solomon or Tornado codes or (1 + )-MDS as for the Mojette codes demonstrated by Parrein [START_REF] Benoît Parrein | Multiple description coding using exact discrete radon transform[END_REF]. The price of the redundancy can be easily paid to gain much more flexibility in the design of the code (which means to choose a pair of shape and set of projections) whereas Solomon codes must be only determined from Galois field GF (2 n ).

For the basic Mojette solution, a rectangular shape P × Q with (P Q) is employed in conjunction with a set of I projections:

S Q = {(p 1 , 1), (p 2 , 1), . . . , (p Q , 1)} .
With these particular directions (∀i, q i = 1) the image is reconstructible with as least Q projections among I according to Katz' criterion Eq. 3 (the height of the smallest non-invertible image equals the number of received projections plus one). The value P × Q gives the amount of data that can be transferred at a time (available bandwidth), I -Q is number of extra projections (descriptions) sent over the channel to protect against data losses. With this simple image support, protection against losses is equal for all data. More precise control of the redundancy can be achieved with a concatenation of image supports of different heights, allowing to fit the protection level to the importance of the data: the less important the data, the higher the image support [START_REF] Benoît Parrein | Multimedia forward error correcting codes for wireless lan[END_REF].

Finally, when the initial data exhibits some intra correlation, phantoms can be used to produce an approximated reconstruction from an iterated algorithm [START_REF] Philippé | Correlation properties of the mojette representation for non-exact image reconstruction[END_REF].

Distributed storage organization and management can follow the same path. Initial data are taken and dispatched after the Mojette operation giving A projections. When a small fixed number B(B A) of different projections issued from the A set are stored on each of the S servers then the Mojette organization is made and gives the following management properties:

if a server is attacked and destroyed by hackers, its content can be restored by any set of other servers as soon as it meets the Mojette reconstruction requirement, -if a server content is stolen, its content can not be useful by itself (several servers have to be attacked at the same time), -from the user side, the ask data will arrive from different sources and routes, which will decrease the network congestion and increase the trust in the data.

Discrete Tomography Tool

The continuous Radon transform described by :

proj(t, θ) = R(t, θ)[f (x, y)] = f (x, y)δ(t -xcosθ + ysinθ)dxdy , (9)
represents the continuous function f (x, y) by an infinite set of projections. The functional projection of f (x, y) onto a spline space {φ(xk), k ∈ Z} leads to an interpolation equation:

f φ (x, y) = k l f (k, l).φ(x -k).φ(y -l) . ( 10
)
When the discrete pixel grid is considered in the tomographic problem, the function f (x, y) in Eq. 9 is replaced by f φ (x, y) of Eq. 10. This leads (after inverting discrete and continuous sum signs) to a definition of the continuous projection from a discrete image and the spline interpolating function.

For tomographic purposes, φ is not taken as φ(x) = δ(x) but rather as a spline function to express the "quality of the physical length" that can be recorded at the detector.

This Mojette sampling has been used both for direct spline FBP algorithm and for iterative conjugate gradient algorithm in the presence of noise with good results. The only remaining problem consists in the rebinning of the initial data recorded from the tomographic device.

Image Watermarking

Watermarking represents the way of transferring invisible data inside data. It can be text in text as in the private correspondence of George Sand with Alfred de Musset or more generally marks informations inside images or videos in order to carry the rights of the author and all the business steps allowing to sell information bytes. The invisibility notion is taken in the human sense, thus the locations to put this information into the carrier information must be given by psychovisual studies. The employed tool that implements the mixture of both information must be flexible enough both to use the given locations and be robust to image attacks as low-pass filtering, small rotation, and image compression for instance.

The Mojette transform in this context has produced different kinds of tools. The first one consists to watermark the projections instead of the image. This technique has also been lately develop by Macq and others in the Radon domain. The clue is to obtain an inconsistent set of informations from the marked projections in order to resist to rotation attacks.

The second tool is related to Normand's inversion algorithm that use pixelbin correspondences at each step. By adding a controlled quantity of noise in the projections (in the proportion of the redundancy Red factor) the reconstruction becomes impossible because of the ill-posed properties of the Radon transform. However, if for each reconstructed pixel the choice can be made (at the coder step) to avoid corrupted bins then a unique reconstructing path can be determined and encoded as the clue to decode the image. In other words, a cryptographic scheme is obtained.

The first unsolved question lies in the possibility to decode the cryptotattooed image thanks to the linearity of the transform even if the modulus addition is employed. As a matter of fact, this problem has to be solve thanks to both the Red factor and the possible decodable errors inside Mojette data : discrete group theory allied with geometry seems the right way to cope with this problem.

Another interesting question arises from the psychovisual locations able to hide a message. A soon as these locations are computed (with a respective strength for each pixel), the best set of projections can only be calculated by testing all the resulting phantoms (very large set) that can match the greater number of positions to assess the robustness of the watermarking. Going from this NP complete problem to a polynomial solution seems feasible using discrete shape descriptions and superimpositions.

We are confident that such a tool will be usefully employed in many other applications because of its intrinsic redundancy properties. However, it should be stressed again that the strong link between projections and discrete directions thus shape is today the only way to get mathematical tools to solve for problems. As soon as this projection structure disapear and is replaced by a simple set of bins, the tools complexity become NP complete instead of linear. This has been the case both for watermarking and for obtaining compression scheme from the right set of bins when each of them reconstructs an entire line of ixels.

Conclusion

This paper has matched the timing for the ten years of a new transform coming from the word of discrete geometry. It allows for reviewing both what has been done in the theory and to demonstrate that discrete geometry can be an important tool for network optimization or distributed storage. These trends must and will be continued and amplified for the next ten years.

Fig. 1 .

 1 Fig. 1. Mojette Transform of a 3 × 3 f(k,l) lattice

Theorem 1 .

 1 Both propositions are equivalent f (k, l) defined on the convex G is reconstructible by {proj pi,qi , i = 1, . . . , I}; -R constructed by I dilations set {{O, (p i , q i )}, i = 1, . . . , I} is not included in G. Corollary 1. Both propositions are equivalent -G is reconstructible by {proj pi,qi , i = 1, . . . , I}; -the erosion of G by R is null.

Fig. 2 .

 2 Fig. 2. Reconstructibility of a convex image

Fig. 4 .

 4 Fig.[START_REF] Dudgeon | Multidimensional Digital Signal Processing[END_REF]. An image support of height Q = 2 with directions (0, 1), (-1, 1) and (1, 1). Any combination of two projections is sufficient to reconstruct the image

Acknowledgements

The authors would like to thank Jean-Marc Chassery for his indefectible help to promote the Mojette transform and Éric Andres and the DGCI board for their kind invitation to write this paper.


Nicolas Normand
email: nicolas.normand@univ-nantes.fr

Andrew Kingston
email: andrew.kingston@univ-nantes.fr

Pierre Évenou
email: pierre.evenou@univ-nantes.fr

A Geometry Driven Reconstruction Algorithm for the Mojette Transform

The Mojette transform is an entirely discrete form of the Radon transform developed in 1995. It is exactly invertible with both the forward and inverse transforms requiring only the addition operation. Over the last 10 years it has found many applications including image watermarking and encryption, tomographic reconstruction, robust data transmission and distributed data storage. This paper presents an elegant and efficient algorithm to directly apply the inverse Mojette transform. The method is derived from the inter-dependance of the "rational" projection vectors (pi, qi) which define the direction of projection over the parallel set of lines b = pilqik. Projection values are acquired by summing the value of image pixels, f (k, l), centered on these lines. The new inversion is up to 5 times faster than previously proposed methods and solves the redundancy issues of these methods.

Introduction

The Mojette transform is a form of Radon transform. It is an entirely discrete mapping which requires only the addition operation and is exactly invertible. It was first proposed by Guédon, Barba and Burger in 1995 [START_REF] Guédon | Psychovisual image coding via an exact discrete Radon transform[END_REF] in the context of psychovisual image coding. It has since been applied in many aspects of image processing such as image analysis [START_REF] Guédon | The Mojette transform: applications in image analysis and coding[END_REF], image watermarking [START_REF] Autrusseau | Image watermarking for copyright protection and data hiding via the Mojette transform[END_REF], image encrytion [START_REF] Autrusseau | Mojette cryptomarking scheme for medical images[END_REF], image compression [START_REF] Autrusseau | Lossless compression based on a discrete and exact Radon transform: A preliminary study[END_REF] and tomographic image reconstruction from projections [START_REF] Servières | A discrete tomography approach to PET reconstruction[END_REF][START_REF] Servières | The Mojette transform: Discrete angles for tomography[END_REF]. The unique properties of the transform have also made it a useful multiple description tool with applications in robust data transmission [START_REF] Parrein | Multimedia forward error correcting codes for wireless LAN[END_REF] and distributed data storage [START_REF] Guédon | Internet distributed image information system[END_REF]. A summary of the evolution and applications of the mojette transform entitled "The Mojette Transform: the Fist Ten Years" [START_REF] Guédon | The Mojette transform: the first ten years[END_REF] was presented at the last DGCI conference.

Since the Mojette transform is pre-dominantly used as a tool, (e.g., for image analysis, to apply a watermark, for channel coding), the transform and inversion procedure should be as efficient as possible especially for real-time applications. This paper presents a inversion algorithm which uses a geometrical approach to streamline the reconstruction process.

Section 2 recalls the definition and some important properties of the Mojette transform as well as the methods for exact inverse utilised to date. Section 3 outlines the proposed geometry driven inversion method initially for simple cases where q i or p i is constant for all I projections and then generalises the results to an arbitrary set of (p i , q i ). A comparison between this method and previously proposed reconstructions is presented in Section 4 followed by a conclusion in Section 5.

The Mojette Transform

The forward transform (Projection)

The linear integration of the discrete 2D function f (k, l) is obtained via the Mojette transform over a set of I pre-defined rational angles, θ i = tan -1 (q i /p i ). The pairs of integers defining the angles, (p i , q i ) must be relatively prime, i.e., gcd(p i , q i ) = 1, and since linear integration is directionally independant, q i is restricted to Z + (except for the case p i = 1, q i = 0) to ensure θ i ∈ [0, π[. Assuming a Dirac pixel model the linear integrations become sums over the pixels centred on the lines b = q i kp i l. The Mojette projection operator is defined as

M{f (k, l)} = Proj(p i , q i , b) = P -1 k=0 Q-1 l=0 f (k, l)δ(b + p i l -q i k), (1)
where δ(η) is the Kronecker function, i.e., δ(η) = 1 if η = 0, otherwise δ(η) = 0. An example of these projections is given in Fig. 1. The number of linesums called "bins" per projection, B, for a P × Q image is found as

B i (P, Q, p i , q i ) = (Q -1)|p| + (P -1)q + 1, (2)
with b ∈ [0, B i -1] for p i ≤ 0 and b ∈ [-(Q -1)p i , (P -1)q i ] otherwise. For a transform with I projections, unique inversion is possible provided the Katz criterion [START_REF] Katz | Questions of uniqueness and resolution in reconstruction from projections[END_REF] is satisfied, i.e.,

P ≤ I-1 i=0 |p i | or Q ≤ I-1 i=0 q i . ( 3
)
This criterion was generalised by Normand, Guédon, Phillipé and Barba [START_REF] Normand | Controlled redundancy for image coding and high-speed transmission[END_REF] for images of arbitrary shape. Their scheme generates the minimum sized ghost functions ,(i.e., functions that exist in the image but disappear in the projections, refer to [START_REF] Katz | Questions of uniqueness and resolution in reconstruction from projections[END_REF] for more detail) as a sequence of 2D convolutions with all two pixel structuring elements formed from the set of projection slopes q i /p i by 1 at (0, 0) and -1 at (p i , q i ). Any array which cannot contain the minimum ghost generated by the projection set therefore has an empty null-space and must have a unique inverse.

The inverse transform (Reconstruction)

When the Mojette transform was first proposed a recursive algebraic method was used for reconstruction. The following year a fast and more direct technique, requiring addition operations only was proposed by Normand, Guédon, Philippé and Barba [START_REF] Normand | Controlled redundancy for image coding and high-speed transmission[END_REF]. It solves for one pixel at a time and subtracts this value from the bins that include this pixel in each of the I projections [START_REF] Normand | Controlled redundancy for image coding and high-speed transmission[END_REF]. The reconstruction propagates from the image corners (where there is only one pixel value in the bins) to the centre. The first step of the inversion for f (k, l), as given in Fig. 1, is shown in Fig. 2a. For each of the P Q pixels there are O(I) operations, so the complexity of this technique is O(IP Q). If the number of projections, I, is chosen to be log(P Q), the Mojette transform has similar complexity to that of the fast Fourier transform [START_REF] Normand | Controlled redundancy for image coding and high-speed transmission[END_REF].

There are two minor problems with this method. First, locating the bins in a projection which can be back projected (i.e., those bins for which only one pixel value remains unknown its corresponding line of projection). Second, determining which one of the pixels, (k, l), in the line of projection, b = q i kp i l, is yet to be reconstructed.

A simple method is utilised to overcome these problems. Two "comptabilité" (or accounting) images are projected with the same projection sets and reconstructed simultaneously with the unknown image. The first of these is a unitary image, i.e., an image where f (k, l) = 1 for all pixels. The second is an index image which labels the pixels according to a raster scan, i.e., f (k, l) = k + lP . The projections of these images assist with the respective problems above. This inversion technique will be referred to as the Comptabilité Mojette Inversion (CMI) method.

In recent years, two back-projection reconstruction methods have been proposed in the context of applying the Mojette transform to reconstruct medical images from continuous projections. The first of these is an exact method which was discovered by Servières, Normand, Guédon and Bizais [START_REF] Servières | The Mojette transform: Discrete angles for tomography[END_REF]. Given all I possible projections in the P × Q array, back-projection (M * ) yields I -1 times the original pixel value plus the sum of the image, f sum (which can be found as b Proj(p i , q i , b) for any projection), i.e.,

M * {Proj(p i , q i , b)} = f (k , l ) = I-1 i=0 Proj(p i , q i , q i k -p i l) = I-1 i=0 f (k, l)δ(q i (k -k ) -p i (l -l )) = (I -1)f (k , l ) + f sum . (4)
The set of projections can be found from (p i , q i ) being all the points visible from the origin, i.e., Farey points, of the P × Q array and all symmetries, (-p i , q i ). Assuming a uniform density of Farey points in the plane, approximately I = 12P Q/π 2 projections are required.

The second back-projection technique uses the conjugate gradient method [START_REF] Servières | Conjugate gradient Mojette reconstruction[END_REF] to minimise ||M * b -M * M f || 2 where f is the reconstructed image. Both of these inversions are relatively stable in the presence of noise and therefore ideal in this context. The exact back-projection however requires a very large number of projections and the conjugate gradient method, while it does give the inverse, is unnecessary in the case of reconstruction from uncorrupted discrete projections.

Since the Mojette transform is often used as a CoDec in data transmission, the most efficient inversion possible is required for real time applications. The following section outlines a very efficient inversion method which is similar to the CMI method but determines the inter-dependance of projections using graph theory to remove the accounting problems. images of infinite size, "on the fly", only the image height Q must be finite, P is not restricted.

When reconstructing an image (according to the above criteria) using the CMI method, the reconstruction can be seen to originate in the image corners retaining a convex region of unknown pixel values and then propagates towards the right. Once the initial trivial section in the corners is completed, it can be noticed that the projections and image rows are linked in that the same projection is utilised to reconstruct the next pixel in the same row. The algorithm proposed here takes advantage of this. To describe how, it is preferable to begin with a simplified case where all projections have a common value for q i of 1.

The case where q

i = 1 for i ∈ Z Q
This case commonly arises when the Mojette transform is utilised for multiple description coding in packet data transmission [START_REF] Parrein | Multimedia forward error correcting codes for wireless LAN[END_REF]. If the set of projections are sorted by p i , (i.e., p 0 < p 1 < . . .), and reconstruction is performed from left to right, then row r of the image, f (k, r) for k ∈ Z P , is reconstructed by the Mojette projection, Proj(p r , 1, b).

Proof. Assume on the contrary that projection, Proj(p

Q-1 , 1, b), is used to re- construct the pixel value f (k, l) on a row other than Q -1, i.e., 0 ≤ l < Q -1.
This implies that the pixel value f (k+p Q-1 , l+1) has already been reconstructed by some projection other than Proj(p Intuitively this can be seen as ordering the projections by the slope (or angle) of their corresponding line of projection. This gives a convex hull to the reconstruction region that ensures the lines of back-projection can cut the hull such that only one pixel on the line lies within the hull; The condition necessary for reconstruction.

Q-1 , 1, b), say Proj(p r , 1, b). Thus the pixel value f (k + p Q-1 -p r ,
Since each projection corresponds to a one row of the image, a dependancy graph can be constructed to show the relationship between the projections in reconstruction. The dependancy graph for the example image of Fig. 1 is given in Fig. 3a. Here vertices correspond to pixels and directed edges represent the dependancies of each pixel on other pixels being reconstructed in the inversion process.

Two simple paths can be found to traverse the graph (as shown in Fig. 3bi and b-iv for the example). The reconstruction process involves beginning to the left of the image, so that only the rightmost of the vertices in the path intersect image pixels in column zero, and reconstructing pixel values according to the path. The distance the path is initially shifted is referred to as the offset. The path is then shifted right 1 pixel and the entire process repeated until the last pixel value in column P -1 of image is reconstructed. This process will be referred to as the Balayage (or sweeping) Mojette Inversion (BMI) method.

The path through the graph termed a reconstruction path does not necessarily have to traverse from one side to the other. Two seperate paths originating from opposite sides of the graph can terminate at a common vertex in the graph (some examples are shown in Fig. 3b-ii and b-iii). To optimise the reconstruction algorithm, it is desired to find the most compact path possible.

The offset for two paths in the graph terminating on row r due to all projection vectors with negative gradient is found as

Offset -(r) = r i=1 max(0, -p i ) + Q-2 i=r max(0, -p i ) = max(0, -p r ) + Q-2 i=1 max(0, -p i ) = max(0, -p r ) + S -, (5)
where S -= Q-2 i=1 max(0, -p i ). Similarly the offset due to all projection vectors with positive gradient is found as:

Offset + (r) = r i=1 max(0, p i ) + Q-2 i=r max(0, p i ) = max(0, p r ) + Q-2 i=1 max(0, p i ) = max(0, p r ) + S + , ( 6
)
where S + = Q-2 i=1 max(0, p i ). The width of the reconstructed path is determined by the maximum of these two offsets. The objective is therefore to find an r which minimises Offset total = max(Offset -(r), Offset + (r)). ( 7)

Let S = S --S + = Q-2
i=1 -p i . Note that if S -> S + then any p r ∈ [0, S] has no effect on Offset total . Similarly, if S -< S + then any p r ∈ [S, 0] has no effect on Offset total . Therefore the optimal p r lies in the range [min(0, S), max(0, S)]. If there is no p r within this range then that which is minimises the following should be selected; p r -min(0, S) + max(0, S)

2 2 = (p r -0.5S) 2 . ( 8
)
Balayage Inversion Algorithm(for

q i = 1 for i ∈ Z I ) Input: Set of
temp ← (p 0 -0.5S) 2 8 r ← 0 9 for i ← 1 to Q -1 do 10 if (p i -0.5S) 2 < temp then 11 temp ← (p i -0.5S) 2 12
r ← i ¯¯ Determine the initial image column offset for each projection 13 offset(r) ← max(max(0, -p r ) + S minus, max(0, p r ) + S plus)

14 for i ← r + 1 to Q -1 do 15 offset(i) ← offset(i -1) + p i-1 16 f or i ← r -1 downto 0 do 17 offset(i) ← offset(i + 1) + p i+1 ¯ Begin reconstructing image, f (k, l), at column k = -offset(r) 18 for k ← -offset(r) to P -1 do 19 for l ← 0 to r -1 do 20 f (k, l) ← Proj l (k -p k l) 21 for i ← 0 to Q -1 do 22 Proj i (k -p i l) ← Proj i (k -p i l) -f (k, l) 23 ¯for l ← Q -1 downto r do 24 k ← k + offset(r) -offset(l) 25 f (k , l) ← Proj l (k -p k l) 26 for i ← 0 to Q -1 do 27 Proj i (k -p i l) ← Proj i (k -p i l) -f (k , l) 28 ¯¯Ē nd
This reconstruction procedure can be trivially generalised to the case where all q i = m for m ∈ Z + . In this instance each projection Proj(p i , m, b) reconstructs m consecutive rows of the image. The reconstruction paths are simlar to that for the above case with q i = 1 but with m passes shifted down a row each time. An example directed graph and reconstruction path for m = 3 is presented in Fig. 4a. Another simple case for the reconstruction occurs when p i is constant for all projections as is discussed in the next section.

The case where

p i = m for i ∈ Z I
This case where m = 1 is the most common type of angle set used for transforming images with minimal redundancy as described in [START_REF] Normand | Controlled redundancy for image coding and high-speed transmission[END_REF]. Since p i is constant, for a proof similar to that given in section 3.1 to apply, the projections must be sorted in order of decreasing q i , (i.e., q 0 > q 1 > . . .), then the r th set of q i rows of the image, i.e., from row R + 1 up to row R + q r where R = r-1 i=0 q i , are reconstructed by projection Proj(m, q r , b).

Since all q i > 0 the set of reconstruction paths all have the same total offset of (I -1)m as shown for the example in Fig. 4b. The reconstruction is similar to that for constant q i , in that it requires multiple passes (q max in this case), however the number of vertices included in each subsequent pass decreases as shown for the example.

3.3

The case where p i ≥ p i-1 and q i ≤ q i-1

The constant q i and constant p i cases from Sections 3 and 3.2 can be amalgamated if, when ordered by slope, the projections have the property that p i is increasing and q i is decreasing, i.e., p i ≥ p i-1 and q i ≤ q i-1 for (0 < i < I). For this case the reconstruction path is straightforward, similar to that for constant p i path with q max passes. The paths through the region with p i ≥ 0 can be constructed independantly to those in the region with p i ≤ 0 similar to the case with constant q i . An example has been presented in Fig. 5.

The general case

The above reconstruction techniques can be generalised to any set of I projections such that I-1 i=0 q i = Q for q i ∈ Z + . As for the previous cases with varying q i , if the set of projections are sorted by slope p i /q i , (i.e., p 0 /q 0 < p 1 /q 0 < . . .), then the r th set of q i rows of the image are reconstructed by the Mojette projection, Proj(p r , q r , b). The proof is again similar to that given in section 3.1.

Since each projection corresponds to q i rows of the image, once again a dependancy graph can be constructed to show the relationship between the projections in reconstruction. However, each vertex of the graph is no longer assured of 2 originating and 2 terminating directed edges. There may only be a single terminating edge and there may be zero or many originating edges depending on the set of (p i , q i ) used to define the projections.

To ensure a reconstruction path in this instance, only the pixels located immediately inside the edge of the convex hull created by the lines of projection ordered by slope are considered. As for the constant q i case, the paths that terminate at the rows with minimum slope are used to generate the most compact convex hull. An example of the selection process is given in Fig. 6a with the dashed line giving the complex hull. The directed graph is then used to determine the reconstruction paths required for these selected vertices as shown in Fig 6b for the example.

Discussion

Although very similar in nature, the balayage inversion algorithm can be shown to be up to 5 times faster than the comptabilité inversion algorithm. Assuming the pixel count and pixel label projections have not been pre-computed (which is possible in the case of image transmission where the incoming array size is known), these must both be determined in 2×O(P QI) operations and during the inversion process, I projection value bins, pixel count bins and pixel label bins must be updated for each of the P × Q pixels (3 × O(P QI) operations), giving a total of approximately 5×O(P QI) operations. In contrast the balayage algorithm requires updating I projection bins for each image pixel once in a single pass across the image in O((P + Offset)QI) operations. This has been demonstrated for three types of angle sets in Table 1 comparing the computation times for the BMI, the CMI with pre-computed unitary and index image projections (CMI-pc) and the complete CMI. In implementation the actual gain in efficiency can be up to an order of magnitude since there is a periodic pattern to the BMI process that can be exploited while this is not the case for the CMI method where the position of the next pixel to be reconstructed is not predetermined at all. Since the proposed method is based on the inter-dependancy of each projection, the algorithm removes the need to search through the projections to find the next candidate bin that can be back-projected. As a pixel is reconstructed the predetermined dependancy graph indicates the pixel that can be reconstructed next and by what projection. This is highlighted by the relative performance of the CMI method for the General case of Table 1 where the projections have a large number of possible reconstruction bins to manage. The CMI algorithm is more robust however, it is more adaptable to any set of projections such as redundant sets where I-1 i=0 q i > Q and sets of partial projection data. The knowledge of which projections are used to reconstruct which rows of the image can be also used to ignore/discard projection bins that are not required for the inversion. This removes unwanted redundancy to optimise Mojette encoding as was investigated by Verbert, Ricordel and Guédon in [START_REF] Verbert | Analysis of the Mojette transform projections for an efficient coding[END_REF].

Conclusion and Future Work

A new inversion algorithm for the Mojette transform has been presented which takes advantage of the knowledge of the interdependancy of projections in reconstruction. The method is more direct and more efficient than previous methods however is less robust in terms of adaptability to any set of projections. This method of reconstruction also automatically enables optimal encoding by the Mojette transform by identifying which projection bins are required for inversion. Developing a BMI method that can be applied to reconstruct a redundant set of projections and can adapt to sets of partial projections as well as determining optimal coding incorporating redundancy are topics of future research.

Fig. 1 .

 1 Fig. 1. Four projections of a 4×4 image f (k, l), Proj(-1, 1, b), Proj(0, 1, b), Proj(1, 1, b) and Proj(2, 1, b).

Fig. 2 .

 2 Fig. 2. Reconstructon via the CMI method. (a) A candidate bin is selected in the projections of the unitary image. (b) The value in the corresponding projection bin of the index image gives the pixel to be reconstructed. (c) The value in the corresponding projection bin of the image gives the pixel value. The projections of all 3 images are then updated simultaneously.

  l) must have been reconstructed and since p Q-1 is the largest in the set of p i , then k + p Q-1p r > k and this pixel is further right than f (k, l). However, it can not be known if reconstructing from left to right; A contradiction.Therefore, only row Q -1 can be reconstructed by Proj(p Q-1 , 1, b). Since reconstruction requires that only one pixel remains unknown in the line of projection, Proj(p Q-1 , 1, b) can not be used to reconstruct any other row. Therefore, this proof can be repeated to show Proj(p Q-2 , 1, b) must reconstruct row Q -2 and so on, down to Proj(p 0 , 1, b).

Fig. 3 .

 3 Fig. 3. (a) The dependancy graph for the example image. (b) The 4 possible reconstruction paths

  Fig 4a Fig 4b 5. Fig 4 -DependancyGraph_p3_q3

Fig. 4 .

 4 Fig. 4. The dependancy graphs (grey) and simplest reconstruction paths (black) for the set of projections (a) Proj(Proj(-2, 3, b), Proj(-1, 3, b), Proj(1, 3, b) and Proj(2, 3, b) which requires 3 passes and (b) Proj(3, 1, b), Proj(3, 2, b), Proj(3, 4, b) and Proj(3, 5, b) which requires qmax = 5 passes.

Fig. 5 .

 5 Fig 4b

Fig. 6 .

 6 Fig 6aFig6b

  projections, Proj(p i , 1, b), ordered with increasing p i Output: Reconstructed image, f(k,l). Begin

	7

2

Compute S -, S + and S 2 S minus ← S plus ← 0

3 for i ← 1 to Q -2 do 4 S minus ← S minus + max(0, -p i ) 5 S plus ← S plus + max(0, p i ) 6 S ← S minus -S plus

Determine the rendezvous row r

Table 1 .

 1 Reconstruction times comparing the BMI with CMI-pc (pre-computed unitary and index image projections) and CMI methods. Times are given as a ratio with respect to the time to perform the forward Mojette transform.

	TEST Const. qi 4096 64 64 (-31,1), (-30,1), ... (31,1), (32,1) 1.08 P BMI Q I {(pi, qi) | i ∈ ZI }	CMI-pc CMI 4.6 6.6
	Const. pi 512 512 9 (1,52),(1,54),(1,55), ... (1,61) 1.05	5.5	7.5
	General	512 512 10 (±95,31), (±63,32), (±31,32),	1.42	8.6	10.6
		(±31,64), (±31,96)

A Geometry Driven ReconstructionFor this method of reconstruction, it is assumed that I-1 i=0 q i = Q. Any redundant projections are ignored. The reconstruction is performed from left to right, (reconstruction from right to left, top to bottom and bottom to top are symmetries of this method). These two properties imply this algorithm can reconstruct

Acknowledgements

AK holds a postdoctoral position at l'Université de Nantes supported by a grant from the Région Pays de la Loire, France.


A Olmos

I C Trelea

F Courtois
email: courtois@ensia.inra.fr

C Bonazzi

G Trystram

DYNAMIC OPTIMAL CONTROL OF BATCH RICE DRYING PROCESS

Keywords: grain breakage, dynamic optimization, optimal control, paddy rice drying, sequential quadratic programming

The drying of paddy rice may result in quality degradation, expressed as a head kernel yield, leading to significant commercial depreciation of the product. A mathematical model of the drying and of the quality degradation process was combined with a dynamic optimization algorithm to determine the drying conditions (air temperature and relative humidity as functions of time) that ensured the highest possible final product quality for a specified drying time and a specified final moisture content. The robustness of the optimal drying strategy with respect to the initial state of the product, to the model parameters and to the initialization of the optimization algorithm was verified. The compromise between the highest achievable final quality and the allowed total drying time was studied. The combination of simulation and optimization yielded a new insight in the rice drying process and in the quality preservation strategies.

INTRODUCTION

The technological processing of paddy rice includes its drying from an initial moisture content contained between 220 and 750 g water / kg d.b. (dry basis) after harvesting down to a final moisture content of 140 g/kg d.b. or less. For some operating conditions the drying process may lead to endosperm breakage due to mechanical stress [START_REF] Kunze | Grain fissuring potentials in harvesting and drying of rice[END_REF]. From an economic point of view, the quality of rice is expressed as a head kernel yield. The term "head rice" denotes milled rice comprised of kernels three-fourths or more the original kernel. The rice Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 [START_REF] Ban | Rice cracking in high rate drying[END_REF], 1319-1345.

price falls with each percentage point lost, so it is critical to maintain optimum conditions during drying, storage and milling operation to obtain a high head rice yield. [START_REF] Siebenmorgen | Role of moisture content in affecting head rice yield[END_REF] Different process strategies in order to improve head rice yield like grain tempering and process interruptions have been proposed. [START_REF] Chen | The study of interrupted drying technique for rough rice[END_REF][START_REF] Inprasit | Effect of drying air temperature and grain temperature of different types of dryer and operation on rice quality[END_REF]. This approach mainly keeps uniform (inner and outer) grain moisture, however a longer operation time is required.

Specific works determined that grain breakage is caused by an unequal shrinking of the endosperm resulting from uneven dehydration of the kernel [START_REF] Kunze | Moisture adsorption in cereal grain technology -a review with emphasis on rice[END_REF] due to high drying temperature [START_REF] Arora | Rice drying cracking versus thermal and mechanical properties[END_REF] and/or high drying rate [START_REF] Ban | Rice cracking in high rate drying[END_REF]. Abud et al. [START_REF] Abud | Processing quality of rough rice during drying -modelling of head rice yield versus moisture gradients and kernel temperature[END_REF] showed that the combination of operating conditions (drying temperature, drying time and air humidity) produces a moisture content gradient inside the grain that can either increase or decrease the mechanical stress.

Fissuring and thus breakage occur if the stress exceeds a certain threshold [START_REF] Kamst | Experimental study of the mechanical properties of rice grains in relation to drying stresses[END_REF]. Abud, et al [START_REF] Abud | Kinetics of mechanical degradation of paddy rice during drying/ Influence of process operations conditions and modelling for control design[END_REF] state an empirical equation for the quality degradation rate as a function of the operating conditions.

The optimal process control is concerned with the selection of the process operating conditions that maximize or minimize a function that relies technical and/or economic parameters called optimization criterion. The optimal operating conditions are usually computed off-line and then applied to the process without change [START_REF] Lee | Dynamic optimization of a continuous polymer reactor using a modified differential evolution algorithm[END_REF]. Thus, control could be successfully accomplished by a very simple controller like the PID [START_REF] Altinten | Tracking performance of control methods[END_REF]. In drying processes, quality improvement and operational cost reduction are two main objectives that, in some cases, are conflicting ones hence both optimization method and criterion have to be carefully chosen.

Boxtel, et al [START_REF] Boxtel | A preliminary study on strategies for optimal fluid-bed drying[END_REF] calculated the optimal control variables trajectories, temperature and flow rate, of a drying process were product quality depends in biological components that could be inactivated by the operation conditions. A deterministic dynamic optimization method is applied. The process model included quality degradation. The optimization criterion is expressed in benefits (quality depended) nets (operational costs are deducted).

Kiranoudis [START_REF] Kiranoudis | Pareto design of conveyor-belt dryers[END_REF] optimized a potato slice drying process with a similar objective).

Quality was defined as product color and the operating conditions were drying air temperature and relative humidity, both kept constant during the process. The mathematical model included the drying state equations, the color evolution kinetics and an estimation of the operation cost. The best operating conditions were determined based on a Pareto graph.

Instead of including quality reduction in the mathematical model, several authors, indirectly avoided quality degradation through the optimization criterion: reducing drying rate

Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 [START_REF] Ban | Rice cracking in high rate drying[END_REF], 1319-1345. [START_REF] Mao | Simulation and optimization of cross flow grain dryers[END_REF] or applying features (e.g. maximal grain temperature and maximal moisture distribution of the product) on the main target as penalty factors [START_REF] Vasconcelos | Development of a supervisory control strategy for the optimal operation of grain dryers[END_REF].

Trelea et al. [START_REF] Trelea | Optimal constrained non-linear control of batch processes: application to corn drying[END_REF] proposed an optimal strategy for corn drying. The optimization criterion was the drying cost taken as a weighted sum of the batch duration and the estimated energy consumption. The optimized operating conditions were the total drying time, the air temperature and relative humidity, subject to constraints on the final product moisture content and the final wet-milling quality. The mathematical model included drying and quality degradation kinetics. The current product state was estimated on-line based on available measurements and the operating conditions were adjusted in real time in order to compensate various disturbances such as higher or lower initial grain moisture content and temporary equipment failures.

In the present work, optimal drying conditions for paddy rice (air temperature and relative humidity trajectories) were determined. The optimization criterion was the final rice quality, expressed as a head kernel yield. The final product moisture content was constrained to a target value in a fixed operation time. The drying time -quality degradation compromise was studied by performing the optimization for several specified batch durations. The dynamic model used for optimization consisted of a compartment drying model for thin layers of rice together with an empirical quality degradation kinetic [START_REF] Abud | Kinetics of mechanical degradation of paddy rice during drying/ Influence of process operations conditions and modelling for control design[END_REF][START_REF] Abud | A compartmental model of thin-layer drying kinetics of rough rice[END_REF]. The dynamic optimization algorithm was based on a sequential quadratic programming technique.

MATERIALS AND METHODS

Product and dryer

The experimental data used in this work was obtained from Abud et al. [START_REF] Abud Archila | Modélisation simultanée des transferts et de l'évolution de la qualité technologique du riz paddy en vue d'optimiser les conditions de séchage[END_REF].

Homogeneous samples of paddy of the Ariete variety from Camargue, France were used. The product was harvested at 250 g/kg d.b. moisture content before vacuum packing at ambient temperature in the dark. The absence of mold was checked. A thin product layer (30 mm, i.e. 600 g of wet paddy) was dried in each experiment. determined using the NF ISO 6646 norm: samples of rough rice were hulled in a husker with two-rubber disk and milled on an abrasive con. Head rice yield is the weight percentage between milled "head rice" and total milled rice. The stated accuracy of the quality determination was 3%.

The automatically controlled grain dryer is schematically shown in Figure 1. It allows time-varying drying conditions ranged in the following bounds: air temperature between 15 and 160°C, air moisture content between the ambient value and 300 g/kg of dry air and air velocity between 0.25 and 2.7 m/s. Drying air is conditioned before entering to the dryer by sensors and PID controllers that assured specific values of air temperature, relative humidity and velocity. Product weight was measured at regular time intervals using electronic scales (air flow was deviated during weighting).

Mathematical model

A compartment model of thin-layer paddy drying kinetics proposed by Abud et al [START_REF] Abud | A compartmental model of thin-layer drying kinetics of rough rice[END_REF] is used in this work.

The process state vector (X) consists of the moisture content in the inner compartment (x 1 ), in the outer compartment (x 2 ), of the grain temperature (T g ) and of the grain quality (Q):

T g Q T x x X 2 1 [1]
The control variables vector (U) depend on control capabilities of the dryer in order to apply the optimal trajectories directly using the dryer controllers. Since air velocities, between 0.26 and 2.12 m/s, have no significant effect in thin layer drying, only air temperature (T a ) and relative air humidity (H R ) were chosen:

T R a H T U [2]
The two compartment model equations are: considered uniform and equal for all the compartments. The significance of the variables is given in Figure 2 and in the nomenclature section.

1 2 1 1 1 x x dt dx g [3]
The rice quality evolution, as described by Abud et al [START_REF] Abud | Kinetics of mechanical degradation of paddy rice during drying/ Influence of process operations conditions and modelling for control design[END_REF], follow a second order kinetic (10):

2 Q K dt dQ [6]
The quality degradation rate coefficient (K) is expressed as a function of the moisture content gradient between the two compartments and of the grain temperature (T g ) via an

Arrhenius-type law:

1 . 273 exp 5 1 2 0 g a T R E x x K K [7]
The average moisture content of the grain takes into account the volume fraction of each compartment ( and ):

2 2 1 1 m x x x [8]
The mass transfer coefficients between the two compartments ( 1 ) and between the outer compartment and the air ( 2 ) are: exp [START_REF] Abud | Kinetics of mechanical degradation of paddy rice during drying/ Influence of process operations conditions and modelling for control design[END_REF] The heat transfer coefficient between the grain surface and the drying air is computed using the empirical relationship proposed by Loncin and Merson (20):

g m T x B B
2 5 v L C [11]
The water activity at the grain surface is computed with a formula similar to the one used by Pfost (21):

4 3 2 2 1 exp exp C T C C x C A g w [12]
The partial vapor pressure at the grain surface is thus:

w sat g g A p p [13]
The numerical values of model parameters are listed in Table 1. For more details on the model derivation and on parameter identification the interested reader is referred to [START_REF] Abud Archila | Modélisation simultanée des transferts et de l'évolution de la qualité technologique du riz paddy en vue d'optimiser les conditions de séchage[END_REF].

Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 [START_REF] Ban | Rice cracking in high rate drying[END_REF], 1319-1345.

Model validation

The mathematical model was validated in the following range of operating conditions: air temperature between 40 and 80°C, relative humidity between 5 and 80% and air velocity between 0.26 and 2.12 m/s. The model accuracy was assessed using the mean absolute error between the simulated and the measured values of the grain moisture content and of the grain quality:

ts measuremen of No. meas m x x x E [14] ts measuremen of No. meas Q Q Q E [15]
A graphical comparison between the simulated and the measured grain moisture content is given in Figure 3a. The computed values are very close to the experimental ones.

The mean absolute error for the 17 drying experiments performed in various operating conditions is 15.1 g/kg of dry rice. In Figure 3b the model accuracy is reported for the various combinations of drying conditions used in the experimental design. The poorest accuracy is obtained for the lowest drying temperature and the highest relative humidity. Except this relatively small area in Figure 3b, the model accuracy is always better than 10 g/kg.

Figure 4a shows the plot of the simulated rice quality against the measured one for 15 drying experiments performed in various operating conditions. The mean modeling error is 3.4%. It can be seen in Figure 4b that the quality is predicted accurately for all combinations of drying conditions, including those with an extreme evaporating capacity: low temperaturehigh relative humidity and vice versa. It should be noted that even the largest modeling error (6.5%) is of the order of the measurement error allowed by the quality determination procedure (3%).

Evaporating capacity is the maximum water quantity than can be evaporated optimization criterion (J) was defined as the negative value of the final quality and minimization of this modified criterion was performed:

J f t t t U 0 ), ( min [16] ) ( f t Q J [17]
Since moisture content reduction and quality preservation are conflicting objectives, a target value for the final moisture content had to be specified as a constraint. Otherwise, too mild-drying conditions would have been selected:

target f m x t x ) ( [ 18
]
The control variables were constrained to lie in the model validity range: [START_REF] Loncin | Food Engineering : principles and selected applications[END_REF] The total drying time reduction is another main economic goal. The tradeoff between fast drying and quality preservation was dealt with by running several optimization calculations with various fixed drying times and constructing a Pareto graph of the highest achievable final quality (Q(t f )) versus the specified batch duration (t f ). The final economic decision was thus left to the user.

max min ) ( U t U U T R a H T U min min min [19] T R a H T U max max max

Initial control profile

The optimization algorithm needs a feasible control profile to start, i.e. temperature and relative humidity profiles that satisfy all the constraints, including the achievement of the final moisture content target in the specified drying time. The calculations were initialized with random control variable profiles so an auxiliary optimization criterion (J aux ) was used to obtain feasible (but generally sub-optimal with respect to the grain quality) controls:

aux t t t U J f 0 ), ( min [21] 2 ) ( target f m aux x t x J [22]
Control profiles obtained from the minimization of the auxiliary criterion (J aux ) were then used as initial profiles for the product quality maximization based on the main criterion (J).

Optimization algorithm

The considered optimal control problem is a dynamic one since the unknown control variables are functions of time. The optimization criterion is either a linear (J) or a quadratic (J aux ) function of the state variables. The constraints are linear functions of the state and control variables. The non-linearity comes from the dynamic model of the process. Taking into account these characteristics of the control problem, a sequential quadratic programming (SQP) algorithm was used for solving it ( 23) coupled with a collocation (time discretization) technique ( 24) and a safeguarded line search based on a "confidence region" method (25).

The main steps of the algorithm are:

Step Step 3: Solution of the local optimization sub-problem. Reliable numeric software exists for the solution of quadratic optimization problems with linear constraints. Commercial software based on an "active constraint set" method was used [START_REF] Press | Numerical recipes in C, Chapter 17: Two point boundary value problems[END_REF].

Step modification of one or more parameters of the reference batch as indicated in each specific case. All unspecified parameters, however, are those of the reference batch.

RESULTS AND DISCUSSION

Convergence of the optimization algorithm

The optimization algorithm is always initialized with random control trajectories as described previously. Robust convergence of the optimization procedure to a unique solution was verified. Figure 5a shows the average and the extreme control profiles obtained after six random initializations of the optimization algorithm for the reference batch. The corresponding state trajectories are given in Figure 5b.

It appears that there is no significant difference between the results of the six calculations. The average standard deviation along the moisture content trajectory equals 1.03 g/kg d.b., which is much less than the model accuracy (15 g/kg) and also less than the allowed experimental error for moisture content measurement (1.5 g/kg). The quality variations seem larger at the figure scale, with an average standard deviation of 0.19% along the trajectory. This variation is, however, insignificant compared to the model accuracy (3.4%) and to the allowed experimental error for quality measurement (3%). The quality uncertainty diminishes at the end of the batch: the final quality equals 76±0.17% for all six runs. Thus, the variability of the optimal control trajectories in Figure 5a does not have any significant effect on the state variables and the convergence of the optimization algorithm appears satisfactory.

Similar studies were performed to check if the convergence was maintained for a different total drying time and various ranges of allowed operating conditions. The observed variations were not significant in all cases. For example, in a 10 hours drying experiment standard deviations were found to be 0.27 g/kg for the moisture content and 0.007% for the product quality, both of them are lower than the reference batch ones.

Based on Figure 5 one can estimate the admissible tolerances for the practical implementation of an optimal drying experiment, without producing significant deviations of the process state. Robustness of the optimal control policy to changes in the initial product state Case 1. The robustness of the optimal control policy was studied by applying without change the optimal control profiles determined for the reference batch to batches with different values of one of the initial product states.

Firstly, the initial moisture content was varied between 230 and 310 g/kg d.b. When the initial moisture content is under the reference value. The drying rate is unnecessarily high since it was calculated for a target final moisture content of 270 g/kg, thus producing a significant quality loss. The opposite is true for the higher initial moisture content (Figure 6a):

where it is impossible to achieve the final moisture target and its final quality is only marginally improved though.

Secondly the initial grain temperature was varied between 20 and 60°C. In all cases the moisture content target is achieved. For 60°C the initial condensation phase (zone A) disappears (Figure 6c) and significant quality degradation occurs during the first 30 minutes (Figure 6d). At 40°C the drying and quality degradation kinetics are relatively similar to the reference conditions (20°C). Finally, the initial grain quality has been varied between 70 and 90%. This has not effect on the drying kinetic (Figure 6e). The quality evolution curves are almost parallel so the differences in the final product quality exclusively reflect the differences in the initial quality (Figure 6f). Case 2. The optimization algorithm was used to re-compute optimal control profiles for each of the initial product state values mentioned previously. For the initial grain temperature and initial product quality variations the results are not significantly different from those of the reference batch (not shown). However, for various initial moisture contents, specifically determined control profiles always ensure a final moisture content equal to the target value (Figure 7a). In order to reach the final moisture content in the specified time, the drying rate at the beginning of the batch was either strongly increased or decreased compared to the reference batch (Figure 7b) by changing the air evaporation capacity (Figure 7c). The quality degradation kinetics (Figure 7d) exhibit the expected behavior: drying conditions are less aggressive in the case of lower initial moisture content. The final product quality was improved compared to the results in Figure 6b (75 instead of 70%).

In summary, appropriate optimal control profiles are worth re-computing only for variations in the initial moisture content of the rice.

Robustness of the optimal control policy to changes in the model parameter values

The consequence of using optimal trajectories computed with a model with either underestimated or overestimated parameters was investigated. The moisture content profiles exhibited significant differences only for variations in the water activity. The grain quality, however, appeared much more sensitive to model parameters. An overestimation of the water activity makes the optimization algorithm to use higher evaporation capacities that cause significant quality degradation (Figure 8a). If the mass transfer coefficient between the inner and the outer grain compartment ( 1 ) is underestimated, the moisture content gradient is higher and the quality degradation is also increased (Figure 8b). A similar effect occurs if the mass transfer coefficient between the outer grain compartment and the air ( 2 ) is higher than expected (Figure 8c): the moisture content of the outer compartment decreases faster, increasing the moisture content gradient between the two compartments. The effects of the mass transfer coefficients on quality degradation are relatively limited, however. Five percent variation of the quality degradation rate coefficient (K 0 ) has a negligible effect (Figure 8d).

In contrast to the case of variations in the initial product state, re-computing optimal control profiles for the various values of the model parameters does not lead to significant performance improvement except for the water activity : the final quality decrease compared to the reference batch becomes 1% instead of 5% in Figure 8a.

Description of the optimal rice drying kinetic

The optimization results of the drying conditions for the reference batch (parameters in Table 2) are reported in Figure 9. The optimization algorithm computes the control variable profiles from which the state and the output variable trajectories are determined using the dynamic process model. Five zones can be distinguished in Figure 9: Zone A. Water condensation on the product surface occurs at the very beginning of the drying process (Figure 9d) because the grain temperature is below the dew point temperature of the air. The water activity of the outer compartment is close to one (Figure 9f) which allows the usage of relatively intense drying conditions without any quality degradation (Figure 9d). Although superficial water activity is equal to 1, the model simulation program calculates the water activity in the whole outer compartment that is slightly less than 1.

Both the air evaporating capacity (Figure 9c) and the drying rate (Figure 9e) are high at the beginning of this zone but they decrease gradually due to the increase of the relative air humidity.

Comment citer ce document : Olmos, A., Trelea Zone B begins when the grain average moisture content falls below its initial value.

The air evaporating capacity is decreased in order to limit internal moisture gradients and thus quality degradation. Two drying strategies are possible: working with conditions of low temperature-low air humidity or high temperature-high air humidity. Both assure weak evaporating capacities, but the drying rate is higher with the second choice.

Since operation duration is one of the constrains, the optimization algorithm uses preferably higher temperature conditions. Therefore, in this zone, the algorithm increases the relative humidity instead of decreasing the air temperature(Figure 9a). Zone C begins when the relative air humidity reaches its maximum allowed value.

Subsequent decrease of the evaporating capacity is only achieved through a temperature decrease (Figure 9a). The drying rate continues to decrease (Figure 9e). Zone D. The water activity at the grain surface reaches the air relative humidity (Figure 9f) and the moisture content in the outer compartment attains its equilibrium value.

The quality degradation rate begins to decrease and the drying rate reaches a minimum (Figure 9e). To continue the drying process, the evaporating capacity has to be increased.

However that increase is dosed in order to cause a succession of mild drying conditions with a quasi-equilibrium state between the superficial water activity and the relative humidity of the air (characteristic in Zone D). The optimization algorithm decreases the relative humidity as much as required in order to slightly decrease the water activity in the outer compartment.

Hence, the intercompartmental moisture content gradient, thus quality degradation rate, is reduced to a minimum. Zone E. Moisture content is far of the target moisture content so the drying rate has to be increased (Figure 9e). Actually, in this zone one third of the total extracted water is removed. This is achieved by a sharp decrease of the relative air humidity (Figure 9a) that causes an increased quality degradation rate, which attains its highest value at the end of the drying process (Figure 9e).

Dependence of the optimal drying conditions on the total drying time

The evolution of the optimal drying strategy when total drying time increases is illustrated in Figure 10. Long drying times allow the usage of mild drying conditions. Lower temperature and higher relative humidity at the beginning of the batch are favored.

Additionally, the control variable profiles become flat for long drying times, approaching constant drying conditions. The average air evaporating capacity decreases as duration increases but its "total" value, i.e. the area under each curve, increases indicating a lower efficiency since the evaporated amount of water is the same in all cases.

Higher drying time allows lower drying rate and leads to reduced quality degradation.

The optimization algorithm has the possibility to extend the drying phases with a low quality damage (zones A and D) and to reduce those with high quality degradation rates (zones B and E). It is well-known that the longer drying time, the better final product quality but also a higher drying cost is thus confirmed. This economic tradeoff is stated explicitly in the next section.

The final quality -total drying time tradeoff A given product, characterized by his initial state (moisture content, quality, grain temperature), dried at the optimal process conditions calculated by the optimization algorithm, yields a particular final quality value. The only way to modify this final value is to change the total drying time.

The highest achievable final product quality is plotted against the allowed drying time in Figure 11. Very fast drying necessarily leads to unacceptable grain breakage. Very slow drying, however, results in minor quality savings. The most convenient drying time has to be selected at the point where the increase of the product price due to higher quality compensates the increase of the process operation cost due to longer drying (principle of zero marginal cost or highest profit). This is essentially an economic decision, which is beyond the scope of this paper. total drying time. The optimal control strategy was very sensitive to the initial moisture content and water activity of the product. In practice, these parameters should be determined carefully in order to apply the most appropriate control profile. The drying time / final quality product tradeoff showed that in shorter operations quality degradation is important, but in too longer ones the quality improvement became negligible. The most appropriate drying conditions should be selected based on economic considerations including product price as a function of its quality and processing cost a function of the total drying time and the drying air parameters.

NOMENCLATURE Symbol Units Significance

A w none Water activity B 10 kg m -3 s -1 Mass transfer coefficient between the two grain compartments at 0°C B 11 kg -1 kg °C-1 Sensitivity coefficient of the mass transfer between the two grain compartments with respect to the mean moisture content and to the grain temperature B 20 kg m -2 Pa -1 s -1 Mass transfer coefficient between the outer grain compartment and the air at 0°C B 21 °C-1 Sensitivity coefficient of the mass transfer between the outer grain compartment and the air with respect to the air temperature

C 1 , C 2 kg kg -1
Sensitivity coefficients of the water activity with respect to the moisture content in the outer grain compartment C 3 , C 4 °C-1 , °C Sensitivity coefficients of the water activity with respect to the grain temperature C 5

Pa °C-1 Constant used in [START_REF] Lee | Dynamic optimization of a continuous polymer reactor using a modified differential evolution algorithm[END_REF] between heat and mass transfer coefficients (kg water) (kg dry matter) -1 Moisture content of the inner grain compartment

C pg J kg -1 °C-
T a H R (b) x m (kg/kg) x m Q (%)
Q (%) (f) Q(0) = 70 % Q(0) = 80 % Q(0) =

  Moisture content was determined by drying in an oven at 130 °C during 2 hours, whit pre-drying for high initial moisture content, at the beginning and at the end of the drying according to the NF V03-707 norm; during the process, moisture content was calculated by mass loss. The measurement accuracy stated by the norm is at least 1.5 g/kg. Rice quality was Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

  the moisture transfer for the inner compartment, Equation 4 for the outer compartment and Equation 5 represents a heat balance, the grain temperature is Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

  adiabatically for one kilogram of drying air. It is calculated as the difference in humidity rations (expressed in g of water per kilogram of dry air) between drying air and the adiabatically saturated one. (22) Constrained optimization problem The selected optimization goal was the maximization of the final rice quality. In order to fit the standard problem formulation required by the optimization algorithm, the Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

  Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

1 :

 1 Discretization of the dynamic problem. The state and control variables are sampled at a finite number of time points. This reduces an infinite-dimensional problem to a finite-dimensional one, which can be handled by a numerical computer. The optimization criterion, the constraints and the dynamic model equations are discretized on the same time grid. Step 2: Construction of a local optimization sub-problem. The non-linear optimization problem is approximated locally (using limited Taylor series expansion) by a quadratic one with linear constraints. The discretized dynamic model equations are treated as additional constraints. The limited Taylor series expansion being valid only locally, a limited search region around the current solution is established using a set of "box" (min-max) constraints.

  Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

  Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

  , I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

  Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

CONCLUSIONA

  sequential quadratic programming technique was used to determine the drying air temperature and relative humidity profiles that maximize the final head kernel yield of paddy rice. The optimization algorithm used a dynamic compartment model of the drying process coupled with a quality degradation kinetic. Constraints on total drying time, final moisture content and allowed range of operating conditions were imposed. The control profiles issued by the optimization method showed robust convergence. Both simulation and optimization yielded a new insight in rice drying process. Quality optimization strategy favored mild drying conditions. Low quality degradation phases (e.g. condensation and quasi equilibrium drying zones)were extended as much as allowed by the Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

x 2 ( 1 ( 1 Mass transfer coefficient between the two grain compartments 2 (

 2112 kg water) (kg dry matter)-1 Moisture content of the outer grain compartment x m (kg water) (kg dry matter)-1 Average grain moisture content x meas (kg water) (kg dry matter)-1 Measured value of the average grain moisture content x target (kg water) (kg dry matter)-1 Target value for the final average grain moisture content W m -2 °C-1Heat transfer coefficient between the grain and the air kg dry matter) m -3 s -kg water) m -2 Pa -1 s -1 Mass transfer coefficient between the outer grain compartment and the drying air g the outer grain compartment Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345. Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

2 FIGURE 2 .

 22 FIGURE 2. Schematic representation of the grain compartments.

FIGURE 3 .FIGURE 4 .

 34 FIGURE 3. Validation of the dynamic drying model. (a) Simulated versus measured moisture content. (b) Model accuracy for various combinations of operating conditions.

FIGURE 5 .

 5 FIGURE 5. Convergence of the optimization algorithm from six random initial control profiles. (a) Control variable profiles: air temperature and relative humidity. (b) Process output trajectories: product moisture content and quality. Solid lines are mean values and dotted lines are extreme values of the six calculations.

90 %

 90 Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

FIGURE 6 .

 6 FIGURE 6. Robustness of the optimal control strategy with respect to changes in the initial product state. The optimal control profiles computed for the reference batch are applied without change. (a, b) Variations in the initial grain moisture content. (c, d) Variations in the initial grain temperature. (e, f) Variations in the initial grain quality. Symbols: ( ) reference batch, ( ) low value of the initial state, ( ) high value of the initial state.

FIGURE 7 .

 7 FIGURE 7. Robustness of the optimal control strategy with respect to changes in the initial moisture content of the grain. The optimal control profiles are re-computed for each value of the initial moisture content. (a) Average product moisture content. (b) Drying rate. (c) Evaporation capacity of the drying air. (d) Grain quality.

FIGURE 8 .

 8 FIGURE 8. Sensitivity of the quality degradation kinetic with respect to changes in the model parameters. The optimal control profiles computed for the reference batch are applied without change. (a) Variations in the water activity at the surface of the grain. (b) Variations in the mass transfer coefficient between the inner and the outer grain compartment. (c) Variations in the mass transfer coefficient between the outer grain compartment and the

  Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

3 FIGURE 9 .

 39 FIGURE 9. Optimal control strategy for the reference batch. (a) Optimal profiles of the control variables: air temperature and relative humidity. (b) Rates of change of the control variables. (c) Air evaporation capacity. (d) Optimal trajectories of the output process variables: grain moisture content and quality. (e) Rates of change of the output variables. (f) Comparison between the water activity at the product surface and the relative humidity of the drying air.

  document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

FIGURE 10 .

 10 FIGURE 10. Optimal control policies for various total drying times. Symbols: ( ) 1 hour, ( ) 2.67 hours, ( ) 4.5 hours, ( ) 7.5 hours. (a) Air temperature. (b) Relative humidity. (c) Evaporation capacity of the drying air. (d) Drying rate. (e) Quality degradation rate. Insert: Quality degradation rate for short drying times. (f) Quality degradation.

TABLE 1 Dynamic

 1 Model Constants and Parameters by Abud 2000 (19) Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345. Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

		Constants
	Symb	Units	Value
	ol
	C 5	Pa °C-1	65
	C pg	J (kg dry matter) -	1300
		1 °C-1
	C pw	J (kg water) -1 °C-	4210
		1
	g	(kg dry matter)	1500
		m -3
	S sg	m 2 m -3	2000
		m 3 m -3	0.6
		m 3 m -3	0.4
	R	J mol -1 K -1	8.32
		Parameters determined from experimental data
	B 10	(kg dry matter)	0.01316
		m -3 s -1
	B 11	(kg water) -1	0.3083
		(kg dry matter) °C-1
	B 20	(kg water) m -2	2.304 10 -9
		Pa -1 s -1
	B 21	°C-1	0.0442
	C 1	(kg water)	0.319
		(kg dry matter) -1
	C 2	(kg water)	0.0493

TABLE 2

 2 Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345. Comment citer ce document : Olmos, A., Trelea, I.-C., Courtois, F., Bonazzi, C., Trystram, G. (2002). Dynamic optimal control of batch rice drying process. Drying Technology, 20 (7), 1319-1345.

	Reference Batch Conditions
	Category	Symbol	Value
	Initial product state	x m (0)	270 g/kg d.b.
		T g (0)	20 °C
		Q(0)	80 %
	Admissible range of the control	T a min	40 °C
	variables	T a max	80 °C
		H R min	5 %
		H R max	80 %
	Imposed conditions	x target	130 g/kg d.b.
		t f	2 hours
	The limit values of control variables and the imposed conditions are
	considered in the optimization strategy as constrains

  drying air. (d) Variations in the quality degradation rate. Symbols: ( ) reference batch, ( ) low value of the model parameter, ( ) high value of the model parameter.

	T a (°C) H R (%)				(a)						x m (kg/kg)	(d)
	80
	60 70										T a						0.28					x m (0)
	40 50										H R						0.23					x m
	30																							Q
	20																	0.18
	0 10	0	A	0.2	0.4 B	0.6	0.8	1	C	1.2	1.4 D	1.6	1 E	0 0.13	A		0.2 0.4 0.6 0.8 B	1 C	1.2 1.4 1.6 1.8 D E
										Time (hours)											Time (hours)
	dT a /dt (°C/h) dH R /dt (%/h)		(b)						dx m /dt (kg/kg/h) (e)	dQ/dt
		50														0.21
		0															0.16
	-100 -50								dH R /dt dT a /dt					0.11				dx m /dt	dQ/
																		0.06
	-150
	-200	0	A	0.2	0.4 B	0.6	0.8		1 C	1.2	1.4 D	1.6	E	0 0.01	A	0.2 0.4 B	0.6 0.8 C	1	1.2 1.4	1.6 1 E
										Time (hours)											Time (hours)
	E v (g/kg)				(c)						H R /100 ( ) A w ( )	(f)
	12 A			B			C		D			E	1
	10																0.8
	8																0.6
	6																0.4
	4																0.2
	2	0			0.2	0.4	0.6	0.8	1	1.2	1.4		1.6	1	0	0	A	0.2	0.4 B	0.6	0.8	1 C	1.2	1.4 D	1.6	E
									Time (hours)											Time (hours)

ACKNOWLEDGEMENTS

The authors gratefully acknowledge to Consejo Nacional de Ciencia y Tecnología (CONACYT -Mexico) for the financial support of this work.