AU695071B2 - Colour conversion method - Google Patents

Description

V

S F Ref: 327990

AUSTRALIA

PATENTS ACT 1990 COMPLETE SPECIFICATION FOR A STANDARD PATENT

ORIGINAL

e C. C

C

Name and Address of Applicant: Actual Inventor(s): Address for Service: Invention Title: ASSOCIATED PROVISIONAL [31] Application No(s) PN1133 Canon Information Systems Research Australia Pty Ltd 1 Thomas Holt Drive North Ryde New South Hales 2113

AUSTRALIA

William Clark Naylor Jnr.

Spruson Ferguson, Patent Attorneys Level 33 St Martins Tower, 31 Market Street Sydney, New South Wales, 2000, Australia Colour Conversion Method APPLICATION DETAILS [33) Country

AU

[32) Application Date 14 February 1995 The following statement is a full description of this invention, including the best method of performing it known to me/us:rni, -1A- COLOUR CONVERSION METHOD Field of the Invention The present invention relates to methods for converting from one colour space to another colour space and to a printing process whereby output is printed on a device having a limited or restricted gamut of colours.

Background Art Each printing device for the printing of colour images has a measurable gamut or range of colours which it is able to reproduce from any primary colours that it may utilise. This gamut will vary from output device to output device, in addition to some variation between output devices of the same kind.

Colour images are often stored and printed in different colour spaces. For example, colour images displayed on a CRT-type device are normally stored in a Red, Green and Blue (RGB) format. This corresponds to the relevant primary colours (RGB) utilised in most CRT-type displays. Colour printing devices normally work in a subtractive colour space that comprises the Cyan, Magenta, Yellow (CMY) or Cyan, Magenta, Yellow and blacK (CMYK) colour spaces. Methods for conversion from one colour space to another are well known, with common methods utilising a sparse lookup table and interpolation between data points.

6 The CIE co-ordinate system defines three standard primaries Y and Z) that 20 represent the amount of spectral energy of a particular colour. Instruments, called colourimeters, can measure the tri-stimulus values of X, Y and Z. Methods are know for mathematically transforming from RGB space or CMYK space to XYZ space and vice versa.

Unfortunately, these methods often require multiple t:)1bi' iookur -in addition to 25 multiplication operations that are required to be performed as part -nterpolation process. Therefore, if it is desired to perform the conversion procc. ithin a standard microcomputer system it is likely that these known methods will take unacceptably long time to complete the colour conversion process of a whole image.

Additionally, as mentioned previously, each output device has its own range or 30 gamut of possible colours that it may print. Therefore, even if the colour conversion process has been carried out from one colour space to another, there is no guarantee that the colour derived in the second colour space will fall within the gamut of the output device. Further processing will be necessary in order to map the unprintable colours to those which are within the gamut of the printing device.

On the one hand, if it is desired to faithfully reproduce an image, extra inks and unique printing processes must be employed. If the goal is to make a pleasing rather than exact reproduction, it has been found that it is unnecessary to reproduce the exact image colour. It is of some importance, however, to attempt to maintain the relative differences between different initial colours.

1( -2- Referring now to Fig. 1, there is shown an example of a conventional CIE chromaticity diagram. The gamut area 1 represents the gamut output by a CRT-type display device. The gamut 2 could represent the gamut of colours printable on a paper printing device such as an ink jet printer.

It can therefore be seen that there is a problem when it is desired to print out an image displayable on a CRT-type device on, say, an ink jet printing device. The ink jet printing device will be unable to print colours e.g. 3 which are outside of its gamut 2 but within the gamut 1 of the CRT-type device.

Summary of the Invention It is an object of the present invention to provide a simple and efficient means for mapping from a first colour space to a second colour space.

In accordance with one aspect of the present invention there is provided a method for displaying colour images, intended to be displayed on a first colour display device, on a second colour display device having a colour gamut different from the first colour display device, said method including the steps of: determining a plurality of input extrema colour values, within the colour gamut of said first colour display device and at the extrema of the gamut of said first colour display, determining a plurality of output extrema colour values corresponding to said input extrema colour values, within the colour gamut of said second colour display and at the extrema of the gamut of said second colour display, determining a mapping for substantially all the values within the colour gamut of the first colour display to substantially all the values within the gamut of the second colour display, said mapping substantially preserving the relative value of differences between colour values relative to the extent of the gamut of said first colour space and differences between corresponding colour values in the gamut of said second colour space relative to the extent of the gamut of said second colour space.

Brief Description of the Drawings A preferred embodiment of the present invention will now be described with reference to the accompanying drawings in which: Fig. 1 illustrates the CIE gamut space of different output devices; Fig. 2 is a CIE diagram illustrating a first methou cf mapping from one printer device colour space to a second printer device colour space; Fig. 3 is a CIE diagram showing a second method of mapping from one printer device colour space to a second printer device colour space; Fig. 4 illustrates, in three dimensions, the different colour spaces of various output devices; Fig. 5 illustrates a barrier function used in the numerical optimisation procedure utilised in the preferred embodiment; and 4 0C 4. 327990 BBLE3 CFP0295AU [o-cisra\bbaefflrhhWA'-,', n -3- Fig. 6 illustrates the problem of "dot gain" experienced by printing devices.

Detailed Description of the Preferred Embodiment In the preferred embodiment of the present invention, an attempt is made to preserve the relative colour differences between any two absolute colours at the expense of a remapping of the colours o be printed.

Referring now to Fig. 2 there is shown a first method of mapping from the CRT type gamut 1 to the ink jet ga-mut space 2. In this method, points outside the ink jet gamut space 2 are mapped to the closest point on the surface of the gamut space 2. For example, the colour 5 which is within the coloured gamut 1 of the CRT type device will be mapped to the colour 6 which lies on the edge of the gamut 2 of the ink jet printing device. Similarly, the colour 7 which is outside the ink jet printer gamut 2 is mapped to the colour 8, which is again on the edge of the gamut 2 of the ink jet type printing device. Those colours e.g. 10 which lie within the gamut space 2 of the ink jet type printing device are mapped to themselves and their values are not altered. The mapping of Fig. 2 is designed to optimise "absolute" colour accuracy. Those colours which are inside the ink jet colour gamut 2 are reproduced perfectly. Those outside the ink jet colour gamut are "clamped" to the nearest points on the boundary.

Referring now to Fig. 3, there is shown the method of colour gamut matching utilised in the preferred embodiment where the strategy adopted is to optimise the mapping to make relative distances between different colours in the image remain the same. The Red R, Green G and Blue B "corners" of the CRT gamut 1 are mapped to the Red Green G' and Blue B' corners of the ink jet gamut 2. This mapping is carried out regardless of the absolute colour values of these corners. Those colours in the interior of the CRT gamut are mapped so as to preserve their relative differences from the R, G, and B corners as well as from each other. For example, in Fig. 3, the point G is mapped to the point the point B is mapped to the point B' and the point R is mapped to the point Additionally, the point 10 will be mapped approximately to the point 11 and the point 12 will be mapped approximately to the point 13.

Referring now to Fig. 4, there is shown the extension of the ideas of Fig. 3 into .i 30 three dimensions in XYZ space. A first colour cube 15 represents the RGB colour cube of displayable colours by a CRT type device. Of course, it will be evident to those skilled in the art, that colour "cube" 15 may not necessarily be an exact cube and can take the form of a parallelepiped or other volume.

Similarly, colour "cube" 16 represents those colours displayable by an ink jet printing device. Colour correction cube 16 can be "skewed" or "rotated" relative to colour correction cube The ink jet printing device normally relies on four colour inks, being Cyan, Magenta, Yellow and blacK (CYMK). The black portion of ink is normally obtained in accordance with the following equations: 327990 BBLE3 CFP0295AU [o:\cirabble.dirbte3]speci.doc:LDP v 4- K =min Y) Equation 1 j r: C K M M' K Equation 2 However, an ink jet printing device is normally a "binary" type device in that normally only two levels of output are possible for each primary colour output comprising either the presence or absence of ink. Therefore, the normal equations 1, 2 are modified to take into account the ink jet printer's binary nature such that if C, M and Y are all to be turned on they are replaced by K. This can be expressed by the following pseudo code: if ON) (M ON) (Y ON)) K <-ON; C <-OFF; M <-OFF;

