function [MunsellSpec MunsellVec Status] = xyYtoMunsell(x, y, Y); % Purpose Convert xyY coordinates to a Munsell specification, by interpolating % over Munsell renotation data. Near the MacAdam limits, extrapolated % Munsell renotation data will be used if necessary. % % Description The Munsell system specifies a local colour by giving its hue (H), value (V), % and chroma(C) in the form HV/C. The value is a number between 0 and 10. % The chroma is a positive number, whose bound depends on hue and value, % as given by the MacAdam limits. The hue specification consists of a letter % designator (B, BG, G, GY, Y, YR, R, RP, P, PB), and a number designator % which is greater than 0, and less than or equal to 10. If chroma is % 0, then the local colour has no hue, and is specified as NV, where N is the % string "N," and V is the value. For example, 5.0R9.0/4.0 is a light pastel % red, while N3 is a dark grey. % % The Commission Internationale de l'Eclairage (CIE) uses a different system % for specifying colours. In their system, standardized in 1931, % a coloured light source is specified by three coordinates: x, y, and Y. % The coordinate Y gives a colour's luminance factor, or relative luminance. % When expressed as a number between 0 and % 100, Y is the intensity of the source as a percentage of some % maximum intensity. This percentage is calculated with regard to the % human photopic luminous efficiency function, which has been established % as part of the CIE 2 degree standard observer. When dealing with % physical samples such as paints, Y is the percentage of a fixed light % source that a paint sample reflects (with regard to the standard % observer). % % The coordinates x and y are chromaticity coordinates. While Y % indicates how light or dark a colour is, x and y % indicate hue and chroma. For example, the colour might be a saturated red, % or a dull green. The CIE chromaticity diagram ([Foley1990, Sect. 13.2.2]) displays % all chromaticities that are possible for physical samples of a fixed % luminance factor. % % There is a conceptual difficulty in converting between Munsell specifications % and CIE coordinates. The difficulty occurs because the Munsell system % applies to local colours, which are defined by spectral reflectance % functions, and CIE coordinates apply to coloured lights, which are defined % by spectral power distributions (SPD). The 1943 Munsell renotation ([Newhall1943]) % bridged this gap by fixing Illuminant C as a light source, % that illuminates the local colour. The local colour is specified in the % Munsell system, and the light reflecting off the local colour is analyzed % in terms of CIE coordinates. The Y value is the percentage of light, % originally of an Illuminant C SPD, that the local colour reflects, % measured in terms of the CIE standard observer. The x and y chromaticity % coordinates can also be found for the reflected light. % % The conversion data for the Munsell renotation was originally given as % a look-up table, expressing a Munsell specification in xyY coordinates % ([Newhall1943, Table I]). The renotation was later updated % in [ASTMD1535-08]. The routine MunsellToxyY % can be used to find xyY coordinates for Munsell specifications that are not % in the original renotation data. % % The interpolation uses two empirical facts. First, the relative luminance % Y is a function solely of the Munsell value, and vice versa. Therefore, % a Munsell section through a fixed value is mapped bijectively to the % chromaticity diagram for the corresponding relative % luminance ([Newhall1943, Figs. 1 to 9]; % reproduced with some modifications in [Agoston1987, Figs. 12.1 to 12.9]). % Second, within a Munsell section of fixed value, the lines of constant % chroma are rings around a central, achromatic point, corresponding to the % Munsell grey of that value. The straight lines radiating from the central % point are roughly lines of constant hue. It is therefore natural to use a % polar coordinate system, (r, theta), about the achromatic point, where r, % the distance from grey, is the chroma, and each angle of theta corresponds % to a particular hue. % % This current routine is an inverse to MunsellToxyY. Given xyY coordinates, we aim to % find the Munsell specification that would generate those xyY coordinates. % We can evaluate the xyY coordinates for any Munsell specification, and must % use that information in the inverse calculation. % % This inversion routine uses the same two empirical facts that the original % function uses. First of all, the Munsell value can be calculated directly % from Y. Secondly, for a fixed value, the level curves of Munsell % hue and chroma are approximately rings and radials about the achromatic % point in xy coordinates. % % The first step in the inversion is therefore a simple Munsell value % calculation, using a table lookup of inputs and outputs of the quintic % polynomial in [ASTMD1535-08, Eq. 2]. The problem then % reduces to a two-dimensional problem: % finding a hue and chroma that give the desired x and y. % % The inversion proceeds iteratively. At each step, a Munsell % specification is selected, whose xyY coordinates are (usually) closer to the % input xyY than the previous step. Closeness is measured by the Euclidean % distance in xy space. When the xyY of the current Munsell specification is % less than some threshold distance, the inversion terminates. % % The starting Munsell specification is chosen with the help of the CIELAB % model. The input xyY are converted into CIELAB coordinates, resulting in % a hue angle and an approximation for chroma. CIELAB coordinates % correspond approximately to Munsell quantities, so they are converted to % a Munsell specification, where the iterative inversion starts. Although % this conversion is not exact, it is close enough to make a useful starting % point. % % At each subsequent step, the xy coordinates are found for the current % Munsell specification. They are converted to (rCur, thetaCur) coordinates about % the achromatic point of value Y, and compared to (rInput, thetaInput), % which correspond to the input x and y. The distance r corresponds to % chroma, and the angle theta corresponds to hue. The goal of the % inversion is to have (rCur, thetaCur) match (rInput, thetaInput) as closely % as possible. % % To achieve this goal, the current chroma and hue will be adjusted. Start by % varying hue while keeping chroma constant. In the CIE chromaticity diagram, % the Munsell hues lie along approximate % radials, whose angles can be read off. The routines MunsellHueToChromDiagHueAngle % and ChromDiagHueAngleToMunsellHue convert between Munsell hues, and the % angles of the radials. The angles thetaCur and thetaInput will probably % differ. Calculate new hues whose angles (theta1, theta2, .., thetan) % bound thetaInput. Then interpolate linearly to find a hue that should % match thetaInput more closely. % % A natural approach to generating theta1 is to choose % (new hue angle) = (current hue angle) + thetaInput - thetaCur, % which roughly corrects for the current angle discrepancy. Further angles % can be generated by % (current hue angle) + 2 * (thetaInput - thetaCur), % (current hue angle) + 3 * (thetaInput - thetaCur), etc. % In practice, two new hue angles are usually sufficient. % % After varying hue without varying chroma, let chroma vary without changing hue. % Each new chroma produces a Munsell specification whose own (r, theta) values % can be calculated. If n chromas are tried, there will be n Munsell % specifications, with n r-values: [r1, r2, ..., rn]. If the chromas are % chosen appropriately, rInput will be bounded by some ri and rj. By linear % interpolation, we find a new chroma value that should be closer to producing % rInput. % % A natural approach to generating r1 is to choose % (new chroma) = (rInput/rCur) * (current chroma), (1) % because this corrects for the current undershoot or overshoot. If rInput % is not between rCur and the r1 generated by (1), then generate r2 from % another chroma given by ((rInput/rCur)^2) * (current chroma). If rInput does % not lie in the interval containing rCur, r1, and r2, then increase the % exponent to 3 to produce an r3, and so on. In practice, it is rarely % necessary to go beyond r1 and r2. % % Proceed iteratively, adjusting the hue and the chroma in alternate % iterations, until the desired accuracy is achieved. % % For a detailed discussion and description of the inversion algorithm, % see [Centore2011]. % % [Agoston1987] G. A. Agoston, Color Theory and Its Application in Art % and Design, 2nd ed., Springer Series in Optical Science, vol. 19, % Springer-Verlag, 1987. % [ASTMD1535-08] ASTM, Standard D 1535-08, "Standard Practice for Specifying Color by the % Munsell System," approved January 1, 2008. % [Centore2011] Paul Centore, "An Open-Source Inversion for the Munsell % Renotation," 2011, unpublished (currently available at centore@99main.com/~centore). % [Foley1990] James D. Foley, Andries van Dam, Steven K. Feiner, & John % F. Hughes, Computer Graphics: Principles and Practice, 2nd ed., % Addison-Wesley Publishing Company, 1990. % [McCamy1992] C. S. McCamy, "Munsell Value as Explicit Functions of CIE % Luminance Factor," COLOR Research and Application, Vol. 17, % 1992, pp. 205-207. % [Newhall1943] S. M. Newhall, D. Nickerson, & D. B. Judd, "Final % Report of the O.S.A. Subcommittee on the Spacing of the Munsell % Colors," Journal of the Optical Society of America, Vol. 33, % Issue 7, pp. 385-418, 1943. % % % Syntax [MunsellSpec MunsellVec Status] = xyYtoMunsell(x, y, Y); % % MunsellSpec A Munsell specification string, in Munsell format, such as 2.3R5.6/10.8. % % MunsellVec A Munsell specification, in ColorLab format, such as [2.3, 5.6, 10.8, 7] . % % [x y Y] CIE coordinates of MunsellSpecString, when % illuminated by Illuminant C. Y is the luminance factor % % Status A return code with two fields. The second field is a list % of possible return messages. The first field is a positive % integer indicating one of the return messages. Two additional fields, % Status.dist and Status.num, might be assigned for debugging. % % Related MunsellToxyY % Functions % % Required MunsellSpecToColorLabFormat, LuminanceFactorToMunsellValue, % Functions MunsellToxyForIntegerMunsellValue, MunsellToxyY, % CIELABtoApproxMunsellSpec, MunsellHueToChromDiagHueAngle, % xyz2lab, lab2perc % % Author Paul Centore (June 13, 2010) % Revision Paul Centore (Jan. 12, 2011) % ---Previouly, Munsell value was calculated in accordance with [Newhall1943]. In the revision, % Munsell value is calculated in accordance with [McCamy1992], which has been % incorporated into [ASTMD1535-08]. % ---Tne Munsell output is now returned as a standard Munsell string, as well as in ColorLab format % Revision Paul Centore (Jan. 22, 2011) % ---Previouly, Munsell value was calculated in accordance with [McCamy1992]. In the revision, % Munsell value is calculated by linear interpolation over a table lookup of evaluations of % [ASTMD1535-08, Eq. 2]. % Revision Paul Centore (May 8, 2012) % ---Changed != to ~= so that code would work in both Matlab and Octave. % Revision Zsolt Kovacs-Vajna (May 14, 2012) % ---Added 'linear' option to interp1 calls when 'extrap' option is used, for compatibility % with Matlab. % Revision Paul Centore (December 26, 2012) % ---Moved from MunsellConversions program to MunsellToolbox. % Revision Paul Centore (August 31, 2013) % ---Moved from MunsellToolbox program to MunsellAndKubelkaMunkToolbox. % Revision Paul Centore (Jan. 1, 2014) % ---Replaced call to IlluminantCWhitePoint with call to roo2xy (from OptProp). % Revision Paul Centore (Feb. 16, 2014) % ---Replaced call to roo2xy (from OptProp) with call to ChromaticityOfWhitePoint % % Copyright 2010-2014 Paul Centore % % This file is part of MunsellAndKubelkaMunkToolbox. % % MunsellAndKubelkaMunkToolbox is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % MunsellAndKubelkaMunkToolbox is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with MunsellAndKubelkaMunkToolbox. If not, see . Status.Messages = {'Success',... [num2str(x),', ',num2str(y),', ',num2str(Y),' beyond gamut of renotation data.'],... 'Exceeded maximum number of iterations.',... [num2str(x),', ',num2str(y),', ',num2str(Y),' outside MacAdam limits.'],... }; % Assign default output values MunsellVec = -99; MunsellSpec = 'ERROR'; Status.ind = -99; % If the colour stimulus is beyond the MacAdam limits, set an error code and return. % Determine the MacAdam limits with regard to illuminant C, which is the Munsell standard. if ~IsWithinMacAdamLimits(x, y, Y, 'C'); Status.ind = 4; return; end TRUE = 1; FALSE = 0; % Set ShowFigure to TRUE if a figure is desired for debugging ShowFigure = false ; if ShowFigure == TRUE figure ; end % The Munsell value is calculated directly from the luminance factor, % using a lookup table of evaluations of [ASTMD1535-08, Eq. 2]. MunsellValues = LuminanceFactorToMunsellValue(Y) ; MunsellValue = MunsellValues.ASTMTableLookup ; % To avoid numerical issues with approximation, set MunsellValue to an % integer if it is very close to an integer. if abs(MunsellValue - round(MunsellValue)) < 0.001 MunsellValue = round(MunsellValue); end % Express xy chromaticity coordinates as polar coordinates around the achromatic % point with Munsell value corresponding to Y. [xcenter ycenter Ytemp StatusCode] = MunsellToxyY([MunsellValue]) ; if StatusCode.ind ~= 1 % disp(['Status 2, line 271']) Status.ind = 2; Status.dist = -99; Status.num = -99; return end % rInput and thetaInput are the values the inversion algorithm will attempt to match [thetaInput rInput] = cart2pol(x-xcenter, y-ycenter); thetaInput = mod((180/pi)*thetaInput, 360) ; % Express in degrees thetaR = [thetaInput rInput] ; %Solution = [rInput thetaInput] % Use the following parameters for the inversion algorithm. ConvergenceThreshold % is the required Euclidean distance in xy coordinates between the input % xy, and the xy corresponding to a Munsell sample. GreyThreshold is the required % distance from the origin for a colour to be considered grey; it is defined % separately from ConvergenceThreshold for numerical reasons. ConvergenceThreshold = 0.0001 ; GreyThreshold = 0.001 ; MaxNumOfTries = 60 ; NumOfTries = 0 ; % Check for grey if rInput < GreyThreshold MunsellVec = [MunsellValue] ; MunsellSpec = ColorLabFormatToMunsellSpec(MunsellVec) ; Status.ind = 1 ; return ; end if ShowFigure == TRUE plot(xcenter, ycenter, '+', 'color', [0.6 0.6 0.6]); hold on plot(x, y, 'kx'); hold on axis('equal'); end % The input xyY are the CIE coordinates of Illuminant C, when reflected % off a Munsell sample [X Y Z] = xyY2XYZ(x, y, Y) ; ReflectedXYZ = [X Y Z] ; % Since Y is a percentage reflection, we want to find the CIE coordinates % for illuminant C, whose Y value is 100. This is the reference % illuminant which is reflected off the Munsell sample. % The chromaticity coordinates for Illuminant C are [xC, yC] = ChromaticityOfWhitePoint('C/2') ; % change made Feb. 16., 2014, by Paul Centore % Convert to XYZ form, and normalize to have luminance 100. This will % be the reference illuminant. [Xr Yr Zr] = xyY2XYZ(xC, yC, Y) ; ReferenceXYZ = [(100/Yr)*Xr, 100, (100/Yr)*Zr]; % ReferenceXYZ are the CIE coordinates for a reference light, whose power spectral distribution % is consistent with Illuminant C, and whose luminance is 100. The relative % luminance of the input sample, when illuminated by Illuminant C, is given by Y, which % varies between 0 and 100, and is the percentage of illuminant C that the sample % reflects. Therefore, RefXYZ can be viewed as the reference illuminant for the % sample, and XYZ as the CIE coordinates for the reflected illuminant. These are the % data needed to apply the CIELAB model. Lab = xyz2lab(ReflectedXYZ, ReferenceXYZ) ; LhC = lab2perc(Lab) ; % Find an initial Munsell specification for the interpolation algorithm, by using a % rough transformation from CIELAB to Munsell coordinates. InitialMunsellSpec is % not an exact inverse for xyY, but should be close. L = LhC(1) ; % Lightness in CIELAB system hab = LhC(2) ; % Hue angle in CIELAB system Cab = LhC(3) ; % C*ab in CIELAB system InitialMunsellSpec = CIELABtoApproxMunsellSpec(L, Cab, hab) ; % Replace value of initial Munsell estimate with known Munsell value TempCLvec = MunsellSpecToColorLabFormat(InitialMunsellSpec) ; TempCLvec(2) = MunsellValue ; % Deliberately underestimate chroma to avoid chromas beyond extrapolated % renotation values TempCLvec(3) = (5/5.5) * TempCLvec(3) ; if ShowFigure == TRUE [xCur yCur YCur StatusCode] = MunsellToxyY(TempCLvec); plot(xCur, yCur, 'gx'); hold on text(xCur, yCur, [num2str(NumOfTries)]); hold on end % Use a loop to iterate over different Munsell specifications, whose xyY coordinates % should progressively approach the input xyY. CurrentCLVec = TempCLvec ; EuclideanDifference = -99 ; while NumOfTries <= MaxNumOfTries % After adjusting chroma without adjusting hue, the next iteration adjusts hue % without adjusting chroma. NumOfTries = NumOfTries + 1; % Extract Munsell quantities from Munsell specification CurrentCLHueNumber = CurrentCLVec(1); CurrentMunsellChroma = CurrentCLVec(3); CurrentCLHueLetterIndex = CurrentCLVec(4); CurrentChromDiagHueAngle = MunsellHueToChromDiagHueAngle(CurrentCLHueNumber,CurrentCLHueLetterIndex); % Check that current chroma is possible for the hue and value in the Munsell specification. If the % chroma is too high, set it to be the maximum chroma. [MaxChroma StatusCode] = MaxChromaForExtrapolatedRenotation(CurrentCLHueNumber, CurrentCLHueLetterIndex, MunsellValue); if CurrentMunsellChroma > MaxChroma CurrentMunsellChroma = MaxChroma ; CurrentCLVec(3) = CurrentMunsellChroma ; end % Find xy coordinates corresponding to the current Munsell specification. [xCur yCur YCur StatusCode] = MunsellToxyY(CurrentCLVec); if StatusCode.ind ~= 1 % disp(['Status 2, line 383']) Status.ind = 2; Status.dist = EuclideanDifference; Status.num = NumOfTries; return end % Hue angles correspond to values of theta. The new hue angle will differ from the % current hue angle approximately as much as the desired theta value (thetaInput) differs from the % current theta value. Call this difference thetaDiff and calculate it. Other % hue angles, with their corresponding thetas and theta differences, will be tried. % Make a list, thetaDiffsVec, of the corresponding theta differences. [thetaCur rCur] = cart2pol(xCur-xcenter, yCur-ycenter); thetaCur = mod((180/pi)*thetaCur, 360) ; % Express in degrees thetaDiff = mod(360 - thetaInput + thetaCur, 360); if thetaDiff > 180 % Adjust for wraparound if necessary thetaDiff = thetaDiff-360 ; end thetaDiffsVec = [thetaDiff] ; % Start a similar list for hue angles that correspond to the theta differences. ChromDiagHueAngles = [CurrentChromDiagHueAngle] ; % Also make a list of how much the new hue angles differ from CurrentChromDiagHueAngle. These % angles will be near zero, and will avoid potential problems with wraparound. ChromDiagHueAngleDiffs = [0]; % Ideally, thetaDiff will be 0 % (thetaInput will agree with the theta corresponding to the new hue angle). % Continue constructing the list of theta differences until it contains both % negative and positive differences; then find the new hue angle by linear interpolation % at the value thetaDiff = 0. ctr = 0; AttemptExtrapolation = FALSE; while sign(min(thetaDiffsVec)) == sign(max(thetaDiffsVec)) && AttemptExtrapolation == FALSE ctr = ctr + 1; if ctr > 10 % Too many attempts. Return with error message Status.ind = 3 ; return; end % Find another hue angle, by increasing the coefficient of (thetaInput-thetaCur). % Construct a trial Munsell specification by using the new hue angle. TempChromDiagHueAngle = mod(CurrentChromDiagHueAngle + ctr*(thetaInput - thetaCur),360); % Record the difference of the trial angle from the current hue angle TempChromDiagHueAngleDiff = mod(ctr*(thetaInput - thetaCur), 360) ; if TempChromDiagHueAngleDiff > 180 TempChromDiagHueAngleDiff = TempChromDiagHueAngleDiff - 360; end [TempHueNumber,TempHueLetterCode] = ChromDiagHueAngleToMunsellHue(TempChromDiagHueAngle); TempCLVec = [TempHueNumber, MunsellValue, CurrentMunsellChroma, TempHueLetterCode]; % Evaluate the trial Munsell specification, convert to polar coordinates, and % calculate the difference of the resulting trial theta and thetaInput [xCurTemp yCurTemp YCurTemp StatusCode] = MunsellToxyY(TempCLVec); if StatusCode.ind ~= 1 % Interpolation is impossible because there are not both positive and negative differences, % but extrapolation is possible if there are at least two data points already if length(thetaDiffsVec) >= 2 AttemptExtrapolation = TRUE ; else % disp(['Status 2, line 442']) Status.ind = 2; Status.dist = EuclideanDifference; Status.num = NumOfTries; return end end if AttemptExtrapolation == FALSE [thetaCurTemp rCurTemp] = cart2pol(xCurTemp-xcenter, yCurTemp-ycenter); thetaCurTemp = mod((180/pi)*thetaCurTemp, 360) ; % Express in degrees thetaDiff = mod(360 - thetaInput + thetaCurTemp, 360) ; if thetaDiff > 180 % Adjust for wraparound if necessary thetaDiff = thetaDiff-360 ; end % Add trial hue angle and theta difference to lists thetaDiffsVec = [thetaDiffsVec, thetaDiff] ; ChromDiagHueAngleDiffs = [ChromDiagHueAngleDiffs, TempChromDiagHueAngleDiff]; ChromDiagHueAngles = [ChromDiagHueAngles, TempChromDiagHueAngle] ; end end % Since the while loop exited successfully, both negative and positive theta % differences have been found, so an extrapolation should % be attempted. Interpolate linearly to estimate the hue % angle that corresponds to a theta difference of 0 [thetaDiffsVecSort I] = sort(thetaDiffsVec) ; ChromDiagHueAnglesSort = ChromDiagHueAngles(I) ; ChromDiagHueAngleDiffs = ChromDiagHueAngleDiffs(I); % The extrapolation option will be used if there is not sufficient data % for interpolation NewChromDiagHueAngleDiff = mod(interp1(thetaDiffsVecSort, ChromDiagHueAngleDiffs, 0, 'linear', 'extrap'),360); NewChromDiagHueAngle = mod(CurrentChromDiagHueAngle + NewChromDiagHueAngleDiff,360); % Adjust the current Munsell specification by replacing the current hue with the % new hue. [NewHueNumber,NewHueLetterCode] = ChromDiagHueAngleToMunsellHue(NewChromDiagHueAngle); CurrentCLVec = [NewHueNumber, MunsellValue, CurrentMunsellChroma, NewHueLetterCode]; if NumOfTries >= 2000 % Use for debugging non-converging cases NumOfTries thetaDiffsVecSort ChromDiagHueAngleDiffs ChromDiagHueAnglesSort NewChromDiagHueAngle [ctr CurrentCLVec] end % Calculate the Euclidean distance between the xy coordinates of the newly % constructed Munsell specification, and the input xy [xCur yCur YCur StatusCode] = MunsellToxyY(CurrentCLVec); if StatusCode.ind ~= 1 % disp(['Status 2, line 494']) % CurrentCLVec % [xCur yCur YCur] % StatusCode Status.ind = 2; Status.dist = EuclideanDifference; Status.num = NumOfTries; return end EuclideanDifference = sqrt(((x-xCur)*(x-xCur)) + ((y-yCur)*(y-yCur))); if ShowFigure == TRUE plot(xCur, yCur, 'gx'); hold on text(xCur, yCur, [num2str(NumOfTries)]); hold on end % If the two xy coordinate pairs are close enough, then exit with a success message if EuclideanDifference < ConvergenceThreshold % Current Munsell spec is close enough MunsellVec = CurrentCLVec; MunsellSpec = ColorLabFormatToMunsellSpec(MunsellVec) ; Status.ind = 1; return end NumOfTries = NumOfTries + 1 ; % Extract Munsell quantities from Munsell specification CurrentCLHueNumber = CurrentCLVec(1); CurrentMunsellChroma = CurrentCLVec(3); CurrentCLHueLetterIndex = CurrentCLVec(4); CurrentChromDiagHueAngle = MunsellHueToChromDiagHueAngle(CurrentCLHueNumber,CurrentCLHueLetterIndex); % Check that current chroma is possible for the hue and value in the Munsell specification. If the % chroma is too high, set it to be the maximum chroma. [MaxChroma StatusCode] = MaxChromaForExtrapolatedRenotation(CurrentCLHueNumber, CurrentCLHueLetterIndex, MunsellValue); if CurrentMunsellChroma > MaxChroma CurrentMunsellChroma = MaxChroma ; CurrentCLVec(3) = CurrentMunsellChroma ; end %[xtmop ytmop Ytmop StstTmop] = MunsellToxyY(CurrentCLVec); %[CurrentCLVec, xtmop ytmop] % Find xy coordinates corresponding to the current Munsell specification. [xCur yCur YCur StatusCode] = MunsellToxyY(CurrentCLVec); if StatusCode.ind ~= 1 % disp(['Status 2, line 542']) Status.ind = 2; Status.dist = EuclideanDifference; Status.num = NumOfTries; return end [thetaCur rCur] = cart2pol(xCur-xcenter, yCur-ycenter) ; thetaCur = mod((180/pi)*thetaCur, 360) ; % Express theta of current point in degrees %TempData = [rCur thetaCur xCur yCur] % For this iteration, keep hue (corresponding to thetaCur) constant, % and let chroma (corresponding to rCur) vary. % Construct a set of chromas, called TempMunsellChromas, and find the (r, theta) % values for Munsell specifications with those chromas. % Make a list of r values, called rTempValues, that correspond to different chromas rTempValues = [rCur] ; TempMunsellChromas = [CurrentMunsellChroma]; if NumOfTries >= 2000 % Use for debugging non-converging cases OneChroma = [NumOfTries rTempValues TempMunsellChromas] end % In order to interpolate, rInput must be within the span of the r values % in rTempValues. Check this condition, and add more chromas and r values % until it is satisfied. Usually, no more than three chromas are necessary. ctr = 0; while rInput < min(rTempValues) || rInput > max(rTempValues) ctr = ctr + 1; if ctr > 10 % Too many attempts to bound rInput. Return with error message Status.ind = 3 ; return; end % Try a new chroma, by increasing the exponent on rInput/rCur TempMunsellChroma = ((rInput./rCur)^ctr) * CurrentMunsellChroma; % Check that current chroma is possible for the hue and value in the Munsell specification. If the % chroma is too high, set it to be the maximum chroma. if TempMunsellChroma > MaxChroma TempMunsellChroma = MaxChroma ; CurrentCLVec(3) = TempMunsellChroma ; end if NumOfTries >= 20000 % Use for debugging non-converging cases [NumOfTries ctr TempMunsellChroma rInput rCur rInput./rCur CurrentMunsellChroma] end % Find (r, theta) values for a Munsell specification which is identical to % the current Munsell specification, except that the chroma is TempMunsellChroma TempCLVec = [CurrentCLHueNumber, MunsellValue, TempMunsellChroma, CurrentCLHueLetterIndex]; [xCurTemp yCurTemp YCurTemp StatusCode] = MunsellToxyY(TempCLVec); if StatusCode.ind ~= 1 %disp(['Status 2, line 597']) Status.ind = 2; Status.dist = EuclideanDifference; Status.num = NumOfTries; return end [thetaCurTemp rCurTemp] = cart2pol(xCurTemp-xcenter, yCurTemp-ycenter); thetaCurTemp = mod((180/pi)*thetaCurTemp, 360) ; % Express in degrees %temprTheta = [thetaCurTemp rCurTemp] % Add r and chroma to lists rTempValues = [rTempValues rCurTemp]; TempMunsellChromas = [TempMunsellChromas TempMunsellChroma]; if NumOfTries >= 2000 % Use for debugging non-converging cases rTempValues TempMunsellChromas end end % Since the while loop has been exited, we have found r values, resulting from % different chroma values, that bound rInput. Linearly interpolate to find a further % chroma, which should be an even better approximation. % Linear interpolation requires the list of r values to be sorted. [rTempSort I] = sort(rTempValues); TempMunsellChromasSort = TempMunsellChromas(I); NewMunsellChroma = interp1(rTempSort, TempMunsellChromasSort, rInput); % Adjust the current Munsell specification, by using the new chroma CurrentCLVec = [CurrentCLHueNumber, MunsellValue, NewMunsellChroma, CurrentCLHueLetterIndex]; if NumOfTries >= 2000 % Use for debugging non-converging cases NumOfTries rTempSort TempMunsellChromasSort NewMunsellChroma CurrentCLVec end % Calculate the Euclidean distance between the xy coordinates of the newly % constructed Munsell specification, and the input xy [xCur yCur YCur StatusCode] = MunsellToxyY(CurrentCLVec); if StatusCode.ind ~= 1 %disp(['Status 2, line 636']) Status.ind = 2; Status.dist = EuclideanDifference; Status.num = NumOfTries; return end EuclideanDifference = sqrt(((x-xCur)*(x-xCur)) + ((y-yCur)*(y-yCur))); if ShowFigure == TRUE plot(xCur, yCur, 'gx'); hold on text(xCur, yCur, [num2str(NumOfTries)]); hold on end % If the two xy coordinate pairs are close enough, then exit with a success message if EuclideanDifference < ConvergenceThreshold MunsellVec = CurrentCLVec ; MunsellSpec = ColorLabFormatToMunsellSpec(MunsellVec) ; Status.ind = 1 ; return end end % If the routine did not exit from the while loop, then the maximum number of % iterations has been tried. Set an error message and return. Status.ind = 3 ; Status.dist = EuclideanDifference ; Status.num = NumOfTries ; return