skcms.cc 85 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394
  1. /*
  2. * Copyright 2018 Google Inc.
  3. *
  4. * Use of this source code is governed by a BSD-style license that can be
  5. * found in the LICENSE file.
  6. */
  7. #include "skcms.h"
  8. #include "skcms_internal.h"
  9. #include <assert.h>
  10. #include <float.h>
  11. #include <limits.h>
  12. #include <stdlib.h>
  13. #include <string.h>
  14. #if defined(__ARM_NEON)
  15. #include <arm_neon.h>
  16. #elif defined(__SSE__)
  17. #include <immintrin.h>
  18. #if defined(__clang__)
  19. // That #include <immintrin.h> is usually enough, but Clang's headers
  20. // "helpfully" skip including the whole kitchen sink when _MSC_VER is
  21. // defined, because lots of programs on Windows would include that and
  22. // it'd be a lot slower. But we want all those headers included so we
  23. // can use their features after runtime checks later.
  24. #include <smmintrin.h>
  25. #include <avxintrin.h>
  26. #include <avx2intrin.h>
  27. #include <avx512fintrin.h>
  28. #include <avx512dqintrin.h>
  29. #endif
  30. #endif
  31. // sizeof(x) will return size_t, which is 32-bit on some machines and 64-bit on others.
  32. // We have better testing on 64-bit machines, so force 32-bit machines to behave like 64-bit.
  33. //
  34. // Please do not use sizeof() directly, and size_t only when required.
  35. // (We have no way of enforcing these requests...)
  36. #define SAFE_SIZEOF(x) ((uint64_t)sizeof(x))
  37. // Same sort of thing for _Layout structs with a variable sized array at the end (named "variable").
  38. #define SAFE_FIXED_SIZE(type) ((uint64_t)offsetof(type, variable))
  39. static const union {
  40. uint32_t bits;
  41. float f;
  42. } inf_ = { 0x7f800000 };
  43. #define INFINITY_ inf_.f
  44. static float fmaxf_(float x, float y) { return x > y ? x : y; }
  45. static float fminf_(float x, float y) { return x < y ? x : y; }
  46. static bool isfinitef_(float x) { return 0 == x*0; }
  47. static float minus_1_ulp(float x) {
  48. int32_t bits;
  49. memcpy(&bits, &x, sizeof(bits));
  50. bits = bits - 1;
  51. memcpy(&x, &bits, sizeof(bits));
  52. return x;
  53. }
  54. static float eval_curve(const skcms_Curve* curve, float x) {
  55. if (curve->table_entries == 0) {
  56. return skcms_TransferFunction_eval(&curve->parametric, x);
  57. }
  58. float ix = fmaxf_(0, fminf_(x, 1)) * (curve->table_entries - 1);
  59. int lo = (int) ix ,
  60. hi = (int)(float)minus_1_ulp(ix + 1.0f);
  61. float t = ix - (float)lo;
  62. float l, h;
  63. if (curve->table_8) {
  64. l = curve->table_8[lo] * (1/255.0f);
  65. h = curve->table_8[hi] * (1/255.0f);
  66. } else {
  67. uint16_t be_l, be_h;
  68. memcpy(&be_l, curve->table_16 + 2*lo, 2);
  69. memcpy(&be_h, curve->table_16 + 2*hi, 2);
  70. uint16_t le_l = ((be_l << 8) | (be_l >> 8)) & 0xffff;
  71. uint16_t le_h = ((be_h << 8) | (be_h >> 8)) & 0xffff;
  72. l = le_l * (1/65535.0f);
  73. h = le_h * (1/65535.0f);
  74. }
  75. return l + (h-l)*t;
  76. }
  77. static float max_roundtrip_error(const skcms_Curve* curve, const skcms_TransferFunction* inv_tf) {
  78. uint32_t N = curve->table_entries > 256 ? curve->table_entries : 256;
  79. const float dx = 1.0f / (N - 1);
  80. float err = 0;
  81. for (uint32_t i = 0; i < N; i++) {
  82. float x = i * dx,
  83. y = eval_curve(curve, x);
  84. err = fmaxf_(err, fabsf_(x - skcms_TransferFunction_eval(inv_tf, y)));
  85. }
  86. return err;
  87. }
  88. bool skcms_AreApproximateInverses(const skcms_Curve* curve, const skcms_TransferFunction* inv_tf) {
  89. return max_roundtrip_error(curve, inv_tf) < (1/512.0f);
  90. }
  91. // Additional ICC signature values that are only used internally
  92. enum {
  93. // File signature
  94. skcms_Signature_acsp = 0x61637370,
  95. // Tag signatures
  96. skcms_Signature_rTRC = 0x72545243,
  97. skcms_Signature_gTRC = 0x67545243,
  98. skcms_Signature_bTRC = 0x62545243,
  99. skcms_Signature_kTRC = 0x6B545243,
  100. skcms_Signature_rXYZ = 0x7258595A,
  101. skcms_Signature_gXYZ = 0x6758595A,
  102. skcms_Signature_bXYZ = 0x6258595A,
  103. skcms_Signature_A2B0 = 0x41324230,
  104. skcms_Signature_A2B1 = 0x41324231,
  105. skcms_Signature_mAB = 0x6D414220,
  106. skcms_Signature_CHAD = 0x63686164,
  107. // Type signatures
  108. skcms_Signature_curv = 0x63757276,
  109. skcms_Signature_mft1 = 0x6D667431,
  110. skcms_Signature_mft2 = 0x6D667432,
  111. skcms_Signature_para = 0x70617261,
  112. skcms_Signature_sf32 = 0x73663332,
  113. // XYZ is also a PCS signature, so it's defined in skcms.h
  114. // skcms_Signature_XYZ = 0x58595A20,
  115. };
  116. static uint16_t read_big_u16(const uint8_t* ptr) {
  117. uint16_t be;
  118. memcpy(&be, ptr, sizeof(be));
  119. #if defined(_MSC_VER)
  120. return _byteswap_ushort(be);
  121. #else
  122. return __builtin_bswap16(be);
  123. #endif
  124. }
  125. static uint32_t read_big_u32(const uint8_t* ptr) {
  126. uint32_t be;
  127. memcpy(&be, ptr, sizeof(be));
  128. #if defined(_MSC_VER)
  129. return _byteswap_ulong(be);
  130. #else
  131. return __builtin_bswap32(be);
  132. #endif
  133. }
  134. static int32_t read_big_i32(const uint8_t* ptr) {
  135. return (int32_t)read_big_u32(ptr);
  136. }
  137. static float read_big_fixed(const uint8_t* ptr) {
  138. return read_big_i32(ptr) * (1.0f / 65536.0f);
  139. }
  140. // Maps to an in-memory profile so that fields line up to the locations specified
  141. // in ICC.1:2010, section 7.2
  142. typedef struct {
  143. uint8_t size [ 4];
  144. uint8_t cmm_type [ 4];
  145. uint8_t version [ 4];
  146. uint8_t profile_class [ 4];
  147. uint8_t data_color_space [ 4];
  148. uint8_t pcs [ 4];
  149. uint8_t creation_date_time [12];
  150. uint8_t signature [ 4];
  151. uint8_t platform [ 4];
  152. uint8_t flags [ 4];
  153. uint8_t device_manufacturer [ 4];
  154. uint8_t device_model [ 4];
  155. uint8_t device_attributes [ 8];
  156. uint8_t rendering_intent [ 4];
  157. uint8_t illuminant_X [ 4];
  158. uint8_t illuminant_Y [ 4];
  159. uint8_t illuminant_Z [ 4];
  160. uint8_t creator [ 4];
  161. uint8_t profile_id [16];
  162. uint8_t reserved [28];
  163. uint8_t tag_count [ 4]; // Technically not part of header, but required
  164. } header_Layout;
  165. typedef struct {
  166. uint8_t signature [4];
  167. uint8_t offset [4];
  168. uint8_t size [4];
  169. } tag_Layout;
  170. static const tag_Layout* get_tag_table(const skcms_ICCProfile* profile) {
  171. return (const tag_Layout*)(profile->buffer + SAFE_SIZEOF(header_Layout));
  172. }
  173. // s15Fixed16ArrayType is technically variable sized, holding N values. However, the only valid
  174. // use of the type is for the CHAD tag that stores exactly nine values.
  175. typedef struct {
  176. uint8_t type [ 4];
  177. uint8_t reserved [ 4];
  178. uint8_t values [36];
  179. } sf32_Layout;
  180. bool skcms_GetCHAD(const skcms_ICCProfile* profile, skcms_Matrix3x3* m) {
  181. skcms_ICCTag tag;
  182. if (!skcms_GetTagBySignature(profile, skcms_Signature_CHAD, &tag)) {
  183. return false;
  184. }
  185. if (tag.type != skcms_Signature_sf32 || tag.size < SAFE_SIZEOF(sf32_Layout)) {
  186. return false;
  187. }
  188. const sf32_Layout* sf32Tag = (const sf32_Layout*)tag.buf;
  189. const uint8_t* values = sf32Tag->values;
  190. for (int r = 0; r < 3; ++r)
  191. for (int c = 0; c < 3; ++c, values += 4) {
  192. m->vals[r][c] = read_big_fixed(values);
  193. }
  194. return true;
  195. }
  196. // XYZType is technically variable sized, holding N XYZ triples. However, the only valid uses of
  197. // the type are for tags/data that store exactly one triple.
  198. typedef struct {
  199. uint8_t type [4];
  200. uint8_t reserved [4];
  201. uint8_t X [4];
  202. uint8_t Y [4];
  203. uint8_t Z [4];
  204. } XYZ_Layout;
  205. static bool read_tag_xyz(const skcms_ICCTag* tag, float* x, float* y, float* z) {
  206. if (tag->type != skcms_Signature_XYZ || tag->size < SAFE_SIZEOF(XYZ_Layout)) {
  207. return false;
  208. }
  209. const XYZ_Layout* xyzTag = (const XYZ_Layout*)tag->buf;
  210. *x = read_big_fixed(xyzTag->X);
  211. *y = read_big_fixed(xyzTag->Y);
  212. *z = read_big_fixed(xyzTag->Z);
  213. return true;
  214. }
  215. static bool read_to_XYZD50(const skcms_ICCTag* rXYZ, const skcms_ICCTag* gXYZ,
  216. const skcms_ICCTag* bXYZ, skcms_Matrix3x3* toXYZ) {
  217. return read_tag_xyz(rXYZ, &toXYZ->vals[0][0], &toXYZ->vals[1][0], &toXYZ->vals[2][0]) &&
  218. read_tag_xyz(gXYZ, &toXYZ->vals[0][1], &toXYZ->vals[1][1], &toXYZ->vals[2][1]) &&
  219. read_tag_xyz(bXYZ, &toXYZ->vals[0][2], &toXYZ->vals[1][2], &toXYZ->vals[2][2]);
  220. }
  221. static bool tf_is_valid(const skcms_TransferFunction* tf) {
  222. // Reject obviously malformed inputs
  223. if (!isfinitef_(tf->a + tf->b + tf->c + tf->d + tf->e + tf->f + tf->g)) {
  224. return false;
  225. }
  226. // All of these parameters should be non-negative
  227. if (tf->a < 0 || tf->c < 0 || tf->d < 0 || tf->g < 0) {
  228. return false;
  229. }
  230. return true;
  231. }
  232. typedef struct {
  233. uint8_t type [4];
  234. uint8_t reserved_a [4];
  235. uint8_t function_type [2];
  236. uint8_t reserved_b [2];
  237. uint8_t variable [1/*variable*/]; // 1, 3, 4, 5, or 7 s15.16, depending on function_type
  238. } para_Layout;
  239. static bool read_curve_para(const uint8_t* buf, uint32_t size,
  240. skcms_Curve* curve, uint32_t* curve_size) {
  241. if (size < SAFE_FIXED_SIZE(para_Layout)) {
  242. return false;
  243. }
  244. const para_Layout* paraTag = (const para_Layout*)buf;
  245. enum { kG = 0, kGAB = 1, kGABC = 2, kGABCD = 3, kGABCDEF = 4 };
  246. uint16_t function_type = read_big_u16(paraTag->function_type);
  247. if (function_type > kGABCDEF) {
  248. return false;
  249. }
  250. static const uint32_t curve_bytes[] = { 4, 12, 16, 20, 28 };
  251. if (size < SAFE_FIXED_SIZE(para_Layout) + curve_bytes[function_type]) {
  252. return false;
  253. }
  254. if (curve_size) {
  255. *curve_size = SAFE_FIXED_SIZE(para_Layout) + curve_bytes[function_type];
  256. }
  257. curve->table_entries = 0;
  258. curve->parametric.a = 1.0f;
  259. curve->parametric.b = 0.0f;
  260. curve->parametric.c = 0.0f;
  261. curve->parametric.d = 0.0f;
  262. curve->parametric.e = 0.0f;
  263. curve->parametric.f = 0.0f;
  264. curve->parametric.g = read_big_fixed(paraTag->variable);
  265. switch (function_type) {
  266. case kGAB:
  267. curve->parametric.a = read_big_fixed(paraTag->variable + 4);
  268. curve->parametric.b = read_big_fixed(paraTag->variable + 8);
  269. if (curve->parametric.a == 0) {
  270. return false;
  271. }
  272. curve->parametric.d = -curve->parametric.b / curve->parametric.a;
  273. break;
  274. case kGABC:
  275. curve->parametric.a = read_big_fixed(paraTag->variable + 4);
  276. curve->parametric.b = read_big_fixed(paraTag->variable + 8);
  277. curve->parametric.e = read_big_fixed(paraTag->variable + 12);
  278. if (curve->parametric.a == 0) {
  279. return false;
  280. }
  281. curve->parametric.d = -curve->parametric.b / curve->parametric.a;
  282. curve->parametric.f = curve->parametric.e;
  283. break;
  284. case kGABCD:
  285. curve->parametric.a = read_big_fixed(paraTag->variable + 4);
  286. curve->parametric.b = read_big_fixed(paraTag->variable + 8);
  287. curve->parametric.c = read_big_fixed(paraTag->variable + 12);
  288. curve->parametric.d = read_big_fixed(paraTag->variable + 16);
  289. break;
  290. case kGABCDEF:
  291. curve->parametric.a = read_big_fixed(paraTag->variable + 4);
  292. curve->parametric.b = read_big_fixed(paraTag->variable + 8);
  293. curve->parametric.c = read_big_fixed(paraTag->variable + 12);
  294. curve->parametric.d = read_big_fixed(paraTag->variable + 16);
  295. curve->parametric.e = read_big_fixed(paraTag->variable + 20);
  296. curve->parametric.f = read_big_fixed(paraTag->variable + 24);
  297. break;
  298. }
  299. return tf_is_valid(&curve->parametric);
  300. }
  301. typedef struct {
  302. uint8_t type [4];
  303. uint8_t reserved [4];
  304. uint8_t value_count [4];
  305. uint8_t variable [1/*variable*/]; // value_count, 8.8 if 1, uint16 (n*65535) if > 1
  306. } curv_Layout;
  307. static bool read_curve_curv(const uint8_t* buf, uint32_t size,
  308. skcms_Curve* curve, uint32_t* curve_size) {
  309. if (size < SAFE_FIXED_SIZE(curv_Layout)) {
  310. return false;
  311. }
  312. const curv_Layout* curvTag = (const curv_Layout*)buf;
  313. uint32_t value_count = read_big_u32(curvTag->value_count);
  314. if (size < SAFE_FIXED_SIZE(curv_Layout) + value_count * SAFE_SIZEOF(uint16_t)) {
  315. return false;
  316. }
  317. if (curve_size) {
  318. *curve_size = SAFE_FIXED_SIZE(curv_Layout) + value_count * SAFE_SIZEOF(uint16_t);
  319. }
  320. if (value_count < 2) {
  321. curve->table_entries = 0;
  322. curve->parametric.a = 1.0f;
  323. curve->parametric.b = 0.0f;
  324. curve->parametric.c = 0.0f;
  325. curve->parametric.d = 0.0f;
  326. curve->parametric.e = 0.0f;
  327. curve->parametric.f = 0.0f;
  328. if (value_count == 0) {
  329. // Empty tables are a shorthand for an identity curve
  330. curve->parametric.g = 1.0f;
  331. } else {
  332. // Single entry tables are a shorthand for simple gamma
  333. curve->parametric.g = read_big_u16(curvTag->variable) * (1.0f / 256.0f);
  334. }
  335. } else {
  336. curve->table_8 = nullptr;
  337. curve->table_16 = curvTag->variable;
  338. curve->table_entries = value_count;
  339. }
  340. return true;
  341. }
  342. // Parses both curveType and parametricCurveType data. Ensures that at most 'size' bytes are read.
  343. // If curve_size is not nullptr, writes the number of bytes used by the curve in (*curve_size).
  344. static bool read_curve(const uint8_t* buf, uint32_t size,
  345. skcms_Curve* curve, uint32_t* curve_size) {
  346. if (!buf || size < 4 || !curve) {
  347. return false;
  348. }
  349. uint32_t type = read_big_u32(buf);
  350. if (type == skcms_Signature_para) {
  351. return read_curve_para(buf, size, curve, curve_size);
  352. } else if (type == skcms_Signature_curv) {
  353. return read_curve_curv(buf, size, curve, curve_size);
  354. }
  355. return false;
  356. }
  357. // mft1 and mft2 share a large chunk of data
  358. typedef struct {
  359. uint8_t type [ 4];
  360. uint8_t reserved_a [ 4];
  361. uint8_t input_channels [ 1];
  362. uint8_t output_channels [ 1];
  363. uint8_t grid_points [ 1];
  364. uint8_t reserved_b [ 1];
  365. uint8_t matrix [36];
  366. } mft_CommonLayout;
  367. typedef struct {
  368. mft_CommonLayout common [1];
  369. uint8_t variable [1/*variable*/];
  370. } mft1_Layout;
  371. typedef struct {
  372. mft_CommonLayout common [1];
  373. uint8_t input_table_entries [2];
  374. uint8_t output_table_entries [2];
  375. uint8_t variable [1/*variable*/];
  376. } mft2_Layout;
  377. static bool read_mft_common(const mft_CommonLayout* mftTag, skcms_A2B* a2b) {
  378. // MFT matrices are applied before the first set of curves, but must be identity unless the
  379. // input is PCSXYZ. We don't support PCSXYZ profiles, so we ignore this matrix. Note that the
  380. // matrix in skcms_A2B is applied later in the pipe, so supporting this would require another
  381. // field/flag.
  382. a2b->matrix_channels = 0;
  383. a2b->input_channels = mftTag->input_channels[0];
  384. a2b->output_channels = mftTag->output_channels[0];
  385. // We require exactly three (ie XYZ/Lab/RGB) output channels
  386. if (a2b->output_channels != ARRAY_COUNT(a2b->output_curves)) {
  387. return false;
  388. }
  389. // We require at least one, and no more than four (ie CMYK) input channels
  390. if (a2b->input_channels < 1 || a2b->input_channels > ARRAY_COUNT(a2b->input_curves)) {
  391. return false;
  392. }
  393. for (uint32_t i = 0; i < a2b->input_channels; ++i) {
  394. a2b->grid_points[i] = mftTag->grid_points[0];
  395. }
  396. // The grid only makes sense with at least two points along each axis
  397. if (a2b->grid_points[0] < 2) {
  398. return false;
  399. }
  400. return true;
  401. }
  402. static bool init_a2b_tables(const uint8_t* table_base, uint64_t max_tables_len, uint32_t byte_width,
  403. uint32_t input_table_entries, uint32_t output_table_entries,
  404. skcms_A2B* a2b) {
  405. // byte_width is 1 or 2, [input|output]_table_entries are in [2, 4096], so no overflow
  406. uint32_t byte_len_per_input_table = input_table_entries * byte_width;
  407. uint32_t byte_len_per_output_table = output_table_entries * byte_width;
  408. // [input|output]_channels are <= 4, so still no overflow
  409. uint32_t byte_len_all_input_tables = a2b->input_channels * byte_len_per_input_table;
  410. uint32_t byte_len_all_output_tables = a2b->output_channels * byte_len_per_output_table;
  411. uint64_t grid_size = a2b->output_channels * byte_width;
  412. for (uint32_t axis = 0; axis < a2b->input_channels; ++axis) {
  413. grid_size *= a2b->grid_points[axis];
  414. }
  415. if (max_tables_len < byte_len_all_input_tables + grid_size + byte_len_all_output_tables) {
  416. return false;
  417. }
  418. for (uint32_t i = 0; i < a2b->input_channels; ++i) {
  419. a2b->input_curves[i].table_entries = input_table_entries;
  420. if (byte_width == 1) {
  421. a2b->input_curves[i].table_8 = table_base + i * byte_len_per_input_table;
  422. a2b->input_curves[i].table_16 = nullptr;
  423. } else {
  424. a2b->input_curves[i].table_8 = nullptr;
  425. a2b->input_curves[i].table_16 = table_base + i * byte_len_per_input_table;
  426. }
  427. }
  428. if (byte_width == 1) {
  429. a2b->grid_8 = table_base + byte_len_all_input_tables;
  430. a2b->grid_16 = nullptr;
  431. } else {
  432. a2b->grid_8 = nullptr;
  433. a2b->grid_16 = table_base + byte_len_all_input_tables;
  434. }
  435. const uint8_t* output_table_base = table_base + byte_len_all_input_tables + grid_size;
  436. for (uint32_t i = 0; i < a2b->output_channels; ++i) {
  437. a2b->output_curves[i].table_entries = output_table_entries;
  438. if (byte_width == 1) {
  439. a2b->output_curves[i].table_8 = output_table_base + i * byte_len_per_output_table;
  440. a2b->output_curves[i].table_16 = nullptr;
  441. } else {
  442. a2b->output_curves[i].table_8 = nullptr;
  443. a2b->output_curves[i].table_16 = output_table_base + i * byte_len_per_output_table;
  444. }
  445. }
  446. return true;
  447. }
  448. static bool read_tag_mft1(const skcms_ICCTag* tag, skcms_A2B* a2b) {
  449. if (tag->size < SAFE_FIXED_SIZE(mft1_Layout)) {
  450. return false;
  451. }
  452. const mft1_Layout* mftTag = (const mft1_Layout*)tag->buf;
  453. if (!read_mft_common(mftTag->common, a2b)) {
  454. return false;
  455. }
  456. uint32_t input_table_entries = 256;
  457. uint32_t output_table_entries = 256;
  458. return init_a2b_tables(mftTag->variable, tag->size - SAFE_FIXED_SIZE(mft1_Layout), 1,
  459. input_table_entries, output_table_entries, a2b);
  460. }
  461. static bool read_tag_mft2(const skcms_ICCTag* tag, skcms_A2B* a2b) {
  462. if (tag->size < SAFE_FIXED_SIZE(mft2_Layout)) {
  463. return false;
  464. }
  465. const mft2_Layout* mftTag = (const mft2_Layout*)tag->buf;
  466. if (!read_mft_common(mftTag->common, a2b)) {
  467. return false;
  468. }
  469. uint32_t input_table_entries = read_big_u16(mftTag->input_table_entries);
  470. uint32_t output_table_entries = read_big_u16(mftTag->output_table_entries);
  471. // ICC spec mandates that 2 <= table_entries <= 4096
  472. if (input_table_entries < 2 || input_table_entries > 4096 ||
  473. output_table_entries < 2 || output_table_entries > 4096) {
  474. return false;
  475. }
  476. return init_a2b_tables(mftTag->variable, tag->size - SAFE_FIXED_SIZE(mft2_Layout), 2,
  477. input_table_entries, output_table_entries, a2b);
  478. }
  479. static bool read_curves(const uint8_t* buf, uint32_t size, uint32_t curve_offset,
  480. uint32_t num_curves, skcms_Curve* curves) {
  481. for (uint32_t i = 0; i < num_curves; ++i) {
  482. if (curve_offset > size) {
  483. return false;
  484. }
  485. uint32_t curve_bytes;
  486. if (!read_curve(buf + curve_offset, size - curve_offset, &curves[i], &curve_bytes)) {
  487. return false;
  488. }
  489. if (curve_bytes > UINT32_MAX - 3) {
  490. return false;
  491. }
  492. curve_bytes = (curve_bytes + 3) & ~3U;
  493. uint64_t new_offset_64 = (uint64_t)curve_offset + curve_bytes;
  494. curve_offset = (uint32_t)new_offset_64;
  495. if (new_offset_64 != curve_offset) {
  496. return false;
  497. }
  498. }
  499. return true;
  500. }
  501. typedef struct {
  502. uint8_t type [ 4];
  503. uint8_t reserved_a [ 4];
  504. uint8_t input_channels [ 1];
  505. uint8_t output_channels [ 1];
  506. uint8_t reserved_b [ 2];
  507. uint8_t b_curve_offset [ 4];
  508. uint8_t matrix_offset [ 4];
  509. uint8_t m_curve_offset [ 4];
  510. uint8_t clut_offset [ 4];
  511. uint8_t a_curve_offset [ 4];
  512. } mAB_Layout;
  513. typedef struct {
  514. uint8_t grid_points [16];
  515. uint8_t grid_byte_width [ 1];
  516. uint8_t reserved [ 3];
  517. uint8_t variable [1/*variable*/];
  518. } mABCLUT_Layout;
  519. static bool read_tag_mab(const skcms_ICCTag* tag, skcms_A2B* a2b, bool pcs_is_xyz) {
  520. if (tag->size < SAFE_SIZEOF(mAB_Layout)) {
  521. return false;
  522. }
  523. const mAB_Layout* mABTag = (const mAB_Layout*)tag->buf;
  524. a2b->input_channels = mABTag->input_channels[0];
  525. a2b->output_channels = mABTag->output_channels[0];
  526. // We require exactly three (ie XYZ/Lab/RGB) output channels
  527. if (a2b->output_channels != ARRAY_COUNT(a2b->output_curves)) {
  528. return false;
  529. }
  530. // We require no more than four (ie CMYK) input channels
  531. if (a2b->input_channels > ARRAY_COUNT(a2b->input_curves)) {
  532. return false;
  533. }
  534. uint32_t b_curve_offset = read_big_u32(mABTag->b_curve_offset);
  535. uint32_t matrix_offset = read_big_u32(mABTag->matrix_offset);
  536. uint32_t m_curve_offset = read_big_u32(mABTag->m_curve_offset);
  537. uint32_t clut_offset = read_big_u32(mABTag->clut_offset);
  538. uint32_t a_curve_offset = read_big_u32(mABTag->a_curve_offset);
  539. // "B" curves must be present
  540. if (0 == b_curve_offset) {
  541. return false;
  542. }
  543. if (!read_curves(tag->buf, tag->size, b_curve_offset, a2b->output_channels,
  544. a2b->output_curves)) {
  545. return false;
  546. }
  547. // "M" curves and Matrix must be used together
  548. if (0 != m_curve_offset) {
  549. if (0 == matrix_offset) {
  550. return false;
  551. }
  552. a2b->matrix_channels = a2b->output_channels;
  553. if (!read_curves(tag->buf, tag->size, m_curve_offset, a2b->matrix_channels,
  554. a2b->matrix_curves)) {
  555. return false;
  556. }
  557. // Read matrix, which is stored as a row-major 3x3, followed by the fourth column
  558. if (tag->size < matrix_offset + 12 * SAFE_SIZEOF(uint32_t)) {
  559. return false;
  560. }
  561. float encoding_factor = pcs_is_xyz ? 65535 / 32768.0f : 1.0f;
  562. const uint8_t* mtx_buf = tag->buf + matrix_offset;
  563. a2b->matrix.vals[0][0] = encoding_factor * read_big_fixed(mtx_buf + 0);
  564. a2b->matrix.vals[0][1] = encoding_factor * read_big_fixed(mtx_buf + 4);
  565. a2b->matrix.vals[0][2] = encoding_factor * read_big_fixed(mtx_buf + 8);
  566. a2b->matrix.vals[1][0] = encoding_factor * read_big_fixed(mtx_buf + 12);
  567. a2b->matrix.vals[1][1] = encoding_factor * read_big_fixed(mtx_buf + 16);
  568. a2b->matrix.vals[1][2] = encoding_factor * read_big_fixed(mtx_buf + 20);
  569. a2b->matrix.vals[2][0] = encoding_factor * read_big_fixed(mtx_buf + 24);
  570. a2b->matrix.vals[2][1] = encoding_factor * read_big_fixed(mtx_buf + 28);
  571. a2b->matrix.vals[2][2] = encoding_factor * read_big_fixed(mtx_buf + 32);
  572. a2b->matrix.vals[0][3] = encoding_factor * read_big_fixed(mtx_buf + 36);
  573. a2b->matrix.vals[1][3] = encoding_factor * read_big_fixed(mtx_buf + 40);
  574. a2b->matrix.vals[2][3] = encoding_factor * read_big_fixed(mtx_buf + 44);
  575. } else {
  576. if (0 != matrix_offset) {
  577. return false;
  578. }
  579. a2b->matrix_channels = 0;
  580. }
  581. // "A" curves and CLUT must be used together
  582. if (0 != a_curve_offset) {
  583. if (0 == clut_offset) {
  584. return false;
  585. }
  586. if (!read_curves(tag->buf, tag->size, a_curve_offset, a2b->input_channels,
  587. a2b->input_curves)) {
  588. return false;
  589. }
  590. if (tag->size < clut_offset + SAFE_FIXED_SIZE(mABCLUT_Layout)) {
  591. return false;
  592. }
  593. const mABCLUT_Layout* clut = (const mABCLUT_Layout*)(tag->buf + clut_offset);
  594. if (clut->grid_byte_width[0] == 1) {
  595. a2b->grid_8 = clut->variable;
  596. a2b->grid_16 = nullptr;
  597. } else if (clut->grid_byte_width[0] == 2) {
  598. a2b->grid_8 = nullptr;
  599. a2b->grid_16 = clut->variable;
  600. } else {
  601. return false;
  602. }
  603. uint64_t grid_size = a2b->output_channels * clut->grid_byte_width[0];
  604. for (uint32_t i = 0; i < a2b->input_channels; ++i) {
  605. a2b->grid_points[i] = clut->grid_points[i];
  606. // The grid only makes sense with at least two points along each axis
  607. if (a2b->grid_points[i] < 2) {
  608. return false;
  609. }
  610. grid_size *= a2b->grid_points[i];
  611. }
  612. if (tag->size < clut_offset + SAFE_FIXED_SIZE(mABCLUT_Layout) + grid_size) {
  613. return false;
  614. }
  615. } else {
  616. if (0 != clut_offset) {
  617. return false;
  618. }
  619. // If there is no CLUT, the number of input and output channels must match
  620. if (a2b->input_channels != a2b->output_channels) {
  621. return false;
  622. }
  623. // Zero out the number of input channels to signal that we're skipping this stage
  624. a2b->input_channels = 0;
  625. }
  626. return true;
  627. }
  628. // If you pass f, we'll fit a possibly-non-zero value for *f.
  629. // If you pass nullptr, we'll assume you want *f to be treated as zero.
  630. static int fit_linear(const skcms_Curve* curve, int N, float tol,
  631. float* c, float* d, float* f = nullptr) {
  632. assert(N > 1);
  633. // We iteratively fit the first points to the TF's linear piece.
  634. // We want the cx + f line to pass through the first and last points we fit exactly.
  635. //
  636. // As we walk along the points we find the minimum and maximum slope of the line before the
  637. // error would exceed our tolerance. We stop when the range [slope_min, slope_max] becomes
  638. // emtpy, when we definitely can't add any more points.
  639. //
  640. // Some points' error intervals may intersect the running interval but not lie fully
  641. // within it. So we keep track of the last point we saw that is a valid end point candidate,
  642. // and once the search is done, back up to build the line through *that* point.
  643. const float dx = 1.0f / (N - 1);
  644. int lin_points = 1;
  645. float f_zero = 0.0f;
  646. if (f) {
  647. *f = eval_curve(curve, 0);
  648. } else {
  649. f = &f_zero;
  650. }
  651. float slope_min = -INFINITY_;
  652. float slope_max = +INFINITY_;
  653. for (int i = 1; i < N; ++i) {
  654. float x = i * dx;
  655. float y = eval_curve(curve, x);
  656. float slope_max_i = (y + tol - *f) / x,
  657. slope_min_i = (y - tol - *f) / x;
  658. if (slope_max_i < slope_min || slope_max < slope_min_i) {
  659. // Slope intervals would no longer overlap.
  660. break;
  661. }
  662. slope_max = fminf_(slope_max, slope_max_i);
  663. slope_min = fmaxf_(slope_min, slope_min_i);
  664. float cur_slope = (y - *f) / x;
  665. if (slope_min <= cur_slope && cur_slope <= slope_max) {
  666. lin_points = i + 1;
  667. *c = cur_slope;
  668. }
  669. }
  670. // Set D to the last point that met our tolerance.
  671. *d = (lin_points - 1) * dx;
  672. return lin_points;
  673. }
  674. static bool read_a2b(const skcms_ICCTag* tag, skcms_A2B* a2b, bool pcs_is_xyz) {
  675. bool ok = false;
  676. if (tag->type == skcms_Signature_mft1) {
  677. ok = read_tag_mft1(tag, a2b);
  678. } else if (tag->type == skcms_Signature_mft2) {
  679. ok = read_tag_mft2(tag, a2b);
  680. } else if (tag->type == skcms_Signature_mAB) {
  681. ok = read_tag_mab(tag, a2b, pcs_is_xyz);
  682. }
  683. if (!ok) {
  684. return false;
  685. }
  686. // Detect and canonicalize identity tables.
  687. skcms_Curve* curves[] = {
  688. a2b->input_channels > 0 ? a2b->input_curves + 0 : nullptr,
  689. a2b->input_channels > 1 ? a2b->input_curves + 1 : nullptr,
  690. a2b->input_channels > 2 ? a2b->input_curves + 2 : nullptr,
  691. a2b->input_channels > 3 ? a2b->input_curves + 3 : nullptr,
  692. a2b->matrix_channels > 0 ? a2b->matrix_curves + 0 : nullptr,
  693. a2b->matrix_channels > 1 ? a2b->matrix_curves + 1 : nullptr,
  694. a2b->matrix_channels > 2 ? a2b->matrix_curves + 2 : nullptr,
  695. a2b->output_channels > 0 ? a2b->output_curves + 0 : nullptr,
  696. a2b->output_channels > 1 ? a2b->output_curves + 1 : nullptr,
  697. a2b->output_channels > 2 ? a2b->output_curves + 2 : nullptr,
  698. };
  699. for (int i = 0; i < ARRAY_COUNT(curves); i++) {
  700. skcms_Curve* curve = curves[i];
  701. if (curve && curve->table_entries && curve->table_entries <= (uint32_t)INT_MAX) {
  702. int N = (int)curve->table_entries;
  703. float c = 0.0f, d = 0.0f, f = 0.0f;
  704. if (N == fit_linear(curve, N, 1.0f/(2*N), &c,&d,&f)
  705. && c == 1.0f
  706. && f == 0.0f) {
  707. curve->table_entries = 0;
  708. curve->table_8 = nullptr;
  709. curve->table_16 = nullptr;
  710. curve->parametric = skcms_TransferFunction{1,1,0,0,0,0,0};
  711. }
  712. }
  713. }
  714. return true;
  715. }
  716. void skcms_GetTagByIndex(const skcms_ICCProfile* profile, uint32_t idx, skcms_ICCTag* tag) {
  717. if (!profile || !profile->buffer || !tag) { return; }
  718. if (idx > profile->tag_count) { return; }
  719. const tag_Layout* tags = get_tag_table(profile);
  720. tag->signature = read_big_u32(tags[idx].signature);
  721. tag->size = read_big_u32(tags[idx].size);
  722. tag->buf = read_big_u32(tags[idx].offset) + profile->buffer;
  723. tag->type = read_big_u32(tag->buf);
  724. }
  725. bool skcms_GetTagBySignature(const skcms_ICCProfile* profile, uint32_t sig, skcms_ICCTag* tag) {
  726. if (!profile || !profile->buffer || !tag) { return false; }
  727. const tag_Layout* tags = get_tag_table(profile);
  728. for (uint32_t i = 0; i < profile->tag_count; ++i) {
  729. if (read_big_u32(tags[i].signature) == sig) {
  730. tag->signature = sig;
  731. tag->size = read_big_u32(tags[i].size);
  732. tag->buf = read_big_u32(tags[i].offset) + profile->buffer;
  733. tag->type = read_big_u32(tag->buf);
  734. return true;
  735. }
  736. }
  737. return false;
  738. }
  739. static bool usable_as_src(const skcms_ICCProfile* profile) {
  740. return profile->has_A2B
  741. || (profile->has_trc && profile->has_toXYZD50);
  742. }
  743. bool skcms_Parse(const void* buf, size_t len, skcms_ICCProfile* profile) {
  744. assert(SAFE_SIZEOF(header_Layout) == 132);
  745. if (!profile) {
  746. return false;
  747. }
  748. memset(profile, 0, SAFE_SIZEOF(*profile));
  749. if (len < SAFE_SIZEOF(header_Layout)) {
  750. return false;
  751. }
  752. // Byte-swap all header fields
  753. const header_Layout* header = (const header_Layout*)buf;
  754. profile->buffer = (const uint8_t*)buf;
  755. profile->size = read_big_u32(header->size);
  756. uint32_t version = read_big_u32(header->version);
  757. profile->data_color_space = read_big_u32(header->data_color_space);
  758. profile->pcs = read_big_u32(header->pcs);
  759. uint32_t signature = read_big_u32(header->signature);
  760. float illuminant_X = read_big_fixed(header->illuminant_X);
  761. float illuminant_Y = read_big_fixed(header->illuminant_Y);
  762. float illuminant_Z = read_big_fixed(header->illuminant_Z);
  763. profile->tag_count = read_big_u32(header->tag_count);
  764. // Validate signature, size (smaller than buffer, large enough to hold tag table),
  765. // and major version
  766. uint64_t tag_table_size = profile->tag_count * SAFE_SIZEOF(tag_Layout);
  767. if (signature != skcms_Signature_acsp ||
  768. profile->size > len ||
  769. profile->size < SAFE_SIZEOF(header_Layout) + tag_table_size ||
  770. (version >> 24) > 4) {
  771. return false;
  772. }
  773. // Validate that illuminant is D50 white
  774. if (fabsf_(illuminant_X - 0.9642f) > 0.0100f ||
  775. fabsf_(illuminant_Y - 1.0000f) > 0.0100f ||
  776. fabsf_(illuminant_Z - 0.8249f) > 0.0100f) {
  777. return false;
  778. }
  779. // Validate that all tag entries have sane offset + size
  780. const tag_Layout* tags = get_tag_table(profile);
  781. for (uint32_t i = 0; i < profile->tag_count; ++i) {
  782. uint32_t tag_offset = read_big_u32(tags[i].offset);
  783. uint32_t tag_size = read_big_u32(tags[i].size);
  784. uint64_t tag_end = (uint64_t)tag_offset + (uint64_t)tag_size;
  785. if (tag_size < 4 || tag_end > profile->size) {
  786. return false;
  787. }
  788. }
  789. if (profile->pcs != skcms_Signature_XYZ && profile->pcs != skcms_Signature_Lab) {
  790. return false;
  791. }
  792. bool pcs_is_xyz = profile->pcs == skcms_Signature_XYZ;
  793. // Pre-parse commonly used tags.
  794. skcms_ICCTag kTRC;
  795. if (profile->data_color_space == skcms_Signature_Gray &&
  796. skcms_GetTagBySignature(profile, skcms_Signature_kTRC, &kTRC)) {
  797. if (!read_curve(kTRC.buf, kTRC.size, &profile->trc[0], nullptr)) {
  798. // Malformed tag
  799. return false;
  800. }
  801. profile->trc[1] = profile->trc[0];
  802. profile->trc[2] = profile->trc[0];
  803. profile->has_trc = true;
  804. if (pcs_is_xyz) {
  805. profile->toXYZD50.vals[0][0] = illuminant_X;
  806. profile->toXYZD50.vals[1][1] = illuminant_Y;
  807. profile->toXYZD50.vals[2][2] = illuminant_Z;
  808. profile->has_toXYZD50 = true;
  809. }
  810. } else {
  811. skcms_ICCTag rTRC, gTRC, bTRC;
  812. if (skcms_GetTagBySignature(profile, skcms_Signature_rTRC, &rTRC) &&
  813. skcms_GetTagBySignature(profile, skcms_Signature_gTRC, &gTRC) &&
  814. skcms_GetTagBySignature(profile, skcms_Signature_bTRC, &bTRC)) {
  815. if (!read_curve(rTRC.buf, rTRC.size, &profile->trc[0], nullptr) ||
  816. !read_curve(gTRC.buf, gTRC.size, &profile->trc[1], nullptr) ||
  817. !read_curve(bTRC.buf, bTRC.size, &profile->trc[2], nullptr)) {
  818. // Malformed TRC tags
  819. return false;
  820. }
  821. profile->has_trc = true;
  822. }
  823. skcms_ICCTag rXYZ, gXYZ, bXYZ;
  824. if (skcms_GetTagBySignature(profile, skcms_Signature_rXYZ, &rXYZ) &&
  825. skcms_GetTagBySignature(profile, skcms_Signature_gXYZ, &gXYZ) &&
  826. skcms_GetTagBySignature(profile, skcms_Signature_bXYZ, &bXYZ)) {
  827. if (!read_to_XYZD50(&rXYZ, &gXYZ, &bXYZ, &profile->toXYZD50)) {
  828. // Malformed XYZ tags
  829. return false;
  830. }
  831. profile->has_toXYZD50 = true;
  832. }
  833. }
  834. skcms_ICCTag a2b_tag;
  835. // For now, we're preferring A2B0, like Skia does and the ICC spec tells us to.
  836. // TODO: prefer A2B1 (relative colormetric) over A2B0 (perceptual)?
  837. // This breaks with the ICC spec, but we think it's a good idea, given that TRC curves
  838. // and all our known users are thinking exclusively in terms of relative colormetric.
  839. const uint32_t sigs[] = { skcms_Signature_A2B0, skcms_Signature_A2B1 };
  840. for (int i = 0; i < ARRAY_COUNT(sigs); i++) {
  841. if (skcms_GetTagBySignature(profile, sigs[i], &a2b_tag)) {
  842. if (!read_a2b(&a2b_tag, &profile->A2B, pcs_is_xyz)) {
  843. // Malformed A2B tag
  844. return false;
  845. }
  846. profile->has_A2B = true;
  847. break;
  848. }
  849. }
  850. return usable_as_src(profile);
  851. }
  852. const skcms_ICCProfile* skcms_sRGB_profile() {
  853. static const skcms_ICCProfile sRGB_profile = {
  854. nullptr, // buffer, moot here
  855. 0, // size, moot here
  856. skcms_Signature_RGB, // data_color_space
  857. skcms_Signature_XYZ, // pcs
  858. 0, // tag count, moot here
  859. // We choose to represent sRGB with its canonical transfer function,
  860. // and with its canonical XYZD50 gamut matrix.
  861. true, // has_trc, followed by the 3 trc curves
  862. {
  863. {{0, {2.4f, (float)(1/1.055), (float)(0.055/1.055), (float)(1/12.92), 0.04045f, 0, 0}}},
  864. {{0, {2.4f, (float)(1/1.055), (float)(0.055/1.055), (float)(1/12.92), 0.04045f, 0, 0}}},
  865. {{0, {2.4f, (float)(1/1.055), (float)(0.055/1.055), (float)(1/12.92), 0.04045f, 0, 0}}},
  866. },
  867. true, // has_toXYZD50, followed by 3x3 toXYZD50 matrix
  868. {{
  869. { 0.436065674f, 0.385147095f, 0.143066406f },
  870. { 0.222488403f, 0.716873169f, 0.060607910f },
  871. { 0.013916016f, 0.097076416f, 0.714096069f },
  872. }},
  873. false, // has_A2B, followed by a2b itself which we don't care about.
  874. {
  875. 0,
  876. {
  877. {{0, {0,0, 0,0,0,0,0}}},
  878. {{0, {0,0, 0,0,0,0,0}}},
  879. {{0, {0,0, 0,0,0,0,0}}},
  880. {{0, {0,0, 0,0,0,0,0}}},
  881. },
  882. {0,0,0,0},
  883. nullptr,
  884. nullptr,
  885. 0,
  886. {
  887. {{0, {0,0, 0,0,0,0,0}}},
  888. {{0, {0,0, 0,0,0,0,0}}},
  889. {{0, {0,0, 0,0,0,0,0}}},
  890. },
  891. {{
  892. { 0,0,0,0 },
  893. { 0,0,0,0 },
  894. { 0,0,0,0 },
  895. }},
  896. 0,
  897. {
  898. {{0, {0,0, 0,0,0,0,0}}},
  899. {{0, {0,0, 0,0,0,0,0}}},
  900. {{0, {0,0, 0,0,0,0,0}}},
  901. },
  902. },
  903. };
  904. return &sRGB_profile;
  905. }
  906. const skcms_ICCProfile* skcms_XYZD50_profile() {
  907. // Just like sRGB above, but with identity transfer functions and toXYZD50 matrix.
  908. static const skcms_ICCProfile XYZD50_profile = {
  909. nullptr, // buffer, moot here
  910. 0, // size, moot here
  911. skcms_Signature_RGB, // data_color_space
  912. skcms_Signature_XYZ, // pcs
  913. 0, // tag count, moot here
  914. true, // has_trc, followed by the 3 trc curves
  915. {
  916. {{0, {1,1, 0,0,0,0,0}}},
  917. {{0, {1,1, 0,0,0,0,0}}},
  918. {{0, {1,1, 0,0,0,0,0}}},
  919. },
  920. true, // has_toXYZD50, followed by 3x3 toXYZD50 matrix
  921. {{
  922. { 1,0,0 },
  923. { 0,1,0 },
  924. { 0,0,1 },
  925. }},
  926. false, // has_A2B, followed by a2b itself which we don't care about.
  927. {
  928. 0,
  929. {
  930. {{0, {0,0, 0,0,0,0,0}}},
  931. {{0, {0,0, 0,0,0,0,0}}},
  932. {{0, {0,0, 0,0,0,0,0}}},
  933. {{0, {0,0, 0,0,0,0,0}}},
  934. },
  935. {0,0,0,0},
  936. nullptr,
  937. nullptr,
  938. 0,
  939. {
  940. {{0, {0,0, 0,0,0,0,0}}},
  941. {{0, {0,0, 0,0,0,0,0}}},
  942. {{0, {0,0, 0,0,0,0,0}}},
  943. },
  944. {{
  945. { 0,0,0,0 },
  946. { 0,0,0,0 },
  947. { 0,0,0,0 },
  948. }},
  949. 0,
  950. {
  951. {{0, {0,0, 0,0,0,0,0}}},
  952. {{0, {0,0, 0,0,0,0,0}}},
  953. {{0, {0,0, 0,0,0,0,0}}},
  954. },
  955. },
  956. };
  957. return &XYZD50_profile;
  958. }
  959. const skcms_TransferFunction* skcms_sRGB_TransferFunction() {
  960. return &skcms_sRGB_profile()->trc[0].parametric;
  961. }
  962. const skcms_TransferFunction* skcms_sRGB_Inverse_TransferFunction() {
  963. static const skcms_TransferFunction sRGB_inv =
  964. {0.416666657f, 1.137283325f, -0.0f, 12.920000076f, 0.003130805f, -0.054969788f, -0.0f};
  965. return &sRGB_inv;
  966. }
  967. const skcms_TransferFunction* skcms_Identity_TransferFunction() {
  968. static const skcms_TransferFunction identity = {1,1,0,0,0,0,0};
  969. return &identity;
  970. }
  971. const uint8_t skcms_252_random_bytes[] = {
  972. 8, 179, 128, 204, 253, 38, 134, 184, 68, 102, 32, 138, 99, 39, 169, 215,
  973. 119, 26, 3, 223, 95, 239, 52, 132, 114, 74, 81, 234, 97, 116, 244, 205, 30,
  974. 154, 173, 12, 51, 159, 122, 153, 61, 226, 236, 178, 229, 55, 181, 220, 191,
  975. 194, 160, 126, 168, 82, 131, 18, 180, 245, 163, 22, 246, 69, 235, 252, 57,
  976. 108, 14, 6, 152, 240, 255, 171, 242, 20, 227, 177, 238, 96, 85, 16, 211,
  977. 70, 200, 149, 155, 146, 127, 145, 100, 151, 109, 19, 165, 208, 195, 164,
  978. 137, 254, 182, 248, 64, 201, 45, 209, 5, 147, 207, 210, 113, 162, 83, 225,
  979. 9, 31, 15, 231, 115, 37, 58, 53, 24, 49, 197, 56, 120, 172, 48, 21, 214,
  980. 129, 111, 11, 50, 187, 196, 34, 60, 103, 71, 144, 47, 203, 77, 80, 232,
  981. 140, 222, 250, 206, 166, 247, 139, 249, 221, 72, 106, 27, 199, 117, 54,
  982. 219, 135, 118, 40, 79, 41, 251, 46, 93, 212, 92, 233, 148, 28, 121, 63,
  983. 123, 158, 105, 59, 29, 42, 143, 23, 0, 107, 176, 87, 104, 183, 156, 193,
  984. 189, 90, 188, 65, 190, 17, 198, 7, 186, 161, 1, 124, 78, 125, 170, 133,
  985. 174, 218, 67, 157, 75, 101, 89, 217, 62, 33, 141, 228, 25, 35, 91, 230, 4,
  986. 2, 13, 73, 86, 167, 237, 84, 243, 44, 185, 66, 130, 110, 150, 142, 216, 88,
  987. 112, 36, 224, 136, 202, 76, 94, 98, 175, 213
  988. };
  989. bool skcms_ApproximatelyEqualProfiles(const skcms_ICCProfile* A, const skcms_ICCProfile* B) {
  990. // Test for exactly equal profiles first.
  991. if (A == B || 0 == memcmp(A,B, sizeof(skcms_ICCProfile))) {
  992. return true;
  993. }
  994. // For now this is the essentially the same strategy we use in test_only.c
  995. // for our skcms_Transform() smoke tests:
  996. // 1) transform A to XYZD50
  997. // 2) transform B to XYZD50
  998. // 3) return true if they're similar enough
  999. // Our current criterion in 3) is maximum 1 bit error per XYZD50 byte.
  1000. // skcms_252_random_bytes are 252 of a random shuffle of all possible bytes.
  1001. // 252 is evenly divisible by 3 and 4. Only 192, 10, 241, and 43 are missing.
  1002. if (A->data_color_space != B->data_color_space) {
  1003. return false;
  1004. }
  1005. // Interpret as RGB_888 if data color space is RGB or GRAY, RGBA_8888 if CMYK.
  1006. // TODO: working with RGBA_8888 either way is probably fastest.
  1007. skcms_PixelFormat fmt = skcms_PixelFormat_RGB_888;
  1008. size_t npixels = 84;
  1009. if (A->data_color_space == skcms_Signature_CMYK) {
  1010. fmt = skcms_PixelFormat_RGBA_8888;
  1011. npixels = 63;
  1012. }
  1013. // TODO: if A or B is a known profile (skcms_sRGB_profile, skcms_XYZD50_profile),
  1014. // use pre-canned results and skip that skcms_Transform() call?
  1015. uint8_t dstA[252],
  1016. dstB[252];
  1017. if (!skcms_Transform(
  1018. skcms_252_random_bytes, fmt, skcms_AlphaFormat_Unpremul, A,
  1019. dstA, skcms_PixelFormat_RGB_888, skcms_AlphaFormat_Unpremul, skcms_XYZD50_profile(),
  1020. npixels)) {
  1021. return false;
  1022. }
  1023. if (!skcms_Transform(
  1024. skcms_252_random_bytes, fmt, skcms_AlphaFormat_Unpremul, B,
  1025. dstB, skcms_PixelFormat_RGB_888, skcms_AlphaFormat_Unpremul, skcms_XYZD50_profile(),
  1026. npixels)) {
  1027. return false;
  1028. }
  1029. // TODO: make sure this final check has reasonable codegen.
  1030. for (size_t i = 0; i < 252; i++) {
  1031. if (abs((int)dstA[i] - (int)dstB[i]) > 1) {
  1032. return false;
  1033. }
  1034. }
  1035. return true;
  1036. }
  1037. bool skcms_TRCs_AreApproximateInverse(const skcms_ICCProfile* profile,
  1038. const skcms_TransferFunction* inv_tf) {
  1039. if (!profile || !profile->has_trc) {
  1040. return false;
  1041. }
  1042. return skcms_AreApproximateInverses(&profile->trc[0], inv_tf) &&
  1043. skcms_AreApproximateInverses(&profile->trc[1], inv_tf) &&
  1044. skcms_AreApproximateInverses(&profile->trc[2], inv_tf);
  1045. }
  1046. static bool is_zero_to_one(float x) {
  1047. return 0 <= x && x <= 1;
  1048. }
  1049. typedef struct { float vals[3]; } skcms_Vector3;
  1050. static skcms_Vector3 mv_mul(const skcms_Matrix3x3* m, const skcms_Vector3* v) {
  1051. skcms_Vector3 dst = {{0,0,0}};
  1052. for (int row = 0; row < 3; ++row) {
  1053. dst.vals[row] = m->vals[row][0] * v->vals[0]
  1054. + m->vals[row][1] * v->vals[1]
  1055. + m->vals[row][2] * v->vals[2];
  1056. }
  1057. return dst;
  1058. }
  1059. bool skcms_PrimariesToXYZD50(float rx, float ry,
  1060. float gx, float gy,
  1061. float bx, float by,
  1062. float wx, float wy,
  1063. skcms_Matrix3x3* toXYZD50) {
  1064. if (!is_zero_to_one(rx) || !is_zero_to_one(ry) ||
  1065. !is_zero_to_one(gx) || !is_zero_to_one(gy) ||
  1066. !is_zero_to_one(bx) || !is_zero_to_one(by) ||
  1067. !is_zero_to_one(wx) || !is_zero_to_one(wy) ||
  1068. !toXYZD50) {
  1069. return false;
  1070. }
  1071. // First, we need to convert xy values (primaries) to XYZ.
  1072. skcms_Matrix3x3 primaries = {{
  1073. { rx, gx, bx },
  1074. { ry, gy, by },
  1075. { 1 - rx - ry, 1 - gx - gy, 1 - bx - by },
  1076. }};
  1077. skcms_Matrix3x3 primaries_inv;
  1078. if (!skcms_Matrix3x3_invert(&primaries, &primaries_inv)) {
  1079. return false;
  1080. }
  1081. // Assumes that Y is 1.0f.
  1082. skcms_Vector3 wXYZ = { { wx / wy, 1, (1 - wx - wy) / wy } };
  1083. skcms_Vector3 XYZ = mv_mul(&primaries_inv, &wXYZ);
  1084. skcms_Matrix3x3 toXYZ = {{
  1085. { XYZ.vals[0], 0, 0 },
  1086. { 0, XYZ.vals[1], 0 },
  1087. { 0, 0, XYZ.vals[2] },
  1088. }};
  1089. toXYZ = skcms_Matrix3x3_concat(&primaries, &toXYZ);
  1090. // Now convert toXYZ matrix to toXYZD50.
  1091. skcms_Vector3 wXYZD50 = { { 0.96422f, 1.0f, 0.82521f } };
  1092. // Calculate the chromatic adaptation matrix. We will use the Bradford method, thus
  1093. // the matrices below. The Bradford method is used by Adobe and is widely considered
  1094. // to be the best.
  1095. skcms_Matrix3x3 xyz_to_lms = {{
  1096. { 0.8951f, 0.2664f, -0.1614f },
  1097. { -0.7502f, 1.7135f, 0.0367f },
  1098. { 0.0389f, -0.0685f, 1.0296f },
  1099. }};
  1100. skcms_Matrix3x3 lms_to_xyz = {{
  1101. { 0.9869929f, -0.1470543f, 0.1599627f },
  1102. { 0.4323053f, 0.5183603f, 0.0492912f },
  1103. { -0.0085287f, 0.0400428f, 0.9684867f },
  1104. }};
  1105. skcms_Vector3 srcCone = mv_mul(&xyz_to_lms, &wXYZ);
  1106. skcms_Vector3 dstCone = mv_mul(&xyz_to_lms, &wXYZD50);
  1107. skcms_Matrix3x3 DXtoD50 = {{
  1108. { dstCone.vals[0] / srcCone.vals[0], 0, 0 },
  1109. { 0, dstCone.vals[1] / srcCone.vals[1], 0 },
  1110. { 0, 0, dstCone.vals[2] / srcCone.vals[2] },
  1111. }};
  1112. DXtoD50 = skcms_Matrix3x3_concat(&DXtoD50, &xyz_to_lms);
  1113. DXtoD50 = skcms_Matrix3x3_concat(&lms_to_xyz, &DXtoD50);
  1114. *toXYZD50 = skcms_Matrix3x3_concat(&DXtoD50, &toXYZ);
  1115. return true;
  1116. }
  1117. bool skcms_Matrix3x3_invert(const skcms_Matrix3x3* src, skcms_Matrix3x3* dst) {
  1118. double a00 = src->vals[0][0],
  1119. a01 = src->vals[1][0],
  1120. a02 = src->vals[2][0],
  1121. a10 = src->vals[0][1],
  1122. a11 = src->vals[1][1],
  1123. a12 = src->vals[2][1],
  1124. a20 = src->vals[0][2],
  1125. a21 = src->vals[1][2],
  1126. a22 = src->vals[2][2];
  1127. double b0 = a00*a11 - a01*a10,
  1128. b1 = a00*a12 - a02*a10,
  1129. b2 = a01*a12 - a02*a11,
  1130. b3 = a20,
  1131. b4 = a21,
  1132. b5 = a22;
  1133. double determinant = b0*b5
  1134. - b1*b4
  1135. + b2*b3;
  1136. if (determinant == 0) {
  1137. return false;
  1138. }
  1139. double invdet = 1.0 / determinant;
  1140. if (invdet > +FLT_MAX || invdet < -FLT_MAX || !isfinitef_((float)invdet)) {
  1141. return false;
  1142. }
  1143. b0 *= invdet;
  1144. b1 *= invdet;
  1145. b2 *= invdet;
  1146. b3 *= invdet;
  1147. b4 *= invdet;
  1148. b5 *= invdet;
  1149. dst->vals[0][0] = (float)( a11*b5 - a12*b4 );
  1150. dst->vals[1][0] = (float)( a02*b4 - a01*b5 );
  1151. dst->vals[2][0] = (float)( + b2 );
  1152. dst->vals[0][1] = (float)( a12*b3 - a10*b5 );
  1153. dst->vals[1][1] = (float)( a00*b5 - a02*b3 );
  1154. dst->vals[2][1] = (float)( - b1 );
  1155. dst->vals[0][2] = (float)( a10*b4 - a11*b3 );
  1156. dst->vals[1][2] = (float)( a01*b3 - a00*b4 );
  1157. dst->vals[2][2] = (float)( + b0 );
  1158. for (int r = 0; r < 3; ++r)
  1159. for (int c = 0; c < 3; ++c) {
  1160. if (!isfinitef_(dst->vals[r][c])) {
  1161. return false;
  1162. }
  1163. }
  1164. return true;
  1165. }
  1166. skcms_Matrix3x3 skcms_Matrix3x3_concat(const skcms_Matrix3x3* A, const skcms_Matrix3x3* B) {
  1167. skcms_Matrix3x3 m = { { { 0,0,0 },{ 0,0,0 },{ 0,0,0 } } };
  1168. for (int r = 0; r < 3; r++)
  1169. for (int c = 0; c < 3; c++) {
  1170. m.vals[r][c] = A->vals[r][0] * B->vals[0][c]
  1171. + A->vals[r][1] * B->vals[1][c]
  1172. + A->vals[r][2] * B->vals[2][c];
  1173. }
  1174. return m;
  1175. }
  1176. #if defined(__clang__) || defined(__GNUC__)
  1177. #define small_memcpy __builtin_memcpy
  1178. #else
  1179. #define small_memcpy memcpy
  1180. #endif
  1181. static float log2f_(float x) {
  1182. // The first approximation of log2(x) is its exponent 'e', minus 127.
  1183. int32_t bits;
  1184. small_memcpy(&bits, &x, sizeof(bits));
  1185. float e = (float)bits * (1.0f / (1<<23));
  1186. // If we use the mantissa too we can refine the error signficantly.
  1187. int32_t m_bits = (bits & 0x007fffff) | 0x3f000000;
  1188. float m;
  1189. small_memcpy(&m, &m_bits, sizeof(m));
  1190. return (e - 124.225514990f
  1191. - 1.498030302f*m
  1192. - 1.725879990f/(0.3520887068f + m));
  1193. }
  1194. static float exp2f_(float x) {
  1195. float fract = x - floorf_(x);
  1196. float fbits = (1.0f * (1<<23)) * (x + 121.274057500f
  1197. - 1.490129070f*fract
  1198. + 27.728023300f/(4.84252568f - fract));
  1199. if (fbits > INT_MAX) {
  1200. return INFINITY_;
  1201. } else if (fbits < INT_MIN) {
  1202. return -INFINITY_;
  1203. }
  1204. int32_t bits = (int32_t)fbits;
  1205. small_memcpy(&x, &bits, sizeof(x));
  1206. return x;
  1207. }
  1208. float powf_(float x, float y) {
  1209. return (x == 0) || (x == 1) ? x
  1210. : exp2f_(log2f_(x) * y);
  1211. }
  1212. float skcms_TransferFunction_eval(const skcms_TransferFunction* tf, float x) {
  1213. float sign = x < 0 ? -1.0f : 1.0f;
  1214. x *= sign;
  1215. return sign * (x < tf->d ? tf->c * x + tf->f
  1216. : powf_(tf->a * x + tf->b, tf->g) + tf->e);
  1217. }
  1218. #if defined(__clang__)
  1219. [[clang::no_sanitize("float-divide-by-zero")]] // Checked for by tf_is_valid() on the way out.
  1220. #endif
  1221. bool skcms_TransferFunction_invert(const skcms_TransferFunction* src, skcms_TransferFunction* dst) {
  1222. if (!tf_is_valid(src)) {
  1223. return false;
  1224. }
  1225. // We're inverting this function, solving for x in terms of y.
  1226. // y = (cx + f) x < d
  1227. // (ax + b)^g + e x ≥ d
  1228. // The inverse of this function can be expressed in the same piecewise form.
  1229. skcms_TransferFunction inv = {0,0,0,0,0,0,0};
  1230. // We'll start by finding the new threshold inv.d.
  1231. // In principle we should be able to find that by solving for y at x=d from either side.
  1232. // (If those two d values aren't the same, it's a discontinuous transfer function.)
  1233. float d_l = src->c * src->d + src->f,
  1234. d_r = powf_(src->a * src->d + src->b, src->g) + src->e;
  1235. if (fabsf_(d_l - d_r) > 1/512.0f) {
  1236. return false;
  1237. }
  1238. inv.d = d_l; // TODO(mtklein): better in practice to choose d_r?
  1239. // When d=0, the linear section collapses to a point. We leave c,d,f all zero in that case.
  1240. if (inv.d > 0) {
  1241. // Inverting the linear section is pretty straightfoward:
  1242. // y = cx + f
  1243. // y - f = cx
  1244. // (1/c)y - f/c = x
  1245. inv.c = 1.0f/src->c;
  1246. inv.f = -src->f/src->c;
  1247. }
  1248. // The interesting part is inverting the nonlinear section:
  1249. // y = (ax + b)^g + e.
  1250. // y - e = (ax + b)^g
  1251. // (y - e)^1/g = ax + b
  1252. // (y - e)^1/g - b = ax
  1253. // (1/a)(y - e)^1/g - b/a = x
  1254. //
  1255. // To make that fit our form, we need to move the (1/a) term inside the exponentiation:
  1256. // let k = (1/a)^g
  1257. // (1/a)( y - e)^1/g - b/a = x
  1258. // (ky - ke)^1/g - b/a = x
  1259. float k = powf_(src->a, -src->g); // (1/a)^g == a^-g
  1260. inv.g = 1.0f / src->g;
  1261. inv.a = k;
  1262. inv.b = -k * src->e;
  1263. inv.e = -src->b / src->a;
  1264. // Now in principle we're done.
  1265. // But to preserve the valuable invariant inv(src(1.0f)) == 1.0f,
  1266. // we'll tweak e. These two values should be close to each other,
  1267. // just down to numerical precision issues, especially from powf_.
  1268. float s = powf_(src->a + src->b, src->g) + src->e;
  1269. inv.e = 1.0f - powf_(inv.a * s + inv.b, inv.g);
  1270. *dst = inv;
  1271. return tf_is_valid(dst);
  1272. }
  1273. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ //
  1274. // From here below we're approximating an skcms_Curve with an skcms_TransferFunction{g,a,b,c,d,e,f}:
  1275. //
  1276. // tf(x) = cx + f x < d
  1277. // tf(x) = (ax + b)^g + e x ≥ d
  1278. //
  1279. // When fitting, we add the additional constraint that both pieces meet at d:
  1280. //
  1281. // cd + f = (ad + b)^g + e
  1282. //
  1283. // Solving for e and folding it through gives an alternate formulation of the non-linear piece:
  1284. //
  1285. // tf(x) = cx + f x < d
  1286. // tf(x) = (ax + b)^g - (ad + b)^g + cd + f x ≥ d
  1287. //
  1288. // Our overall strategy is then:
  1289. // For a couple tolerances,
  1290. // - fit_linear(): fit c,d,f iteratively to as many points as our tolerance allows
  1291. // - invert c,d,f
  1292. // - fit_nonlinear(): fit g,a,b using Gauss-Newton given those inverted c,d,f
  1293. // (and by constraint, inverted e) to the inverse of the table.
  1294. // Return the parameters with least maximum error.
  1295. //
  1296. // To run Gauss-Newton to find g,a,b, we'll also need the gradient of the residuals
  1297. // of round-trip f_inv(x), the inverse of the non-linear piece of f(x).
  1298. //
  1299. // let y = Table(x)
  1300. // r(x) = x - f_inv(y)
  1301. //
  1302. // ∂r/∂g = ln(ay + b)*(ay + b)^g
  1303. // - ln(ad + b)*(ad + b)^g
  1304. // ∂r/∂a = yg(ay + b)^(g-1)
  1305. // - dg(ad + b)^(g-1)
  1306. // ∂r/∂b = g(ay + b)^(g-1)
  1307. // - g(ad + b)^(g-1)
  1308. // Return the residual of roundtripping skcms_Curve(x) through f_inv(y) with parameters P,
  1309. // and fill out the gradient of the residual into dfdP.
  1310. static float rg_nonlinear(float x,
  1311. const skcms_Curve* curve,
  1312. const skcms_TransferFunction* tf,
  1313. const float P[3],
  1314. float dfdP[3]) {
  1315. const float y = eval_curve(curve, x);
  1316. const float g = P[0], a = P[1], b = P[2],
  1317. c = tf->c, d = tf->d, f = tf->f;
  1318. const float Y = fmaxf_(a*y + b, 0.0f),
  1319. D = a*d + b;
  1320. assert (D >= 0);
  1321. // The gradient.
  1322. dfdP[0] = 0.69314718f*log2f_(Y)*powf_(Y, g)
  1323. - 0.69314718f*log2f_(D)*powf_(D, g);
  1324. dfdP[1] = y*g*powf_(Y, g-1)
  1325. - d*g*powf_(D, g-1);
  1326. dfdP[2] = g*powf_(Y, g-1)
  1327. - g*powf_(D, g-1);
  1328. // The residual.
  1329. const float f_inv = powf_(Y, g)
  1330. - powf_(D, g)
  1331. + c*d + f;
  1332. return x - f_inv;
  1333. }
  1334. static bool gauss_newton_step(const skcms_Curve* curve,
  1335. const skcms_TransferFunction* tf,
  1336. float P[3],
  1337. float x0, float dx, int N) {
  1338. // We'll sample x from the range [x0,x1] (both inclusive) N times with even spacing.
  1339. //
  1340. // We want to do P' = P + (Jf^T Jf)^-1 Jf^T r(P),
  1341. // where r(P) is the residual vector
  1342. // and Jf is the Jacobian matrix of f(), ∂r/∂P.
  1343. //
  1344. // Let's review the shape of each of these expressions:
  1345. // r(P) is [N x 1], a column vector with one entry per value of x tested
  1346. // Jf is [N x 3], a matrix with an entry for each (x,P) pair
  1347. // Jf^T is [3 x N], the transpose of Jf
  1348. //
  1349. // Jf^T Jf is [3 x N] * [N x 3] == [3 x 3], a 3x3 matrix,
  1350. // and so is its inverse (Jf^T Jf)^-1
  1351. // Jf^T r(P) is [3 x N] * [N x 1] == [3 x 1], a column vector with the same shape as P
  1352. //
  1353. // Our implementation strategy to get to the final ∆P is
  1354. // 1) evaluate Jf^T Jf, call that lhs
  1355. // 2) evaluate Jf^T r(P), call that rhs
  1356. // 3) invert lhs
  1357. // 4) multiply inverse lhs by rhs
  1358. //
  1359. // This is a friendly implementation strategy because we don't have to have any
  1360. // buffers that scale with N, and equally nice don't have to perform any matrix
  1361. // operations that are variable size.
  1362. //
  1363. // Other implementation strategies could trade this off, e.g. evaluating the
  1364. // pseudoinverse of Jf ( (Jf^T Jf)^-1 Jf^T ) directly, then multiplying that by
  1365. // the residuals. That would probably require implementing singular value
  1366. // decomposition, and would create a [3 x N] matrix to be multiplied by the
  1367. // [N x 1] residual vector, but on the upside I think that'd eliminate the
  1368. // possibility of this gauss_newton_step() function ever failing.
  1369. // 0) start off with lhs and rhs safely zeroed.
  1370. skcms_Matrix3x3 lhs = {{ {0,0,0}, {0,0,0}, {0,0,0} }};
  1371. skcms_Vector3 rhs = { {0,0,0} };
  1372. // 1,2) evaluate lhs and evaluate rhs
  1373. // We want to evaluate Jf only once, but both lhs and rhs involve Jf^T,
  1374. // so we'll have to update lhs and rhs at the same time.
  1375. for (int i = 0; i < N; i++) {
  1376. float x = x0 + i*dx;
  1377. float dfdP[3] = {0,0,0};
  1378. float resid = rg_nonlinear(x,curve,tf,P, dfdP);
  1379. for (int r = 0; r < 3; r++) {
  1380. for (int c = 0; c < 3; c++) {
  1381. lhs.vals[r][c] += dfdP[r] * dfdP[c];
  1382. }
  1383. rhs.vals[r] += dfdP[r] * resid;
  1384. }
  1385. }
  1386. // If any of the 3 P parameters are unused, this matrix will be singular.
  1387. // Detect those cases and fix them up to indentity instead, so we can invert.
  1388. for (int k = 0; k < 3; k++) {
  1389. if (lhs.vals[0][k]==0 && lhs.vals[1][k]==0 && lhs.vals[2][k]==0 &&
  1390. lhs.vals[k][0]==0 && lhs.vals[k][1]==0 && lhs.vals[k][2]==0) {
  1391. lhs.vals[k][k] = 1;
  1392. }
  1393. }
  1394. // 3) invert lhs
  1395. skcms_Matrix3x3 lhs_inv;
  1396. if (!skcms_Matrix3x3_invert(&lhs, &lhs_inv)) {
  1397. return false;
  1398. }
  1399. // 4) multiply inverse lhs by rhs
  1400. skcms_Vector3 dP = mv_mul(&lhs_inv, &rhs);
  1401. P[0] += dP.vals[0];
  1402. P[1] += dP.vals[1];
  1403. P[2] += dP.vals[2];
  1404. return isfinitef_(P[0]) && isfinitef_(P[1]) && isfinitef_(P[2]);
  1405. }
  1406. // Fit the points in [L,N) to the non-linear piece of tf, or return false if we can't.
  1407. static bool fit_nonlinear(const skcms_Curve* curve, int L, int N, skcms_TransferFunction* tf) {
  1408. float P[3] = { tf->g, tf->a, tf->b };
  1409. // No matter where we start, dx should always represent N even steps from 0 to 1.
  1410. const float dx = 1.0f / (N-1);
  1411. for (int j = 0; j < 3/*TODO: tune*/; j++) {
  1412. // These extra constraints a >= 0 and ad+b >= 0 are not modeled in the optimization.
  1413. // We don't really know how to fix up a if it goes negative.
  1414. if (P[1] < 0) {
  1415. return false;
  1416. }
  1417. // If ad+b goes negative, we feel just barely not uneasy enough to tweak b so ad+b is zero.
  1418. if (P[1] * tf->d + P[2] < 0) {
  1419. P[2] = -P[1] * tf->d;
  1420. }
  1421. assert (P[1] >= 0 &&
  1422. P[1] * tf->d + P[2] >= 0);
  1423. if (!gauss_newton_step(curve, tf,
  1424. P,
  1425. L*dx, dx, N-L)) {
  1426. return false;
  1427. }
  1428. }
  1429. // We need to apply our fixups one last time
  1430. if (P[1] < 0) {
  1431. return false;
  1432. }
  1433. if (P[1] * tf->d + P[2] < 0) {
  1434. P[2] = -P[1] * tf->d;
  1435. }
  1436. tf->g = P[0];
  1437. tf->a = P[1];
  1438. tf->b = P[2];
  1439. tf->e = tf->c*tf->d + tf->f
  1440. - powf_(tf->a*tf->d + tf->b, tf->g);
  1441. return true;
  1442. }
  1443. bool skcms_ApproximateCurve(const skcms_Curve* curve,
  1444. skcms_TransferFunction* approx,
  1445. float* max_error) {
  1446. if (!curve || !approx || !max_error) {
  1447. return false;
  1448. }
  1449. if (curve->table_entries == 0) {
  1450. // No point approximating an skcms_TransferFunction with an skcms_TransferFunction!
  1451. return false;
  1452. }
  1453. if (curve->table_entries == 1 || curve->table_entries > (uint32_t)INT_MAX) {
  1454. // We need at least two points, and must put some reasonable cap on the maximum number.
  1455. return false;
  1456. }
  1457. int N = (int)curve->table_entries;
  1458. const float dx = 1.0f / (N - 1);
  1459. *max_error = INFINITY_;
  1460. const float kTolerances[] = { 1.5f / 65535.0f, 1.0f / 512.0f };
  1461. for (int t = 0; t < ARRAY_COUNT(kTolerances); t++) {
  1462. skcms_TransferFunction tf,
  1463. tf_inv;
  1464. // It's problematic to fit curves with non-zero f, so always force it to zero explicitly.
  1465. tf.f = 0.0f;
  1466. int L = fit_linear(curve, N, kTolerances[t], &tf.c, &tf.d);
  1467. if (L == N) {
  1468. // If the entire data set was linear, move the coefficients to the nonlinear portion
  1469. // with G == 1. This lets use a canonical representation with d == 0.
  1470. tf.g = 1;
  1471. tf.a = tf.c;
  1472. tf.b = tf.f;
  1473. tf.c = tf.d = tf.e = tf.f = 0;
  1474. } else if (L == N - 1) {
  1475. // Degenerate case with only two points in the nonlinear segment. Solve directly.
  1476. tf.g = 1;
  1477. tf.a = (eval_curve(curve, (N-1)*dx) -
  1478. eval_curve(curve, (N-2)*dx))
  1479. / dx;
  1480. tf.b = eval_curve(curve, (N-2)*dx)
  1481. - tf.a * (N-2)*dx;
  1482. tf.e = 0;
  1483. } else {
  1484. // Start by guessing a gamma-only curve through the midpoint.
  1485. int mid = (L + N) / 2;
  1486. float mid_x = mid / (N - 1.0f);
  1487. float mid_y = eval_curve(curve, mid_x);
  1488. tf.g = log2f_(mid_y) / log2f_(mid_x);
  1489. tf.a = 1;
  1490. tf.b = 0;
  1491. tf.e = tf.c*tf.d + tf.f
  1492. - powf_(tf.a*tf.d + tf.b, tf.g);
  1493. if (!skcms_TransferFunction_invert(&tf, &tf_inv) ||
  1494. !fit_nonlinear(curve, L,N, &tf_inv)) {
  1495. continue;
  1496. }
  1497. // We fit tf_inv, so calculate tf to keep in sync.
  1498. if (!skcms_TransferFunction_invert(&tf_inv, &tf)) {
  1499. continue;
  1500. }
  1501. }
  1502. // We find our error by roundtripping the table through tf_inv.
  1503. //
  1504. // (The most likely use case for this approximation is to be inverted and
  1505. // used as the transfer function for a destination color space.)
  1506. //
  1507. // We've kept tf and tf_inv in sync above, but we can't guarantee that tf is
  1508. // invertible, so re-verify that here (and use the new inverse for testing).
  1509. if (!skcms_TransferFunction_invert(&tf, &tf_inv)) {
  1510. continue;
  1511. }
  1512. float err = max_roundtrip_error(curve, &tf_inv);
  1513. if (*max_error > err) {
  1514. *max_error = err;
  1515. *approx = tf;
  1516. }
  1517. }
  1518. return isfinitef_(*max_error);
  1519. }
  1520. // ~~~~ Impl. of skcms_Transform() ~~~~
  1521. typedef enum {
  1522. Op_load_a8,
  1523. Op_load_g8,
  1524. Op_load_8888_palette8,
  1525. Op_load_4444,
  1526. Op_load_565,
  1527. Op_load_888,
  1528. Op_load_8888,
  1529. Op_load_1010102,
  1530. Op_load_161616LE,
  1531. Op_load_16161616LE,
  1532. Op_load_161616BE,
  1533. Op_load_16161616BE,
  1534. Op_load_hhh,
  1535. Op_load_hhhh,
  1536. Op_load_fff,
  1537. Op_load_ffff,
  1538. Op_swap_rb,
  1539. Op_clamp,
  1540. Op_invert,
  1541. Op_force_opaque,
  1542. Op_premul,
  1543. Op_unpremul,
  1544. Op_matrix_3x3,
  1545. Op_matrix_3x4,
  1546. Op_lab_to_xyz,
  1547. Op_tf_r,
  1548. Op_tf_g,
  1549. Op_tf_b,
  1550. Op_tf_a,
  1551. Op_table_r,
  1552. Op_table_g,
  1553. Op_table_b,
  1554. Op_table_a,
  1555. Op_clut,
  1556. Op_store_a8,
  1557. Op_store_g8,
  1558. Op_store_4444,
  1559. Op_store_565,
  1560. Op_store_888,
  1561. Op_store_8888,
  1562. Op_store_1010102,
  1563. Op_store_161616LE,
  1564. Op_store_16161616LE,
  1565. Op_store_161616BE,
  1566. Op_store_16161616BE,
  1567. Op_store_hhh,
  1568. Op_store_hhhh,
  1569. Op_store_fff,
  1570. Op_store_ffff,
  1571. } Op;
  1572. #if defined(__clang__)
  1573. template <int N, typename T> using Vec = T __attribute__((ext_vector_type(N)));
  1574. #elif defined(__GNUC__)
  1575. // For some reason GCC accepts this nonsense, but not the more straightforward version,
  1576. // template <int N, typename T> using Vec = T __attribute__((vector_size(N*sizeof(T))));
  1577. template <int N, typename T>
  1578. struct VecHelper { typedef T __attribute__((vector_size(N*sizeof(T)))) V; };
  1579. template <int N, typename T> using Vec = typename VecHelper<N,T>::V;
  1580. #endif
  1581. // First, instantiate our default exec_ops() implementation using the default compiliation target.
  1582. namespace baseline {
  1583. #if defined(SKCMS_PORTABLE) || !(defined(__clang__) || defined(__GNUC__)) \
  1584. || (defined(__EMSCRIPTEN_major__) && !defined(__wasm_simd128__))
  1585. #define N 1
  1586. template <typename T> using V = T;
  1587. using Color = float;
  1588. #elif defined(__AVX512F__)
  1589. #define N 16
  1590. template <typename T> using V = Vec<N,T>;
  1591. using Color = float;
  1592. #elif defined(__AVX__)
  1593. #define N 8
  1594. template <typename T> using V = Vec<N,T>;
  1595. using Color = float;
  1596. #elif defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && defined(SKCMS_OPT_INTO_NEON_FP16)
  1597. #define N 8
  1598. template <typename T> using V = Vec<N,T>;
  1599. using Color = _Float16;
  1600. #else
  1601. #define N 4
  1602. template <typename T> using V = Vec<N,T>;
  1603. using Color = float;
  1604. #endif
  1605. #include "src/Transform_inl.h"
  1606. #undef N
  1607. }
  1608. // Now, instantiate any other versions of run_program() we may want for runtime detection.
  1609. #if !defined(SKCMS_PORTABLE) && \
  1610. !defined(SKCMS_NO_RUNTIME_CPU_DETECTION) && \
  1611. (( defined(__clang__) && __clang_major__ >= 5) || \
  1612. (!defined(__clang__) && defined(__GNUC__))) \
  1613. && defined(__x86_64__)
  1614. #if !defined(__AVX2__)
  1615. #if defined(__clang__)
  1616. #pragma clang attribute push(__attribute__((target("avx2,f16c"))), apply_to=function)
  1617. #elif defined(__GNUC__)
  1618. #pragma GCC push_options
  1619. #pragma GCC target("avx2,f16c")
  1620. #endif
  1621. namespace hsw {
  1622. #define USING_AVX
  1623. #define USING_AVX_F16C
  1624. #define USING_AVX2
  1625. #define N 8
  1626. template <typename T> using V = Vec<N,T>;
  1627. using Color = float;
  1628. #include "src/Transform_inl.h"
  1629. // src/Transform_inl.h will undefine USING_* for us.
  1630. #undef N
  1631. }
  1632. #if defined(__clang__)
  1633. #pragma clang attribute pop
  1634. #elif defined(__GNUC__)
  1635. #pragma GCC pop_options
  1636. #endif
  1637. #define TEST_FOR_HSW
  1638. #endif
  1639. #if !defined(__AVX512F__)
  1640. #if defined(__clang__)
  1641. #pragma clang attribute push(__attribute__((target("avx512f,avx512dq,avx512cd,avx512bw,avx512vl"))), apply_to=function)
  1642. #elif defined(__GNUC__)
  1643. #pragma GCC push_options
  1644. #pragma GCC target("avx512f,avx512dq,avx512cd,avx512bw,avx512vl")
  1645. #endif
  1646. namespace skx {
  1647. #define USING_AVX512F
  1648. #define N 16
  1649. template <typename T> using V = Vec<N,T>;
  1650. using Color = float;
  1651. #include "src/Transform_inl.h"
  1652. // src/Transform_inl.h will undefine USING_* for us.
  1653. #undef N
  1654. }
  1655. #if defined(__clang__)
  1656. #pragma clang attribute pop
  1657. #elif defined(__GNUC__)
  1658. #pragma GCC pop_options
  1659. #endif
  1660. #define TEST_FOR_SKX
  1661. #endif
  1662. #if defined(TEST_FOR_HSW) || defined(TEST_FOR_SKX)
  1663. enum class CpuType { None, HSW, SKX };
  1664. static CpuType cpu_type() {
  1665. static const CpuType type = []{
  1666. // See http://www.sandpile.org/x86/cpuid.htm
  1667. // First, a basic cpuid(1) lets us check prerequisites for HSW, SKX.
  1668. uint32_t eax, ebx, ecx, edx;
  1669. __asm__ __volatile__("cpuid" : "=a"(eax), "=b"(ebx), "=c"(ecx), "=d"(edx)
  1670. : "0"(1), "2"(0));
  1671. if ((edx & (1u<<25)) && // SSE
  1672. (edx & (1u<<26)) && // SSE2
  1673. (ecx & (1u<< 0)) && // SSE3
  1674. (ecx & (1u<< 9)) && // SSSE3
  1675. (ecx & (1u<<12)) && // FMA (N.B. not used, avoided even)
  1676. (ecx & (1u<<19)) && // SSE4.1
  1677. (ecx & (1u<<20)) && // SSE4.2
  1678. (ecx & (1u<<26)) && // XSAVE
  1679. (ecx & (1u<<27)) && // OSXSAVE
  1680. (ecx & (1u<<28)) && // AVX
  1681. (ecx & (1u<<29))) { // F16C
  1682. // Call cpuid(7) to check for AVX2 and AVX-512 bits.
  1683. __asm__ __volatile__("cpuid" : "=a"(eax), "=b"(ebx), "=c"(ecx), "=d"(edx)
  1684. : "0"(7), "2"(0));
  1685. // eax from xgetbv(0) will tell us whether XMM, YMM, and ZMM state is saved.
  1686. uint32_t xcr0, dont_need_edx;
  1687. __asm__ __volatile__("xgetbv" : "=a"(xcr0), "=d"(dont_need_edx) : "c"(0));
  1688. if ((xcr0 & (1u<<1)) && // XMM register state saved?
  1689. (xcr0 & (1u<<2)) && // YMM register state saved?
  1690. (ebx & (1u<<5))) { // AVX2
  1691. // At this point we're at least HSW. Continue checking for SKX.
  1692. if ((xcr0 & (1u<< 5)) && // Opmasks state saved?
  1693. (xcr0 & (1u<< 6)) && // First 16 ZMM registers saved?
  1694. (xcr0 & (1u<< 7)) && // High 16 ZMM registers saved?
  1695. (ebx & (1u<<16)) && // AVX512F
  1696. (ebx & (1u<<17)) && // AVX512DQ
  1697. (ebx & (1u<<28)) && // AVX512CD
  1698. (ebx & (1u<<30)) && // AVX512BW
  1699. (ebx & (1u<<31))) { // AVX512VL
  1700. return CpuType::SKX;
  1701. }
  1702. return CpuType::HSW;
  1703. }
  1704. }
  1705. return CpuType::None;
  1706. }();
  1707. return type;
  1708. }
  1709. #endif
  1710. #endif
  1711. static bool is_identity_tf(const skcms_TransferFunction* tf) {
  1712. return tf->g == 1 && tf->a == 1
  1713. && tf->b == 0 && tf->c == 0 && tf->d == 0 && tf->e == 0 && tf->f == 0;
  1714. }
  1715. typedef struct {
  1716. Op op;
  1717. const void* arg;
  1718. } OpAndArg;
  1719. static OpAndArg select_curve_op(const skcms_Curve* curve, int channel) {
  1720. static const struct { Op parametric, table; } ops[] = {
  1721. { Op_tf_r, Op_table_r },
  1722. { Op_tf_g, Op_table_g },
  1723. { Op_tf_b, Op_table_b },
  1724. { Op_tf_a, Op_table_a },
  1725. };
  1726. const OpAndArg noop = { Op_load_a8/*doesn't matter*/, nullptr };
  1727. if (curve->table_entries == 0) {
  1728. return is_identity_tf(&curve->parametric)
  1729. ? noop
  1730. : OpAndArg{ ops[channel].parametric, &curve->parametric };
  1731. }
  1732. return OpAndArg{ ops[channel].table, curve };
  1733. }
  1734. static size_t bytes_per_pixel(skcms_PixelFormat fmt) {
  1735. switch (fmt >> 1) { // ignore rgb/bgr
  1736. case skcms_PixelFormat_A_8 >> 1: return 1;
  1737. case skcms_PixelFormat_G_8 >> 1: return 1;
  1738. case skcms_PixelFormat_RGBA_8888_Palette8 >> 1: return 1;
  1739. case skcms_PixelFormat_ABGR_4444 >> 1: return 2;
  1740. case skcms_PixelFormat_RGB_565 >> 1: return 2;
  1741. case skcms_PixelFormat_RGB_888 >> 1: return 3;
  1742. case skcms_PixelFormat_RGBA_8888 >> 1: return 4;
  1743. case skcms_PixelFormat_RGBA_1010102 >> 1: return 4;
  1744. case skcms_PixelFormat_RGB_161616LE >> 1: return 6;
  1745. case skcms_PixelFormat_RGBA_16161616LE >> 1: return 8;
  1746. case skcms_PixelFormat_RGB_161616BE >> 1: return 6;
  1747. case skcms_PixelFormat_RGBA_16161616BE >> 1: return 8;
  1748. case skcms_PixelFormat_RGB_hhh_Norm >> 1: return 6;
  1749. case skcms_PixelFormat_RGBA_hhhh_Norm >> 1: return 8;
  1750. case skcms_PixelFormat_RGB_hhh >> 1: return 6;
  1751. case skcms_PixelFormat_RGBA_hhhh >> 1: return 8;
  1752. case skcms_PixelFormat_RGB_fff >> 1: return 12;
  1753. case skcms_PixelFormat_RGBA_ffff >> 1: return 16;
  1754. }
  1755. assert(false);
  1756. return 0;
  1757. }
  1758. static bool prep_for_destination(const skcms_ICCProfile* profile,
  1759. skcms_Matrix3x3* fromXYZD50,
  1760. skcms_TransferFunction* invR,
  1761. skcms_TransferFunction* invG,
  1762. skcms_TransferFunction* invB) {
  1763. // We only support destinations with parametric transfer functions
  1764. // and with gamuts that can be transformed from XYZD50.
  1765. return profile->has_trc
  1766. && profile->has_toXYZD50
  1767. && profile->trc[0].table_entries == 0
  1768. && profile->trc[1].table_entries == 0
  1769. && profile->trc[2].table_entries == 0
  1770. && skcms_TransferFunction_invert(&profile->trc[0].parametric, invR)
  1771. && skcms_TransferFunction_invert(&profile->trc[1].parametric, invG)
  1772. && skcms_TransferFunction_invert(&profile->trc[2].parametric, invB)
  1773. && skcms_Matrix3x3_invert(&profile->toXYZD50, fromXYZD50);
  1774. }
  1775. bool skcms_Transform(const void* src,
  1776. skcms_PixelFormat srcFmt,
  1777. skcms_AlphaFormat srcAlpha,
  1778. const skcms_ICCProfile* srcProfile,
  1779. void* dst,
  1780. skcms_PixelFormat dstFmt,
  1781. skcms_AlphaFormat dstAlpha,
  1782. const skcms_ICCProfile* dstProfile,
  1783. size_t npixels) {
  1784. return skcms_TransformWithPalette(src, srcFmt, srcAlpha, srcProfile,
  1785. dst, dstFmt, dstAlpha, dstProfile,
  1786. npixels, nullptr);
  1787. }
  1788. bool skcms_TransformWithPalette(const void* src,
  1789. skcms_PixelFormat srcFmt,
  1790. skcms_AlphaFormat srcAlpha,
  1791. const skcms_ICCProfile* srcProfile,
  1792. void* dst,
  1793. skcms_PixelFormat dstFmt,
  1794. skcms_AlphaFormat dstAlpha,
  1795. const skcms_ICCProfile* dstProfile,
  1796. size_t nz,
  1797. const void* palette) {
  1798. const size_t dst_bpp = bytes_per_pixel(dstFmt),
  1799. src_bpp = bytes_per_pixel(srcFmt);
  1800. // Let's just refuse if the request is absurdly big.
  1801. if (nz * dst_bpp > INT_MAX || nz * src_bpp > INT_MAX) {
  1802. return false;
  1803. }
  1804. int n = (int)nz;
  1805. // Null profiles default to sRGB. Passing null for both is handy when doing format conversion.
  1806. if (!srcProfile) {
  1807. srcProfile = skcms_sRGB_profile();
  1808. }
  1809. if (!dstProfile) {
  1810. dstProfile = skcms_sRGB_profile();
  1811. }
  1812. // We can't transform in place unless the PixelFormats are the same size.
  1813. if (dst == src && dst_bpp != src_bpp) {
  1814. return false;
  1815. }
  1816. // TODO: more careful alias rejection (like, dst == src + 1)?
  1817. if (needs_palette(srcFmt) && !palette) {
  1818. return false;
  1819. }
  1820. Op program [32];
  1821. const void* arguments[32];
  1822. Op* ops = program;
  1823. const void** args = arguments;
  1824. skcms_TransferFunction inv_dst_tf_r, inv_dst_tf_g, inv_dst_tf_b;
  1825. skcms_Matrix3x3 from_xyz;
  1826. switch (srcFmt >> 1) {
  1827. default: return false;
  1828. case skcms_PixelFormat_A_8 >> 1: *ops++ = Op_load_a8; break;
  1829. case skcms_PixelFormat_G_8 >> 1: *ops++ = Op_load_g8; break;
  1830. case skcms_PixelFormat_ABGR_4444 >> 1: *ops++ = Op_load_4444; break;
  1831. case skcms_PixelFormat_RGB_565 >> 1: *ops++ = Op_load_565; break;
  1832. case skcms_PixelFormat_RGB_888 >> 1: *ops++ = Op_load_888; break;
  1833. case skcms_PixelFormat_RGBA_8888 >> 1: *ops++ = Op_load_8888; break;
  1834. case skcms_PixelFormat_RGBA_1010102 >> 1: *ops++ = Op_load_1010102; break;
  1835. case skcms_PixelFormat_RGB_161616LE >> 1: *ops++ = Op_load_161616LE; break;
  1836. case skcms_PixelFormat_RGBA_16161616LE >> 1: *ops++ = Op_load_16161616LE; break;
  1837. case skcms_PixelFormat_RGB_161616BE >> 1: *ops++ = Op_load_161616BE; break;
  1838. case skcms_PixelFormat_RGBA_16161616BE >> 1: *ops++ = Op_load_16161616BE; break;
  1839. case skcms_PixelFormat_RGB_hhh_Norm >> 1: *ops++ = Op_load_hhh; break;
  1840. case skcms_PixelFormat_RGBA_hhhh_Norm >> 1: *ops++ = Op_load_hhhh; break;
  1841. case skcms_PixelFormat_RGB_hhh >> 1: *ops++ = Op_load_hhh; break;
  1842. case skcms_PixelFormat_RGBA_hhhh >> 1: *ops++ = Op_load_hhhh; break;
  1843. case skcms_PixelFormat_RGB_fff >> 1: *ops++ = Op_load_fff; break;
  1844. case skcms_PixelFormat_RGBA_ffff >> 1: *ops++ = Op_load_ffff; break;
  1845. case skcms_PixelFormat_RGBA_8888_Palette8 >> 1: *ops++ = Op_load_8888_palette8;
  1846. *args++ = palette;
  1847. break;
  1848. }
  1849. if (srcFmt == skcms_PixelFormat_RGB_hhh_Norm ||
  1850. srcFmt == skcms_PixelFormat_RGBA_hhhh_Norm) {
  1851. *ops++ = Op_clamp;
  1852. }
  1853. if (srcFmt & 1) {
  1854. *ops++ = Op_swap_rb;
  1855. }
  1856. skcms_ICCProfile gray_dst_profile;
  1857. if ((dstFmt >> 1) == (skcms_PixelFormat_G_8 >> 1)) {
  1858. // When transforming to gray, stop at XYZ (by setting toXYZ to identity), then transform
  1859. // luminance (Y) by the destination transfer function.
  1860. gray_dst_profile = *dstProfile;
  1861. skcms_SetXYZD50(&gray_dst_profile, &skcms_XYZD50_profile()->toXYZD50);
  1862. dstProfile = &gray_dst_profile;
  1863. }
  1864. if (srcProfile->data_color_space == skcms_Signature_CMYK) {
  1865. // Photoshop creates CMYK images as inverse CMYK.
  1866. // These happen to be the only ones we've _ever_ seen.
  1867. *ops++ = Op_invert;
  1868. // With CMYK, ignore the alpha type, to avoid changing K or conflating CMY with K.
  1869. srcAlpha = skcms_AlphaFormat_Unpremul;
  1870. }
  1871. if (srcAlpha == skcms_AlphaFormat_Opaque) {
  1872. *ops++ = Op_force_opaque;
  1873. } else if (srcAlpha == skcms_AlphaFormat_PremulAsEncoded) {
  1874. *ops++ = Op_unpremul;
  1875. }
  1876. if (dstProfile != srcProfile) {
  1877. if (!prep_for_destination(dstProfile,
  1878. &from_xyz, &inv_dst_tf_r, &inv_dst_tf_b, &inv_dst_tf_g)) {
  1879. return false;
  1880. }
  1881. if (srcProfile->has_A2B) {
  1882. if (srcProfile->A2B.input_channels) {
  1883. for (int i = 0; i < (int)srcProfile->A2B.input_channels; i++) {
  1884. OpAndArg oa = select_curve_op(&srcProfile->A2B.input_curves[i], i);
  1885. if (oa.arg) {
  1886. *ops++ = oa.op;
  1887. *args++ = oa.arg;
  1888. }
  1889. }
  1890. *ops++ = Op_clamp;
  1891. *ops++ = Op_clut;
  1892. *args++ = &srcProfile->A2B;
  1893. }
  1894. if (srcProfile->A2B.matrix_channels == 3) {
  1895. for (int i = 0; i < 3; i++) {
  1896. OpAndArg oa = select_curve_op(&srcProfile->A2B.matrix_curves[i], i);
  1897. if (oa.arg) {
  1898. *ops++ = oa.op;
  1899. *args++ = oa.arg;
  1900. }
  1901. }
  1902. static const skcms_Matrix3x4 I = {{
  1903. {1,0,0,0},
  1904. {0,1,0,0},
  1905. {0,0,1,0},
  1906. }};
  1907. if (0 != memcmp(&I, &srcProfile->A2B.matrix, sizeof(I))) {
  1908. *ops++ = Op_matrix_3x4;
  1909. *args++ = &srcProfile->A2B.matrix;
  1910. }
  1911. }
  1912. if (srcProfile->A2B.output_channels == 3) {
  1913. for (int i = 0; i < 3; i++) {
  1914. OpAndArg oa = select_curve_op(&srcProfile->A2B.output_curves[i], i);
  1915. if (oa.arg) {
  1916. *ops++ = oa.op;
  1917. *args++ = oa.arg;
  1918. }
  1919. }
  1920. }
  1921. if (srcProfile->pcs == skcms_Signature_Lab) {
  1922. *ops++ = Op_lab_to_xyz;
  1923. }
  1924. } else if (srcProfile->has_trc && srcProfile->has_toXYZD50) {
  1925. for (int i = 0; i < 3; i++) {
  1926. OpAndArg oa = select_curve_op(&srcProfile->trc[i], i);
  1927. if (oa.arg) {
  1928. *ops++ = oa.op;
  1929. *args++ = oa.arg;
  1930. }
  1931. }
  1932. } else {
  1933. return false;
  1934. }
  1935. // A2B sources should already be in XYZD50 at this point.
  1936. // Others still need to be transformed using their toXYZD50 matrix.
  1937. // N.B. There are profiles that contain both A2B tags and toXYZD50 matrices.
  1938. // If we use the A2B tags, we need to ignore the XYZD50 matrix entirely.
  1939. assert (srcProfile->has_A2B || srcProfile->has_toXYZD50);
  1940. static const skcms_Matrix3x3 I = {{
  1941. { 1.0f, 0.0f, 0.0f },
  1942. { 0.0f, 1.0f, 0.0f },
  1943. { 0.0f, 0.0f, 1.0f },
  1944. }};
  1945. const skcms_Matrix3x3* to_xyz = srcProfile->has_A2B ? &I : &srcProfile->toXYZD50;
  1946. // There's a chance the source and destination gamuts are identical,
  1947. // in which case we can skip the gamut transform.
  1948. if (0 != memcmp(&dstProfile->toXYZD50, to_xyz, sizeof(skcms_Matrix3x3))) {
  1949. // Concat the entire gamut transform into from_xyz,
  1950. // now slightly misnamed but it's a handy spot to stash the result.
  1951. from_xyz = skcms_Matrix3x3_concat(&from_xyz, to_xyz);
  1952. *ops++ = Op_matrix_3x3;
  1953. *args++ = &from_xyz;
  1954. }
  1955. // Encode back to dst RGB using its parametric transfer functions.
  1956. if (!is_identity_tf(&inv_dst_tf_r)) { *ops++ = Op_tf_r; *args++ = &inv_dst_tf_r; }
  1957. if (!is_identity_tf(&inv_dst_tf_g)) { *ops++ = Op_tf_g; *args++ = &inv_dst_tf_g; }
  1958. if (!is_identity_tf(&inv_dst_tf_b)) { *ops++ = Op_tf_b; *args++ = &inv_dst_tf_b; }
  1959. }
  1960. // Clamp here before premul to make sure we're clamping to normalized values _and_ gamut,
  1961. // not just to values that fit in [0,1].
  1962. //
  1963. // E.g. r = 1.1, a = 0.5 would fit fine in fixed point after premul (ra=0.55,a=0.5),
  1964. // but would be carrying r > 1, which is really unexpected for downstream consumers.
  1965. if (dstFmt < skcms_PixelFormat_RGB_hhh) {
  1966. *ops++ = Op_clamp;
  1967. }
  1968. if (dstAlpha == skcms_AlphaFormat_Opaque) {
  1969. *ops++ = Op_force_opaque;
  1970. } else if (dstAlpha == skcms_AlphaFormat_PremulAsEncoded) {
  1971. *ops++ = Op_premul;
  1972. }
  1973. if (dstFmt & 1) {
  1974. *ops++ = Op_swap_rb;
  1975. }
  1976. switch (dstFmt >> 1) {
  1977. default: return false;
  1978. case skcms_PixelFormat_A_8 >> 1: *ops++ = Op_store_a8; break;
  1979. case skcms_PixelFormat_G_8 >> 1: *ops++ = Op_store_g8; break;
  1980. case skcms_PixelFormat_ABGR_4444 >> 1: *ops++ = Op_store_4444; break;
  1981. case skcms_PixelFormat_RGB_565 >> 1: *ops++ = Op_store_565; break;
  1982. case skcms_PixelFormat_RGB_888 >> 1: *ops++ = Op_store_888; break;
  1983. case skcms_PixelFormat_RGBA_8888 >> 1: *ops++ = Op_store_8888; break;
  1984. case skcms_PixelFormat_RGBA_1010102 >> 1: *ops++ = Op_store_1010102; break;
  1985. case skcms_PixelFormat_RGB_161616LE >> 1: *ops++ = Op_store_161616LE; break;
  1986. case skcms_PixelFormat_RGBA_16161616LE >> 1: *ops++ = Op_store_16161616LE; break;
  1987. case skcms_PixelFormat_RGB_161616BE >> 1: *ops++ = Op_store_161616BE; break;
  1988. case skcms_PixelFormat_RGBA_16161616BE >> 1: *ops++ = Op_store_16161616BE; break;
  1989. case skcms_PixelFormat_RGB_hhh_Norm >> 1: *ops++ = Op_store_hhh; break;
  1990. case skcms_PixelFormat_RGBA_hhhh_Norm >> 1: *ops++ = Op_store_hhhh; break;
  1991. case skcms_PixelFormat_RGB_hhh >> 1: *ops++ = Op_store_hhh; break;
  1992. case skcms_PixelFormat_RGBA_hhhh >> 1: *ops++ = Op_store_hhhh; break;
  1993. case skcms_PixelFormat_RGB_fff >> 1: *ops++ = Op_store_fff; break;
  1994. case skcms_PixelFormat_RGBA_ffff >> 1: *ops++ = Op_store_ffff; break;
  1995. }
  1996. auto run = baseline::run_program;
  1997. #if defined(TEST_FOR_HSW)
  1998. switch (cpu_type()) {
  1999. case CpuType::None: break;
  2000. case CpuType::HSW: run = hsw::run_program; break;
  2001. case CpuType::SKX: run = hsw::run_program; break;
  2002. }
  2003. #endif
  2004. #if defined(TEST_FOR_SKX)
  2005. switch (cpu_type()) {
  2006. case CpuType::None: break;
  2007. case CpuType::HSW: break;
  2008. case CpuType::SKX: run = skx::run_program; break;
  2009. }
  2010. #endif
  2011. run(program, arguments, (const char*)src, (char*)dst, n, src_bpp,dst_bpp);
  2012. return true;
  2013. }
  2014. static void assert_usable_as_destination(const skcms_ICCProfile* profile) {
  2015. #if defined(NDEBUG)
  2016. (void)profile;
  2017. #else
  2018. skcms_Matrix3x3 fromXYZD50;
  2019. skcms_TransferFunction invR, invG, invB;
  2020. assert(prep_for_destination(profile, &fromXYZD50, &invR, &invG, &invB));
  2021. #endif
  2022. }
  2023. bool skcms_MakeUsableAsDestination(skcms_ICCProfile* profile) {
  2024. skcms_Matrix3x3 fromXYZD50;
  2025. if (!profile->has_trc || !profile->has_toXYZD50
  2026. || !skcms_Matrix3x3_invert(&profile->toXYZD50, &fromXYZD50)) {
  2027. return false;
  2028. }
  2029. skcms_TransferFunction tf[3];
  2030. for (int i = 0; i < 3; i++) {
  2031. skcms_TransferFunction inv;
  2032. if (profile->trc[i].table_entries == 0
  2033. && skcms_TransferFunction_invert(&profile->trc[i].parametric, &inv)) {
  2034. tf[i] = profile->trc[i].parametric;
  2035. continue;
  2036. }
  2037. float max_error;
  2038. // Parametric curves from skcms_ApproximateCurve() are guaranteed to be invertible.
  2039. if (!skcms_ApproximateCurve(&profile->trc[i], &tf[i], &max_error)) {
  2040. return false;
  2041. }
  2042. }
  2043. for (int i = 0; i < 3; ++i) {
  2044. profile->trc[i].table_entries = 0;
  2045. profile->trc[i].parametric = tf[i];
  2046. }
  2047. assert_usable_as_destination(profile);
  2048. return true;
  2049. }
  2050. bool skcms_MakeUsableAsDestinationWithSingleCurve(skcms_ICCProfile* profile) {
  2051. // Operate on a copy of profile, so we can choose the best TF for the original curves
  2052. skcms_ICCProfile result = *profile;
  2053. if (!skcms_MakeUsableAsDestination(&result)) {
  2054. return false;
  2055. }
  2056. int best_tf = 0;
  2057. float min_max_error = INFINITY_;
  2058. for (int i = 0; i < 3; i++) {
  2059. skcms_TransferFunction inv;
  2060. if (!skcms_TransferFunction_invert(&result.trc[i].parametric, &inv)) {
  2061. return false;
  2062. }
  2063. float err = 0;
  2064. for (int j = 0; j < 3; ++j) {
  2065. err = fmaxf_(err, max_roundtrip_error(&profile->trc[j], &inv));
  2066. }
  2067. if (min_max_error > err) {
  2068. min_max_error = err;
  2069. best_tf = i;
  2070. }
  2071. }
  2072. for (int i = 0; i < 3; i++) {
  2073. result.trc[i].parametric = result.trc[best_tf].parametric;
  2074. }
  2075. *profile = result;
  2076. assert_usable_as_destination(profile);
  2077. return true;
  2078. }