Y<-OFF;

else K <-OFF; leave CMY unchanged A number of assumptions are made in the implementation of the colour mapping of the preferred embodiment. Firstly, it is assumed that all points in the interior of both colour cubes can be physically achieved by mixing various amounts of the corner colours and, secondly, that no point exterior to each cube can be achieved physically by mixing colours from its corners. These assumptions are not often valid as measurements with a colourimeter can easily show that the assumptions are violated in practice. Therefore, methods will be discussed to allow the minimisation of errors caused by physical deviations from these theoretical assumptions.

Further, it will be assumed that there is a one to one, non-linear mapping function which models the distortion caused by the half toning and/or printing by the ink jet printing device. As such a function exists, an inverse must also exist and this inverse is the colour correction function which is required. If such an inverse can be found, then, theoretically, it can be used to correct all the colour distortions of the ink jet printing device resulting in a printing system which will preserve relative colour. In practice, it is only possible to approximate this inverse function. The method of constructing an inverse approximation consists of the following four steps: 1. A functional form of the inverse is assumed to exist. This functional form will map the point from the three dimensional RGB space of input values to points on a mathematically identical RGB space. The INPUT points to this function will represent the OUTPUT of the ink jet printe,; system. The OUTPUT points of this function will represent the INPUT to the ink jet printing system. It is further assumed that the eight corner points (Fig. 4) will map to themselves. The function will further have a number of free parameters, which are set to optimise the fit between the function and measured experimental data. The function can therefore be written as: (Rin, Gin, Bin) colour_correction (Rprint, Gprint, Bprint; cO, cl, Equation 3 where cO, cl etc represent the free parameters, Rprint represents the red portion of the colour space value of the ink jet printer, and Rin presents the red portion of the value "input to the ink jet system.

2. The various RGB colour sample points are printed utilising the target ink jet printing system and their corresponding colour values are measured using a 20 colourimeter. The points printed and measured should include the eight corners of the RGB colour cube 16 of Fig. 4. The measurements taken by the colourimeter can be utilised in fitting the colour correction function. One suitable method is to print out a 5x5x5 grid of equally spaced colour points in the input colour space of the ink jet printing system.

3. The sample points printed out are then measured in a linear colour space Si" such as XYZ. If measurements are required to be formed in a non-linear form of colour space, such as CIE or Lab, the measured numbers must first be transformed into a linear space.

In order to fit the colour correction function (Equation it is necessary to 30 express these XYZ numbers as RGB coordinates. As set out previously, it is assumed that the corner points of the colour cube 15 are mapped to the corner points of the RGB colour cube 16. The other printed colours are assigned coordinates in RGB space by comparing their measured colour values with the measured colours of the corners of the RGB cube The sample data points printed by the above process and measured by the colourimeter will consist of pairs of Rin, Gin, Bin input values and Xprint, Yprint, Zprint output values as follows: -6- (Rin, Gin, Bin) (Xprint, Yprint, Zprint) [i] Equation 4 Of course, this assumes that the points are measured in a linear space such as XYZ, however, methods are known to transform from a nonlinear space such as CIE or Lab to XYZ.

In order to fit the colour correction function, the output colours must be expressed in RGB coordinates. Therefore, it is required that for all data points, the pairs as set out in "quation 4 be transformed into the following form: (Rin, Gin, Bin) (Rprint, Gprint, Bprint) [i] Equation As mentioned previously, it is assumed that the corner points of colour cube are correctly mapped to the corner points of colour cube 16. Thus, as a result of this assumption we have the following relationship: r r Ir r 1 r r ri r r s r r r r i rr r D r crfcs~ r r (Rin=0, (Rin=0, (Rin=0, (Rin=0, (Rin 1, (Rin=l, (Rin= 1, (Rin 1, Gin=0, Bin=0) Gin=0, Bin=1) Gin= 1, Bin=0) Gin= 1, Bin= 1) Gin=0, Bin=0) Gin=0, Bin=1) Gin=1, Bin=0) Gin=1, Bin=1) (Rprint=0, (Rprint=0, (Rprint=0, (Rprint= 0, (Rprint= 1, (Rprint= 1, (Rprint 1, (Rprint 1, Gprint=0, Bprint Gprint=0, Bprint Gprint= 1, Bprint Gprint= 1, Bprint Gprint=0, Bprint Gprint=0, Bprint Gprint= 1, Bprint Gprint= 1, Bprint Equation 6 As colourimeter measurements have been made of these corner points, the colourimeter measurements of the other measured sample points can be expressed as 25 convex linear combinations of the above corner point measurements. The convex linear combination weightings or coefficients can then lead to an assignment of RGB print values to these measured points.

t

A

327990 BBE3 cFP0295AU [o:\cisra\bble.'iir\bble3]speci.doc:IAD 7- If the colourimeter measurements for the eight corner points are as follows: (Xprint, Yprint, Zprint) [0,0,0] (Xprint, Yprint, Zprint) [0,0,1] (Xprint, Yprint, Zprint) [0,1,0] (Xprint, Yprint, Zprint) 1] (Xprint, Yprint, Zprint) [1,0,0] (Xprint, Yprint, Zprint) [1,0,1] (Xprint, Yprint, Zprint) [1,1,0] (Xprint, Yprint, Zprint) [1,1,1] Equation 7 and if the printed RGB gamut is assumed to be of a cubic nature, and the measurement space is assumed to be linear, it is possible to calculate the XYZ measurements for other RGB print values by the use of trilinear interpolation which can be determined as follows: (Xpredict, Ypredict, Zpredict) XYZpredict (Rprint, Gprint, Bprint)

S..

or

L

4i i (1-Rprint) (1-Rprint)* (1 -Rprint) (1-Rprint)* Rprint Rprint Rprint Rprint (1-Gprint) (1-Gprint) Gprint Gprint (1-Gprint) (1 -Gprint) Gprint Gprint (1-Bprint)* Bprint (1-Bprint)* Bprint (l-Bprint)* Bprint (1-Bprint)* Bprint (Xprint, (Xprint, (Xprint, (Xprint, (Xprint, (Xprint, (Xprint, (Xprint, Yprint, Zprint) Yprint, Zprint) Yprint, Zprint) Yprint, Zprint) Yprint, Zprint) Yprint, Zprint) Yprint, Zprint) Yprint, Zprint) [0,0,0] [0,0,1] [0,1,0] [0,1,1] [1,0,0] [1,0,1] [1,1,0] [1,1,1] Equation 8 Since measurements have been taken for various (Xprint, Yprint, Zprint) it can be assumed that Equation 8 can be utilised to estimate values for the variables Rprint, Gprint, Bprint by minimisation of the error in the predicted (Xpredict, Ypredict, Zpredict) output value. The error between predicted and actual values can be expressed as an error function of the following form: error (Rprint,Gprint, Bprint) j(Xprint, Yprint, Zprint) XYZpredict (Rprint, Gprint, Bprint) 1 2 Equation 9 Conjugate direction methods can then be utilised with respl, to Equation 9 to choose (Rprint, Gprint, Bprint) so as to minimise the error over all the data samples.

CCC.

.C.

-8- Conjugate direction methods are numerical methods utilised for minimising non-linear functions in multiple variables. Conjugate direction methods are well kno-'n to those skilled in the art of numerical analysis and appear in a number of standard texts. For example, reference is made to the text "Linear and Nonlinear Programming" by David G Luenberger, Second Edition, 1984, Addison-Wesley Publishing Company, Reading, Massachusetts and to the standard text "Numerical Recipes in C, The Art of Scientific Computing" by Press et al, published by Cambridge University Press 1992.

As the colour cube 16 is not a perfect cube and as a result of a number of the assumptions made, and because of noise in the measured data, the application of well known conjugate direction methods directly to the above formulas can yield values which are well outside the admissible range of colour cube 16. Simply moving these values to the closest points on colour cube 16 was found to not yield the best RGB values in the admissible range. A better method is to utilise a barrier function in the conjugate gradient direction method optimisation procedure. The utilisation of a barrier functions is more fully described in the above texts. The barrier function can be used to confine the RGB values to the admissible range during optimisation. A barrier function of the following form is suitable: barrier(x) if (x<0) result 0 if (0<x<1) 1 l-x 2 result co Equation The barrier function takes the form as shown in Fig. 5 with an asymptotic value of plus infinity at x 1.

By adding the barrier function to the error objective, the objective function is caused to go to positive infinity for small deviations outside the admissible colour range. The new objective function is therefore: objective (Rprint, Gprint, Bprint) error (Rprint, Gprint, Bprint) barrier (s (Rprint barrier Rprint) barrier (s (Gprint barrier Gprint) barrier (s (Bprint barrier Bprint) Equation 11 S C t V t C L i.

CrI( where is a scaling factor which determines how far outside the admissible colour range the solution values are allowed.

Setting s to positive infinity will, in theory, prevent any deviation from the admissible range. However, this also has the known effect of driving the condition number of the corresponding Hessian of the optimisation problem to infinity, which will dramatically slow or stop convergence of the conjugate direction methods.

Experiment has found that a value of s =300 is a good compromise value to use.

The objective of equation 9 is easily minimised using a conjugate direction method. In the present case, Powell's conjugate direction method was found to be adequate to optimise equation 11. The optimisation is performed for every measured sample point, to give a (Rprint, Gprint, Bprint) vector. These vector values are then used to fit the colour correction function.

4. Fitting the colour correction function.

The fitting of the colour correction function to the measured data is a nonlinear process and, as such, a nonlinear fitting method must be utilised. To fit the colour correction function, an error function is again constructed, however this error function measures the average deviation of predicted data from the measured data. Conjugate direction methods can then be used on this error function in order to find its minimum.

The error function utilised was a simple sum of the squares error function of the 20 following form: error(c, cl, SI(Rini, Gini, Bini) colour_ correction (Rprinti, Gprinti, Bprinti; cO, cl, 2 all measured data points i Equation 12 The free parameters cj of the above formula can be adjusted to give the minimum error using a conjugate direction method. The colour correction function maps points from the three dimensional RGB space 1 0<r 1 0<b_ 1 0 b 1 Equation 13 to points in a mathematically identical RGB space. As noted previously, the INPUT points to the colour correction function represent the OUTPUT of the ink jet printing system, and the OUTPUT points of the colour correction function represent the INPUT to the ink jet printing system. It is further assumed that the eight corner points map to 327990 BBLES CFP0295AU [o:\cisra\bble.dir\bble3]speci.doc:IAI) ii .c r themselves. Further, as noted previously the function has a number of free parameters cj which are required to be optimised.

The colour correction function is the inverse of a print distortion function which can be written in the following form: (Rprint, Gprint, Bprint) printdistortion (Rin, Gin, Bin) Equation 14 In the present embodiment, the form of the print distortion function is first derived and then inverted to get the colour correction function.

Several considerations enter into the choice of the colour correction function's actual form. These considerations include: 1. The function should be able to accurately approximate c'.cal realit 2. The function should be able to be efficiently prgrammed for incorporation into driver software for a microcomputer or the like so that an inordinate delay does not occur in the printing out of images.

3. The function should have as few free parameters, cj, as possible so that fewer measurements are needed to obtain a good fit to the sample data without overfitting occurring.

One popular colour correction function, utilised with CRT type devices, is the "gamma correction" function which takes the following form: 1 Rin Rprint Y r Gin Gprint' Bin Bprint' b Equation However, this function was found to give poor quality results when applied to printing devices such as ink jet printers. This is mor rrobably as a result of its intended purposes, which is in use for CRT type displays and hence there is no physical reason why it should produce good results for an ink jet printing system. A more thorough analysis can therefore be undertaken as to the physicai operation of a printer such as an ink jet type printer. These physical principles of operation can then be used to construct a more suitable print distortion function which can then be inverted to get the colour correction function.

The most significant distortion in devices such as ink jet printing devices is that caused by "dot gain". Inkjet printing devices, such as the Canon BJ600, work on the i

I

-r-7M BBLE3 CFM0295AU 3o:.isralbble.dirbble3lieci.doc:LDP -11principal of squirting ink onto paper in the form of a large array of dots or pixels.

Commonly, these devices can squirt 300 or 400 dots per inch. Methods, commonly known as halftoning, are utilised so that the collation of the dots printed approximates a desired full coloured image. In a first, oversimplified representation, each dot can be considered uniform and modelled for example as a uniformed sized round dot 40 of Fig. o. Unfortunately, in practice, when a single dot is ejected onto paper or the like, it is absorbed by the paper and "spreads out" such that it has a non-uniform boundary and each dot has a variation in size as shown by the dot 41 of Fig. 6. Also, when a series of dots 42-44 are printed close together, the moisture from a first dot 42 can effect subsequent dots e.g. 43 such that their combined area is again substantially increased. These problems are well known to those in the printing art, and fall under the general heading of "dot gain".

The problems of dot gain appear to be physically unavoidable and the amount of dot gain depends on the details of the arrangement of the dots on a page. A very simple and robust assumption is that the dots of ink are projected onto the page in a random distribution. The random distribution assumption leads to the following differential equation for the amount of colour produced by a given amount of ink squirted on a page: dC GdI 1-C Equation 16 In this equation C is the amount of colour, I is the amount of ink, and G is the amount that each dot "spreads out" when printed and can be considered a measure of the dot gain. This differential equation states that the fraction of uncoloured or uncovered paper which gets coloured by a small amount of ink is given by the dot gain 25 G. The differential equation of Equation 16 can be easily integrated to yield C 1 e-GI Equation 17 with the boundary condition that C 0 when I 0. Equation 17 in its present form implies that the colour coverage never reaches 100%. Equation 17 can be renormalised such that colour coverage at I 1 is equal to 100%. This renormalisation takes the following form: *0 t

I

1-e-G Equation 18 Equation 18 -12- Thus Equation 18 represents a form of the print distortion function. In order to get the colour correction function, Equation 18 must be inverted. This inversion takes the following form: In e- C)

G

Equation 19 Equation 19 applies only to one colour of ink. Inkjet printing devices, printing on, for example white paper, utilise a "subtractive" colour printing process. One form of well known full colour subtractive colour printing process utilises three colours of ink, being Cyan, Magenta and Yellow with methods for transferring from one colour printing process to the other being well known.

The utilisation of three subtractive colour inks yields the following series of equations: -In e G) cyanprint) cyanil G -In magentaprnt) magentai, -In yellowprint) yellowi =G SEquation The cyan, magenta and yellow CMY values can be easily converted to and from corresponding RGB values utilising the following relationship.

1 R M' 1 G Y' 1 I B Equation 21 S. Reasonably good results are produced when utilizing Equation 20 as the colour correction function is free of the parameter G ad the results produced were found to be reasonably good. Equation 20 has the advantage of being easy to fit using conjugate gradient methods as it only has one free parameter G. Further, it would be easy to 25 implement in inkjet printer driver software. However, Equation 20 has a number of significant disadvantages. These include: The conversion of CMY to K is totally ignored in the above equations.

327990 BBLE3 CFP0295AU [o:\cisra\bble.dir\bble3]speci.doc:IAD

-I

V$

II

-13- (ii) The dot gain is not constant but rather appears to be a function of the total amount of ink or moisture on the page. More ink on a page means that more moisture is present and this normally tends to change the tendency for a particular ink to interfuse into the paper. This unwanted effect tends to destroy the possibility of independently correcting for each colour.

These effects can be modelled by making the dot gain a function of the amount of ink on the page. The amount of ink on the page can be roughly modelled as follows: kprint min (cyanprint, magentaprint, yellowprint) ink (cyanprint kprint) (magentaprint kprint) (yellowprint kprint) kprint Equation 22 The functional form of the relationship between the dot gain G and the total amount of ink is not known. However, assuming that the functional form is reasonably well behaved, it will have a Taylor's series expansion of the following form: G(ink) go gi ink g 2 ink 2 Equation 23 where go, gl etc are unknown constants. This Taylor's series expansion can be truncated and the parameters fitted to the certain measured colourimeter data. S ubstantial improvement can result from using the second linear term, but no substantial improvement normally results from using higher order terms. Therefore, the Taylor's 20 series expansion of Equation 23 is truncated to take the following form: G(ink) go g, ink Equation 24 a 9 i 327990 BBLE3 CFP029U [o:\isra\bble.dirbble3]speci.do:IAD i -14- The final form for the colour correction function is therefore: (Rin, Gin, Bin) colour_ correction (Rprint, Gprint, Bprint; go, g Internals: cyanprint 1 Rprint magentaprint 1 Gprint yellowprint 1- Bprint kprint min (cyanprint, magentaprint, yellowprin) ink (cyanprint kplin) (magentaprint kprint) (yellowprint- kprint) kprint G(ink) go g* ink cyanin single_color_correct (cyanprint) magentain single_color_correct (magentaprint) yellowin singlecolor_correct (yellowprint) where -In

I)

single_ color_ correct

G

Rin 1- Cyanin .Gin 1 magentain Bin 1 yellowin S: Equation 20 Above function is a reasonably simple one and has only two free parameters, go and gl. Fitting go and gl to a set uf 125 measured colourimeter values can be done in reasonable computer time using Powell's conjugate direction method. Tests have ri shown that the results yield excellent colour reproduction on ink jet printing devices such as Canon Inc's BJ600. Different types of paper will produce different dot gains and hence lead to different values for go and gl.

The foregoing describes only one main embodiment of the present invention, with minor modifications. Further modifications obvious to those skilled in the art, can be made thereto without departing from the scope of the invention.

SAppendix 1 sets out a C code program for performing the first Powell method S 30 optimisation, Appendix 2 sets out C code source for performing the second Powell ,1 optimisation, Appendix 3 sets out various non-standard library functions utilised by Appendix 1 and 2.

327990 TBL3 CFP0M2SAU [o:\cisra\bble.dir\bble3speci.doc:IAD Appendix 1 #include <stdio.h> #include "wnlib.h 1 #include "wnswap.h" #define NUMDIVISIONS #define NUMSAMPLES (NUMDJVTIONS*NUM DIVISIONS*NUMDfVISIONS) local double **i1n,**out; local double cub e[2] local void triliniear-initerpolate(out-vect, in-vect) double out-vect[3],in-vect[3];

S

*555 4555

S

555.

'S.

4* S. S S 5* 55*5 e *45.

0*

S.

*5

S

int il,i2,i3; 25 double wl,w2,w3; wn-zero vw,-ct(out-vect,3); for(il =0;il <22; +il) 30 if(il 1 wi in vect[0];, 35 else wi i.0-in vect[0];

A,

-16for(i2=O;i2 ++i2) if(i2 1) w2 =in-vect[1]; else w2 for(i3 =O;i3 2; ++i3) if(i3 1) w3 in-vectr2]; eIlse w3 1 .0-in-vect[2];, wn-add-scaled vect(out-vect,cube[i1] [12] [i3] ,w1 *w2*w3 local void make-cubeO iut il,i2,i3,index; int r-i,g-i,b-i; 35 for(il =O;il for(i2=0;12<2;++12) for(i3 =0;i3 r i 12;

S

Sb ten..

ft V 327990 BBLE3 CFP0295AU 32790 BBE3 FPO25AU[o:\cisra\bble.dir\bble3I.neri dn,.4TAT) -17bi ii; index (NUMDLVISIONS-1)* (NUMDIVTSIONS*NUMDIVISIONS*b i +NUMDIVISIONS*r-i g-i); printfQ"(%d, il, i2, B3,index); wn copyvect(cube[r i] [b_i],out~index] #define SCALE (300.0) local double barrier ev.al(vect) double vect[3]; extern double barriero; m i; double sum;

W..

sum 0.0; for(i=0;i<3;++i) ~0 :sum =barrier(-SCALEIect[i]); sum return(sum); local double *save out sample; local double function(vect) double -vectli3]; extern double wn dist2O; double actual([3], dis t, bar; trilinear iiterpolate(actual ,vect); dist wn dist2(actual ,save out sample, 3); bar barrier-eval~vect); return(dist+ bar); local void bound num(pn) double *pfl; jf(*pn 0.0) =p 0. 0; else jf(*pn 1 .0) *pn CC~cco9 local void bound-out sample(out sample) double out sample[3]; int i; for(i=0;i bound num(&(out sample[i])); local void renormalize-samnple(out samnple) double out sample[3]; int i; int code; double vect[3]; double val; for(i=O;i vect[i] save out-sample out-sample; conj_direction -method(&code, &val ,vect, 3, (function) ,WNIHUGE); printf("code %d,val %lg\n",code,val); wnprint-vect~vect,3); *66: wn-copy vect(out sample,vect,3); bound-out-samplc(out sample); 35 local void renormalizeO 327990 BBLE3 CFP0295AU 32790 BLE CFO29AU[o:\cisra~bbk,..,;r-lbhe3]speci.doc:IAD -20-1 ilt i; make-cubeQ; for(i=0;i <NUMSAMPLES; printf("i renormalize-Sample(out[il); local void read data_ine(px,py,pz) double *px,*py,*pz; ilt linenum; scanf(" %d X %lg Y %lg Z %lg\n",&linenum,px,py,pz);

I.

I

I*0~

I..

0*C 4 9* C local void read-dataO 25 double in-red ini-areen in-blue,x,y,z; ilt count; count 0; C. a

C.

eC S

C.

*0 C~ C. o e for(inblue in-blue 1. 0000 1; inblue 0- for~in.red=.O;in -red <1.000 1; n-red 01(5. 0-1 for(in,_green 0; ingreen 1. 0000 1; in,_green+ 01(5. 0-1. 35 read-data line(&x,&y,&z); in[count]J[0] =in-red; in[count]I[1] in-green; in[countl[2] in,_blue; 20700f) 1DnTrP~ ('T1;ATT 21out[count][O] X out[count][11] Y out[count][2) z; 1 out[count] [0] out[count][1] out[couW] [2] 100.0*in-redU', 20. 0 *111blue; prin.tfQ'in %lg, out %lg, in[courit] ,in[count][1] ,in[count] out[count] [0],out[count][1] ,out~count] +count; wn-assert(count =NUM SAMPLES); renormalizeo; t 4 4~R~ p..

K

4**

C..

0*a te *a.

4 a wn-swap (in, out, double for(count count< NUMSAMPLES; -F +count) {rnf"(I,%f %f in [count] in[count] [1],in[count] out[count] [0],out[count] out[count] void maino e=r double wn-random-double-betweenO,liypoto,atan2O,wn-fmodO; nnet~type -nnet; int i; 22wn-gpmake(".no free"); wn-make mat(&inNUMSAMPLES,3); wn-make-mat(&outNUMSAMPLES,3); read-dataO; fit-wn-model(out inNUMSAMPLESWNIHUGE); 4~ S .4 4

S~JJ

.,r C C' 23- Appendix 2 physics fitting code #include <mth.h> #include malloc.h #include "wnlib li" #include "wnabs.i #include "wnmax .h" #define NUMPARAMS 2 *4 0 *400 0* 0 local int iteration-Count; static double **save_1n; static double **save-out; static int save-num samples; local void compute extra-ink(p~extra-ink,froin) 0 04 00*0 0 9 000000 0 30 double *pextra -ink; double from[3]; f double c,m,y,k,ink,max; c .O-from[0]; mn .O-from[l]; y 1.0-from[2]; 24max =wn-max(c,m); max =wn-max(max,y); if(max 0. 0) *pextra-ink 1. 0; return; k =wn min(c,r); k =wn min(k,y); ,k OC ink k (rn-k) *pextra ink =ink/max; *pextmra-ink =ink; local void physics colormap(to. from, corig,k,vect) double frorn[3] ,to[3],c_orig,k,vect[NUMPARAMS]; Q externdouble expo,logo; 9999double c,n,y,x extra-ink; *compute extra ink(&extra ink,from);, c =c-orig*(1I.0+k*(extra-ink-1.0)); c =c-orig; c vect[0] +vect[1]*extra ink+vect[2]*extra ink*extra ink; 327990 l3BLE3 CFM *1 printf("c =%lf extra-ink extra-inky; n for(i=0;i< y x -1og(1.O-n*y)/c; to[i] printf(" (%lIf,%lIf,%lIf) (%lIf,%lIf, from[1, from[1,from[2], to[1, to[t, to[2]); local void eval-func(out, in,vect) :double out[3,in[3,vect[NUMPARAMS]; 30 cs corig vect[0]; k vect[l]; physics colormap(out, in,c orig,k,vect); I A-TA 327990 ]3BLE3 CFP029SAU #if 0 local void eval-func(out,in ,vect) double out[3], in[33 ,vect[NUMPARAMS]; extern double expo; double c,k,x,y; int i; c vect[O]; x l.0-in[i]; y outri] =l.O-y; #endif *ftft, ft...

ft...

25 local double function(vect) double vectLU; W t.

extern double wnidist2o,sqrto; int i; double sum; double actual-out[3]; sum 0.0; for(i i save nums samples; i) evalfunc(actual-out save in[i] ,vect); 11ri7QQ TRIE r.PO,295ATT 27sum =wn-dist2(actual-out save out[i],3); printf("c =%lg,k %lg\n",vect[0],vect[1]); printfQ'ob %lg,ave error gn" sum, sqrt(sum/(save num samples return(sum); #if 0 local void gradient(grad,vect) double grad[],vect[]; 0* 0 *0a0 a 0~* Ot 0* a 00 imt ij; double actual-out[3],partials[NUMPARAMS]; wn-zero vect(grad for(i=0;i <save-num-samples; 25 eval-func(actual-out,save-in[i]); foro =0;j actual-outoj] 2 .0*(actual outj] -save out[i][j); backwardpropagate(save-nnet,actual-out); printfQ'grad.------\") printf("iteration count %dn" 1 ,iteration-count); 28fflush(stdout); iteration-count;

I

#endif

I

II.

I *4 o* R

V

void fit-wn model (out, in, num samples, num-iterations) double **out; double **in; int num-samples; int num-iterations; extern double wni-normal-distributionQo; double vect[NUMPARAMS]; int i, code; double val-mmn; wngpmake(" no-free") save -in in; save-out out; save-num samples num-samples; vect[0] vect[l] 1 0; vect [0] vect[1] *4 ~4~4 t t iteration-count =0; 1 coni gradient-metliod(&code,&val -min,vect, size, (function), (gradient),num-iterations); qnconjgradient-method(&code,&val -min,vect, size, (function), (gradient) ,num-iterations); 29powell-methiod(&code, &val min, vect ,NUMPARAMS, (functioin),ium-iterationls); printf("code %dn" ,code); printfQ'val-min %lffn",val-mii); wngpfreeo; o Oa.* a S a 5% {ttc(

C

Appendix 3

NAME

conjdir conjugate directions function minimization package

SYNOPSIS

void conj_directioii_method(&code,&val_min, vect,len,pfunction,max_iterations) int code; double valmin; double vect[]; int len; double (*pfunction)(/* vect int max_iterations; void conj_gradient_method(&code,&val_min, vect,len,pfunction,pgradient,maxiterations) int code; double val_min; double vect[]; int len; double (*pfunction)(/* vect void (*pgradient)(/* grad,vect int max_iterations; double barrier(x) double x; double dbarrier(x) 30 double x; double penalty(x) double x; double dpenalty(x) Sdouble x; a

DESCRIPTION

This package is useful for minimizing continuous, differentiable Sfunctions of many variables. The routines assume that the function to minimize is well approximated by a quadratic form with a positive definite Hessian.

For functions which are well approximated this way, convergence rates are reasonably predictable. If all goes well, these routines converge to a local minimum; this local minimum is not necessarily a global minimum or even a good local minimum. If the function is not differentiable, the algorithms often get stuck in very -31sub-optimal solutions, and no warning is given that this has happened.

"conj_direction_method" conducts a conjugate directions search without the need for derivatives of the function to minimize.

Consequently, it is the easiest to use, but it runs very slowly for medium-sized and large problems. "code" indicates the status of the completed search. "val_min" is the objective function value at the final solution. "vect" should be loaded with a solution used to begin the search; upon return it contains the final solution.

"len" is the number of variables in "vect". "pfunction" is a pointer to the function to minimize. "max_iterations" is the maximum number of iterations allowed; set to WNIHUGE for no limit.

"conj_gradientmethod" conducts a conjugate directions search using derivatives of the function to minimize. Using derivatives usually results in a dramatic speedup, but it makes programming much more complicated. Furthermore, if the derivatives are in error, no warning is given; the result is either very slow convergence or convergence to an incorrect result. The arguments mean the same as the arguments to "conj directionmethod".

The additional argument "pgradient" is a pointer to a function i which computes the gradient of the function to minimize; that is S 25 grad[i]= d function d vect[i] Because of the difficulty in ensuring that the derivatives are S: correct, it is usually best to start by experimenting on a small problem using "conjdirection_method", even if you ultimately 30 intend to use "conj gradientmethod". The solutions obtained from "conj_gradientmethod" using derivatives should be similar or the same as solutions obtained from "conj_directionsmethod" not using derivatives. If this is not so, you.r derivatives are 35 probably faulty.

"barrier", "dbarrier", "penalty", and "dpenalty" are useful for implementing barrier and penalty methods of non-linear programming.

barrier(x) returns 0 for x 0 and infinity for x 1. Between 0 and 1, it moves up from 0 to infinity. It is continuous and differentible for allx 1.

dbarrier(x) returns the derivative of barrier(x).

penalty(x) returns 0 for x 0 and x^2 for x 0. It is continuous and differentiable for all x.

dpenalty(x) returns the derivative of penalty(x).

-32-

RESOURCES

conjugatedirectionmethod runs with AVERAGE CASE: time len^2 (len^2 <time to evaluate function>) stack memory 1 dynamic memory len*len conjugate_gradient_method runs with AVERAGE CASE: time len <time to evaluate function> time to evaluate gradien>) stack memory 1 dynamic memory len Run time depends heavily on the problem being solved. If the function is badly behaved, convergence can take much longer than these times.

If the function is not too badly behaved and the Hessian has distinct eigenvalues, the times given here usually apply. If the eigenvalues of the Hessian are not distinct, much faster convergence times are possible.

S. DIAGNOSTICS ai code WN SUCCESS means successful completion, optimum found t, code WN_SUBOPTIMAL means termination before optimum found 30 code WN_UNBOUNDED means function is unbounded

BUGS

SEE ALSO

AUTHOR

Will Naylor ci

C*

327990 BBLE3 CFP0295AU 327990 I3BLEI-AUno:\cisr\hhie;hhb 4 ,TTh -33- 1 void force optimize stopo void minimize ld(pvalmin,px_min,xstart,function) void conj directionmethod(pcode,pvalmin,vectlenpfunctionmax iterations) #include <stdio.h> #include <math.h> #include "wnlib.h" #include "wnswap.h" #include "wnmat.h"

S

o a.

9*

S..

*Sa.

5 C *r #define TOLERANCE #define ROOTTOLERANCE #define TOLERANCE #define ROOTTOLERANCE 25 #define TOLERANCE #define ROOTTOLERANCE #define TOLERANCE #define ROOTTOLERANCE #define TOLERANCE 30 #define ROOTTOLERANCE (1.0e-13) (1.0e-7) (1.Oe-i') (1.0e-7) (1.0e-16) (1.0e-8) (1.0e-20) (1.Oe-10) (1.0e-15) (3.3e-8) #define GOLDENRATIO 1.618034 #define GOLDENSECTION 0.3819660 #define MAXMAGNIFICATION 100.0 #define MAXBRENTITERATIONS 200 #define SHIFT3Qxl,x2,-x3) Cx2); x2)= Qx3);} #define SHIFT4Cx x2_-x3,_x4) {SHIFT3Cx1,_x2,_x3);

I

327990 BBE3 Cr-.nn-?-AU [o:\cisra\bble.dir~bble3]speci.doc:IAD 34bool force optimize stop flag FALSE; void force optimize stop() force optimize stop flag TRUE; local double max(x,y) double x,y; if(x =y) return(x); else return(y); 4.p**t .4 4 e* S a. ta 44 *4 4 .4 4.544.

local double min(x,y) 30 double x,y; if(x y) return(x); else return(y); local double inherit-sign~a,b) double a,b; if(b 0. 0) return(fabs(a)); else return(-fabs(a)); /assume *pfa,*pfb,*pfc already computed local void bracket-minimum(pax,pbx,pcx,pfa,pfb ,pfc,pfunction) double *pa,*pbx,*pcx, *pfa,*pfb, *pfc; double (*pfunction) (/*double 4t 4 4I4* 4.4, 4040 Ode 44 4* 4

S*

o 44 double ulim,u,r,q,fu; 30 if(*pfc *pfa) {nsa(pa *fdul) wn-swap (*pfa, *pfcx,double); while((*ptb *pfc)&&not(force op~timize stop -flag)) r -(*pbx-*pax)*(*pfb-*pfc); q (*pbx-*pcx)*(*pfb-*pfa); 327990 DDLE3 CF110295AU 32790 BBE3 FPO2SAU[o.\k isra\bble.dir\bbks3]speci.doc:JAD 36u (*pbx) *pax)*r)/ O*inlierit sign(max(fabs(q-r) ,TOLERANCE) ulim (*pbx) +MAXMAGNJF1CATION*(*pcx-*pbx); jf((*pbx-u)*(u-*pcx) 0.0) fu (*pfunctjon)(u); if(fu *pfc) *Paxj *pbx; *pbx u; *pfa *pfb; *pfb fu; return; else if(fu *pfb) *p fc =fu; return; u (*pcx)±GOLDBNRATIO*(*pcx-*pbx); fu =(*Pfunction)(U); 25 else if((*pcx u)*(u-ulim) 0.0) fu (*pfunction)(u); if(fu *pfc) SHIFT4(*pbx, *pcx,u, *pcx+ GOLDENRATIO*(*pcx *pbx)); SHIFT4(*pfb, *pfc,fj, (*pfunctjon) else if((u-ulim)*(ullin-*pcx) 0.0) u ulim; fu ~-(*pfunction)(u); else L~4~ ftdtO ~I 4* 5

S

ft 327990 BBLE3 CFP029SAU 3279O BBB3 FPO2SAU[o:\cisra\bble.dir\bble3]speci.doc:IAD 37u =(*pcx GOLDEN RA fIO*(*pcx-*pbx); fu =(*pfunctjon)(u); SHIFT4(*pax, *pbx, *pcx,u); SHIFT4(*pfa, *pffi, *pfc, fu); 4 1* assume *pfa,*pfh,*pfc already computed local void brent-method(pax,pbx,pcx,pfa, pfh,pf c,pfuniction) 4 double *pax,*pbx,*pcx, *pfa,*pfb, *pfc; double (*pfulnctio-i) (/*double int iterations; double d, fu, fv,fw,p, q,r, toll, tol2,u,v,w,xm, e,etemp; iff orce-op timize-stop flag) force op~timize stop flag =FALSE; return; jf(*pajc *PCX) 30 *wn-swap (*pax, *pc.-.double); wn-swap(*pfa, *pfc, double); wn assert(*pax *pbx); wn-assert(*pbx *pc); wn-assert(*pfb *pfa); 38e e=d WNFHUGE; V =*pax; fv *p fa; w =*PCX; fiv *pfc; for(iterations (iterations< <i,/AX BRENTITERAT'IONS)&&inot(force__optimize-sto p flag); d- iterations) xm =0.5*(*pax "pcx), toll =ROOTTOLERANCE*fabs(*pbx)+ROOTTOLERANCE; tol2 if(fabs(*pbx xm) (tol2.-0.5*(*pcx *pax))) return; if(fabs(e) toll) r -(*pbx w)*(*pfb f) q (*pbx v)*(*pfb fw); p (*pbx v)*q (*pbx- )r 25 q if(q 0. 0) p p C: 30 q fabs(q); etemp =e: e d .:if( (fabs(p) =fabs(0 11 (p =q*(*pax *pbx)) (p q*(*pcx *pbx)) -111-1 39e (*pbx =xm)?(*pax *pbx):(*pcx *pbx); d =GOLDENSECTJON*e; else d p/q; u =*pbx d; if(((u *pax) tol2)1 I ((*pcx u) tol2)) d inherit sign(toll1,xm *pbx); else e =(*pbx (*pax *pbx): (*pcx *pbx); d GOLDENSECTION*e; 69 ft ft...

ft.

546 ft 6* ft **CBftft u (fabs(d) =toll)?(*pbx d):(*pbx inherit -Sign(toll fu =(*pfunction)(u); iffu pfb) iffu *pbx) *pcx *pbx; SHJFT4(v,w, *pbx,u); SHIFT4(fv,fw, *pfb,fij); else if(u *pbx) *pax u; else *PX= U; if(fu fw)~ (w HF3v~~) SHIFT3(v,fw,u); else if(fu I I(v IfI(v w)) v U; fv fu; force_optimize stop flag FALSE; fprintf(stderr, "brent-method: too many iterations. ~tess.

LYA

U~

1~ It..

assume *pfa, *pfb ,*pfc already computed local void minimize 1 d raw(pax,pbx,pcx,pfa,pfb ,pfc,pfunctioni) double *pax,*pbx, *pcx,*pfa, *pfh,*pfc; double (*pfunction) (/*double x*I); bracket minimum~pax,pbx,pcx,pfa,pfb ,pfc,pfunction brent_method(pax,pbx,pcx,pfa,pfb ,pfc,pfunction); void minimize 1 d(pvalmin,px _min,xstrt,pfunction) double *pvalmin,*px,_min,x start; double (*pfunction) (/*double x*I); I- 41double ax,bx,cx,fe,fb,fc; force optimize stop flag FALSE; ax x-start; bx cx x start+ fa =(*pfunction)(ax); fb (*pfunction)(bx); fc (*pfunction)(cx); minimize 1 d raw(&ax, &bx ,&cx &fb ,&fc ,pfunction); *pval-ii =f *pxmin bx; 4 4 6 of*.

04#6 4 4 6 46.6 t m6 local void initializepowell directions (xi ,len) double **xi; mnt len; int i; wn-zero mat(xi,len,len); for(i=O;i <len; xi[iffi]

I

local mnt save-len; local double *buffer-vect, *save-vect, *save-direction; O'rflnn nOTI 01 nOAIOC ATT local double (*savepfunction) (/*double vecto]*/); local double powell-line-function(x) t C double x; wn-copyvect(save-vect, buffer-vect,save-len); wn-add-scaled-vect(save vect save-direction,x,save-leni); return((*savepfunction) (save-vect));

I

local void powell-line-minimize(vect, directioin,len,pval-min,pfunctioni) double vect[], directionfl; mnt len; double *pval-mim; double (*pfunction)(/*double vectfl*I); double ax,bx,cx,fa,fb,fc; wn-copyvect(buffer-vect,vect,lein); save-vect =vect; save len len; save direction =direction; savepfunction =pfunction; ax 0.0; bx cx 1.0; fa *pval-min; fb =powell-line function(bx); fc powell-line-function(cx);

I

V

minimize_1 d-raw(&ax, &bx, &cx, &fa, &fb, &fc, (powell-line function)); -43printf("l %lf\n",bx); wn -copyvect(vect, buffer -vect,len); wn-add-scaled-vect~vect, direction, bx, len); *pval min fb; local double square(x) double x;

I

return(x*x); void conj direction method~pcode,pva_min vect len,pfunction, max-iterationF) hit *pcode; double *pval-mmn; double vecto; int len; double (*pfunctjon) (/*double vectoU*/); ilt max-iterations; int i, ibig,j ,iteration; double t,fptt,f ,deI; double *pt, *ptt, *xit; double **xj; force_optimize stop flag FALSE; wn ,gpmake('no-free"); 327990 BBLE3 CF.P0295KU 32790 BBE3 FPO25AUlo:.\cisra~bble.dir\bble3]speci.doc:IAD 44wn-.make vect(&buffer-vect,len); wn -make mat(&xi,len, len); wni-make vect(&pt,len); wn-make vect(&ptt, len); wn-make vect(&xit,len); initial izepowell directions (xi, len); *pval min =(*pfunctjon)(yect); wn copy vect(pt,vect,len); iteration =0; for(iteration=O;; iteration) fp =*pval mill; del 0.0 for(i=O;i <len; 20 fptt *p-val min; powell-linle-minimize~vect,xi[i], len,pval rn *In,pfunctioni); if(fabs(fptt-(*pva1 minl)) del) del fabs(fptt-(*pval mini)); *~.ibig i :if(2 .0*fabs(fp-(*pval mmn)) TOLERANCE*(fabs(fp) +fabs(*pval-min))) 30 -*pcode WNSUCCESS; force,_optimize stopjfag FALSE; w_gpfreO return; if( force optimize-stop flag ((iteration max iterations) &&(max,_iterations WNIHUGE)) *pcode WN_-SUBOPTIMAL; force _optimize-stop flag FALSE; wn-gpfreeo; return;

I

foroj= Oj len; pttlj] 2 .0*vecto]-ptUj]; xitlj] vectlj]-ptojl; wn-copyvect(pt,vect, len); fptt (*pfunction)(ptt); if(fptt fp) t= 2. O*(fp-2 .0*(*pval min) +fptt) *square(fp-(*pval rnini-del) de1*square(fp-fptt); if(t 0. 0) powell-line-minimize(vect, xit, len,pval -Mii,pfu: tion); wn-Copyvect(xi[ibig] ,xit,len); Lv.

4 4 -4,1199n RRUThT3 CFP0295AIT 46void conj gradient-rnethod(pcode,pval-min,vect, Ieni,pfunction,pgradient,max-iterations) #include <sti.h> #include "wnlib.h" #include "wnmat.h" #include wnswap.h" #include "wnabs .h" Pt mE.

S.

#define GOLDENRATIO 1.618034 #define GOLDENSECTION 0.38 19660 #define MAXSEARCHITERATIONS 100 extern bool force_optimize stop flag; local double *buffervyect; Fit parabola to 3 points (x Parabola is of the form y a*(x-x0)A2 b Return a, xO, b.

local void fitparabola 3p(pa,px0,pb,xl ,x2,x3,yl,y2,y3) p A 327990 BBLE3 -,!:IP029SAU 32790 BLE CFO29AU[o:\cisra\bble.dir\bble31sneci.dot,-TADl 47double *pa, *pxO, *pb; double xi,x2,x3, yl,y2,y3; double x12,x23,dyl2,dy23,ddy,diff; I* sort the X's *I if(not(xl x3)) wn-swap(xl 1,x3, double); wn-swap (y1, y3, double); if(x2 xl1) wn-swap(xl 1,x2, double); wn swap (y l,y2, double); else if(x3 x2) wn-swap (x2, x3, double); wn-swap(y2,y3, double); w I er~l x) wn,_assert(xl x3); x12 0.5 *(xl+x2); x23 0.5*(x2+x3);

S

CS'.

dy12 (y2-yl)I(x2-xl); dy23 (y3-y2)I(x'3-x2); ddy =dy23-dyl2; *pa =0.5*ddyI(x23-x12); if(ddy 0) *pxO 0.0; 48- *pb =y 2

I

else f *pxQ (dy23*x12-dyt2*x23)/ddy; diff x2-(*pxO); *pb y2-(*pa)*diff*diff; Fit parabola to 2 points and a derivative at one point.

Parabola is of the form y a(x-xO)^2 b Return a, AO, b.

20 local void fitparabola 2pd~pa,px0,pb,x1 ,x2,y 1,dy l,y2) double *pa, *pxO, *pb; double xl,x2, yl,dyl,y2; double x 12,dyl12, ddy, diff; 30 wnassert(xl x2); .t* 0 4 SC

(CCC

t S4 *054*5 x12 .5*(xl+x2); dy12 (y2-yl)/(x2-xl); ddy dy12-dyI; *pa 0.5*ddy/(x12.x1); if(ddy 0. 0) III-nr DflT72 r(IOCAT A*-%ihi-,2i A--.TAn 49- *pxO 0. 0; *pb =yl1; else *px0 (dy12*lx1-dy1*xl2)/ddy; diff =x1-(*pxO); *pb =y1-(*pa)*diff*diff; local bool vects-equal(vl,v2,len) double vl[],v2fl; int len; 20 -for(i=0;i<len;++i) if(vl[i] v2[i]) return(FALSE); return(TRUE); local -vfid eval-fuinction~pval,vect,direction,x,1en ,PLfunction) double *pval; double vectfl,directionol; double x; int len; double (*pfunctioil) (/*double wn-add. vect-and-scaled-vect(buffer vect ,vect,direction,x,len); *pval (*pfunction) (buffer vect); local void line-Minimize(pval,pjuinpjlen, vect, direction, last-val ,last_g ,lastjumplei,len,pfuctioni) double *pval, *pjumplen; double vect[]; double direction[]; double last-yal; double last 20 double lastjumpjlen; mnt 1ren; dcuble (*pfunction) (/*ouble extern double wn-dot-vects wn-random-double-betweenO; double a,xO,yO,b,dyl,y2,x3,y,last x3,lasty3 t iteration count; a:dyl wn dot vects (direction, last gjen); 30 lastjump-len wn-sign(dyl)*wnabs(lastjumpjen); .:eval_functioni(&y2,vect,direction,lastjump-len,len,pfunction); fitparabola_2pd(&a, &xC,&b 0,lastjump len last-val dy l,y2); if(a 0. 0) eval-function(&yO,vect, direction,x0,len,pfunction); 327990 RULE3 CFP0295AU -51if((yO last-val)&&(yO y2)) 1* parabolic fit succeeded *pvaI yO; *pjumplen =xA; wnviadd-scaled vect(vect,clirection,xO ,len); return; parabolic fit failed x3 lastjump_len; y 2 if(y2 lastval) 1* project inwar, 20 iteration-count 0; do x3 GLDN-ETIN x3 (=GOLDENSECTIONw admdul btenO911) 4. 4 _-function(&y3 ,vect, direction,x3 ,len,pfunction-); Cal, wn-add-vect-and-scaled vect(buffer-vect vect,direction,x3 ,len); C, ~if(vects-equalfkbuffer-vect,vect,len)) break; y3 (*pfunctjon) (buffer-vect); 52iteration-Count; while((y3 =last (iteration count MAXSEARCHITERATIONS)); if(y3 last-vat) *pval =y3 *pjumplen x3; wn-add.-scaled vect~vect, direction, x3, len); return; if(y2 last vat) *pvaO last-val; *pjulplen z,0.0; a. 4 o at ~a Ga., at return; projeui outward x3 lastjur-npjen; y3 =y2 iteration-count 0; do f last-x3 lasty3 x3; y 3 x3 *1'

GOLDENRATIO;

x3 (GOLDENRATIO*wn random-double-between(O.9, eval-function(&y3 ,vect,direction,x3 ,len,pfunctLion); t

I

V.

IJ

53iteration-count; while((y3 last y3)&&(y3 last val)&& (iteration-count MAXSEARCHITERATIONS)); *pval last y3; *pjuirnpleni last -x3; wn-add-scaled-vect(vect, direction, last-x3, len);

CC

t tC t t r tr V return; void conj gradient method(pcode,pval-Min vect, len,pfunction,pgradient,max-iterationis) mnt *pcode; double *pval-min; double vectfl; int len; double (*pfunction) (/*double vect[1 1) void (*pgradient) (/*double grad[] double vect[] ilt max-iterations; extern double wn-normal-distributionO; extern double wn -norm2O,wn normo,wn-dot-vectso; intj, ,iteration, no -move-count; double norm2_g ,norm2._lastg,val,,last-val numerator,beta, jump__len, lastjumplen; double *last g, *direction; wngpmake("no free"); force-optimize-stop flag FALSE; wn-make-vect(&buffer-vect,len); wn-make vect(&g ,1en); g =buffer-vect; share memory with buffer-vect wn-make vect(&last g, len); wn make vect(&direction, len); jump_len no move-count 0; val (*pfunction) (vect); (*pgradient)(g,vect); wn-copy_vei(direction, g, len); f'or (iteration=0; iteration< <max iterations; +iteration) norm2g wn-norm2(g,len); if(norm2_g 0) unlikely wngpfreeo; *pval min val; *pcode =WNSUCCESS: return; lastjumplen =jump len; last-val val; wn -copyvect~lastg, g,1len); 30 norm2lastg norm2g9; C line minimize (&val,&jump_len, vect, directionjlast val last g ,astjumplen,len,pfunction); ifforce optimize-stop_flag) Iforce _optimnize stop flag =FALSE; wngpfreeo; -9.

I

I

r *pval-minf= vat; *pcode WNSUBOPTIMAL; return; wn-assert(val last val); if(val ==last val) if(no move-count =2) wn-gpfreeo; *pval-min vat; *pcode WNSUCCESS; return; else +no-move-count; jump len lastjumplen; tie .5.

e 4,5*

S

S

4554 a 4**S C C C 4~ 20 else no-move-count 0; gradient) (g,vect); numerator 0.0; foroj Oj len; +j) numerator gj*gj; numerator g~j]*(g~j]4astg9j]); beta numerator/norm2last-g; -1 11 56wn-copy vect(direction,g ,len); wn-add-vect-and-scaled-vect(direction,g ,direction,beta,leni); i t v t t It prinitfQ'ob %lg,beta %lg,numerator %g,denom =%lg,norm2(direction) val, beta, numerator, norrn2_last g, wn-orn2 (direction)); printf("iteration %d,ob ,iteration,val); wngpfreeo; *pval -mini val; *pcode WN-SUBOPTIMAL; 327990 BBLE3 CFP0295AU 32790 BBE3 FPO25AU[o:\cisra\bble.dir\bblc3]speci.doc:IAD double barrier(x) double dbarrier(x) double penalty(x) double dpenalty(x) #include <stdio.h> #include <mathi.h> #include "wnlib .h" double barrier(x) double x; C t if(x 0.0) return(0.0); else if(x return(WN-FHJGE); else double dbarrier(x) 37,7990 BBLE3 CFP029SAU 32790 BB133 o:\cisra\bble.dir\bblc3lshneci.doc:IAD 58double x; iffx 0.0) return(0O); else if(x t0 return(WNFHUGE); else

LI

double diff; diff return(2 .0*x/(diff*diff)); double penalty(x) double x; if(x 0.0) S 30 return(0. 0); else return(x*x); double dpenalty(x) double x; if(x 0.0) return(0.0); else return(2. 0*x); t', 327990 BBLE3 CFP0295AU 32790BBE3 FP29AU[o:\cis ra\bbl e. dir\bbl e3 Isp eci. doc:IAD

Claims (9)

1. A method for displaying colour images, intended to be displayed on a first colour display device, on a second colour display device having a colour gamut different from the first colour display device, said method including the steps of: determining a plurality of input extrema colour values, within the colour gamut of said first colour display device and at the extrema of the gamut of said first colour display; determining a plurality of output extrema colour values corresponding to said input extrema colour values, within the colour gamut of said second colour display and at the extrema of the gamut of said second colour display; and determining a mapping for substantially all the values within the colour gamut of the first colour display to substantially all the values within the gamut of the second colour display, said mapping substantially preserving the relative value of differences between colour values relative to the extent of the gamut of said first colour space and differences between corresponding colour values in the gamut of said second colour space relative to the extent of the gamut of said second colour space.

2. A method for displaying colour images as claimed in claim 1 wherein said mapping includes determining a colour correction function mapping colour gamut values within the gamut of said second colour display to colour gamut values within the j colour gamut of said first colour display device. tt said colour correction function includes a number of free parameters, the values of which are determined by measurement of output colour values of said second colour display device.

4. A method for displaying colour images as claimed in claim 3 wherein A method for displaying colour images as claimed in claim 3 wherein Ssaid measurements are made in an alternative colour space and converted to the colour space of said colour gamut of said second colour display device.

6. A method for displaying colour images as claimed in claim 5 wherein said conversion utilising a predetermined number of said output extrema values to determine said conversion. ''7nn nlT r VI" t Mrm-"nmoATT ro:\cisra\bble.dir\bble3lsveci.doc:LDP -61-

7. A method for displaying colour images as claimed in claim 6 wherein said conversion includes utilising a barrier function to ensure that said converted measurements are within the colour gamut of said second colour display device.

8. A method for displaying colour images as claimed in any one of claims 2 to 7 wherein said colour correction function is substantially of the form: -In (1-(l-e G C) I= G where C is a colour value within the colour gamut of said second colour display, G is a predetermined value or a predetermined fuction, and I is the corresponding colour within the gamut of said first colour display device.

9. A method for displaying colour images as claimed in claim 8 wherein said colour correction function is applied separately for each colour component of said second colour display device. A method of displaying colour images as claimed in claim 9 wherein said second colour display device is an ink printing device and G is a function of the amount of ink currently printed by said second colour display device.

11. A method of displaying colour images as claimed in claim 10 wherein G is of the form: G(ink) go glx ink where go and gl are constants determined from measurement of output colour values of said second colour display and "ink" is an approximation of the amount of ink printed at a present pixel position of said image.

12. A method for displaying colour images substantially as described herein with reference to Figs 3 to 5 of the drawings. 6" q 66 t ~c~ DATED this TWELFTH day of FEBRUARY 19)6 Canon Information Systems Research Australia Pty Ltd Patent Attorneys for the Applicant SPRUSON FERGUSON ABSTRACT COLOUR CONVERSION METHOD A method for displaying colour images, intended to be displayed on a first colour display device, on a second colour display device having a colour gamut different from a colour gamut of the first colour display device is disclosed. The method includes the a first step of determining a plurality of input extrema colour values (RGB), within the colour gamut of the first colour display device and at the extrema of the gamut of "aid first colour display. Next, the method determines a plurality of output extrema colour values corresponding to the input extrema colour values (RGB), wit'hin the colour gamut of the second colour display and at the extrema of the ,o gamut o( said second colour display. Then, the method determines a mapping for 15 substantially all the values within the colour gamut of the first colour display to substantially all the values within the gamut of the second colour display, the mapping substantially preserving the relative value of differences between colour values relative to the extent of the gamut of the first colour space and differences between corresponding colour values in the gamut of the second colour space relative to the extent of the gamut of the second colour space. 14