查看“︁Module:Hct/HctSolver”︁的源代码

--[[
  Copyright 2021 Google LLC
 
  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at
 
      http://www.apache.org/licenses/LICENSE-2.0
 
  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
]]
 
--[[
  This file has been modified. The original version is at
      https://github.com/material-foundation/material-color-utilities
]]

local colorUtils = require('Module:Hct/ColorUtils');
local mathUtils = require('Module:Hct/MathUtils');
local Cam16 = require('Module:Hct/Cam16');


--[[ A class that solves the HCT equation. ]]
local HctSolver = {
	SCALED_DISCOUNT_FROM_LINRGB = {
		{0.001200833568784504, 0.002389694492170889, 0.0002795742885861124},
		{0.0005891086651375999, 0.0029785502573438758, 0.0003270666104008398},
		{0.00010146692491640572, 0.0005364214359186694, 0.0032979401770712076},
	};
	
	LINRGB_FROM_SCALED_DISCOUNT = {
		{1373.2198709594231, -1100.4251190754821, -7.278681089101213},
		{-271.815969077903, 559.6580465940733, -32.46047482791194},
		{1.9622899599665666, -57.173814538844006, 308.7233197812385},
	};
	
	Y_FROM_LINRGB = {0.2126, 0.7152, 0.0722};

	CRITICAL_PLANES = {
		0.015176349177441876, 0.045529047532325624, 0.07588174588720938,
		0.10623444424209313,  0.13658714259697685,  0.16693984095186062,
		0.19729253930674434,  0.2276452376616281,   0.2579979360165119,
		0.28835063437139563,  0.3188300904430532,   0.350925934958123,
		0.3848314933096426,   0.42057480301049466,  0.458183274052838,
		0.4976837250274023,   0.5391024159806381,   0.5824650784040898,
		0.6277969426914107,   0.6751227633498623,   0.7244668422128921,
		0.775853049866786,    0.829304845476233,    0.8848452951698498,
		0.942497089126609,    1.0022825574869039,   1.0642236851973577,
		1.1283421258858297,   1.1946592148522128,   1.2631959812511864,
		1.3339731595349034,   1.407011200216447,    1.4823302800086415,
		1.5599503113873272,   1.6398909516233677,   1.7221716113234105,
		1.8068114625156377,   1.8938294463134073,   1.9832442801866852,
		2.075074464868551,    2.1693382909216234,   2.2660538449872063,
		2.36523901573795,     2.4669114995532007,   2.5710888059345764,
		2.6777882626779785,   2.7870270208169257,   2.898822059350997,
		3.0131901897720907,   3.1301480604002863,   3.2497121605402226,
		3.3718988244681087,   3.4967242352587946,   3.624204428461639,
		3.754355295633311,    3.887192587735158,    4.022731918402185,
		4.160988767090289,    4.301978482107941,    4.445716283538092,
		4.592217266055746,    4.741496401646282,    4.893568542229298,
		5.048448422192488,    5.20615066083972,     5.3666897647573375,
		5.5300801301023865,   5.696336044816294,    5.865471690767354,
		6.037501145825082,    6.212438385869475,    6.390297286737924,
		6.571091626112461,    6.7548350853498045,   6.941541251256611,
		7.131223617812143,    7.323895587840543,    7.5195704746346665,
		7.7182615035334345,   7.919981813454504,    8.124744458384042,
		8.332562408825165,    8.543448553206703,    8.757415699253682,
		8.974476575321063,    9.194643831691977,    9.417930041841839,
		9.644347703669503,    9.873909240696694,    10.106627003236781,
		10.342513269534024,   10.58158024687427,    10.8238400726681,
		11.069304815507364,   11.317986476196008,   11.569896988756009,
		11.825048221409341,   12.083451977536606,   12.345119996613247,
		12.610063955123938,   12.878295467455942,   13.149826086772048,
		13.42466730586372,    13.702830557985108,   13.984327217668513,
		14.269168601521828,   14.55736596900856,    14.848930523210871,
		15.143873411576273,   15.44220572664832,    15.743938506781891,
		16.04908273684337,    16.35764934889634,    16.66964922287304,
		16.985093187232053,   17.30399201960269,    17.62635644741625,
		17.95219714852476,    18.281524751807332,   18.614349837764564,
		18.95068293910138,    19.290534541298456,   19.633915083172692,
		19.98083495742689,    20.331304511189067,   20.685334046541502,
		21.042933821039977,   21.404114048223256,   21.76888489811322,
		22.137256497705877,   22.50923893145328,    22.884842241736916,
		23.264076429332462,   23.6469514538663,     24.033477234264016,
		24.42366364919083,    24.817520537484558,   25.21505769858089,
		25.61628489293138,    26.021211842414342,   26.429848230738664,
		26.842203703840827,   27.258287870275353,   27.678110301598522,
		28.10168053274597,    28.529008062403893,   28.96010235337422,
		29.39497283293396,    29.83362889318845,    30.276079891419332,
		30.722335150426627,   31.172403958865512,   31.62629557157785,
		32.08401920991837,    32.54558406207592,    33.010999283389665,
		33.4802739966603,     33.953417292456834,   34.430438229418264,
		34.911345834551085,   35.39614910352207,    35.88485700094671,
		36.37747846067349,    36.87402238606382,    37.37449765026789,
		37.87891309649659,    38.38727753828926,    38.89959975977785,
		39.41588851594697,    39.93615253289054,    40.460400508064545,
		40.98864111053629,    41.520882981230194,   42.05713473317016,
		42.597404951718396,   43.141702194811224,   43.6900349931913,
		44.24241185063697,    44.798841244188324,   45.35933162437017,
		45.92389141541209,    46.49252901546552,    47.065252796817916,
		47.64207110610409,    48.22299226451468,    48.808024568002054,
		49.3971762874833,     49.9904556690408,     50.587870934119984,
		51.189430279724725,   51.79514187861014,    52.40501387947288,
		53.0190544071392,     53.637271562750364,   54.259673423945976,
		54.88626804504493,    55.517063457223934,   56.15206766869424,
		56.79128866487574,    57.43473440856916,    58.08241284012621,
		58.734331877617365,   59.39049941699807,    60.05092333227251,
		60.715611475655585,   61.38457167773311,    62.057811747619894,
		62.7353394731159,     63.417162620860914,   64.10328893648692,
		64.79372614476921,    65.48848194977529,    66.18756403501224,
		66.89098006357258,    67.59873767827808,    68.31084450182222,
		69.02730813691093,    69.74813616640164,    70.47333615344107,
		71.20291564160104,    71.93688215501312,    72.67524319850172,
		73.41800625771542,    74.16517879925733,    74.9167682708136,
		75.67278210128072,    76.43322770089146,    77.1981124613393,
		77.96744375590167,    78.74122893956174,    79.51947534912904,
		80.30219030335869,    81.08938110306934,    81.88105503125999,
		82.67721935322541,    83.4778813166706,     84.28304815182372,
		85.09272707154808,    85.90692527145302,    86.72564993000343,
		87.54890820862819,    88.3767072518277,     89.2090541872801,
		90.04595612594655,    90.88742016217518,    91.73345337380438,
		92.58406282226491,    93.43925555268066,    94.29903859396902,
		95.16341895893969,    96.03240364439274,    96.9059996312159,
		97.78421388448044,    98.6670533535366,     99.55452497210776,
	};
};

