function MunsellValues = LuminanceFactorToMunsellValue(LumFac);
% Purpose		Convert the luminance factor of an object colour, into a
%				Munsell value
%
% Description	A non-luminous sample of a colour reflects a certain percentage
%				of a light that illuminates it.  For human colour vision, this
%				percentage is calculated with respect to the human luminous
%				efficiency function, given by the photopic data for the CIE 1931
%				standard observer.  The Munsell value is a similar concept,
%				quantifying how light or dark a non-luminous colour sample is.
%				The reflected percentage of light is called the luminance factor.
%
%				The Munsell renotation ([Newhall1943]) redefined Munsell value as
%				a function of the percentage of light a sample reflects when
%				illuminated by Illuminant C.  Later, colour appearance models such as
%				CIELAB also expressed Munsell value as similar, but not identical,
%				functions of colorimetric measurements.  [McCamy1992] gives another
%				expression in closed form.
%
%				This routine calculates multiple possible Munsell values.  Currently, the first
%				value comes from a quintic polynomial suggested in [Newhall1943, p.
%				417].  The second value is taken from the CIELAB colour appearance
%				model, as presented in ([Fairchild2005, Sect. 10.3]).  The third
%				value uses the expression in [McCamy1992].  The fourth value is obtained by
%				linear interpolation over a precalculated table of luminance factors and
%				values from the ASTM quintic ([ASTMD1535-08, Eq. 2]).  The possible Munsell values
%				are returned in a data structure, which can be easily extended to include
%				other conversion expressions.
%
%				[ASTMD1535-08] ASTM, Standard D 1535-08, "Standard Practice for Specifying Color by the
%					Munsell System," approved January 1, 2008.
%				[Fairchild2005] Mark D. Fairchild, Color Appearance Models, 2nd ed.,
%					John Wiley & Sons, Ltd., 2005.
%				[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		MunsellValues = LuminanceFactorToMunsellValue(LumFac);
%
%				LumFac			Percentage of light reflected by a non-luminous
%								colour sample, measured with respect to the photopic
%								luminous efficiency function.  This variable should be between
%								0 and 100.
%
%				MunsellValues	A data structure.  The first field is the input luminance
%								factor that is being converted.  The remaining fields are
%								the Munsell values when calculated by different methods.
%
% Related		MunsellValueToLuminanceFactor
% Functions
%
% Required		None
% Functions
%
% Author		Paul Centore (May 11, 2010)
% Revised by	Paul Centore (January 22, 2011)
%	Revisions:  Changed name from ReflectancePercentageToMunsellValue to LuminanceFactorToMunsellValue
%				Changed output format from vector to data structure, for easy extensibility
%				Added McCamy conversion expression
%				Added ASTM Table Lookup option
% Revision		Paul Centore (December 26, 2012)  
%				 ---Moved from MunsellConversions program to MunsellToolbox.
% Revision		Paul Centore (August 31, 2013)  
%				 ---Moved from MunsellToolbox program to MunsellAndKubelkaMunkToolbox.
%
% Copyright 2010, 2011, 2012 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 <http://www.gnu.org/licenses/>.

% Calculate look-up table of luminance factors from quintic polynomial in 
% [Newhall1943].  Once the look-up table is made, interpolate linearly to find the 
% Munsell value that corresponds to a given luminance factor.
AllMunsVals = [0:0.1:10.0]						;
% The first option of the if-statement calculates a list of luminance
% factors, and outputs them to the command screen.  The list can then be
% copied and pasted into the second option of the if-statement.  In this way,
% the calculation is done only once.
if 0
   Term1 =  1.2219    * AllMunsVals				;
   Term2 = -0.23111   * (AllMunsVals.^2)			;
   Term3 =  0.23951   * (AllMunsVals.^3)			;
   Term4 = -0.021009  * (AllMunsVals.^4)			;
   Term5 =  0.0008404 * (AllMunsVals.^5)			;
   AllLumFaccs = Term1 + Term2 + Term3 + Term4 + Term5		;
   for i = 1:length(AllLumFaccs)
      disp([sprintf('      %7.4f,...', AllLumFaccs(i))]);
   end
else					% Values generated by running code immediately above
   AllLumFaccs = ...
      [0.0000,...
       0.1201,...
       0.2370,...
       0.3521,...
       0.4666,...
       0.5818,...
       0.6990,...
       0.8193,...
       0.9439,...
       1.0738,...
       1.2101,...
       1.3538,...
       1.5059,...
       1.6672,...
       1.8387,...
       2.0212,...
       2.2156,...
       2.4225,...
       2.6428,...
       2.8771,...
       3.1262,...
       3.3906,...
       3.6711,...
       3.9681,...
       4.2822,...
       4.6141,...
       4.9641,...
       5.3327,...
       5.7204,...
       6.1277,...
       6.5550,...
       7.0025,...
       7.4708,...
       7.9601,...
       8.4709,...
       9.0033,...
       9.5577,...
      10.1344,...
      10.7336,...
      11.3557,...
      12.0007,...
      12.6691,...
      13.3610,...
      14.0766,...
      14.8161,...
      15.5797,...
      16.3676,...
      17.1801,...
      18.0172,...
      18.8791,...
      19.7661,...
      20.6783,...
      21.6159,...
      22.5791,...
      23.5680,...
      24.5829,...
      25.6239,...
      26.6912,...
      27.7850,...
      28.9055,...
      30.0529,...
      31.2275,...
      32.4294,...
      33.6590,...
      34.9165,...
      36.2021,...
      37.5161,...
      38.8589,...
      40.2307,...
      41.6319,...
      43.0628,...
      44.5238,...
      46.0153,...
      47.5377,...
      49.0914,...
      50.6769,...
      52.2946,...
      53.9451,...
      55.6289,...
      57.3466,...
      59.0986,...
      60.8858,...
      62.7086,...
      64.5678,...
      66.4641,...
      68.3982,...
      70.3710,...
      72.3831,...
      74.4355,...
      76.5291,...
      78.6647,...
      80.8434,...
      83.0661,...
      85.3339,...
      87.6479,...
      90.0092,...
      92.4190,...
      94.8784,...
      97.3889,...
      99.9516,...
      102.5680]				;
end
NewhallMunsellValue = interp1(AllLumFaccs, AllMunsVals, LumFac);

% Calculate lightness (L*) from CIELAB model ([Fairchild2005, Sect. 10.3]),
% and divide by 10 to get the Munsell value.
LumFacAsFrac = LumFac/100	;
if LumFacAsFrac > 0.008856
   fofomega = LumFacAsFrac^(1/3)	;
else
   fofomega = 7.787*LumFacAsFrac + (16/116)	;
end
Lstar = 116*fofomega - 16		;
CIELABMunsellValue = Lstar/10	;

% Calculate Munsell Value using the closed-form approximation from [McCamy1992].
% This approximation has been incorporated in [ASTMD1535-08]
Y = LumFac	;
if Y <= 0.9
   McCamyMunsellValue = 0.87445 * (Y^0.9967);
else
   McCamyMunsellValue = 2.49268*(Y^(1/3)) - 1.5614 ...
                        - (0.985/(((0.1073*Y-3.084)^2) + 7.54)) ...
						+ (0.0133/(Y^2.3)) ...
						+ 0.0084*sin(4.1*(Y^(1/3))+1) ...
						+ (0.0221/Y)*sin(0.39*(Y-2)) ...
						- (0.0037/(0.44*Y)) * sin(1.28*(Y-0.53));
end

AllMunsVals = [0:0.02:10.0]						;
if 0
% The first option of the if-statement calculates a list of luminance
% factors, using the expression in ([ASTMD1535-08, Eq. 2]),
% for the values in AllMunsVals,
% and outputs them to the command screen.  The list can then be
% copied and pasted into the second option of the if-statement.  In this way,
% the calculation is done only once.
   Term1 =  1.1914    * AllMunsVals				;
   Term2 = -0.22533   * (AllMunsVals.^2)			;
   Term3 =  0.23352   * (AllMunsVals.^3)			;
   Term4 = -0.020484  * (AllMunsVals.^4)			;
   Term5 =  0.00081939 * (AllMunsVals.^5)		;
   AllASTMLuminanceFactors = Term1 + Term2 + Term3 + Term4 + Term5	;
   for i = 1:length(AllASTMLuminanceFactors)
      disp([sprintf('      %7.6f,...', AllASTMLuminanceFactors(i))]);
   end
else
   AllASTMLuminanceFactors = ...
	 [0.000000,...
      0.023740,...
      0.047310,...
      0.070723,...
      0.093989,...
      0.117118,...
      0.140123,...
      0.163012,...
      0.185799,...
      0.208492,...
      0.231102,...
      0.253641,...
      0.276118,...
      0.298543,...
      0.320928,...
      0.343281,...
      0.365614,...
      0.387936,...
      0.410257,...
      0.432587,...
      0.454936,...
      0.477314,...
      0.499730,...
      0.522194,...
      0.544715,...
      0.567303,...
      0.589967,...
      0.612717,...
      0.635561,...
      0.658509,...
      0.681571,...
      0.704754,...
      0.728068,...
      0.751522,...
      0.775125,...
      0.798885,...
      0.822812,...
      0.846913,...
      0.871197,...
      0.895673,...
      0.920349,...
      0.945234,...
      0.970336,...
      0.995662,...
      1.021222,...
      1.047023,...
      1.073073,...
      1.099381,...
      1.125954,...
      1.152799,...
      1.179925,...
      1.207340,...
      1.235051,...
      1.263065,...
      1.291391,...
      1.320035,...
      1.349005,...
      1.378308,...
      1.407952,...
      1.437944,...
      1.468291,...
      1.498999,...
      1.530077,...
      1.561531,...
      1.593367,...
      1.625594,...
      1.658217,...
      1.691243,...
      1.724679,...
      1.758532,...
      1.792808,...
      1.827513,...
      1.862655,...
      1.898239,...
      1.934272,...
      1.970759,...
      2.007709,...
      2.045125,...
      2.083016,...
      2.121386,...
      2.160241,...
      2.199588,...
      2.239433,...
      2.279781,...
      2.320638,...
      2.362010,...
      2.403902,...
      2.446321,...
      2.489271,...
      2.532759,...
      2.576790,...
      2.621368,...
      2.666500,...
      2.712192,...
      2.758447,...
      2.805272,...
      2.852671,...
      2.900650,...
      2.949214,...
      2.998368,...
      3.048116,...
      3.098465,...
      3.149418,...
      3.200980,...
      3.253157,...
      3.305953,...
      3.359373,...
      3.413421,...
      3.468102,...
      3.523421,...
      3.579382,...
      3.635990,...
      3.693248,...
      3.751163,...
      3.809737,...
      3.868975,...
      3.928881,...
      3.989460,...
      4.050716,...
      4.112653,...
      4.175275,...
      4.238586,...
      4.302590,...
      4.367291,...
      4.432693,...
      4.498800,...
      4.565616,...
      4.633144,...
      4.701389,...
      4.770353,...
      4.840042,...
      4.910458,...
      4.981605,...
      5.053487,...
      5.126107,...
      5.199468,...
      5.273575,...
      5.348431,...
      5.424038,...
      5.500401,...
      5.577523,...
      5.655406,...
      5.734055,...
      5.813473,...
      5.893662,...
      5.974626,...
      6.056368,...
      6.138891,...
      6.222198,...
      6.306293,...
      6.391178,...
      6.476856,...
      6.563330,...
      6.650603,...
      6.738677,...
      6.827557,...
      6.917244,...
      7.007741,...
      7.099052,...
      7.191178,...
      7.284123,...
      7.377889,...
      7.472478,...
      7.567894,...
      7.664139,...
      7.761215,...
      7.859126,...
      7.957873,...
      8.057459,...
      8.157887,...
      8.259158,...
      8.361276,...
      8.464242,...
      8.568059,...
      8.672730,...
      8.778256,...
      8.884640,...
      8.991885,...
      9.099991,...
      9.208963,...
      9.318801,...
      9.429508,...
      9.541086,...
      9.653537,...
      9.766863,...
      9.881067,...
      9.996150,...
      10.112115,...
      10.228963,...
      10.346696,...
      10.465317,...
      10.584827,...
      10.705228,...
      10.826523,...
      10.948712,...
      11.071799,...
      11.195785,...
      11.320671,...
      11.446459,...
      11.573152,...
      11.700751,...
      11.829258,...
      11.958675,...
      12.089003,...
      12.220244,...
      12.352400,...
      12.485472,...
      12.619463,...
      12.754374,...
      12.890206,...
      13.026961,...
      13.164642,...
      13.303248,...
      13.442783,...
      13.583248,...
      13.724643,...
      13.866972,...
      14.010235,...
      14.154434,...
      14.299571,...
      14.445646,...
      14.592662,...
      14.740620,...
      14.889522,...
      15.039369,...
      15.190162,...
      15.341904,...
      15.494595,...
      15.648237,...
      15.802831,...
      15.958379,...
      16.114883,...
      16.272343,...
      16.430761,...
      16.590140,...
      16.750479,...
      16.911780,...
      17.074046,...
      17.237277,...
      17.401474,...
      17.566640,...
      17.732774,...
      17.899880,...
      18.067958,...
      18.237010,...
      18.407037,...
      18.578040,...
      18.750021,...
      18.922981,...
      19.096921,...
      19.271844,...
      19.447749,...
      19.624640,...
      19.802516,...
      19.981380,...
      20.161233,...
      20.342076,...
      20.523910,...
      20.706737,...
      20.890559,...
      21.075377,...
      21.261191,...
      21.448004,...
      21.635817,...
      21.824631,...
      22.014448,...
      22.205269,...
      22.397096,...
      22.589929,...
      22.783771,...
      22.978623,...
      23.174486,...
      23.371361,...
      23.569251,...
      23.768157,...
      23.968079,...
      24.169020,...
      24.370981,...
      24.573964,...
      24.777969,...
      24.982999,...
      25.189055,...
      25.396139,...
      25.604251,...
      25.813394,...
      26.023570,...
      26.234779,...
      26.447023,...
      26.660305,...
      26.874624,...
      27.089984,...
      27.306386,...
      27.523831,...
      27.742321,...
      27.961858,...
      28.182443,...
      28.404078,...
      28.626765,...
      28.850505,...
      29.075300,...
      29.301153,...
      29.528064,...
      29.756035,...
      29.985069,...
      30.215167,...
      30.446331,...
      30.678562,...
      30.911863,...
      31.146236,...
      31.381682,...
      31.618204,...
      31.855802,...
      32.094480,...
      32.334239,...
      32.575081,...
      32.817008,...
      33.060023,...
      33.304126,...
      33.549321,...
      33.795610,...
      34.042994,...
      34.291476,...
      34.541057,...
      34.791741,...
      35.043528,...
      35.296423,...
      35.550425,...
      35.805539,...
      36.061766,...
      36.319109,...
      36.577570,...
      36.837151,...
      37.097854,...
      37.359683,...
      37.622640,...
      37.886726,...
      38.151945,...
      38.418299,...
      38.685791,...
      38.954422,...
      39.224197,...
      39.495117,...
      39.767186,...
      40.040405,...
      40.314777,...
      40.590306,...
      40.866994,...
      41.144844,...
      41.423859,...
      41.704041,...
      41.985394,...
      42.267920,...
      42.551623,...
      42.836505,...
      43.122569,...
      43.409819,...
      43.698257,...
      43.987888,...
      44.278713,...
      44.570736,...
      44.863961,...
      45.158390,...
      45.454027,...
      45.750875,...
      46.048938,...
      46.348219,...
      46.648721,...
      46.950448,...
      47.253403,...
      47.557591,...
      47.863014,...
      48.169676,...
      48.477581,...
      48.786732,...
      49.097134,...
      49.408790,...
      49.721704,...
      50.035879,...
      50.351320,...
      50.668031,...
      50.986015,...
      51.305276,...
      51.625820,...
      51.947649,...
      52.270767,...
      52.595180,...
      52.920891,...
      53.247904,...
      53.576225,...
      53.905856,...
      54.236803,...
      54.569070,...
      54.902662,...
      55.237582,...
      55.573836,...
      55.911428,...
      56.250364,...
      56.590646,...
      56.932281,...
      57.275273,...
      57.619628,...
      57.965348,...
      58.312441,...
      58.660911,...
      59.010762,...
      59.362000,...
      59.714631,...
      60.068658,...
      60.424088,...
      60.780926,...
      61.139177,...
      61.498846,...
      61.859939,...
      62.222462,...
      62.586419,...
      62.951817,...
      63.318661,...
      63.686957,...
      64.056710,...
      64.427926,...
      64.800612,...
      65.174773,...
      65.550415,...
      65.927544,...
      66.306166,...
      66.686287,...
      67.067913,...
      67.451051,...
      67.835707,...
      68.221887,...
      68.609598,...
      68.998846,...
      69.389637,...
      69.781979,...
      70.175877,...
      70.571338,...
      70.968370,...
      71.366978,...
      71.767171,...
      72.168954,...
      72.572334,...
      72.977319,...
      73.383917,...
      73.792132,...
      74.201974,...
      74.613450,...
      75.026566,...
      75.441330,...
      75.857750,...
      76.275833,...
      76.695586,...
      77.117018,...
      77.540136,...
      77.964947,...
      78.391460,...
      78.819682,...
      79.249622,...
      79.681288,...
      80.114687,...
      80.549827,...
      80.986718,...
      81.425366,...
      81.865781,...
      82.307971,...
      82.751945,...
      83.197711,...
      83.645277,...
      84.094652,...
      84.545845,...
      84.998866,...
      85.453722,...
      85.910422,...
      86.368977,...
      86.829394,...
      87.291683,...
      87.755853,...
      88.221914,...
      88.689875,...
      89.159745,...
      89.631535,...
      90.105252,...
      90.580908,...
      91.058511,...
      91.538072,...
      92.019601,...
      92.503107,...
      92.988601,...
      93.476093,...
      93.965592,...
      94.457109,...
      94.950655,...
      95.446239,...
      95.943873,...
      96.443567,...
      96.945331,...
      97.449176,...
      97.955113,...
      98.463153,...
      98.973307,...
      99.485586,...
      100.000000];
end
ASTMTableLookupMunsellValue = interp1(AllASTMLuminanceFactors, AllMunsVals, LumFac);

MunsellValues.OriginalLuminanceFactor = LumFac						;
MunsellValues.Newhall1943             = NewhallMunsellValue			;
MunsellValues.CIELAB                  = CIELABMunsellValue			;
MunsellValues.McCamy1992              = McCamyMunsellValue			;
MunsellValues.ASTMTableLookup         = ASTMTableLookupMunsellValue	;
return