Full Width [alt+shift+f] Shortcuts [alt+shift+k] TRY SIMPLE MODE
Sign Up [alt+shift+s] Log In [alt+shift+l]
111
As I'm exploring the fantastic world of indie game development lately, I end up watching a large number of video tutorials on the subject. Even though the quality of the content is pretty variable, I'm very grateful to the creators for it. That being said, I couldn't help noticing this particular bit times and times again: a = lerp(a, B, delta * RATE) Behind this apparent banal call hides a terrible curse, forever perpetrated by innocent souls on the Internet. In this article we will study what it's trying to achieve, how it works, why it's wrong, and then we'll come up with a good solution to the initial problem. The usual warning: I don't have a mathematics or academic background, so the article is addressed at other neanderthals like myself, who managed to understand that pressing keys on a keyboard make pixels turn on and off. What is it? Let's start from the beginning. We're in a game engine main loop callback called at a regular interval (roughly), passing down the time...
a year ago

Improve your reading experience

Logged in users get linked directly to articles resulting in a better reading experience. Please login for free, it takes less than 1 minute.

More from A small freedom area RSS

Hacking window titles to help OBS

This write-up is meant to present the rationale and technical details behind a tiny project I wrote the other day, WTH, or WindowTitleHack, which is meant to force a constant window name for apps that keep changing it (I'm looking specifically at Firefox and Krita, but there are probably many others). Why tho? I've been streaming on Twitch from Linux (X11) with a barebone OBS Studio setup for a while now, and while most of the experience has been relatively smooth, one particularly striking frustration has been dealing with windows detection. If you don't want to capture the whole desktop for privacy reasons or simply to have control over the scene layout depending on the currently focused app, you need to rely on the Window Capture (XComposite) source. This works mostly fine, and it is actually able to track windows even when their title bar is renamed. But obviously, upon restart it can't find it again because both the window title and the window ID changed, meaning you have to redo your setup by reselecting the window again. It would have been acceptable if that was the only issue I had, but one of the more advanced feature I'm extensively using is the Advanced Scene Switcher (the builtin one, available through the Tools menu). This tool is a basic window title pattern matching system that allows automatic scene switches depending on the current window. Since it doesn't even seem to support regex, it's troublesome to have it reliably work with apps constantly changing their title (and even if it had, there is no guarantee that the app would leave a recognizable matchable pattern in its title). Hacking Windows One unreliable hack would be to spam xdotool commands to correct the window title. This could be a resource hog, and it would create quite a bunch of races. One slight improvement over this would be to use xprop -spy, but that wouldn't address the race conditions (since we would adjust the title after it's been already changed). So how do we deal with that properly? Well, on X11 with the reference library (Xlib) there are actually various (actually a lot of) ways of changing the title bar. It took me a while to identify which call(s) to target, but ended up with the following call graph, where each function is actually exposed publicly: From this we can easily see that we only need to hook the deepest function XChangeProperty, and check if the property is XA_WM_NAME (or its "modern" sibling, _NET_WM_NAME). How do we do that? With the help of the LD_PRELOAD environment variable and a dynamic library that implements a custom XChangeProperty. First, we grab the original function: #include <dlfcn.h> /* A type matching the prototype of the target function */ typedef int (*XChangeProperty_func_type)( Display *display, Window w, Atom property, Atom type, int format, int mode, const unsigned char *data, int nelements ); /* [...] */ XChangeProperty_func_type XChangeProperty_orig = dlsym(RTLD_NEXT, "XChangeProperty"); We also need to craft a custom _NET_WM_NAME atom: _NET_WM_NAME = XInternAtom(display, "_NET_WM_NAME", 0); With this we are now able to intercept all the WM_NAME events and override them with our own: if (property == XA_WM_NAME || property == _NET_WM_NAME) { data = (const unsigned char *)new_title; nelements = (int)strlen(new_title); } return XChangeProperty_orig(display, w, property, type, format, mode, data, nelements); We wrap all of this into our own redefinition of XChangeProperty and… that's pretty much it. Now due to a long history of development, Xlib has been "deprecated" and superseded by libxcb. Both are widely used, but fortunately the APIs are more or less similar. The function to hook is xcb_change_property, and defining _NET_WM_NAME is slightly more cumbered but not exactly challenging: const xcb_intern_atom_cookie_t cookie = xcb_intern_atom(conn, 0, strlen("_NET_WM_NAME"), "_NET_WM_NAME"); xcb_intern_atom_reply_t *reply = xcb_intern_atom_reply(conn, cookie, NULL); if (reply) _NET_WM_NAME = reply->atom; free(reply); Aside from that, the code is pretty much the same. Configuration To pass down the custom title to override, I've been relying on an environment variable WTH_TITLE. From a user point of view, it looks like this: LD_PRELOAD="builddir/libwth.so" WTH_TITLE="Krita4ever" krita We could probably improve the usability by creating a wrapping tool (so that we could have something such as ./wth --title=Krita4ever krita). Unfortunately I wasn't yet able to make a self-referencing executable accepted by LD_PRELOAD, so for now the manual LD_PRELOAD and WTH_TITLE environment will do just fine. Thread safety To avoid a bunch of redundant function roundtrips we need to globally cache a few things: the new title (to avoid fetching it in the environment all the time), the original functions (to save the dlsym call), and _NET_WM_NAME. Those are loaded lazily at the first function call, but we have no guarantee with regards to concurrent calls on that hooked function so we must create our own lock. I initially though about using pthread_once but unfortunately the initialization callback mechanism doesn't allow any custom argument. Again, this is merely a slight annoyance since we can implement our own in a few lines of code: /* The "once" API is similar to pthread_once but allows a custom function argument */ struct wth_once { pthread_mutex_t lock; int initialized; }; #define WTH_ONCE_INITIALIZER {.lock=PTHREAD_MUTEX_INITIALIZER} typedef void (*init_func_type)(void *user_arg); void wth_init_once(struct wth_once *once, init_func_type init_func, void *user_arg) { pthread_mutex_lock(&once->lock); if (!once->initialized) { init_func(user_arg); once->initialized = 1; } pthread_mutex_unlock(&once->lock); } Which we use like this: static struct wth_once once = WTH_ONCE_INITIALIZER; static void init_once(void *user_arg) { Display *display = user_arg; /* [...] */ } /* [...] */ wth_init_once(&once, init_once, display); The End? I've been delaying doing this project for weeks because it felt complex at first glance, but it actually just took me a few hours. Probably the same amount of time it took me to write this article. While the project is admittedly really small, it still feel like a nice accomplishment. I hope it's useful to other people. Now, the Wayland support is probably the most obvious improvement the project can receive, but I don't have such a setup locally to test yet, so this is postponed for an undetermined amount of time. The code is released with a permissive license (MIT); if you want to contribute you can open a pull request but getting in touch with me first is appreciated to avoid unnecessary and overlapping efforts.

over a year ago 99 votes
Improving color quantization heuristics

In 2015, I wrote an article about how the palette color quantization was improved in FFmpeg in order to make nice animated GIF files. For some reason, to this day this is one of my most popular article. As time passed, my experience with colors grew and I ended up being quite ashamed and frustrated with the state of these filters. A lot of the code was naive (when not terribly wrong), despite the apparent good results. One of the major change I wanted to do was to evaluate the color distances using a perceptually uniform colorspace, instead of using a naive euclidean distance of RGB triplets. As usual it felt like a week-end long project; after all, all I have to do is change the distance function to work in a different space, right? Well, if you're following my blog you might have noticed I've add numerous adventures that stacked up on each others: I had to work out the colorspace with integer arithmetic first ...which forced me to look into integer division more deeply ...which confronted me to all sort of undefined behaviours in the process And when I finally reached the point where I could make the switch to OkLab (the perceptual colorspace), a few experiments showed that the flavour of the core algorithm I was using might contain some fundamental flaws, or at least was not implementing optimal heuristics. So here we go again, quickly enough I find myself starting a new research study in the pursuit of understanding how to put pixels on the screen. This write-up is the story of yet another self-inflicted struggle. Palette quantization But what is palette quantization? It essentially refers to the process of reducing the number of available colors of an image down to a smaller subset. In sRGB, an image can have up to 16.7 million colors. In practice though it's generally much less, to the surprise of no one. Still, it's not rare to have a few hundreds of thousands different colors in a single picture. Our goal is to reduce that to something like 256 colors that represent them best, and use these colors to create a new picture. Why you may ask? There are multiple reasons, here are some: Improve size compression (this is a lossy operation of course, and using dithering on top might actually defeat the original purpose) Some codecs might not support anything else than limited palettes (GIF or subtitles codecs are examples) Various artistic purposes Following is an example of a picture quantized at different levels: Original (26125 colors) Quantized to 8bpp (256 colors) Quantized to 2bpp (4 colors) This color quantization process can be roughly summarized in a 4-steps based process: Sample the input image: we build an histogram of all the colors in the picture (basically a simple statistical analysis) Design a colormap: we build the palette through various means using the histograms Create a pixel mapping which associates a color (one that can be found in the input image) with another (one that can be found in the newly created palette) Image quantizing: we use the color mapping to build our new image. This step may also involve some dithering. The study here will focus on step 2 (which itself relies on step 1). Colormap design algorithms A palette is simply a set of colors. It can be represented in various ways, for example here in 2D and 3D: To generate such a palette, all sort of algorithms exists. They are usually classified into 2 large categories: Dividing/splitting algorithms (such as Median-Cut and its various flavors) Clustering algorithms (such as K-means, maximin distance, (E)LBG or pairwise clustering) The former are faster but non-optimal while the latter are slower but better. The problem is NP-complete, meaning it's possible to find the optimal solution but it can be extremely costly. On the other hand, it's possible to find "local optimums" at minimal cost. Since I'm working within FFmpeg, speed has always been a priority. This was the reason that motivated me to initially implement the Median-Cut over a more expensive algorithm. The rough picture of the algorithm is relatively easy to grasp. Assuming we want a palette of K colors: A set S of all the colors in the input picture is constructed, along with a respective set W of the weight of each color (how much they appear) Since the colors are expressed as RGB triplets, they can be encapsulated in one big cuboid, or box The box is cut in two along one of the axis (R, G or B) on the median (hence the name of the algorithm) If we don't have a total K boxes yet, pick one of them and go back to previous step All the colors in each of the K boxes are then averaged to form the color palette entries Here is how the process looks like visually: Median-Cut algorithm targeting 16 boxes You may have spotted in this video that the colors are not expressed in RGB but in Lab: this is because instead of representing the colors in a traditional RGB colorspace, we are instead using the OkLab colorspace which has the property of being perceptually uniform. It doesn't really change the Median Cut algorithm, but it definitely has an impact on the resulting palette. One striking limitation of this algorithm is that we are working exclusively with cuboids: the cuts are limited to an axis, we are not cutting along an arbitrary plane or a more complex shape. Think of it like working with voxels instead of more free-form geometries. The main benefit is that the algorithm is pretty simple to implement. Now the description provided earlier conveniently avoided describing two important aspects happening in step 3 and 4: How do we choose the next box to split? How do we choose along which axis of the box we make the cut? I pondered about that for a quite a long time. An overview of the possible heuristics In bulk, some of the heuristics I started thinking of: should we take the box that has the longest axis across all boxes? should we take the box that has the largest volume? should we take the box that has the biggest Mean Squared Error when compared to its average color? should we take the box that has the axis with the biggest MSE? assuming we choose to go with the MSE, should it be normalized across all boxes? should we even account for the weight of each color or consider them equal? what about the axis? Is it better to pick the longest or the one with the higher MSE? I tried to formalize these questions mathematically to the best of my limited abilities. So let's start by saying that all the colors c of given box are stored in a N×M 2D-array following the matrix notation: L₁L₂L₃…Lₘ a₁a₂a₃…aₘ b₁b₂b₃…bₘ N is the number of components (3 in our case, whether it's RGB or Lab), and M the number of colors in that box. You can visualize this as a list of vectors as well, where c_{i,j} is the color at row i and column j. With that in mind we can sketch the following diagram representing the tree of heuristic possibilities to implement: Mathematicians are going to kill me for doodling random notes all over this perfectly understandable symbols gibberish, but I believe this is required for the human beings reading this article. In summary, we end up with a total of 24 combinations to try out: 2 axis selection heuristics: cut the axis with the maximum error squared cut the axis with the maximum length 3 operators: maximum measurement out of all the channels product of the measurements of all the channels sum of the measurements of all the channels 4 measurements: error squared, honoring weights error squared, not honoring weights error squared, honoring weights, normalized length of the axis If we start to intuitively think about which ones are likely going to perform the best, we quickly realize that we haven't actually formalized what we are trying to achieve. Such a rookie mistake. Clarifying this will help us getting a better feeling about the likely outcome. I chose to target an output that minimizes the MSE against the reference image, in a perceptual way. Said differently, trying to make the perceptual distance between an input and output color pixel as minimal as possible. This is an arbitrary and debatable target, but it's relatively simple and objective to evaluate if we have faith in the selected perceptual model. Another appropriate metric could have been to find the ideal palette through another algorithm, and compare against that instead. Doing that unfortunately implied that I would trust that other algorithm, its implementation, and that I have enough computing power. So to summarize, we want to minimize the MSE between the input and output, evaluated in the OkLab colorspace. This can be expressed with the following formula: Where: P is a partition (which we constrain to a box in our implementation) C the set of colors in the partition P w the weight of a color c a single color µ the average color of the set C Special thanks to criver for helping me a ton on the math area, this last formula is from them. Looking at the formula, we can see how similar it is to certain branches of the heuristics tree, so we can start getting an intuition about the result of the experiment. Experiment language Short deviation from the main topic (feel free to skip to the next section): working in C within FFmpeg quickly became a hurdle more than anything. Aside from the lack of flexibility, the implicit casts destroying the precision deceitfully, and the undefined behaviours, all kind of C quirks went in the way several times, which made me question my sanity. This one typically severly messed me up while trying to average the colors: #include <stdio.h> #include <stdint.h> int main (void) { const int32_t x = -30; const uint32_t y = 10; const uint32_t a = 30; const int32_t b = -10; printf("%d×%u=%d\n", x, y, x * y); printf("%u×%d=%d\n", a, b, a * b); printf("%d/%u=%d\n", x, y, x / y); printf("%u/%d=%d\n", a, b, a / b); return 0; } % cc -Wall -Wextra -fsanitize=undefined test.c -o test && ./test -30×10=-300 30×-10=-300 -30/10=429496726 30/-10=0 Anyway, I know this is obvious but if you aren't already doing that I suggest you build your experiments in another language, Python or whatever, and rewrite them in C later when you figured out your expected output. Re-implementing what I needed in Python didn't take me long. It was, and still is obviously much slower at runtime, but that's fine. There is a lot of room for speed improvement, typically by relying on numpy (which I didn't bother with). Experiment results I created a research repository for the occasion. The code to reproduce and the results can be found in the color quantization README. In short, based on the results, we can conclude that: Overall, the box that has the axis with the largest non-normalized weighted sum of squared error is the best candidate in the box selection algorithm Overall, cutting the axis with the largest weighted sum of squared error is the best axis cut selection algorithm To my surprise, normalizing the weights per box is not a good idea. I initially observed that by trial and error, which was actually one of the main motivator for this research. I initially thought normalizing each box was necessary in order to compare them against each others (such that they are compared on a common ground). My loose explanation of the phenomenon was that not normalizing causes a bias towards boxes with many colors, but that's actually exactly what we want. I believe it can also be explained by our evaluation function: we want to minimize the error across the whole set of colors, so small partitions (in color counts) must not be made stronger. At least not in the context of the target we chose. It's also interesting to see how the max() seems to perform better than the sum() of the variance of each component most of the time. Admittedly my picture samples set is not that big, which may imply that more experiments to confirm that tendency are required. In retrospective, this might have been quickly predictable to someone with a mathematical background. But since I don't have that, nor do I trust my abstract thinking much, I'm kind of forced to try things out often. This is likely one of the many instances where I spent way too much energy on something obvious from the beginning, but I have the hope it will actually provide some useful information for other lost souls out there. Known limitations There are two main limitations I want to discuss before closing this article. The first one is related to minimizing the MSE even more. K-means refinement We know the Median-Cut actually provides a rough estimate of the optimal palette. One thing we could do is use it as a first step before refinement, for example by running a few K-means iterations as post-processing (how much refinement/iterations could be a user control). The general idea of K-means is to progressively move each colors individually to a more appropriate box, that is a box for which the color distance to the average color of that box is smaller. I started implementing that in a very naive way, so it's extremely slow, but that's something to investigate further because it definitely improves the results. Most of the academic literature seems to suggest the use of the K-means clustering, but all of them require some startup step. Some come up with various heuristics, some use PCA, but I've yet to see one that rely on Median-Cut as first pass; maybe that's not such a good idea, but who knows. Bias toward perceived lightness Another more annoying problem for which I have no solution for is with regards to the human perception being much more sensitive to light changes than hue. If you look at the first demo with the parrot, you may have observed the boxes are kind of thin. This is because the a and b components (respectively how green/red and blue/yellow the color is) have a much smaller amplitude compared to the L (perceived lightness). Here is a side by side comparison of the spread of colors between a stretched and normalized view: You may rightfully question whether this is a problem or not. In practice, this means that when K is low (let's say smaller than 8 or even 16), cuts along L will almost always be preferred, causing the picture to be heavily desaturated. This is because it tries to preserve the most significant attribute in human perception: the lightness. That particular picture is actually a pathological study case: 4 colors 8 colors 12 colors 16 colors We can see the hue timidly appearing around K=16 (specifically it starts being more strongly noticeable starting the cut K=13). Conclusion For now, I'm mostly done with this "week-end long project" into which I actually poured 2 or 3 months of lifetime. The FFmpeg patchset will likely be upstreamed soon so everyone should hopefully be able to benefit from it in the next release. It will also come with additional dithering methods, which implementation actually was a relaxing distraction from all this hardship. There are still many ways of improving this work, but it's the end of the line for me, so I'll trust the Internet with it.

over a year ago 58 votes
Porting OkLab colorspace to integer arithmetic

For reasons I'll explain in a futur write-up, I needed to make use of a perceptually uniform colorspace in some computer vision code. OkLab from Björn Ottosson was a great candidate given how simple the implementation is. But there is a plot twist: I needed the code to be deterministic for the tests to be portable across a large variety of architecture, systems and configurations. Several solutions were offered to me, including reworking the test framework to support a difference mechanism with threshold, but having done that in another project I can confidently say that it's not trivial (when not downright impossible in certain cases). Another approach would have been to hardcode the libc math functions, but even then I wasn't confident the floating point arithmetic would determinism would be guaranteed in all cases. So I ended up choosing to port the code to integer arithmetic. I'm sure many people would disagree with that approach, but: code determinism is guaranteed not all FPU are that efficient, typically on embedded it can now be used in the kernel; while this is far-fetched for OkLab (though maybe someone needs some color management in v4l2 or something), sRGB transforms might have their use cases it's a learning experience which can be re-used in other circumstances working on the integer arithmetic versions unlocked various optimizations for the normal case Note: I'm following Björn Ottosson will to have OkLab code in the public domain as well as under MIT license, so this "dual licensing" applies to all the code presented in this article. Warning: The integer arithmetics in this write-up can only work if your language behaves the same as C99 (or more recent) with regard to integer division. See this previous article on integer division for more information. Quick summary of uniform colorspaces For those unfamiliar with color management, one of the main benefit of a uniform colorspace like OkLab is that the euclidean distance between two colors is directly correlated with the human perception of these colors. More concretely, if we want to evaluate the distance between the RGB triplets (R₀,G₀,B₀) and (R₁,G₁,B₁), one may naively compute the euclidean distance √((R₀-R₁)²+(G₀-G₁)²+(B₀-B₁)²). Unfortunately, even if the RGB is gamma expanded into linear values, the computed distance will actually be pretty far from reflecting how the human eye perceive this difference. It typically isn't going to be consistent when applied to another pair of colors. With OkLab (and many other uniform colorspaces), the colors are also identified with 3D coordinates, but instead of (R,G,B) we call them (L,a,b) (which is an entirely different 3D space). In that space √((L₀-L₁)²+(a₀-a₁)²+(b₀-b₁)² (called ΔE, or Delta-E) is expected to be aligned with human perception of color differences. Of course, this is just one model, and it doesn't take into account many parameters. For instance, the perception of a color depends a lot on the surrounding colors. Still, these models are much better than working with RGB triplets, which don't make much sense visually speaking. Reference code / diagram In this study case, We will be focusing on the transform that goes from sRGB to OkLab, and back again into sRGB. Only the first part is interesting if we want the color distance, but sometimes we also want to alter a color uniformly and thus we need the 2nd part as well to reconstruct an sRGB color from it. We are only considering the sRGB input and output for simplicity, which means we will be inlining the sRGB color transfer in the pipeline. If you're not familiar with gamma compression, there are many resources about it on the Internet which you may want to look into first. Here is a diagram of the complete pipeline: And the corresponding code (of the 4 circles in the diagram) we will be porting: struct Lab { float L, a, b; } uint8_t linear_f32_to_srgb_u8(float x) { if (x <= 0.0) { return 0; } else if (x >= 1.0) { return 0xff; } else { const float v = x < 0.0031308f ? x*12.92f : 1.055f*powf(x, 1.f/2.4f) - 0.055f; return lrintf(v * 255.f); } } float srgb_u8_to_linear_f32(uint8_t x) { const float v = x / 255.f; return v < 0.04045f ? v/12.92f : powf((v+0.055f)/1.055f, 2.4f); } struct Lab srgb_u8_to_oklab_f32(uint32_t srgb) { const float r = srgb_u8_to_linear_f32(srgb >> 16 & 0xff); const float g = srgb_u8_to_linear_f32(srgb >> 8 & 0xff); const float b = srgb_u8_to_linear_f32(srgb & 0xff); const float l = 0.4122214708f * r + 0.5363325363f * g + 0.0514459929f * b; const float m = 0.2119034982f * r + 0.6806995451f * g + 0.1073969566f * b; const float s = 0.0883024619f * r + 0.2817188376f * g + 0.6299787005f * b; const float l_ = cbrtf(l); const float m_ = cbrtf(m); const float s_ = cbrtf(s); const struct Lab ret = { .L = 0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_, .a = 1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_, .b = 0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_, }; return ret; } uint32_t oklab_f32_to_srgb_u8(struct Lab c) { const float l_ = c.L + 0.3963377774f * c.a + 0.2158037573f * c.b; const float m_ = c.L - 0.1055613458f * c.a - 0.0638541728f * c.b; const float s_ = c.L - 0.0894841775f * c.a - 1.2914855480f * c.b; const float l = l_*l_*l_; const float m = m_*m_*m_; const float s = s_*s_*s_; const uint8_t r = linear_f32_to_srgb_u8(+4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s); const uint8_t g = linear_f32_to_srgb_u8(-1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s); const uint8_t b = linear_f32_to_srgb_u8(-0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s); return r<<16 | g<<8 | b; } sRGB to Linear The first step is converting the sRGB color to linear values. That sRGB transfer function can be intimidating, but it's pretty much a simple power function: The input is 8-bit ([0x00;0xff] for each of the 3 channels) which means we can use a simple 256 values lookup table containing the precomputed resulting linear values. Note that we can already do that with the reference code with a table remapping the 8-bit index into a float value. For our integer version we need to pick an arbitrary precision for the linear representation. 8-bit is not going to be enough precision, so we're going to pick the next power of two to be space efficient: 16-bit. We will be using the constant K=(1<<16)-1=0xffff to refer to this scale. Alternatively we could rely on a fixed point mapping (an integer for the decimal part and another integer for the fractional part), but in our case pretty much everything is normalized so the decimal part doesn't really matter. /** * Table mapping formula: * f(x) = x < 0.04045 ? x/12.92 : ((x+0.055)/1.055)^2.4 (sRGB EOTF) * Where x is the normalized index in the table and f(x) the value in the table. * f(x) is remapped to [0;K] and rounded. */ static const uint16_t srgb2linear[256] = { 0x0000, 0x0014, 0x0028, 0x003c, 0x0050, 0x0063, 0x0077, 0x008b, 0x009f, 0x00b3, 0x00c7, 0x00db, 0x00f1, 0x0108, 0x0120, 0x0139, 0x0154, 0x016f, 0x018c, 0x01ab, 0x01ca, 0x01eb, 0x020e, 0x0232, 0x0257, 0x027d, 0x02a5, 0x02ce, 0x02f9, 0x0325, 0x0353, 0x0382, 0x03b3, 0x03e5, 0x0418, 0x044d, 0x0484, 0x04bc, 0x04f6, 0x0532, 0x056f, 0x05ad, 0x05ed, 0x062f, 0x0673, 0x06b8, 0x06fe, 0x0747, 0x0791, 0x07dd, 0x082a, 0x087a, 0x08ca, 0x091d, 0x0972, 0x09c8, 0x0a20, 0x0a79, 0x0ad5, 0x0b32, 0x0b91, 0x0bf2, 0x0c55, 0x0cba, 0x0d20, 0x0d88, 0x0df2, 0x0e5e, 0x0ecc, 0x0f3c, 0x0fae, 0x1021, 0x1097, 0x110e, 0x1188, 0x1203, 0x1280, 0x1300, 0x1381, 0x1404, 0x1489, 0x1510, 0x159a, 0x1625, 0x16b2, 0x1741, 0x17d3, 0x1866, 0x18fb, 0x1993, 0x1a2c, 0x1ac8, 0x1b66, 0x1c06, 0x1ca7, 0x1d4c, 0x1df2, 0x1e9a, 0x1f44, 0x1ff1, 0x20a0, 0x2150, 0x2204, 0x22b9, 0x2370, 0x242a, 0x24e5, 0x25a3, 0x2664, 0x2726, 0x27eb, 0x28b1, 0x297b, 0x2a46, 0x2b14, 0x2be3, 0x2cb6, 0x2d8a, 0x2e61, 0x2f3a, 0x3015, 0x30f2, 0x31d2, 0x32b4, 0x3399, 0x3480, 0x3569, 0x3655, 0x3742, 0x3833, 0x3925, 0x3a1a, 0x3b12, 0x3c0b, 0x3d07, 0x3e06, 0x3f07, 0x400a, 0x4110, 0x4218, 0x4323, 0x4430, 0x453f, 0x4651, 0x4765, 0x487c, 0x4995, 0x4ab1, 0x4bcf, 0x4cf0, 0x4e13, 0x4f39, 0x5061, 0x518c, 0x52b9, 0x53e9, 0x551b, 0x5650, 0x5787, 0x58c1, 0x59fe, 0x5b3d, 0x5c7e, 0x5dc2, 0x5f09, 0x6052, 0x619e, 0x62ed, 0x643e, 0x6591, 0x66e8, 0x6840, 0x699c, 0x6afa, 0x6c5b, 0x6dbe, 0x6f24, 0x708d, 0x71f8, 0x7366, 0x74d7, 0x764a, 0x77c0, 0x7939, 0x7ab4, 0x7c32, 0x7db3, 0x7f37, 0x80bd, 0x8246, 0x83d1, 0x855f, 0x86f0, 0x8884, 0x8a1b, 0x8bb4, 0x8d50, 0x8eef, 0x9090, 0x9235, 0x93dc, 0x9586, 0x9732, 0x98e2, 0x9a94, 0x9c49, 0x9e01, 0x9fbb, 0xa179, 0xa339, 0xa4fc, 0xa6c2, 0xa88b, 0xaa56, 0xac25, 0xadf6, 0xafca, 0xb1a1, 0xb37b, 0xb557, 0xb737, 0xb919, 0xbaff, 0xbce7, 0xbed2, 0xc0c0, 0xc2b1, 0xc4a5, 0xc69c, 0xc895, 0xca92, 0xcc91, 0xce94, 0xd099, 0xd2a1, 0xd4ad, 0xd6bb, 0xd8cc, 0xdae0, 0xdcf7, 0xdf11, 0xe12e, 0xe34e, 0xe571, 0xe797, 0xe9c0, 0xebec, 0xee1b, 0xf04d, 0xf282, 0xf4ba, 0xf6f5, 0xf933, 0xfb74, 0xfdb8, 0xffff, }; int32_t srgb_u8_to_linear_int(uint8_t x) { return (int32_t)srgb2linear[x]; } You may have noticed that we are returning the value in a i32: this is to ease arithmetic operations (preserving the 16-bit unsigned precision would have overflow warping implications when working with the value). Linear to OkLab The OkLab is expressed in a virtually continuous space (floats). If we feed all 16.7 millions sRGB colors to the OkLab transform we get the following ranges in output: min Lab: 0.000000 -0.233887 -0.311528 max Lab: 1.000000 0.276216 0.198570 We observe that L is always positive and neatly within [0;1] while a and b are in a more restricted and signed range. Multiple choices are offered to us with regard to the integer representation we pick. Since we chose 16-bit for the input linear value, it makes sense to preserve that precision for Lab. For the L component, this fits neatly ([0;1] in the ref maps to [0;0xffff] in the integer version), but for the a and b component, not so much. We could pick a signed 16-bit, but that would imply a 15-bit precision for the arithmetic and 1-bit for the sign, which is going to be troublesome: we want to preserve the same precision for L, a and b since the whole point of this operation is to have a uniform space. Instead, I decided to go with 16-bits of precision, with one extra bit for the sign (which will be used for a and b), and thus storing Lab in 3 signed i32. Alternatively, we could decide to have a 15-bit precision with an extra bit for the sign by using 3 i16. This should work mostly fine but having the values fit exactly the boundaries of the storage can be problematic in various situations, typically anything that involves boundary checks and overflows. Picking a larger storage simplifies a bunch of things. Looking at srgb_u8_to_oklab_f32 we quickly see that for most of the function it's simple arithmetic, but we have a cube root (cbrt()), so let's study that first. Cube root All the cbrt inputs are driven by this: const float l = 0.4122214708f * r + 0.5363325363f * g + 0.0514459929f * b; const float m = 0.2119034982f * r + 0.6806995451f * g + 0.1073969566f * b; const float s = 0.0883024619f * r + 0.2817188376f * g + 0.6299787005f * b; This might not be obvious at first glance but here l, m and s all are in [0;1] range (the sum of the coefficients of each row is 1), so we will only need to deal with this range in our cbrt implementation. This simplifies greatly the problem! Now, what does it look like? This function is simply the inverse of f(x)=x³, which is a more convenient function to work with. And I have some great news: not long ago, I wrote an article on how to inverse a function, so that's exactly what we are going to do here: inverse f(x)=x³. What we first need though is a good approximation of the curve. A straight line is probably fine but we could try to use some symbolic regression in order to get some sort of rough polynomial approximation. PySR can do that in a few lines of code: import numpy as np from pysr import PySRRegressor # 25 points of ³√x within [0;1] x = np.linspace(0, 1, 25).reshape(-1, 1) y = x ** (1/3) model = PySRRegressor(model_selection="accuracy", binary_operators=["+", "-", "*"], niterations=200) r = model.fit(x, y, variable_names=["x"]) print(r) The output is not deterministic for some reason (which is quite annoying) and the expressions provided usually follows a wonky form. Still, in my run it seemed to take a liking to the following polynomial: u₀ = x³ - 2.19893x² + 2.01593x + 0.219407 (reformatted in a sane polynomial form thanks to WolframAlpha). Note that increasing the number of data points is not really a good idea because we quickly start being confronted to Runge's phenomenon. No need to overthink it, 25 points is just fine. Now we can make a few Newton iterations. For that, we need the derivative of f(uₙ)=uₙ³-x, so f'(uₙ)=3uₙ² and thus the iteration expressions can be obtained easily: uₙ₊₁ = uₙ - (f(uₙ)-v)/f'(uₙ) = uₙ - (uₙ³-v)/(3uₙ²) = (2uₙ³+v)/(3uₙ²) If you don't understand what the hell is going on here, check the article referred to earlier, we're simply following the recipe here. Now I had a look into how most libc compute cbrt, and despite sometimes referring to Newton iterations, they were actually using Halley iterations. So we're going to do the same (not lying, just the Halley part). To get the Halley iteration instead of Newton, we need the first but also the second derivative of f(uₙ)=uₙ³-x (f'(uₙ)=3uₙ² and f"(uₙ)=6uₙ) from which we deduce a relatively simple expression: uₙ₊₁ = uₙ-2f(uₙ)f'(uₙ)/(2f'(uₙ)²-f(uₙ)f"(uₙ)) = uₙ(2x+uₙ³)/(x+2uₙ³) We have everything we need to approximate a cube root of a real between 0 and 1. In Python a complete implementation would be as simple as this snippet: b, c, d = -2.19893, 2.01593, 0.219407 def cbrt01(x): # We only support [0;1] if x <= 0: return 0 if x >= 1: return 1 # Initial approximation u = x**3 + b*x**2 + c*x + d # 2 Halley iterations u = u * (2*x+u**3) / (x+2*u**3) u = u * (2*x+u**3) / (x+2*u**3) return u But now we need to scale the floating values up into 16-bit integers. First of all, in the integer version our x is actually in K scale, which means we want to express u according to X=x·K. Similarly, we want to use B=b·K, C=c·K and D=d·K instead of b, c and d because we have no way of expressing the former as integer otherwise. Finally, we're not actually going to compute u₀ but u₀·K because we're preserving the scale through the function. We have: u₀·K = K·(x³ + bx² + cx + d) = K·((x·K)³/K³ + b(x·K)²/K² + c(x·K)/K + d) = K·(X³/K³ + bX²/K² + cX/K + d) = X³·K/K³ + bX²·K/K² + cX·K/K + d·K = X³/K² + BX²/K² + CX/K + D = X³/K² + BX²/K² + CX/K + D = (X³ + BX²)/K² + CX/K + D = ((X³ + BX²)/K + CX)/K + D = (X(X² + BX)/K + CX)/K + D U₀ = (X(X(X + B)/K + CX)/K + D With this we have a relatively cheap expression where the K divisions would still preserve enough precision even if evaluated as integer division. We can do the same for the Halley iteration. I spare you the algebra, the expression u(2x+u³) / (x+2u³) becomes (U(2X+U³/K²)) / (X+2U³/K²). Looking at this expression you may start to worry about overflows, and that would fair since even K² is getting dangerously close to the sun (it's actually already larger than INT32_MAX). For this reason we're going to cheat and simply use 64-bit arithmetic in this function. I believe we could reduce the risk of overflow, but I don't think there is a way to remain in 32-bit without nasty compromises anyway. This is also why in the code below you'll notice the constants are suffixed with LL (to force long-long/64-bit arithmetic). Beware that overflows are a terrible predicament to get into as they will lead to undefined behaviour. Do not underestimate this risk. You might not detect them early enough, and missing them may mislead you when interpreting the results. For this reason, I strongly suggest to always build with -fsanitize=undefined during test and development. I don't do that often, but for this kind of research, I also highly recommend to first write tests that cover all possible integers input (when applicable) so that overflows are detected as soon as possible. Before we write the integer version of our function, we need to address rounding. In the case of the initial approximation I don't think we need to bother, but for our Halley iteration we're going to need as much precision as we can get. Since we know U is positive (remember we're evaluating cbrt(x) where x is in [0;1]), we can use the (a+b/2)/b rounding formula. Our function finally just looks like: #define K2 ((int64_t)K*K) int32_t cbrt01_int(int32_t x) { int64_t u; /* Approximation curve is for the [0;1] range */ if (x <= 0) return 0; if (x >= K) return K; /* Initial approximation: x³ - 2.19893x² + 2.01593x + 0.219407 */ u = x*(x*(x - 144107LL) / K + 132114LL) / K + 14379LL; /* Refine with 2 Halley iterations. */ for (int i = 0; i < 2; i++) { const int64_t u3 = u*u*u; const int64_t den = x + (2*u3 + K2/2) / K2; u = (u * (2*x + (u3 + K2/2) / K2) + den/2) / den; } return u; } Cute, isn't it? If we test the accuracy of this function by calling it for all the possible values we actually get extremely good results. Here is a test code: int main(void) { float max_diff = 0; float total_diff = 0; for (int i = 0; i <= K; i++) { const float ref = cbrtf(i / (float)K); const float out = cbrt01_int(i) / (float)K; const float d = fabs(ref - out); if (d > max_diff) max_diff = d; total_diff += d; } printf("max_diff=%f total_diff=%f avg_diff=%f\n", max_diff, total_diff, total_diff / (K + 1)); return 0; } Output: max_diff=0.030831 total_diff=0.816078 avg_diff=0.000012 If we want to trade precision for speed, we could adjust the function to use Newton iterations, and maybe remove the rounding. Back to the core Going back to our sRGB-to-OkLab function, everything should look straightforward to implement now. There is one thing though, while lms computation (at the beginning of the function) is exclusively working with positive values, the output Lab value expression is signed. For this reason we will need a more involved rounded division, so referring again to my last article we will use: static int64_t div_round64(int64_t a, int64_t b) { return (a^b)<0 ? (a-b/2)/b : (a+b/2)/b; } And thus, we have: struct LabInt { int32_t L, a, b; }; struct LabInt srgb_u8_to_oklab_int(uint32_t srgb) { const int32_t r = (int32_t)srgb2linear[srgb >> 16 & 0xff]; const int32_t g = (int32_t)srgb2linear[srgb >> 8 & 0xff]; const int32_t b = (int32_t)srgb2linear[srgb & 0xff]; // Note: lms can actually be slightly over K due to rounded coefficients const int32_t l = (27015LL*r + 35149LL*g + 3372LL*b + K/2) / K; const int32_t m = (13887LL*r + 44610LL*g + 7038LL*b + K/2) / K; const int32_t s = (5787LL*r + 18462LL*g + 41286LL*b + K/2) / K; const int32_t l_ = cbrt01_int(l); const int32_t m_ = cbrt01_int(m); const int32_t s_ = cbrt01_int(s); const struct LabInt ret = { .L = div_round64( 13792LL*l_ + 52010LL*m_ - 267LL*s_, K), .a = div_round64(129628LL*l_ - 159158LL*m_ + 29530LL*s_, K), .b = div_round64( 1698LL*l_ + 51299LL*m_ - 52997LL*s_, K), }; return ret; } The note in this code is here to remind us that we have to saturate lms to a maximum of K (corresponding to 1.0 with floats), which is what we're doing in cbrt01_int(). At this point we can already work within the OkLab space but we're only half-way through the pain. Fortunately, things are going to be easier from now on. OkLab to sRGB Our OkLab-to-sRGB function relies on the Linear-to-sRGB function (at the end), so we're going to deal with it first. Linear to sRGB Contrary to sRGB-to-Linear it's going to be tricky to rely on a table because it would be way too large to hold all possible values (since it would require K entries). I initially considered computing powf(x, 1.f/2.4f) with integer arithmetic somehow, but this is much more involved than how we managed to implement cbrt. So instead I thought about approximating the curve with a bunch of points (stored in a table), and then approximate any intermediate value with a linear interpolation, that is as if the point were joined through small segments. We gave 256 16-bit entries to srgb2linear, so if we were to give as much storage to linear2srgb we could have a table of 512 8-bit entries (our output is 8-bit). Here it is: /** * Table mapping formula: * f(x) = x < 0.0031308 ? x*12.92 : (1.055)*x^(1/2.4)-0.055 (sRGB OETF) * Where x is the normalized index in the table and f(x) the value in the table. * f(x) is remapped to [0;0xff] and rounded. * * Since a 16-bit table is too large, we reduce its precision to 9-bit. */ static const uint8_t linear2srgb[P + 1] = { 0x00, 0x06, 0x0d, 0x12, 0x16, 0x19, 0x1c, 0x1f, 0x22, 0x24, 0x26, 0x28, 0x2a, 0x2c, 0x2e, 0x30, 0x32, 0x33, 0x35, 0x36, 0x38, 0x39, 0x3b, 0x3c, 0x3d, 0x3e, 0x40, 0x41, 0x42, 0x43, 0x45, 0x46, 0x47, 0x48, 0x49, 0x4a, 0x4b, 0x4c, 0x4d, 0x4e, 0x4f, 0x50, 0x51, 0x52, 0x53, 0x54, 0x55, 0x56, 0x56, 0x57, 0x58, 0x59, 0x5a, 0x5b, 0x5b, 0x5c, 0x5d, 0x5e, 0x5f, 0x5f, 0x60, 0x61, 0x62, 0x62, 0x63, 0x64, 0x65, 0x65, 0x66, 0x67, 0x67, 0x68, 0x69, 0x6a, 0x6a, 0x6b, 0x6c, 0x6c, 0x6d, 0x6e, 0x6e, 0x6f, 0x6f, 0x70, 0x71, 0x71, 0x72, 0x73, 0x73, 0x74, 0x74, 0x75, 0x76, 0x76, 0x77, 0x77, 0x78, 0x79, 0x79, 0x7a, 0x7a, 0x7b, 0x7b, 0x7c, 0x7d, 0x7d, 0x7e, 0x7e, 0x7f, 0x7f, 0x80, 0x80, 0x81, 0x81, 0x82, 0x82, 0x83, 0x84, 0x84, 0x85, 0x85, 0x86, 0x86, 0x87, 0x87, 0x88, 0x88, 0x89, 0x89, 0x8a, 0x8a, 0x8b, 0x8b, 0x8c, 0x8c, 0x8c, 0x8d, 0x8d, 0x8e, 0x8e, 0x8f, 0x8f, 0x90, 0x90, 0x91, 0x91, 0x92, 0x92, 0x93, 0x93, 0x93, 0x94, 0x94, 0x95, 0x95, 0x96, 0x96, 0x97, 0x97, 0x97, 0x98, 0x98, 0x99, 0x99, 0x9a, 0x9a, 0x9a, 0x9b, 0x9b, 0x9c, 0x9c, 0x9c, 0x9d, 0x9d, 0x9e, 0x9e, 0x9f, 0x9f, 0x9f, 0xa0, 0xa0, 0xa1, 0xa1, 0xa1, 0xa2, 0xa2, 0xa3, 0xa3, 0xa3, 0xa4, 0xa4, 0xa5, 0xa5, 0xa5, 0xa6, 0xa6, 0xa6, 0xa7, 0xa7, 0xa8, 0xa8, 0xa8, 0xa9, 0xa9, 0xa9, 0xaa, 0xaa, 0xab, 0xab, 0xab, 0xac, 0xac, 0xac, 0xad, 0xad, 0xae, 0xae, 0xae, 0xaf, 0xaf, 0xaf, 0xb0, 0xb0, 0xb0, 0xb1, 0xb1, 0xb1, 0xb2, 0xb2, 0xb3, 0xb3, 0xb3, 0xb4, 0xb4, 0xb4, 0xb5, 0xb5, 0xb5, 0xb6, 0xb6, 0xb6, 0xb7, 0xb7, 0xb7, 0xb8, 0xb8, 0xb8, 0xb9, 0xb9, 0xb9, 0xba, 0xba, 0xba, 0xbb, 0xbb, 0xbb, 0xbc, 0xbc, 0xbc, 0xbd, 0xbd, 0xbd, 0xbe, 0xbe, 0xbe, 0xbf, 0xbf, 0xbf, 0xc0, 0xc0, 0xc0, 0xc1, 0xc1, 0xc1, 0xc1, 0xc2, 0xc2, 0xc2, 0xc3, 0xc3, 0xc3, 0xc4, 0xc4, 0xc4, 0xc5, 0xc5, 0xc5, 0xc6, 0xc6, 0xc6, 0xc6, 0xc7, 0xc7, 0xc7, 0xc8, 0xc8, 0xc8, 0xc9, 0xc9, 0xc9, 0xc9, 0xca, 0xca, 0xca, 0xcb, 0xcb, 0xcb, 0xcc, 0xcc, 0xcc, 0xcc, 0xcd, 0xcd, 0xcd, 0xce, 0xce, 0xce, 0xce, 0xcf, 0xcf, 0xcf, 0xd0, 0xd0, 0xd0, 0xd0, 0xd1, 0xd1, 0xd1, 0xd2, 0xd2, 0xd2, 0xd2, 0xd3, 0xd3, 0xd3, 0xd4, 0xd4, 0xd4, 0xd4, 0xd5, 0xd5, 0xd5, 0xd6, 0xd6, 0xd6, 0xd6, 0xd7, 0xd7, 0xd7, 0xd7, 0xd8, 0xd8, 0xd8, 0xd9, 0xd9, 0xd9, 0xd9, 0xda, 0xda, 0xda, 0xda, 0xdb, 0xdb, 0xdb, 0xdc, 0xdc, 0xdc, 0xdc, 0xdd, 0xdd, 0xdd, 0xdd, 0xde, 0xde, 0xde, 0xde, 0xdf, 0xdf, 0xdf, 0xe0, 0xe0, 0xe0, 0xe0, 0xe1, 0xe1, 0xe1, 0xe1, 0xe2, 0xe2, 0xe2, 0xe2, 0xe3, 0xe3, 0xe3, 0xe3, 0xe4, 0xe4, 0xe4, 0xe4, 0xe5, 0xe5, 0xe5, 0xe5, 0xe6, 0xe6, 0xe6, 0xe6, 0xe7, 0xe7, 0xe7, 0xe7, 0xe8, 0xe8, 0xe8, 0xe8, 0xe9, 0xe9, 0xe9, 0xe9, 0xea, 0xea, 0xea, 0xea, 0xeb, 0xeb, 0xeb, 0xeb, 0xec, 0xec, 0xec, 0xec, 0xed, 0xed, 0xed, 0xed, 0xee, 0xee, 0xee, 0xee, 0xef, 0xef, 0xef, 0xef, 0xef, 0xf0, 0xf0, 0xf0, 0xf0, 0xf1, 0xf1, 0xf1, 0xf1, 0xf2, 0xf2, 0xf2, 0xf2, 0xf3, 0xf3, 0xf3, 0xf3, 0xf3, 0xf4, 0xf4, 0xf4, 0xf4, 0xf5, 0xf5, 0xf5, 0xf5, 0xf6, 0xf6, 0xf6, 0xf6, 0xf6, 0xf7, 0xf7, 0xf7, 0xf7, 0xf8, 0xf8, 0xf8, 0xf8, 0xf9, 0xf9, 0xf9, 0xf9, 0xf9, 0xfa, 0xfa, 0xfa, 0xfa, 0xfb, 0xfb, 0xfb, 0xfb, 0xfb, 0xfc, 0xfc, 0xfc, 0xfc, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfe, 0xfe, 0xfe, 0xfe, 0xff, 0xff, 0xff, }; Again we're going to start with the floating point version as it's easier to reason with. We have a precision P of 9-bits: P = (1<<9)-1 = 511 = 0x1ff. But for the sake of understanding the math, the following diagram will assume a P of 3 so that we can clearly see the segment divisions: The input of our table is an integer index which needs to be calculated according to our input x. But as stated earlier, we won't need one but two indices in order to interpolate a point between 2 discrete values from our table. We will refer to these indices as iₚ and iₙ, which can be computed like this: i = x·P iₚ = ⌊i⌋ iₙ = iₚ + 1 (⌊a⌋ means floor(a)) In order to get an approximation of y according to i, we simply need a linear remapping: the ratio of i between iₚ and iₙ is the same ratio as y between yₚ and yₙ. So yet again we're going to rely on the most useful maths formulas: remap(iₚ,iₙ,yₚ,yₙ,i) = mix(yₚ,yₙ,linear(iₚ,iₙ,i)). The ratio r we're computing as an input to the y-mix can be simplified a bit: r = linear(iₚ,iₙ,i) = (i-iₚ) / (iₙ-iₚ) = i-iₚ = x·P - ⌊x·P⌋ = fract(x·P) So in the end our formula is simply: y = mix(yₚ,yₙ,fract(x·P)) Translated into C we can write it like this: uint8_t linear_f32_to_srgb_u8_fast(float x) { if (x <= 0.f) { return 0; } else if (x >= 1.f) { return 0xff; } else { const float i = x * P; const int32_t idx = (int32_t)floorf(i); const float y0 = linear2srgb[idx]; const float y1 = linear2srgb[idx + 1]; const float r = i - idx; return lrintf(mix(y0, y1, r)); } } Note: in case you are concerned about idx+1 overflowing, floorf((1.0-FLT_EPSILON)*P) is P-1, so this is safe. Linear to sRGB, integer version In the integer version, our function input x is within [0;K], so we need to make a few adjustments. The first issue we have is that with integer arithmetic our i and idx are the same. We have X=x·K as input, so i = idx = X·P/K because we are using an integer division, which in this case is equivalent to the floor() expression in the float version. So while it's a simple and fast way to get yₚ and yₙ, we have an issue figuring out the ratio r. One tool we have is the modulo operator: the integer division is destructive of the fractional part, but fortunately the modulo (the rest of the division) gives this information back. It can also be obtained for free most of the time because CPU division instructions tend to also provide that modulo as well without extra computation. If we give m = (X·P) % K, we have the fractional part of the division expressed in the K scale, which means we can derivate our ratio r from it: r = m / K. Slipping the K division in our mix() expression we end up with the following code: uint8_t linear_int_to_srgb_u8(int32_t x) { if (x <= 0) { return 0; } else if (x >= K) { return 0xff; } else { const int32_t xP = x * P; const int32_t i = xP / K; const int32_t m = xP % K; const int32_t y0 = linear2srgb[i]; const int32_t y1 = linear2srgb[i + 1]; return (m * (y1 - y0) + K/2) / K + y0; } } Testing this function for all the possible input of x, the biggest inaccuracy is a off-by-one, which concerns 6280 of the 65536 possible values (less than 10%): 2886 "off by -1" and 3394 "off by +1". It matches exactly the inaccuracy of the float version of this function, so I think we can be pretty happy with it. Given how good this approach is, we could also consider applying the same strategy for cbrt, so this is left as an exercise to the reader. Back to the core We're finally in our last function. Using everything we've learned so far, it can be trivially converted to integer arithmetic: uint32_t oklab_int_to_srgb_u8(struct LabInt c) { const int64_t l_ = c.L + div_round64(25974LL * c.a, K) + div_round64( 14143LL * c.b, K); const int64_t m_ = c.L + div_round64(-6918LL * c.a, K) + div_round64( -4185LL * c.b, K); const int64_t s_ = c.L + div_round64(-5864LL * c.a, K) + div_round64(-84638LL * c.b, K); const int32_t l = l_*l_*l_ / K2; const int32_t m = m_*m_*m_ / K2; const int32_t s = s_*s_*s_ / K2; const uint8_t r = linear_int_to_srgb_u8((267169LL * l - 216771LL * m + 15137LL * s + K/2) / K); const uint8_t g = linear_int_to_srgb_u8((-83127LL * l + 171030LL * m - 22368LL * s + K/2) / K); const uint8_t b = linear_int_to_srgb_u8(( -275LL * l - 46099LL * m + 111909LL * s + K/2) / K); return r<<16 | g<<8 | b; } Important things to notice: we're storing l_, m_ and s_ in 64-bits values so that the following cubic do not overflow we're using div_round64 for part of the expressions of l_, m_ and s_ because they are using signed sub-expressions we're using a naive integer division in r, g and b because the value is expected to be positive Evaluation We're finally there. In the end the complete code is less than 200 lines of code and even less for the optimized float one (assuming we don't implement our own cbrt). The complete code, test functions and benchmarks tools can be found on Github. Accuracy Comparing the integer version to the reference float gives use the following results: sRGB to OkLab: max_diff=0.000883 total_diff=0.051189 OkLab to sRGB: max_diff_r=2 max_diff_g=1 max_diff_b=1 I find these results pretty decent for an integer version, but you're free to disagree and improve them. Speed The benchmarks are also interesting: on my main workstation (Intel® Core™ i7-12700, glibc 2.36, GCC 12.2.0), the integer arithmetic is slightly slower that the optimized float version: Command Mean [s] Min [s] Max [s] Relative Reference 1.425 ± 0.008 1.414 1.439 1.59 ± 0.01 Fast float 0.897 ± 0.005 0.888 0.902 1.00 Integer arithmetic 0.937 ± 0.006 0.926 0.947 1.04 ± 0.01 Observations: The FPU is definitely fast in modern CPUs Both integer and optimized float versions are destroying the reference code (note that this only because of the transfer functions optimizations, as we have no change in the OkLab functions themselves in the optimized float version) On the other hand, on one of my random ARM board (NanoPI NEO 2 with a Cortex A53, glibc 2.35, GCC 12.1.0), I get different results: Command Mean [s] Min [s] Max [s] Relative Reference 27.678 ± 0.009 27.673 27.703 2.04 ± 0.00 Fast float 15.769 ± 0.001 15.767 15.772 1.16 ± 0.00 Integer arithmetic 13.551 ± 0.001 13.550 13.553 1.00 Not that much faster proportionally speaking, but the integer version is still significantly faster overall on such low-end device. Conclusion This took me ages to complete, way longer than I expected but I'm pretty happy with the end results and with everything I learned in the process. Also, you may have noticed how much I referred to previous work; this has been particularly satisfying from my point of view (re-using previous toolboxes means they were actually useful). This write-up won't be an exception to the rule: in a later article, I will make use of OkLab for another project I've been working on for a while now. See you soon!

over a year ago 58 votes
GCC undefined behaviors are getting wild

Happy with my recent breakthrough in understanding C integer divisions after weeks of struggle, I was minding my own business having fun writing integer arithmetic code. Life was good, when suddenly… zsh: segmentation fault (core dumped). That code wasn't messing with memory much so it was more likely to be a side effect of an arithmetic overflow or something. Using -fsanitize=undefined quickly identified the issue, which confirmed the presence of an integer overflow. The fix was easy but something felt off. I was under the impression my code was robust enough against that kind of honest mistake. Turns out, the protecting condition I had in place should indeed have been enough, so I tried to extract a minimal reproducible case: #include <stdint.h> #include <stdio.h> #include <stdlib.h> uint8_t tab[0x1ff + 1]; uint8_t f(int32_t x) { if (x < 0) return 0; int32_t i = x * 0x1ff / 0xffff; if (i >= 0 && i < sizeof(tab)) { printf("tab[%d] looks safe because %d is between [0;%d[\n", i, i, (int)sizeof(tab)); return tab[i]; } return 0; } int main(int ac, char **av) { return f(atoi(av[1])); } The overflow can happen on x * 0x1ff. Since an integer overflow is undefined, GCC makes the assumption that it cannot happen, ever. In practice in this case it does, but the i >= 0 && i < sizeof(tab) condition should be enough to take care of it, whatever crazy value it becomes, right? Well, I have bad news: % cc -Wall -O2 overflow.c -o overflow && ./overflow 50000000 tab[62183] looks safe because 62183 is between [0;512[ zsh: segmentation fault (core dumped) ./overflow 50000000 Note: this is GCC 12.2.0 on x86-64. We have i=62183 as the result of the overflow, and nevertheless the execution violates the gate condition, spout a non-sense lie, go straight into dereferencing tab, and die miserably. Let's study what GCC is doing here. Firing up Ghidra we observe the following decompiled code: uint8_t f(int x) { int tmp; if (-1 < x) { tmp = x * 0x1ff; if (tmp < 0x1fffe00) { printf("tab[%d] looks safe because %d is between [0;%d[\n",(ulong)(uint)tmp / 0xffff, (ulong)(uint)tmp / 0xffff,0x200); return tab[(int)((uint)tmp / 0xffff)]; } } return '\0'; } When I said GCC makes the assumption that it cannot happen this is what I meant: tmp is not supposed to overflow so part of the condition I had in place was simply removed. More specifically since x can not be lesser than 0, and since GCC assumes a multiplication cannot overflow into a random value (that could be negative) because it is undefined behaviour, it then decides to drop the "redundant" i >= 0 condition because "it cannot happen". I reported that exact issue to GCC to make sure it wasn't a bug, and it was indeed confirmed to me that the undefined behaviour of an integer overflow is not limited in scope to whatever insane value it could take: it is apparently perfectly acceptable to mess up the code flow entirely. While I understand how attractive it can be from an optimization point of view, the paranoid developer in me is straight up terrified by the perspective of a single integer overflow removing security protection and causing such havoc. I've worked several years in a project where the integer overflows were (and probably still are) legion. Identifying and fixing of all them is likely a lifetime mission of several opinionated individuals. I'm expecting this article to make the rust crew go in a crusade again, and I think I might be with them this time. Edit: it was made clear to me while reading Predrag's blog that the key to my misunderstanding boils down to this: "Undefined behavior is not the same as implementation-defined behavior". While I was indeed talking about undefined behaviour, subconsciously I was thinking that the behaviour of an overflow on a multiplication would be "implementation-defined behaviour". This is not the case, it is indeed an undefined behaviour, and yes the compiler is free to do whatever it wants to because it is compliant with the specifications. It's my mistake of course, but to my defense, despite the arrogant comments I read, this confusion happens a lot. This happens I believe because it's violating the Principle of least astonishment. To illustrate this I'll take this interesting old OpenBSD developer blog post being concerned about the result of the multiplication rather than the invalidation of any guarantee with regard to what's going to happen to the execution flow (before and after). This is not uncommon and in my opinion perfectly understandable.

over a year ago 54 votes

More in programming

The Framework Desktop is a beast

I've been running the Framework Desktop for a few months here in Copenhagen now. It's an incredible machine. It's completely quiet, even under heavy, stress-all-cores load. It's tiny too, at just 4.5L of volume, especially compared to my old beautiful but bulky North tower running the 7950X — yet it's faster! And finally, it's simply funky, quirky, and fun! In some ways, the Framework Desktop is a curious machine. Desktop PCs are already very user-repairable! So why is Framework even bringing their talents to this domain? In the laptop realm, they're basically alone with that concept, but in the desktop space, it's rather crowded already. Yet it somehow still makes sense. Partly because Framework has gone with the AMD Ryzen AI Max 395+, which is technically a laptop CPU. You can find it in the ASUS ROG Flow Z13 and the HP ZBook Ultra. Which means it'll fit in a tiny footprint, and Framework apparently just wanted to see what they could do in that form factor. They clearly had fun with it. Look at mine: There are 21 little tiles on the front that you can get in a bunch of different colors or with logos from Framework. Or you can 3D print your own! It's a welcome change in aesthetic from the brushed aluminum or gamer-focused RGBs approach that most of the competition is taking. But let's cut to the benchmarks. That's really why you'd buy a machine like the Framework Desktop. There are significantly cheaper mini PCs available from Beelink and others, but so far, Framework has the only AMD 395+ unit on sale that's completely silent (the GMKTec very much is not, nor is the Z3 Flow). And for me, that's just a dealbreaker. I can't listen to roaring fans anymore. Here's the key benchmark for me: That's the only type of multi-core workload I really sit around waiting on these days, and the Framework Desktop absolutely crushes it. It's almost twice as fast as the Beelink SER8 and still a solid third faster than the Beelink SER9 too. Of course, it's also a lot more expensive, but you're clearly getting some multi-core bang for your buck here! It's even a more dramatic difference to the Macs. It's a solid 40% faster than the M4 Max and 50% faster than the M4 Pro! Now some will say "that's just because Docker is faster on Linux," and they're not entirely wrong. Docker runs natively on Linux, so for this test, where the MySQL/Redis/ElasticSearch data stores run in Docker while Ruby and the app code runs natively, that's part of the answer. Last I checked, it was about 25% of the difference. But so what? Docker is an integral part of the workflow for tons of developers. We use it to be able to run different versions of MySQL, Redis, and ElasticSearch for different applications on the same machine at the same time. You can't really do that without Docker. So this is what Real World benchmarks reveal. It's not just about having a Docker advantage, though. The AMD 395+ is also incredibly potent in RAW CPU performance. Those 16 Zen5 cores are running at 5.1GHz, and in Geekbench 6 multicore, this is how they stack up: Basically matching the M4 Max! And a good chunk faster than the M4 Pro (as well as other AMDs and Intel's 14900K!). No wonder that it's crazy quick with a full-core stress test like running 30,000 assertions for our HEY test suite. To be fair, the M4s are faster in single-core performance. Apple holds the crown there. It's about 20%. And you'll see that in benchmarks like Speedometer, which mostly measures JavaScript single-core performance. The Framework Desktop puts out 670 vs 744 on the M4 Pro on Speedometer 2.1. On SP 3.1, it's an even bigger difference with 35 vs 50. But I've found that all these computers feel fast enough in single-core performance these days. I can't actually feel the difference browsing on a machine that does 670 vs 744 on SP2.1. Hell, I can barely feel the difference between the SER8, which does 506, and the M4 Pro! The only time I actually feel like I'm waiting on anything is in multi-core workloads like the HEY test suite, and here the AMD 395+ is very near the fastest you can get for a consumer desktop machine today at any price. It gets even better when you bring price into the equation, though. The Framework Desktop with 64GB RAM + 2TB NVMe is $1,876. To get a Mac Studio with similar specs — M4 Max, 64GB RAM, 2TB NVMe — you'll literally spend nearly twice as much at $3,299! If you go for 128GB RAM, you'll spend $2,276 on the Framework, but $4,099 on the Mac. And it'll still be way slower for development work using Docker! The Framework Desktop is simply a great deal. Speaking of 64GB vs 128GB, I've been running the 64GB version, and I almost never get anywhere close to the limits. I think the highest I've seen in regular use is about 20GB of RAM in action. Linux is really efficient. Especially when you're using a window manager like Hyprland, as we do in Omarchy. The only reason you really want to go for the full 128GB RAM is to run local LLM models. The AMD 395+ uses unified memory, like Apple, so nearly all of it is addressable to be used by the GPU. That means you can run monster models, like the new 120b gpt-oss from OpenAI. Framework has a video showing them pushing out 40 tokens/second doing just that. That seems about in range of the numbers I've seen from the M4 Max, which also seem in the 40-50 token/second range, but I'll defer to folks who benchmark local LLMs for the exact details on that. I tried running the new gpt-oss-20b on my 64GB machine, though, and I wasn't exactly blown away by the accuracy. In fact, I'd say it was pretty bad. I mean, exceptionally cool that it's doable, but very far off the frontier models we have access to as SaaS. So personally, this isn't yet something I actually use all that much in day-to-day development. I want the best models running at full speed, and right now that means SaaS. So if you just want the best, small computer that runs Linux superbly well out of the box, you should buy the Framework Desktop. It's completely quiet, fantastically fast, and super fun to look at. But I think it's also fair to mention that you can get something like a Beelink SER9 for half the price! Yes, it's also only 2/3 the performance in multi-core, but it's just as fast in single-core. Most developers could totally get away with the SER9, and barely notice what they were missing. But there are just as many people for whom the extra $1,000 is worth the price to run the test suite 40 seconds quicker! You know who you are. Oh, before I close, I also need to mention that this thing is a gaming powerhouse. It basically punches about as hard as an RTX 4060! With an iGPU! That's kinda crazy. Totally new territory on the PC side for integrated graphics. ETA Prime has a video showing the same chip in the GMK Tech running premier games at 1440p High Settings at great frame rates. You can run most games under Linux these days too (thanks Valve and Steam Deck!), but if you need to dual boot with Windows, the dual NVMe slots in the Framework Desktop come very handy. Framework did good with this one. AMD really blew it out of the water with the 395+. We're spoiled to have such incredible hardware available for Linux at such appealing discounts over similar stuff from Cupertino. What a great time to love open source software and tinker-friendly hardware!

21 hours ago 4 votes
Writing: Blog Posts and Songs

I was listening to a podcast interview with the Jackson Browne (American singer/songwriter, political activist, and inductee into the Rock and Roll Hall of Fame) and the interviewer asks him how he approaches writing songs with social commentaries and critiques — something along the lines of: “How do you get from the New York Times headline on a social subject to the emotional heart of a song that matters to each individual?” Browne discusses how if you’re too subtle, people won’t know what you’re talking about. And if you’re too direct, you run the risk of making people feel like they’re being scolded. Here’s what he says about his songwriting: I want this to sound like you and I were drinking in a bar and we’re just talking about what’s going on in the world. Not as if you’re at some elevated place and lecturing people about something they should know about but don’t but [you think] they should care. You have to get to people where [they are, where] they do care and where they do know. I think that’s a great insight for anyone looking to have a connecting, effective voice. I know for me, it’s really easily to slide into a lecturing voice — you “should” do this and you “shouldn’t” do that. But I like Browne’s framing of trying to have an informal, conversational tone that meets people where they are. Like you’re discussing an issue in the bar, rather than listening to a sermon. Chris Coyier is the canonical example of this that comes to mind. I still think of this post from CSS Tricks where Chris talks about how to have submit buttons that go to different URLs: When you submit that form, it’s going to go to the URL /submit. Say you need another submit button that submits to a different URL. It doesn’t matter why. There is always a reason for things. The web is a big place and all that. He doesn’t conjure up some universally-applicable, justified rationale for why he’s sharing this method. Nor is there any pontificating on why this is “good” or “bad”. Instead, like most of Chris’ stuff, I read it as a humble acknowledgement of the practicalities at hand — “Hey, the world is a big place. People have to do crafty things to make their stuff work. And if you’re in that situation, here’s something that might help what ails ya.” I want to work on developing that kind of a voice because I love reading voices like that. Email · Mastodon · Bluesky

2 days ago 4 votes
Doing versus Delegating

A staff+ skill

2 days ago 7 votes
p-fast trie, but smaller

Previously, I wrote some sketchy ideas for what I call a p-fast trie, which is basically a wide fan-out variant of an x-fast trie. It allows you to find the longest matching prefix or nearest predecessor or successor of a query string in a set of names in O(log k) time, where k is the key length. My initial sketch was more complicated and greedy for space than necessary, so here’s a simplified revision. (“p” now stands for prefix.) layout A p-fast trie stores a lexicographically ordered set of names. A name is a sequence of characters from some small-ish character set. For example, DNS names can be represented as a set of about 50 letters, digits, punctuation and escape characters, usually one per byte of name. Names that are arbitrary bit strings can be split into chunks of 6 bits to make a set of 64 characters. Every unique prefix of every name is added to a hash table. An entry in the hash table contains: A shared reference to the closest name lexicographically greater than or equal to the prefix. Multiple hash table entries will refer to the same name. A reference to a name might instead be a reference to a leaf object containing the name. The length of the prefix. To save space, each prefix is not stored separately, but implied by the combination of the closest name and prefix length. A bitmap with one bit per possible character, corresponding to the next character after this prefix. For every other prefix that matches this prefix and is one character longer than this prefix, a bit is set in the bitmap corresponding to the last character of the longer prefix. search The basic algorithm is a longest-prefix match. Look up the query string in the hash table. If there’s a match, great, done. Otherwise proceed by binary chop on the length of the query string. If the prefix isn’t in the hash table, reduce the prefix length and search again. (If the empty prefix isn’t in the hash table then there are no names to find.) If the prefix is in the hash table, check the next character of the query string in the bitmap. If its bit is set, increase the prefix length and search again. Otherwise, this prefix is the answer. predecessor Instead of putting leaf objects in a linked list, we can use a more complicated search algorithm to find names lexicographically closest to the query string. It’s tricky because a longest-prefix match can land in the wrong branch of the implicit trie. Here’s an outline of a predecessor search; successor requires more thought. During the binary chop, when we find a prefix in the hash table, compare the complete query string against the complete name that the hash table entry refers to (the closest name greater than or equal to the common prefix). If the name is greater than the query string we’re in the wrong branch of the trie, so reduce the length of the prefix and search again. Otherwise search the set bits in the bitmap for one corresponding to the greatest character less than the query string’s next character; if there is one remember it and the prefix length. This will be the top of the sub-trie containing the predecessor, unless we find a longer match. If the next character’s bit is set in the bitmap, continue searching with a longer prefix, else stop. When the binary chop has finished, we need to walk down the predecessor sub-trie to find its greatest leaf. This must be done one character at a time – there’s no shortcut. thoughts In my previous note I wondered how the number of search steps in a p-fast trie compares to a qp-trie. I have some old numbers measuring the average depth of binary, 4-bit, 5-bit, 6-bit and 4-bit, 5-bit, dns qp-trie variants. A DNS-trie varies between 7 and 15 deep on average, depending on the data set. The number of steps for a search matches the depth for exact-match lookups, and is up to twice the depth for predecessor searches. A p-fast trie is at most 9 hash table probes for DNS names, and unlikely to be more than 7. I didn’t record the average length of names in my benchmark data sets, but I guess they would be 8–32 characters, meaning 3–5 probes. Which is far fewer than a qp-trie, though I suspect a hash table probe takes more time than chasing a qp-trie pointer. (But this kind of guesstimate is notoriously likely to be wrong!) However, a predecessor search might need 30 probes to walk down the p-fast trie, which I think suggests a linked list of leaf objects is a better option.

2 days ago 4 votes
Software books I wish I could read

New Logic for Programmers Release! v0.11 is now available! This is over 20% longer than v0.10, with a new chapter on code proofs, three chapter overhauls, and more! Full release notes here. Software books I wish I could read I'm writing Logic for Programmers because it's a book I wanted to have ten years ago. I had to learn everything in it the hard way, which is why I'm ensuring that everybody else can learn it the easy way. Books occupy a sort of weird niche in software. We're great at sharing information via blogs and git repos and entire websites. These have many benefits over books: they're free, they're easily accessible, they can be updated quickly, they can even be interactive. But no blog post has influenced me as profoundly as Data and Reality or Making Software. There is no blog or talk about debugging as good as the Debugging book. It might not be anything deeper than "people spend more time per word on writing books than blog posts". I dunno. So here are some other books I wish I could read. I don't think any of them exist yet but it's a big world out there. Also while they're probably best as books, a website or a series of blog posts would be ok too. Everything about Configurations The whole topic of how we configure software, whether by CLI flags, environmental vars, or JSON/YAML/XML/Dhall files. What causes the configuration complexity clock? How do we distinguish between basic, advanced, and developer-only configuration options? When should we disallow configuration? How do we test all possible configurations for correctness? Why do so many widespread outages trace back to misconfiguration, and how do we prevent them? I also want the same for plugin systems. Manifests, permissions, common APIs and architectures, etc. Configuration management is more universal, though, since everybody either uses software with configuration or has made software with configuration. The Big Book of Complicated Data Schemas I guess this would kind of be like Schema.org, except with a lot more on the "why" and not the what. Why is important for the Volcano model to have a "smokingAllowed" field?1 I'd see this less as "here's your guide to putting Volcanos in your database" and more "here's recurring motifs in modeling interesting domains", to help a person see sources of complexity in their own domain. Does something crop up if the references can form a cycle? If a relationship needs to be strictly temporary, or a reference can change type? Bonus: path dependence in data models, where an additional requirement leads to a vastly different ideal data model that a company couldn't do because they made the old model. (This has got to exist, right? Business modeling is a big enough domain that this must exist. Maybe The Essence of Software touches on this? Man I feel bad I haven't read that yet.) Computer Science for Software Engineers Yes, I checked, this book does not exist (though maybe this is the same thing). I don't have any formal software education; everything I know was either self-taught or learned on the job. But it's way easier to learn software engineering that way than computer science. And I bet there's a lot of other engineers in the same boat. This book wouldn't have to be comprehensive or instructive: just enough about each topic to understand why it's an area of study and appreciate how research in it eventually finds its way into practice. MISU Patterns MISU, or "Make Illegal States Unrepresentable", is the idea of designing system invariants in the structure of your data. For example, if a Contact needs at least one of email or phone to be non-null, make it a sum type over EmailContact, PhoneContact, EmailPhoneContact (from this post). MISU is great. Most MISU in the wild look very different than that, though, because the concept of MISU is so broad there's lots of different ways to achieve it. And that means there are "patterns": smart constructors, product types, properly using sets, newtypes to some degree, etc. Some of them are specific to typed FP, while others can be used in even untyped languages. Someone oughta make a pattern book. My one request would be to not give them cutesy names. Do something like the Aarne–Thompson–Uther Index, where items are given names like "Recognition by manner of throwing cakes of different weights into faces of old uncles". Names can come later. The Tools of '25 Not something I'd read, but something to recommend to junior engineers. Starting out it's easy to think the only bit that matters is the language or framework and not realize the enormous amount of surrounding tooling you'll have to learn. This book would cover the basics of tools that enough developers will probably use at some point: git, VSCode, very basic Unix and bash, curl. Maybe the general concepts of tools that appear in every ecosystem, like package managers, build tools, task runners. That might be easier if we specialize this to one particular domain, like webdev or data science. Ideally the book would only have to be updated every five years or so. No LLM stuff because I don't expect the tooling will be stable through 2026, to say nothing of 2030. A History of Obsolete Optimizations Probably better as a really long blog series. Each chapter would be broken up into two parts: A deep dive into a brilliant, elegant, insightful historical optimization designed to work within the constraints of that era's computing technology What we started doing instead, once we had more compute/network/storage available. c.f. A Spellchecker Used to Be a Major Feat of Software Engineering. Bonus topics would be brilliance obsoleted by standardization (like what people did before git and json were universal), optimizations we do today that may not stand the test of time, and optimizations from the past that did. Sphinx Internals I need this. I've spent so much goddamn time digging around in Sphinx and docutils source code I'm gonna throw up. Systems Distributed Talk Today! Online premier's at noon central / 5 PM UTC, here! I'll be hanging out to answer questions and be awkward. You ever watch a recording of your own talk? It's real uncomfortable! In this case because it's a field on one of Volcano's supertypes. I guess schemas gotta follow LSP too ↩

2 days ago 9 votes