--[[
  Sanitizes a small enough angle in radians.

  @param angle An angle in radians; must not deviate too much from 0.
  @return A coterminal angle between 0 and 2pi.
]]
function HctSolver.sanitizeRadians(angle)
	return (angle + math.pi * 8) % (math.pi * 2);
end

--[[
  Delinearizes an RGB component, returning a floating-point number.

  @param rgbComponent 0.0 <= rgb_component <= 100.0, represents linear R/G/B channel
  @return 0.0 <= output <= 255.0, color channel converted to regular RGB space
]]
function HctSolver.trueDelinearized(rgbComponent)
	local normalized = rgbComponent / 100.0;
	local delinearized = 0.0;
	if normalized <= 0.0031308 then
		delinearized = normalized * 12.92;
	else
		delinearized = 1.055 * math.pow(normalized, 1.0 / 2.4) - 0.055;
	end
	return delinearized * 255.0;
end

function HctSolver.chromaticAdaptation(component)
	local af = math.pow(math.abs(component), 0.42);
	return mathUtils.signum(component) * 400.0 * af / (af + 27.13);
end

--[[
  Returns the hue of a linear RGB color in CAM16.

  @param linrgb The linear RGB coordinates of a color.
  @return The hue of the color in CAM16, in radians.
]]
function HctSolver.hueOf(linrgb)
	local scaledDiscount =
		mathUtils.matrixMultiply(linrgb, HctSolver.SCALED_DISCOUNT_FROM_LINRGB);
	local rA = HctSolver.chromaticAdaptation(scaledDiscount[1]);
	local gA = HctSolver.chromaticAdaptation(scaledDiscount[2]);
	local bA = HctSolver.chromaticAdaptation(scaledDiscount[3]);
	-- redness-greenness
	local a = (11.0 * rA + -12.0 * gA + bA) / 11.0;
	-- yellowness-blueness
	local b = (rA + gA - 2.0 * bA) / 9.0;
	return math.atan2(b, a);
