Method and device for measuring the colour of an object
09671329 ยท 2017-06-06
Assignee
Inventors
Cpc classification
G01N21/4738
PHYSICS
G01N21/255
PHYSICS
International classification
G01J3/46
PHYSICS
Abstract
The present invention relates to a method for measuring the uniform diffuse reflectance R.sup.OBJ() at least at one point on an object (30) using a device (10) comprising a means (11) capable of emitting color illuminants expressed in the form of luminous flux and an electronic color image sensor (12). The present invention also relates to a device (10) comprising a means (11) for emitting color illuminants expressed as luminous flux of colors and an electronic color image sensor (12), for measuring the uniform diffuse reflectance R.sup.OBJ() at least at one point on an object (30) placed in a zone located opposite and substantially perpendicular to the said means (11) capable of emitting colors and located in the field of vision of the said electronic color image sensor (12) and being subjected to an external illuminant expressed as a constant and unknown external environmental luminous flux (40) denoted I.sup.ext().
Claims
1. A method for measuring the uniform diffuse reflectance R.sup.OBJ() at least at one point of an object (30) by using a device (10) comprising a means (11) capable of emitting coloured illuminants expressed in the form of luminous flux and an electronic colour image sensor (12), characterized in that it comprises the following steps: placing of the said object (30) in a zone located opposite and substantially perpendicular to the said means (11) capable of emitting coloured illuminants in the form of luminous fluxes of colours and located in the field of vision of the said electronic colour image sensor (12), the said object (30) also being subjected to an external illuminant in the form of a constant and unknown surrounding external luminous flux (40) I.sup.ext(), where denotes the wavelength; emission by the said means (11) of a series of N illuminants S.sup.source().sub.i (with N being a natural number greater than one, i varying from 1 to N and being the wavelength), S.sup.source().sub.i being known as a function of the input parameters of the said means (11) capable of emitting luminous fluxes of colours, capture by the said electronic colour image sensor (12) of the luminous flux reflected at least at one point of the said object (30) and entering in the sensor, the said luminous flux being denoted as E.sup.capteur().sub.i, with N being a natural number strictly greater than two, i varying from 1 to N and being the wavelength; and obtaining of N equations E.sub.i: E.sup.capteur().sub.i=R.sup.OBJ()*(I.sup.ext()+S.sup.source().sub.i) due to the additive nature of the wave light and by definition of the uniform diffuse reflectance R.sup.OBJ() at least at one point of the object (30); and determination by the said device (10) of the two unknown continuous functions R.sup.OBJ() and I.sup.ext() by solving the system of N equations E.sub.i: by integrating each equation E.sub.i on the intersection of the source and sensor spectra, by denoting x, y and z the sensitivities in the colorimetric base selected, each equation E.sub.i then generating three E.sub.i integrated equations:
E.sup.capteur().sub.i*x()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.SOURCE().sub.i)*x()*d
E.sup.capteur().sub.i*y()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.SOURCE().sub.i)*x()*d
E.sup.capteur().sub.i*z()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.SOURCE().sub.i)*x()*d by calculating the numerical value corresponding to the left-hand side of the Ei integrated equations with the use of the output parameters of the digital image sensor; and by expressing the two unknown continuous functions R.sup.OBJ() and I.sup.ext() with the use of a finite number of interpolation points (.sub.j, y.sub.j) connected by at least one interpolation function s() for maintaining the continuous nature of the said unknown continuous functions R.sup.OBJ() and I.sup.ext(), the .sub.j being wavelengths selected in the intersection of the source and sensor spectra and being input parameters of the method, chosen to minimize the number of interpolation points for a given precision; and by finding the parameters y.sub.j of the functions R.sup.OBJ() and I.sup.ext() that minimize the least squares system A*XB.sub.2 resulting from the E.sub.i integrated equations.
2. A method according to claim 1, characterized in that it further includes a step of determining the value of the external illuminant I.sup.ext().
3. A method according to claim 1, characterized in that it includes in addition, a step of transcription of the function R.sup.OBJ() of uniform diffuse reflectance at least at one point of the object (30) into the CIE XYZ coordinates for a given illuminant.
4. A method according to claim 1, characterized in that the series of N illuminants is of the same order of magnitude as the number of interpolation points for determining the values of the uniform diffuse reflectance R.sup.OBJ() at least at one point of the object (30) and of the external illuminant I.sup.ext().
5. A method according to claim 1, characterized in that it includes a step of determining the values of the uniform diffuse reflectance R.sup.OBJ() at least at one point of the object (30) and the external illuminant I.sup.ext() in several spectral bands.
6. A method according to claim 1, characterised in that the said device (10) makes use of a screen for emitting flashes of colour and an electronic image sensor for sensing and capturing the light reflected by the target object.
7. A method according to claim 1, characterised in that the said device (10) is a camera unit or a camera with in built or removable flash.
8. A method according to claim 1, characterized in that the said device (10) implements waveguides for ensuring effective transiting of the emission and reception of flashes of colours.
9. A method according to claim 1, characterized in that it is implemented in order to take spectrometric photographs of objects and to make chromatic adjustments (balancing of whites) at will.
10. A method according to claim 1, characterized in that it is implemented in order to measure the colour of an element included in the following group: materials, solids, liquids, gases, paintings, tapestries, graphics, textiles, plastics, woods, metals, soils, minerals, plants and foods.
11. A method according to claim 1, characterized in that it is implemented for the measurement of colours for medical or cosmetic purposes on human beings and living organisms of at least one element included in the following group: skin, pimples, moles, hair, fur, makeup, and teeth.
12. A method according to claim 1, characterized in that it is implemented for the use of colour barcodes, of one or more dimensions.
13. A method according to claim 1, characterized in that it is implemented with a view to assisting people having colour blindness and/or who are blind.
14. A device (10) comprising the means (11) capable of emitting colour illuminants in the form of luminous flux of colours and an electronic colour image sensor (12), for measuring the uniform diffuse reflectance R.sup.OBJ() at least at one point of an object (30) placed in a zone located opposite and substantially perpendicular to the said means (11) capable of emitting colours and located in the field of vision of the said electronic colour image sensor (12), and also being subjected to an external illuminant in the form of a constant and unknown surrounding external luminous flux (40) denoted as I.sup.ext(), characterized in that it comprises the means for: emitting a series of N illuminants S.sup.source().sub.i (with N being a natural number strictly greater than two, i varying from 1 to N and being the wavelength), S.sup.source().sub.i being known as a function of the input parameters of the said means (11) capable of emitting luminous fluxes of colours, capture by the said electronic colour image sensor (12) of the luminous flux reflected at least at one point of the said object (30) and entering in the sensor, the said luminous flux being denoted as E.sup.capteur().sub.i, with N being a natural number greater than one, i varying from 1 to N and being the wavelength; and obtaining of N equations (E.sub.i): E.sup.capteur().sub.i=R.sup.OBJ()*(I.sup.ext()+S.sup.source().sub.i) due to the additive nature of the wave light and by definition of the uniform diffuse reflectance R.sup.OBJ() at least at one point of the object (30); and determining the two unknown continuous functions R.sup.OBJ() and I.sup.ext() by solving the system of N equations E.sub.1: by integrating each equation E.sub.i on the intersection of the source and sensor spectra, by denoting x, y and z the sensitivities in the colorimetric base selected, each equation E.sub.i then generating three E.sub.i integrated equations:
E.sup.capteur().sub.i*x()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.SOURCE().sub.i)*x()*d
E.sup.capteur().sub.i*y()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.SOURCE().sub.i)*x()*d
E.sup.capteur().sub.i*z()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.SOURCE().sub.i)*x()*d by calculating the numerical value corresponding to the left-hand side of the Ei integrated equations with the use of the output parameters of the digital image sensor; and by expressing the two unknown continuous functions R.sup.OBJ() and I.sup.ext() with the use of a finite number of interpolation points (.sub.j, y.sub.j) connected by at least one interpolation function s() for maintaining the continuous nature of the said unknown continuous functions R.sup.OBJ() and I.sup.ext(), the .sub.j being selected wavelengths in the intersection of the source and sensor spectra and being input parameters of the method, chosen for minimizing the number of interpolation points for a given precision; and by finding the parameters y.sub.j of the curves R.sup.OBJ() and I.sup.ext() that minimize the least squares system A*XB.sub.2 resulting from the E.sub.i integrated equations.
15. A method according to claim 2, characterized in that the series of N illuminants is of the same order of magnitude as the number of interpolation points for determining the values of the uniform diffuse reflectance R.sup.OBJ() at least at one point of the object (30) and of the external illuminant I.sup.ext().
16. A method according to claim 2, characterized in that it includes a step of determining the values of the uniform diffuse reflectance R.sup.OBJ() at least at one point of the object (30) and the external illuminant I.sup.ext() in several spectral bands.
Description
BRIEF DESCRIPTION OF THE DRAWINGS
(1) The invention will be better understood with the aid of the description, provided here below for purely explanatory purposes, of an embodiment of the invention, with reference to the Figures in which:
(2)
(3)
(4)
DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION
(5) By way of a preliminary remark, it should be noted that the notations in this section are as follows: sensor=electronic colour image sensor (the video function is not used in the context of calibration and benchmarking, only the taking of static images is used), source=source of illuminants/flashes of colours (screen, diode(s), laser(s), etc). (R, G, B).sup.sourcei=non-linear input chrominance parameters of the source of colour illuminants (RGBW.sup.source i or RGBY.sup.source i for the devices with 4 and more primaries), BL.sup.source i=non-linear input luminance parameter of the source of colour illuminants (for example the back lighting of a LCD screen); (R, G, B).sup.capteur i=non linear output chrominance parameters of the luminous flux captured by the colour image sensor, BV.sup.capteur i=non linear output luminance parameter of the luminous flux captured by the colour image sensor (for example Brightness Value of Japanese standardEXIFExchangeable Image File Format); OBJ=colour object to be measured; EXT=external environmental illuminant; R/G/B=equation valid for each primaryred (red) or green (green) or blue (blue); OBJ/EXT=equation valid for the coloured objector foroutside external illuminant; x/y/z()=equation valid for each spectralsensitivity X(), y() and z(); spectralsensitivities of the CIE 1931 standard observer=x.sup.EC.sup._.sup.CIE.sup._.sup.1931 (), y.sup.EC.sup._.sup.CIE.sup._.sup.1931 (), z.sup.EC.sup._.sup.CIE.sup._.sup.1931 () for [380 nm; 780 nm]; spectral sensitivities of the electronic colour image sensor: x.sup.EC.sup._.sup.capteur () y.sup.EC.sup._.sup.capteur (), z.sup.EC.sup._.sup.capteur (), The superscript notations EC_XXX signify that it is situated in the XXX colorimetric space.
(6) In the figures and the exemplary embodiment described here below, the means 11 capable of emitting colours is an emissive display screen. It is understood that this means 11 capable of emitting colours may also be one or more multicolour diode(s), one or more multicolour laser(s), one or more coloured flash(es) or any other means capable of emitting colours (a colour=a spectral energy function in the range of wavelength considered). Moreover, in the figures and the exemplary embodiment described here below, the colour of the object is reduced to the uniform diffuse reflectance. It is understood that the said method is capable of capturing the other components of reflectance, in particular specular reflectance for objects with high shine and satin like shine via the display on the screen of a white figure on a black background and the detection of the resultant specular reflection lobe by a gradient analysis in the image focal plane of the colour image sensor.
(7)
(8) In one embodiment, the computer application 13 is partially or totally distributed in the network or cloud (cloud in Anglo-Saxon terminology).
(9)
E.sup.capteur().sub.i*x()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.source().sub.i)*x()*d
E.sup.capteur().sub.i*x()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.source().sub.i)*y()*d
E.sup.capteur().sub.i*x()*d=R.sup.OBJ()*(I.sup.ext()+S.sup.source().sub.i)*z()*d, by calculating the numerical value corresponding to the left-hand side of the Ei integrated equations with the use of the output parameters of the digital image sensor; and by expressing the two unknown continuous functions R.sup.OBJ() and I.sup.ext() with the use of a finite number of interpolation points (.sub.j, y.sub.j) connected by at least one interpolation function s() for maintaining the continuous nature of the said unknown continuous functions R.sup.OBJ() and I.sup.ext(), the .sub.j being selected wavelengths in the intersection of the source and sensor spectra and being input parameters of the method, chosen to minimize the number of interpolation points for a given precision; and by finding the parameters y.sub.j of the functions R.sup.OBJ() and I.sup.ext() that minimize the least squares system A*XB.sub.2 resulting from the E.sub.i integrated equations.
(10)
(11) A consequence of
(12)
(13) Now a description will be provided of the treatment process for determining the unknowns R.sup.OBJ() and I.sup.EXT() from the N flashes of colour.
(14) In the embodiment shown in
(15) Then the source 11 successively emits a succession of flashes of colours based on the flashing algorithm described in the Annexe 1 which seeks to optimise the conditioning of the matrix A, while remaining within the screen and camera gamuts and ensuring that the external illuminant remains constant over the course of the flashes.
(16) This algorithm takes as input a valid number of flashes N required. N depends on the number n of interpolation points required to determine the functions R.sup.OBJ() and I.sup.ext(). At a minimum, 3*N(2n+2) in order for the system A*X=B to be property determined (see the end of the Detailed Description), thus N2/3*(n+1).
(17) The flashing algorithm returns as output a sequence of N valid flashes, with for each flash, the input parameters of the source and the output parameters of the sensor. The index of the sequence of flashes is denoted as i.
(18) Note: If the accuracy required in order to determine the functions R.sup.OBJ() and I.sup.ext() is of the order of 10 nm, Nn30. Taking into account the reactivity of present day common electronic components, the entire flashing cycle takes about ten seconds at a maximum. For a non-uniform target (for example a texture), the method uses an algorithm for image stabilisation with at least one point of contrast.
(19) Note: If the source of the flashes of colours 11 is an emissive display screen, it is capable of generating a large number of flashes of colours, typically N being able to reach 2.sup.24=16.7 million for 38 bit RGB systems.
(20) For each flash i, the light emitted by the source 11 is reflected on the coloured target 30 and then it enters in the image sensor 12 which gives the equation (Ei): E.sup.capteur()i=R.sup.OBJ ()*(S.sup.source()i+I.sup.ext()).
Developing the equation (E.sub.i): E.sup.capteur().sub.i=R.sup.OBJ()*[S.sup.source().sub.i+I.sup.ext()]=R.sup.OBJ()*S.sup.source().sub.i+R.sup.OBJ()*I.sup.ext();
Let's write I.sup.EXT REF()=R.sup.OBJ()*I.sup.ext()
The equation (E.sub.i) becomes: E.sup.capteur().sub.i=R.sup.OBJ()*S.sup.source().sub.i+I.sup.EXT REF()
(21) Firstly, the method will interpolate the functions R.sup.OBJ() and I.sup.EXT REF () with the use of a sealed cubic spline function and (n.sup.OBJ/EXT REF+1) interpolation points of coordinates [x.sub.k=.sup.OBJ/EXT REF.sub.k: y.sub.k=y.sup.OBJ/EXT REF.sub.k] for k=0 to n.sup.OBJ/EXT REF such that: all the .sup.OBJ/EXT REF.sub.k are included in the intersection of the spectra of the sensor and the source [1.sup.mesure; 2.sup.mesure] (see
(22) the zero slope of the splines is zero at the ends: p0=p n.sup.OBJ/EXT REF=0
(23) The y.sup.OBJ/EXT REF.sub.k are the unknowns that the method will determine.
(24) As the functions R.sup.OBJ() and I.sup.EXT REF() are sealed spline functions with zero slope at the ends, R.sup.OBJ() and I.sup.EXT REF() can be written in a linear form:
(25)
(26) with: l between l and n.sup.OBJ/EXT REF such that: .sub.l1<.sub.l and l=0 if =1.sup.mesure; .sup.EXT REF(l, k, )=a.sub.k.sup.OBJ/EXT REF+b.sub.k.sup.OBJ/EXT REF*(.sub.l1)+c.sub.k.sup.OBJ/EXT REF*(.sub.l1).sup.2*(.sub.l)+d.sub.k.sup.OBJ/EXT REF*(.sub.l1)*(.sub.l).sup.2, for k=0 to n.sup.OBJ/EXT REF; a.sub.k.sup.OBJ/EXT REF=0 for k=0 to l1; a.sub.l1.sup.OBJ/EXT REF=1; a.sub.l.sup.OBJ/EXT REF=0; a.sub.k.sup.OBJ/EXT REF=0 for k=l+1 to n.sup.OBJ/EXT REF; b.sub.k.sup.OBJ/EXT REF=0 for k=0 to l1;
(27)
(28)
(29)
(30)
(31) These .sub.l,k.sup.OBJ/EXT REF are calculated by the inversion (by pivot) of the following linear system, characteristic of sealed cubic spline functions; they are only a function of .sub.l (l=0 to n.sup.OBJEXT REF):
(32)
(33) Secondly, the method knows the transfer function of the source of flashes of colours f.sup.source that gives the functions S.sup.source ().sub.l based on input parameters of the source of chrominance (R, G, B).sup.source (C.sup.source R/G/B i) and the source of luminance BL.sup.source i=C.sup.source BL i. This transfer function is determined based on the factory theoretical output values of the electronic device and/or from a calibration done prior to the measurement. It is unuseful to necessarily redo this calibration prior to each measurement of colour. This transfer function is also available for devices with 4 and more primary colours (red-green-blue-white, red-green-blue-yellow, etc).
(34) By way of Illustration, one form of the transfer function of the source is the following it being clarified that general consumer electronic devices usually seek to comply with the sRGB standard: f.sup.source(C.sup.source R/G/B/BLi)()=S.sup.source().sub.i f.sup.source(C.sup.source R/G/B/BLi)()=C.sup.source BLlinaire i*(S.sup.source().sub.i.sup.R+S.sup.source().sub.i.sup.G+S.sup.source().sub.i.sup.B); C.sup.source BLlinaire i=(a.sup.source BL*C.sup.source BLi+b.sup.source BL).sup..sup.
(35) Thirdly, the method knows the transfer function of the electronic colour image sensor f.sup.capteur that gives the chromaticity coordinates (X, Y, Z).sup.EC.sup._.sup.mesure i in the colorimetric space of measurement of the luminous flux E.sup.capteur()i entering in the sensor, as a function of its output parameters of chrominance (R, G, B).sup.capteur i=(C.sup.capteur R/G/B i) and of luminance BV.sup.capteur i=(C.sup.capteur BV i). This transfer function is determined based on the factory theoretical output values of the electronic device and/or from a calibration done prior to the measurement. It is unuseful to necessarily redo this calibration prior to each measurement of colour.
(36) By way of Illustration, one form of the transfer function of the electronic colour image sensor is the following it being clarified that general consumer electronic devices usually seek to comply with the sRGB standard: Luminance Y: The method extracts from the EXIF the Brightness Value Bv.sup.capteur for calculating the luminance B.sup.capteur of the entering luminous flux (Bv=Log.sub.2 (B/N/K) in cd/cm^2) with N=1/3.125 and K=10.7, then the method determines the luminance Y=K*B (K is a calibration parameter resulting from various losses: light scattering from the display screen, the absorption of the lens, etc). Chrominance (x, y).sup.EC.sup._.sup.mesure: firstly, the method linearises the three (RGB).sup.capteur coordinates with the use of a parametric function gamma (g, a, b) (f(x)=(a*x+b)^g) and it obtains the three (RGB_linaire).sup.EC.sup._.sup.capteur coordinates. Secondly, the method converts the 3 (RGB_linaire).sup.EC.sup._.sup.capteur coordinates into 3 (RGB_raw).sup.EC.sup._.sup.capteur coordinates by multiplication with a 33 [WB.sup.capteur] matrix which corresponds to the values of the white balance. The white balance consists of performing a chromatic adjustment to go from a D65 white (sRGB reference) to an estimated white. Thirdly, the method converts the 3 (RGB_raw).sup.EC.sup._.sup.capteur coordinates into 3 (X, Y, Z).sup.EC.sup._.sup.capteur coordinates by multiplication with the transition matrix 33 [p.sup.EC.sup._.sup.capteur>EC.sup._.sup.mesure] which corresponds to a change in the vector base to get from the colorimetric space of the sensor to a vector subs pace which is the colorimetric measurement space. Fourthly, the method converts the three (X, Y, Z).sup.EC.sup._.sup.capteur coordinates into (x, y).sup.EC.sup._.sup.capteur coordinates.
(37) Now developing the equations (E) by integrating them on the intersection of the spectra of the sensor and the source [1.sup.mesure; 2.sup.mesure] (see
(38)
(39) By interchanging the order of summations, the equations (Ei integrated) become:
(40)
(41) Let's write n.sup.OBJ+EXT REF=(n.sup.OBJ+1)+(n.sup.EXT+1);
(42) Let's write the vector X.sup.OBJ+EXT REF of dimension n.sup.OBJ+EXT REF such that: (X.sup.OBJ+EXT REF).sup.T=(y.sup.OBJ 0, . . . , y.sup.OBJ n.sup.OBJ, y.sup.EXT REF 0, . . . , y.sup.EXT REF n.sup.EXT REF); (X.sup.OBJ+EXT REF).sup.T=X.sup.OBJ+EXT REF 1, . . . , X.sup.OBJ+EXT REF n.sup.OBJ+EXT REF);
(43) Let's write the variables .sup.OBJ+EXT REF (i,k*, X/Y/Z.sup.EC.sup._.sup.mesure) such that:
(44) if 1kn.sup.OBJ+1:
(45)
.sup.OBJ+EXT REF(i,k,X/Y/Z.sup.EC.sup._.sup.mesure)=.sub.l=1.sup.n.sup.
(46) The equations (Ei integrated) may be written in the form of the following 3*N equations with the only unknowns being the X.sup.OBJ+EXT REF k:
f.sup.capteur(C.sup.capteur R/G/B/BVi)X/Y/Z.sup.EC.sup._.sup.mesure=.sub.k=1.sup.n.sup.
(47) Let's write the Jacobian matrix A of dimension (3*N, n.sup.OBJ+EXT REF):
(48)
(49) Let the vector B equal to f.sup.capteur (C.sup.capteur R/V/B/BV i) X/Y/Z.sup.EC.sup._.sup.mesure of dimension 3*N (i comprised between 1 and N).
(50) The equations (Ei integrated) then form the linear system A*X=B.
(51) The method will use the linear least squares algorithm for minimizing A.Math.XB.sub.2. The minimum is reached for: Xmin=(A.sup.T.Math.A).sup.1.Math.A.sup.T.Math.B, hence the values of the interpolation points of R.sup.OBJ() and I.sup.EXT REF () and thus I.sup.EXT()=R.sup.OBJ()/I.sup.EXT REF ().
(52) There are three conditions to be met: it is shown that the matrix A.sup.T.Math.A is invertible if and only if A is injective, which is true with
(53) .sub.l=1.sup.NS.sup.source().sub.i>0 over the interval [1.sup.mesure; 2.sup.mesure] and with h=max (.sub.k+1.sub.k) sufficiently small. The mathematical demonstration is described in Annexe 2. In addition, X.sup.OBJ+EXT REF k0 for all the k because they are the ordinates of energy flow.
(54) The first condition is satisfied by the construction of [1.sup.mesure, 2.sup.mesure] as the intersection of the source and sensor spectra (see
(55) As regards the third condition, the search to find Xmin is performed under the constraining condition X.sup.OBJ+EXT REF k0 for all the k, by using, among others, the NLLS (Non Linear Least Squares) algorithm of Lawson and Hanson (Jet Propulsion Laboratory of the National Aeronautics and Space AdministrationNASA, Solving Least Squares Problems; SIAM Editions) or, more generally, a quadratic optimisation algorithm.
(56) Note: The process is possible with other interpolation functions in linear or non-linear mode. The non-linear mode is also to be considered for moving from a system of 3*N equations based on the X/Y/Z.sup.EC.sup._.sup.mesure components to a system of N equations in the form of Euclidean type standards ( . . . .sub.2) or of deltaE type of coordinates (X, Y, Z).sup.EC.sup._.sup.mesure.
(57) Note: The method works on the one hand, with sensors that provide access to the coordinates (R,G,B_raw) and, on the other hand, with the sensors that do not provide access to the coordinates (R, G, B_raw), that do not provide access to the values of the white balance matrix, but provide the ability to lock the white balance values during the flashing process. In this second case, the values of the white balance matrix become additional unknowns to be determined (9 at the maximum). In order to solve this AX=B system with increased white balance unknowns, the technique consists of increasing the number of flashes so as to have available an over-determined system, and then to solve either in the non-linear mode or in the linear iterative mode in the following manner (inter alia, when the possible values of parameters of the white balance are comprised within a finite set of discrete values): QR decomposition of the system with 2 blocks (X | White Balance), determination of X with the first block by developing a hypothesis with respect to the parameters of white balance, injection of values of X into the second block, determination of the parameters of white balance, and then reinjection into the first block in order to iterate X, and so on.
(58) Two embodiments for viewing are possible for the invention in colorimeter mode in order to enrich the user experience and to improve the accuracy of the colour chart mode of the invention: First Embodiment: display screen 11 pointing in a direction away from the eyes of the user Second Embodiment: display screen 11 pointing in the direction of the eyes of the user
First Embodiment: Display Screen Pointing in a Direction Away from the Eyes of the User
(59) The screen 11 displays a white figure on a black background (round/bands/square . . . ) so as to materialise the place to be viewed, then the user presses on the shutter release button normally used to take a photo (see button on the side for an iPhone S or at the bottom of the screen away from the position of the video conference camera).
(60) A particular point of interest of this first embodiment is to be insensitive to the light environment which ensures the operation of the invention in virtually all the usual places, even those very brightly lit (interiors, exteriors). This quality is based on the very short distance between the target, the screen 11 and the image sensor 12, on the high brightness of the field emission display screens and on the fact that there are several barriers indeed serving as blocks to the surrounding light rays: the support for the display screen and the image sensor, hand/finger(s) of the user on the sides.
Second Embodiment: Display Screen Pointing in the Direction of the Eyes of the User
(61) The display screen 11 is placed facing towards the external illuminant, typically facing the sky. The user then comes to place on the top of the screen the coloured object to be measured (part close to the image sensor). In order to facilitate the manoeuvre, the screen 11 is divided into two parts, the upper part close to the image sensor 12 is used to emit the luminous flux for the measurement, the lower part serves to provide a feedback loop to the user on the place that is pointed (wysiwyg).
(62) A particular point of interest of this second embodiment is the ability to measure the reflectance of the screen R.sup.ecran ().
(63) The device 10 according to the present invention makes it possible to transform into a spectrometer any electronic equipment unit having a field emission display screen 11 and an image sensor 12 located on the side, including amongst others, phones, tablets, PDAs (personal digital assistants), computers and monitors/TVs equipped with a video conferencing camera. It is also intended for example for cameras and camcorders having a swivel/removable screen.
(64) The device 10 according to the present invention offers the possibility of opening up new avenues for mobile or hand held uses for accurate measurement of absolute colour in a daily context, without requiring a dedicated device. The following list is not exhaustive: the measurement of colours on objects of all kinds: materials, liquids, gases, paints, paintings, tapestries, graphics, textiles, plastics, wood, metals, soils, minerals, plants and foods, etc; the use of colour bar codes (in one or more dimensions); the measurement of colours, for medical or cosmetic purposes with respect to humans and living beings/organisms: skin, pimples, moles, hair, fur/hair coat, makeup, teeth, etc; the measurement of colours with a view to assisting people having colour blindness and/or who are blind; the measurement of the colour of object(s) appearing in a photograph in order to make chromatic adjustments and guide the white balance algorithm instead of using a supplementary colour chart of calibrated colours.
(65) The invention is described in the foregoing sections solely by way of example. It is understood that the person skilled in the art is capable of developing different variants of the invention without however in any way departing from the scope of the patent.
Annexe 1: Flashing Algorithm
(66) Introduction
(67) The present annexe defines a flashing algorithm that is aimed at optimising the conditioning of the matrix A by seeking to exclude to the maximum possible degree the said flashes and to distribute them in as homogeneous a manner as possible in the colorimetric measurement space, while remaining within the gamuts of the source and the electronic image sensor and while also ensuring that the exterior illuminant remains constant during flashes.
(68) Description of the Flashing Algorithm
(69) The said algorithm takes as input a number of required valid flashes N.
(70) The said algorithm returns as output a sequence of N valid flashes (with for each flash, the input parameters of the source and the output parameters of the sensor. The index of this sequence of flashes is denoted as i (flash_i).
(71) The said algorithm interacts with the method and the device according to the invention that emit flashes of colours towards a target object and that capture the reflected light with an electronic image sensor (including among others, a smartphone with the screen that emits flashes and the video conferencing camera which captures light reflected by the target object).
(72) The said algorithm works in the presence of an external illuminant I.sup.ext ().
(73) The algorithm is broken down into 4 main phases, and each phase into steps: Phase 1: emission of a flash black 1 Phase 2: finding the 3 vertices that maximise the flashing triangle Phase 3: homogeneous completion of the number of required valid flashes Phase 4: emission of a flash black 2
(74) Notation: In the sections that follow, the coordinates (dispR %, dispG %, dispB %) denote the chrominance input parameters of the source of the flashes of colours as a percentage ([0, 1]). DispBL % denotes the luminance input parameter of the source of the flashes of colours as a percentage ([0;1]). In practice, they are binary values, generally ranging between 0 and 255 (dispECMY, dispECMx, dispECMy) denote the coordinates of the luminous flux emitted by the source in the colourimtrique space Yxy.sup.EC.sup._.sup.mesure.
(75) Notation: In the sections that follow, the coordinates (camR %_k, camG %_k, camB %_k) denote the chrominance output parameters of the electronic image sensor as a percentage ([0;1]). In practice, they are binary values, generally ranging between 0 and 255. camBv denotes the Brightness Value which is the luminance output parameter of the electronic image sensor (see standard EXIF). (camECMX, camECMY, camECMZ) and (camECMY, camECMx, camECMy) denote the coordinates of the luminous flux emitted by the source respectively in the colorimetric space XYZ.sup.EC.sup._.sup.mesure and in the chromaticity diagram thereof.
(76) Notation: In the sections that follow, the values dispEps and camEps signify that the values of the threshold below which the input values of the source and the output values of the sensor shall be considered to be nil.
(77) Phase 1: Emission of a Flash black 1
(78) The objective of Phase 1 is twofold: to check whether the external illuminant is compatible with the physical limitations of the device and to check to ensure that external illuminant is constant during the flashing.
(79) The algorithm outputs a black flash, that is to say, it performs a measurement with the source switched off, that is dispR %=dispG %=dispB %=0 (and dispBL=0 if possible).
(80) The algorithm is faced with two possible cases:
(81) Case 1: The value camECMY is strictly higher than a limit camECMY_IEmax defined as a function of the maximum luminance of the source of the flashes of colours, which signifies that the luminance of the external illuminant is too strong compared to the physical limitations of the source of the flashes of colours. The measurement is not possible. The algorithm interrupts the flashing process by generating an error code; Case 2: The value camECMY is less than or equal to a limit camECMY_IEmax defined as a function of the maximum luminance of the source of the flashes of colours, which signifies that the luminance of the external illuminant is compatible with the physical limitations of the source of the flashes of colours. The algorithm stores in the memory the values of chrominance and luminance measured by the colour image sensor in order to ensure that the illuminant remains constant during the flashing, and then the algorithm proceeds to step 1 of the phase 2.
(82) Note: In this present Annexe 1, there is mention of only two black flashes for checking the constancy of the external illuminant, a first one at the beginning of the flashing and a second one at the end. Depending on the conditions of use of the device according to the invention, it is quite possible to reduce or increase the number of black flashes to check the constancy of the external illuminant.
(83) Phase 2: Finding the 3 Vertices that Maximise the Flashing Triangle
(84) The objective of Phase 2 is to build a flashing sub-gamut of maximum size for a given precision, while remaining compatible with the colour gamut of the source of light flashes and that of the electronic image sensor.
(85) Let's define the function f: (dispECMx, dispECMy).fwdarw.(camECMx, camECMy) of the gamut of the source of the flashes of colours to the gamut of electronic image sensor.
(86) The source-sensor sub-gamut is defined as being the intersection of the gamut of the sensor and the image by f of the gamut of the source. The source sub-gamut is defined as being the inverse image by f of the sensor sub-gamut.
(87) In other words, the objective of Phase 2 is to determine a triangle included in the sensor-source sub-gamut of maximum size for a given precision. This triangle is called flashing triangle.
(88) Step 1: Finding a Pivot
(89) The purpose of Step 1 is to define a flash pivot within the sensor gamut.
(90) The algorithm goes to step 1.1 by using a pre-defined flash pivot value for the iteration k=1, for example a white flash of maximum luminance (dispR %_1=dispG %_1=dispB %_1=dispBL %_1=1).
(91) Step 1.1: Testing of a Flash Pivot of Index k
(92) The purpose of step 1.1 is to test whether the proposed flash pivot of index k is within the sensor gamut.
(93) The algorithm asks the source to generate the flash pivot k and it retrieves the output coordinates of the sensor (camR %_k, camG %_k, camB %_k, camBv_k).
(94) The algorithm is faced with two possible cases: Case 1.1-1: At least one of the components (camR %_k, camG %_k, camB %_k) is less than or equal to camEps, which signifies that the flash is outside the gamut of the sensor. The algorithm proceeds to step #1.2 Case 1.1-2: All the components (camR %_k, camG %_k, camB %_k) are strictly greater than camEps, which signifies that the flash is within the gamut of the sensor. The algorithm goes to Step 2.
(95) Step 1.2: Generation of a New Flash-Pivot of Index k+1
(96) The purpose of Step 1.2 is to generate a new flash pivot of index k+1 that is within the sensor gamut unlike the previous flash pivot of index k.
(97) The algorithm is faced with two possible cases depending upon (camR %_k, camG %_k, camB %_k) Case 1.2-1: Only one of the three components (camR %_k, camG %_k, camB %_k) is zero with the flash pivot of index k. The algorithm will then try a new flash pivot (k+1) with a greater proportion of the component that is cancelled. By way of an illustration, if camR %.sub.k<camEps:
(98)
(99) The algorithm proceeds to step 1.1 with this flash pivot k+1 Case 1.2-2: Two of the three components are zero (camR %_k, camG %_k, camB %_k). The algorithm will then try a new flash pivot (k+1) with a higher proportion of both the two components that were cancelled. By way of an illustration, if camR %_k<camEps and camG %_k<camEps.
dispR %.sub.k+1=*dispR %.sub.k
dispG %.sub.k+1=*dispG %.sub.k
dispB %.sub.k+1=1*dispR %.sub.k*dispG %.sub.k
with >1, >1 et *dispR %.sub.k+*dispG %.sub.k<1
(100) The algorithm goes to step 1.1 with this flash pivot k+1. Case 1.2-3: the algorithm is unable to find a new flash pivot k+1. Then the algorithm declares the impossibility of performing the measurement and exits the algorithm by returning an error code.
(101) Step 2: Shortcut for Going from the Flash Pivot to a First Flashing Triangle of a Lane Size
(102) The purpose of Step 2 is to save on flashes so as to cause the increase in size of the flashing triangle around the flash pivot.
(103) To do this, the algorithm has available a database with the input of the input and output values of reference flashes and with the outputting of the coordinates of the associated flash triangles. This database is thus a shortcut for causing the flashing triangle to grow.
(104) More specifically, the algorithm frames the flash pivot with the 3 closest reference flashes (in the sense of a Euclidean norm in the measurement space). Then it mixes the coordinates of the reference flashing triangles by using the barycentric coordinates of the flash pivot in the triangle of the 3 reference flashes. The algorithm thus obtains a flashing triangle of large size around the flash pivot. It then performs a homothetic transformation or homogeneous dilation of ratio K<1 on this triangle with respect to the flash pivot so as to take a safety margin and it asks the source to emit the 3 flashes of colour corresponding to the three vertices of the said triangle. The algorithm then goes to the step 4.1 for analysing the results.
(105) If the algorithm is not able to frame the flash pivot with 3 reference flashes, it selects the closest reference flash from the database (if this latter is sufficiently close in relation to an adjustable threshold) then it causes the 3 flashes to be emitted as before and it goes to the step 4.
(106) If there is no reference flash sufficiently close in the database, the algorithm proceeds to step 3.
(107) Step 3: Passage from the Flash Pivot to a First Flashing Triangle of Small Size
(108) The purpose of the Step 3 is to create a first flashing triangle of small size by using three points, the first being the sensor output for the flash pivot (camR %_k, camG %_k, camB %_k).
(109) The algorithm orders the values (camR %_k, camG %_k, camB %_k), and then generates two new flashes, the first having more (in proportion) of the lowest component, the second having more of the two lowest components.
(110) By way of illustration, if camG %_k<camR %_k<camB %_k, the algorithm outputs a first flash having more (in proportion) of the green component, and then a second flash having more of the green component and more of the red component.
(111) Step 3.1: Testing of the Flashing Triangle of Small Size
(112) The purpose of step 3.1 is to test whether the small sized flashing triangle is included in the sensor gamut.
(113) The algorithm is faced with 3 possible cases:
(114) Case 3.1-1: at least one flash exits out of the sensor gamut. The algorithm proceeds to step 3.2. Case 3.1-2: the sensor outputs for these flashes are aligned with the sensor output for the flash pivot. The algorithm declares the impossibility of performing the measurement and exits the algorithm by returning an error code. Case 3.1-3: the small sized flashing triangle is included in the sensor gamut and the points are not aligned. The algorithm proceeds to the step 4.
(115) Step 3.2: Generation of a New Flashing Triangle of Small Size
(116) The purpose of step 3.2 is to generate a new small sized flashing triangle by replacing the flashes whose sensor output has exited out of the sensor gamut.
(117) The algorithm replaces the flash or flashes that have failed by adding less of the components added. The algorithm proceeds to step 3.1 with this new base.
(118) Step 4: Maximisation of the Size of the Flashing Triangle
(119) The purpose of step 4 is to construct the flashing triangle of the maximum size for a given precision.
(120) The transformation f is assumed to be linear and the algorithm determines this transformation by using the data from the three points of the base of the sensor gamut. The algorithm infers therefrom the source-sensor sub-gamut and then determines the flashing triangle with maximum surface area included in this sub-gamut. It then performs a homothetic transformation or homogeneous dilation of ratio K<1 on this triangle with respect to the flash pivot so as to take a safety margin and it asks the source to emit the 3 flashes of colour corresponding to the three vertices of the said triangle. The algorithm then goes to the step 4.1 for analysing the results.
(121) Step 4.1: Testing of a Flashing Triangle of Order k
(122) The purpose of step 4.1 is to test whether the flashing triangle is included in the sensor gamut.
(123) The algorithm is faced with two possible cases:
(124) Case 4.1-1: a vertex of the flashing triangle exits out of the sensor gamut. The algorithm proceeds to step 4.2. Case 4.1-2: the flashing triangle is considered to be not sufficiently large, because at least one of its vertices is too far from the vertex of the closest sensor gamut (in the sense of the Euclidean norm within the measurement space), and the flash corresponding to this vertex is also too far from the vertex of the closest source gamut. The algorithm proceeds to the step 4.3. Case 4.1-3: the flashing triangle is satisfactory. The algorithm proceeds to the Phase 3.
(125) Step 4.2: Generation of a New Vertex of Index k+1 of the Flashing Triangle Following the Exit from the Sensor Gamut of the Vertex of Index k
(126) The purpose of step 4.2 is to generate a new vertex of index k+1 of the flashing triangle within the sensor gamut unlike the vertex of index k and with coordinates (camR %_k, camG %_k, camB %_k).
(127) The algorithm is faced with two possible cases depending upon (camR %_k, camG %_k, camB %_k):
(128) Case 4.2-1: Only one of the three components (camR %_k, camG %_k, camB %_k) is zero. The algorithm will then try a new flash pivot (k+1) with a greater proportion of the component that is cancelled. By way of an illustration, if camR %_k<camEps,
(129)
(130) The algorithm proceeds to the step 1.1 with this flash pivot of k+1. Case 4.2-2: Two of the three components are zero (camR %_k, camG %_k, camB %_k). The algorithm will then try a new flash pivot (k+1) with a higher proportion of both the two components that were cancelled. By way of an illustration, if camR %_k<camEps and camG %_k<camEps, dispR %_k+1=*dispR %_k,
(131)
(132) The algorithm proceeds to the step 1.1 with this flash pivot k+1. Case 4.2-3: the algorithm is unable to find a new flash pivot k+1. Then the algorithm declares the impossibility of performing the measurement and exits the algorithm by returning an error code.
(133) Step 4.3: Generation of a New Vertex for the Flashing Triangle Because of the Insufficient Size of the Flashing Triangle
(134) The purpose of step 4.3 is to enlarge the flashing triangle because at least one of its vertices is too far from the vertex of the closest sensor gamut and the flash flash_k_1 corresponding to this vertex is also too far from the vertex of the closest source gamut.
(135) To be noted are flash_k_2 and flash_k_3 the two flashes whose sensor outputs are the other two vertices of the flashing triangle.
(136) The algorithm generates two flashes respectively barycentre of the system {(flash flash_k_1, a), (flash flash_k_2, 1)} and barycentre of the system {(flash flash_k_1, ), (flash flash_k_3, 1)} (for example: =0.2). These two flashes and the flash flash_k_1 form a triangle on which the transformation f is assumed affine. The algorithm determines this transformation by using the data from the three points and infers therefrom the source-sensor sub-gamut and then determines the vertex that provides the ability to obtain the flashing triangle with maximum surface area included in this sub-gamut. It then performs a homothetic transformation or homogeneous dilation of ratio K<1 on this point with respect to the preceding vertex so as to take a safety margin and it asks the source to emit the flash of colour corresponding to the said point. If the sensor output exits out of the sensor gamut, the algorithm repeats the operation with a higher margin (smaller K), otherwise the sensor output replaces the preceding vertex and forms along with the other vertices a new flashing triangle. The algorithm proceeds to the step 4.1 with this new flashing triangle.
(137) Phase 3: Homogeneous Completion of the Number of Valid Flashes Required
(138) The objective of Phase 3 is to generate a sequence of N flashes whose sensor outputs are distributed in a homogenous manner within the flashing triangle defined during the Phase 2.
(139) The algorithm thus has available the coordinates of the vertices of the flashing triangle as well as a certain number of intermediate flashes that were used to determine the flashing triangle during the phase 2 (at least one: the flash pivot). These points are denoted as (camECMx_k, camECMy_k).
(140) Step 1: Generation of a Grid of N Points in the Flashing Triangle
(141) The purpose of step 1 is to generate a grid of N points distributed in a homogeneous manner within the flashing triangle.
(142) The algorithm generates a grid of N points distributed as described in
(143) The algorithm proceeds to the Step 2 with this grid of N points within the flashing triangle, denoted as (camECMx_g_k, camECMy_g_k).
(144) Step 2: Generation of a Flash
(145) The purpose of this step 2 is to generate flashes whose sensor output approximates ideal points located on the grid defined during the step 1.
(146) For each point (camECMx_g_k, camECMy_g_k), the algorithm determines the triangle with vertices of three unaligned points(camECMx_k, camECMy_k) which minimizes the sum of squares of the distances from the point (camECMx_g_k, camECMy_g_k) to the vertices of the triangle. This sum is referred to as the distance to the triangle.
(147) The algorithm chooses the point (camECMx_g_k, camECMy_g_k) whose distance to the triangle is the smallest.
(148) The transformation f is assumed to be linear on this triangle and its neighbouring triangles, and it is determined using the data from the 3 vertices of the triangle. The algorithm thus determines the flash (dispECMx_g_k, dispECMy_g_k) the sensor output of which will be (camECMx_g_k, camECMy_g_k).
(149) If the sensor output for the flash (dispECMx_g_k, dispECMy_g_k) exits out of the sensor gamut, the flash is not retained.
(150) The point (camECMx_g_k, camECMy_g_k) is removed from the list of points to be brought close.
(151) The algorithm is faced with two possible cases: Case 2-1: there no longer remain any points (camECMx_g_k, camECMy_g_k) to be brought close. The algorithm proceeds to the step #3 with the set of flashes which have been generated over the course of phases 2 and 3. Case 2-2: there still remains at least one point (camECMx_g_k, camECMy_g_k) to be brought close. The algorithm proceeds to the step 2 with the new list of points (camECMx_g_k, camECMy_g_k).
(152) Step 3: Selection of Flashes to be Returned
(153) The algorithm selects only the flashes whose sensor output is sufficiently close to a point of the grid generated in step 1 (in the sense of the Euclidean distance in the colorimetric space of measurement). If two points are sufficiently close to a point in the grid, then only the closest point is selected. The algorithm then passes on to the phase 4.
(154) Note: if at least one flash whose sensor output is strictly included within the sensor gamut is selected, then the condition .sub.i (i=1).sup.TN(S.sup.source ().sub.i)>0 over the interval [1.sup.mesure; 2.sup.mesure] used in the demonstration of the invertibility of A.sup.T*A in Annexe 2 is satisfied.
(155) Phase 4: Emission of a Flash black 2
(156) The objective of this Phase 4 is to check to ensure that external illuminant is constant during the flashing.
(157) The algorithm outputs a black flash, that is to say that it proceeds to perform a measurement with the source switched off, that is dispR %=dispG %=dispB %=0 (and dispBL=0 if possible).
(158) The algorithm is faced with two possible cases: Case 1: The difference between the values of chrominance and luminance measured by the colour image sensor between the current black flash and the preceding black flash is strictly greater than a given threshold. This signifies that external illuminant has changed over the course of the flashes. The algorithm interrupts the flashing process by generating an error code; Case 2: The difference between the values of chrominance and luminance measured by the colour image sensor between the current black flash and the preceding black flash is less than or equal to a given threshold. This signifies that external illuminant has not changed over the course of the flashes and thus that the flash process is valid. The algorithm returns the sequence of N valid flashes created with for each flash, the input parameters of the source and the output parameters of the sensor.
Appendix 2: Demonstration of the Invertibility of AT*A
(159) It is shown that A.sup.T*A is invertible in the case where .sub.l=1.sup.NS.sup.source().sub.i>0 over the interval [1.sup.mesure; 2.sup.mesure] (the flashes cover all of the wavelengths considered for determining R.sup.OBJ() and I.sup.ext()), and for h=max(.sub.i+1.sub.i) that is sufficiently small.
(160) Let x be an eigen vector associated with the eigen value
(161) There is on the one hand (A.sup.T*A.Math.x, x)=(.Math.x, x)=.Math.(x, x)=x.sup.2
(162) And on the other hand (A.sup.T*A.Math.x, x)=(A.Math.x, A.Math.x)=A.Math.x.sup.2
(163) Thus =A.Math.x.sup.2/x.sup.2
(164) Thus is positive and it is so strictly if A is injective.
(165) Thus if A is injective, all the eigen values of A.sup.T*A are strictly positive, therefore A.sup.T*A is invertible.
(166) Conversely, if A is not injective, there exists non zero x such that A.Math.x=0. Thus one then gets A.sup.T*Ax=0 therefore A.sup.T*A is not invertible.
(167) Finally, A.sup.T*A is invertible if and only if A is injective.
(168) Lemma 1: The Integral of a Sealed Cubic Spline S() with Zero Slopes at the Ends [S()*d is Always Positive when all the y.sub.i are Positive
(169) Case No 1: If all the y.sub.i=0,
(170) Then s()=0 for all the due to the zero slopes at the edges, hence Int{S()*d}=0
(171) Case 2: If there Exists One y>0 and if all the Other y.sub.i are Zero:
(172) y_i creates a positive belly over [.sub.i1; .sub.i][.sub.i; .sub.i+1] (area >0) and negative bellies over [.sub.i2; .sub.i1] and [.sub.i+1; .sub.i+2] (area <0), and so on up to the edges. As the spline function minimizes the energy of the tongue, the surface area of the positive belly [.sub.i1; .sub.i][.sub.i; .sub.i+1] is greater than those of the negative bellies [.sub.i2; .sub.i1] and [.sub.i+1; .sub.i+2]. This is due to the fact that |p.sub.i+1|>|p.sub.i+2|> . . . >|p.sub.n|=0 and |p.sub.i1|>|p.sub.i2|> . . . >|p.sub.0|=0. As a result, the surface area of the positive bellies is greater than the negative bellies and therefore S()*d>0
(173) Case 3: If there Exist Two y.sub.i>0:
(174) Case no 3.1: y.sub.i and y.sub.i+1>0 (contiguous): This creates a large positive belly. The reasoning given for the case no 2 is applicable.
(175) Case no 3.2: y.sub.i and y.sub.i+2>0 (y.sub.i+1=0): there are 2 positive contiguous bellies. The reasoning given for the case no 2 is applicable.
(176) Case no 3.3: y.sub.i and y.sub.i+3>0 (y.sub.i+1=y.sub.i+2): there is 1 positive belly, 1 negative belly, 1 positive belly. The reasoning given for the case no 2 is applicable.
(177) Case no 3.4: y.sub.i and y.sub.i+k>0 with k>3. Same as previously given for Case no 3.3.
(178) Case 4: If there Exist at Least 3 y.sub.i>0 (General Case)
(179) Thus the reasoning given for the cases under Case no 3 is applicable. QED
(180) Lemma 2: The Integral of a Sealed Cubic Spline S() with Zero Slopes at the Ends Multiplied by a Function K()>0 ([S()*K()*d) is Always Positive when all the y.sub.i are Positive
(181) Lemma 1 (case no 2) is not directly applicable because the term K() can greatly reduce the area of the positive bellies and increase that of the negative bellies.
(182) The trick is to increase the number of interpolation points in order to reduce the surface area of the negative bellies.
(183) Let us be positioned in the scenario of case no 2 of the Lemma 1. The error of a sealed spline is bounded by the following formula: |f(x)s (x|a*h^4 with a=5/384*max[a; b]{|f(4)(E)|}>0 (constant value), and h=max{|x.sub.ix.sub.i1|)}.
(184) That is to say that: f(x)a*hs(x)f(x)+a*h
(185) Since f(x)0 (R.sup.OBJ() and I.sup.EXT() are the energy flows), s(x)a*h, thus the maximum surface area of the negative bellies [.sub.i2; .sub.i1][.sub.i+2] is equal to 2*a*h^2.
(186) The surface area of the positive bellies remains constant because the increase in the number of interpolation points will create multiple contiguous y.sub.i>0.
(187) As a consequence, there is an h for which the surface area of the positive bellies is strictly greater than that of the negative bellies.
(188) Returning to the Demonstration
(189) Let x be such that A*x=0 with x formed of the ordinates of the sealed cubic spline with zero slopes at the ends representing R.sup.OBJ() and I.sup.ext().
(190) Showing that x=0
(191) We get for all i=1, . . . , N and for any sensitivity x/y/z():
R.sup.OBJ()*S.sup.source().sub.i*x/y/z()*d+I.sup.ext()*x/y/z()*d=0
where R.sup.OBJ() and I.sup.ext() are the sealed cubic spline functions with zero slopes at the ends.
Summing for the 3*N equations, we obtain:
(192)
With x()+y()+z()>0 and .sub.i=1.sup.NS.sup.source().sub.i>0 for all , it is deduced from Lemma 2 that there exists a sufficiently small h such that:
(193)
and
I.sup.ext()*(x()+y()+z())*d=0
Since R.sup.OBJ()=0 and I.sup.ext()=0 for all ,
ie x=0. A is injective. Therefore A.sup.T*A is invertible. QED