end

function HctSolver.areInCyclicOrder(a, b, c)
	local deltaAB = HctSolver.sanitizeRadians(b - a);
	local deltaAC = HctSolver.sanitizeRadians(c - a);
	return deltaAB < deltaAC;
end

--[[
  Solves the lerp equation.

  @param source The starting number.
  @param mid The number in the middle.
  @param target The ending number.
  @return A number t such that lerp(source, target, t) = mid.
]]
function HctSolver.intercept(source, mid, target)
	return (mid - source) / (target - source);
end

function HctSolver.lerpPoint(source, t, target)
	return {
		source[1] + (target[1] - source[1]) * t,
		source[2] + (target[2] - source[2]) * t,
		source[3] + (target[3] - source[3]) * t,
	};
end

--[[
  Intersects a segment with a plane.

  @param source The coordinates of point A.
  @param coordinate The R-, G-, or B-coordinate of the plane.
  @param target The coordinates of point B.
  @param axis The axis the plane is perpendicular with. (1: R, 2: G, 3: B)
  @return The intersection point of the segment AB with the plane
      R=coordinate, G=coordinate, or B=coordinate
]]
function HctSolver.setCoordinate(source, coordinate, target, axis)
	local t = HctSolver.intercept(source[axis], coordinate, target[axis]);
	return HctSolver.lerpPoint(source, t, target);
end

function HctSolver.isBounded(x)
	return 0.0 <= x and x <= 100.0;
end

--[[
  Returns the nth possible vertex of the polygonal intersection.

  @param y The Y value of the plane.
  @param n The zero-based index of the point. 0 <= n <= 11.
  @return The nth possible vertex of the polygonal intersection of the y plane
      and the RGB cube, in linear RGB coordinates, if it exists. If this
      possible vertex lies outside of the cube, {-1.0, -1.0, -1.0} is returned.
]]
function HctSolver.nthVertex(y, n)
	local kR = HctSolver.Y_FROM_LINRGB[1];
	local kG = HctSolver.Y_FROM_LINRGB[2];
	local kB = HctSolver.Y_FROM_LINRGB[3];
	local coordA = n % 4 <= 1 and 0.0 or 100.0;
	local coordB = n % 2 == 0 and 0.0 or 100.0;
	if n < 4 then
		local g = coordA;
		local b = coordB;
		local r = (y - g * kG - b * kB) / kR;
		if HctSolver.isBounded(r) then
			return {r, g, b};
		else
			return {-1.0, -1.0, -1.0};
		end
	elseif n < 8 then
		local b = coordA;
		local r = coordB;
		local g = (y - r * kR - b * kB) / kG;
		if HctSolver.isBounded(g) then
			return {r, g, b};
		else
			return {-1.0, -1.0, -1.0};
		end
	else
		local r = coordA;
		local g = coordB;
		local b = (y - r * kR - g * kG) / kB;
		if HctSolver.isBounded(b) then
			return {r, g, b};
		else
			return {-1.0, -1.0, -1.0};
		end
	end
end

--[[
  Finds the segment containing the desired color.

  @param y The Y value of the color.
  @param targetHue The hue of the color.
  @return Two sets of linear RGB coordinates, each corresponding to an endpoint
      of the segment containing the desired color.
]]
function HctSolver.bisectToSegment(y, targetHue)
	local left = {-1.0, -1.0, -1.0};
	local right = left;
	local leftHue = 0.0;
	local rightHue = 0.0;
	local initialized = false;
	local uncut = true;
	for n = 0, 11 do
		local mid = HctSolver.nthVertex(y, n);
		if mid[1] >= 0 then
			local midHue = HctSolver.hueOf(mid);
			if initialized then
				if uncut or HctSolver.areInCyclicOrder(leftHue, midHue, rightHue) then
					uncut = false;
					if HctSolver.areInCyclicOrder(leftHue, targetHue, midHue) then
						right = mid;
						rightHue = midHue;
					else
						left = mid;
						leftHue = midHue;
					end
				end
			else
				left = mid;
				right = mid;
				leftHue = midHue;
				rightHue = midHue;
				initialized = true;
			end
		end
	end
	return left, right;
end

function HctSolver.midpoint(a, b)
	return {
		(a[1] + b[1]) / 2,
		(a[2] + b[2]) / 2,
		(a[3] + b[3]) / 2,
	};
end

function HctSolver.criticalPlaneBelow(x)
	return math.floor(x - 0.5);
end

function HctSolver.criticalPlaneAbove(x)
	return math.ceil(x - 0.5);
end

--[[
  Finds a color with the given Y and hue on the boundary of the
  cube.

  @param y The Y value of the color.
  @param targetHue The hue of the color.
  @return The desired color, in linear RGB coordinates.
]]
function HctSolver.bisectToLimit(y, targetHue)
	local left, right = HctSolver.bisectToSegment(y, targetHue);
	local leftHue = HctSolver.hueOf(left);
	for axis = 1, 3 do
		if left[axis] ~= right[axis] then
			local lPlane = -1;
			local rPlane = 255;
			if left[axis] < right[axis] then
				lPlane = HctSolver.criticalPlaneBelow(
					HctSolver.trueDelinearized(left[axis]));
				rPlane = HctSolver.criticalPlaneAbove(
					HctSolver.trueDelinearized(right[axis]));
			else
				lPlane = HctSolver.criticalPlaneAbove(
					HctSolver.trueDelinearized(left[axis]));
				rPlane = HctSolver.criticalPlaneBelow(
					HctSolver.trueDelinearized(right[axis]));
			end
			for i = 1, 8 do
				if math.abs(rPlane - lPlane) <= 1 then
					break;
				end
				local mPlane = math.floor((lPlane + rPlane) / 2.0);
				local midPlaneCoordinate = HctSolver.CRITICAL_PLANES[mPlane + 1];
				local mid = HctSolver.setCoordinate(left, midPlaneCoordinate, right, axis);
				local midHue = HctSolver.hueOf(mid);
				if HctSolver.areInCyclicOrder(leftHue, targetHue, midHue) then
					right = mid;
					rPlane = mPlane;
				else
					left = mid;
					leftHue = midHue;
					lPlane = mPlane;
				end
			end
		end
	end
	return HctSolver.midpoint(left, right);
end

function HctSolver.inverseChromaticAdaptation(adapted)
	local adaptedAbs = math.abs(adapted);
	local base = math.max(0, 27.13 * adaptedAbs / (400.0 - adaptedAbs));
	return mathUtils.signum(adapted) * math.pow(base, 1.0 / 0.42);
end

--[[
  Finds a color with the given hue, chroma, and Y.

  @param hueRadians The desired hue in radians.
  @param chroma The desired chroma.
  @param y The desired Y.
  @return The desired color as a hexadecimal integer, if found; 0 otherwise.
]]
function HctSolver.findResultByJ(hueRadians, chroma, y)
	-- Initial estimate of j.
	local j = math.sqrt(y) * 11.0;
	-- ===========================================================
	-- Operations inlined from Cam16 to avoid repeated calculation
	-- ===========================================================
	local viewingConditions = require('Module:Hct/ViewingConditions');
	local tInnerCoeff =
		1 / math.pow(1.64 - math.pow(0.29, viewingConditions.n), 0.73);
	local eHue = 0.25 * (math.cos(hueRadians + 2.0) + 3.8);
	local p1 =
		eHue * (50000.0 / 13.0) * viewingConditions.nc * viewingConditions.ncb;
	local hSin = math.sin(hueRadians);
	local hCos = math.cos(hueRadians);
	for iterationRound = 1, 5 do
		-- ===========================================================
		-- Operations inlined from Cam16 to avoid repeated calculation
		-- ===========================================================
		local jNormalized = j / 100.0;
		local alpha =
			(chroma == 0.0 or j == 0.0) and 0.0 or chroma / math.sqrt(jNormalized);
		local t = math.pow(alpha * tInnerCoeff, 1.0 / 0.9);
		local ac = viewingConditions.aw *
			math.pow(jNormalized, 1.0 / viewingConditions.c / viewingConditions.z);
		local p2 = ac / viewingConditions.nbb;
		local gamma = 23.0 * (p2 + 0.305) * t /
				(23.0 * p1 + 11 * t * hCos + 108.0 * t * hSin);
		local a = gamma * hCos;
		local b = gamma * hSin;
		local rA = (460.0 * p2 + 451.0 * a + 288.0 * b) / 1403.0;
		local gA = (460.0 * p2 - 891.0 * a - 261.0 * b) / 1403.0;
		local bA = (460.0 * p2 - 220.0 * a - 6300.0 * b) / 1403.0;
		local rCScaled = HctSolver.inverseChromaticAdaptation(rA);
		local gCScaled = HctSolver.inverseChromaticAdaptation(gA);
		local bCScaled = HctSolver.inverseChromaticAdaptation(bA);
		local linrgb = mathUtils.matrixMultiply(
			{rCScaled, gCScaled, bCScaled},
			HctSolver.LINRGB_FROM_SCALED_DISCOUNT
		);
		-- ===========================================================
		-- Operations inlined from Cam16 to avoid repeated calculation
		-- ===========================================================
		if linrgb[1] < 0 or linrgb[2] < 0 or linrgb[3] < 0 then
			return 0;
		end
		local kR = HctSolver.Y_FROM_LINRGB[1];
		local kG = HctSolver.Y_FROM_LINRGB[2];
		local kB = HctSolver.Y_FROM_LINRGB[3];
		local fnj = kR * linrgb[1] + kG * linrgb[2] + kB * linrgb[3];
		if fnj <= 0 then
			return 0;
		end
		if iterationRound == 5 or math.abs(fnj - y) < 0.002 then
			if linrgb[1] > 100.01 or linrgb[2] > 100.01 or linrgb[3] > 100.01 then
				return 0;
			end
			return colorUtils.argbFromLinrgb(linrgb);
		end
		-- Iterates with Newton method,
		-- Using 2 * fn(j) / j as the approximation of fn'(j)
		j = j - (fnj - y) * j / (2 * fnj);
	end
	return 0;
end

--[[
  Finds an sRGB color with the given hue, chroma, and L*, if possible.

  @param hueDegrees The desired hue, in degrees.
  @param chroma The desired chroma.
  @param lstar The desired L*.
  @return A hexadecimal representing the sRGB color. The color has sufficiently
      close hue, chroma, and L* to the desired values, if possible; otherwise,
      the hue and L* will be sufficiently close, and chroma will be maximized.
]]
function HctSolver.solveToInt(hueDegrees, chroma, lstar)
	if chroma < 0.0001 or lstar < 0.0001 or lstar > 99.9999 then
		return colorUtils.argbFromLstar(lstar);
	end
	hueDegrees = mathUtils.sanitizeDegreesDouble(hueDegrees);
	local hueRadians = math.rad(hueDegrees);
	local y = colorUtils.yFromLstar(lstar);
	local exactAnswer = HctSolver.findResultByJ(hueRadians, chroma, y);
	if exactAnswer ~= 0 then
		return exactAnswer;
	end
	local linrgb = HctSolver.bisectToLimit(y, hueRadians);
	return colorUtils.argbFromLinrgb(linrgb);
end

--[[
  Finds an sRGB color with the given hue, chroma, and L*, if possible.

  @param hueDegrees The desired hue, in degrees.
  @param chroma The desired chroma.
  @param lstar The desired L*.
  @return An CAM16 object representing the sRGB color. The color has
      sufficiently close hue, chroma, and L* to the desired values, if possible;
      otherwise, the hue and L* will be sufficiently close, and chroma will be
      maximized.
]]
function HctSolver.solveToCam(hueDegrees, chroma, lstar)
	return Cam16.fromInt(HctSolver.solveToInt(hueDegrees, chroma, lstar));
end

return HctSolver;