Fast HSV to RGB Conversion, Hacker News

Calculation complexity

The conversion between color spaces on small 8-bit CPUs should preferably be done using 8-bit math. The native size is the fastest way for those CPUs to do calculations. The goal is to devise an algorithm that is both as fast and as accurate as possible within the constraints of a small 8-bit CPU.

Using 8-bit math (and a little 16 – bit math) will result is small errors. The first source of error consist of rounding and truncation errors due to finite sizes. The second source can be found in calculation shortcuts. The algorithm must be optimized for small CPUs and that means that some calculations are more difficult than others. The difficult ones can often be replaced if some error is accepted. The trick is to reduce the error to a minimal level.

Mapping bottom level (equation 1)

Equation 1 with scaling can be written as follows:

(level_ {bottom}={255 cdot Big [ {V over 255} cdot Big( 1.0 – {S over 255} Big) Big]}={{V cdot (255 – S)} over {255}} )

The equation can be partially rewritten because (255 – S= overline {S} ) which is the one’s complement of (S ) in 8-bit integer calculation. The multiplication is ( text {8-bit} times text {8-bit}= text {16 – bit} ), with the complication that the maximum values ​​are 255 and the result reaches only a maximum of 0xfe 01. The correct scaling factor in the calculation is (1 / 255 ), but that requires a rather awkward division. A much better division would use a power of two and result in a division by 256: (level_ {bottom}={{V cdot (255 – S)} over {256}} ). A division by 256 can be implemented as a simple shift operation. Using 256 in the division also has the advantage of being a multiple of 8 bits, which in turn is a complete byte. This enables us to ignore the LSByte of the multiplication result completely and just go with the MSByte. However, the result isnot exactand only reaches a level of (255 / 256 ) from the target.

A compensation for the wrong division is to lift the numerator with a (1 / 256 ) part of the original calculation. The wrong division introduces an error factor of (255 / 256 approx 0. 9961 ). The effective result with compensation becomes (255 / 256 255 / (256 cdot 256) approx 0. 9999955 ). A final correction is to lift the rounding error of the division by adding 1 to the nominator. The result has no residual error and maps perfectly (see also appendix aboutfloating point calculation).

(level_ {bottom}={{V cdot ( (- S) {{V cdot) 255 – S)} over {256}} 1} over {256}} )

The equation is easily implemented because the correction is self-similar. The advantage of a self-similar correction is the ability to use an intermediate calculation result.

Equation 1 can be written in C as:

uint8_tv;uint8_ts;  (uint)  _tnumerator;uint8_tlvl_bot;  numerator=v * (uint8_t) (~ s); numerator =numerator>>8; numerator =1; lvl_bot=numerator>>8;

The order of execution of the 1 corrective addition and the self-similar correction turns out to be insignificant. Normally it is important to perform these types of correction in the right order, but there is no difference in this case. The order of the 1 addition is only significant if a) a carry occurs between the two bytes of the 16 – bit multiplication result, and b) this simultaneously results in a secondary byte-level carry which would otherwise not occur. However, this is not the case because the addition that would cause the carry will do so anyway, whether it is done as first or secondary carry. There is code-wise no performance advantage in executing the 1 addition first or second.

Monochromatic values ​​

Monochromatic values ​​in the HSV cylinder are those where the saturation is zero. This is a special case, where a perfect mapping can be achieved by bypassing the calculation and assigning the RGB values ​​to be equal to the V parameter. Even though equation 1 maps perfectly, the other equations might not. Any set of HSV values ​​where (S=0 ) are trivially detected and can simply be handled separately.

(Mapping top level) 2)

Equation 2 is simply a one-to-one mapping of the (V ) parameter scaled at 8-bit unsigned integer level. Equation 2 can be written in C as:

uint8_tv;uint8_tlvl_top;  lvl_top=v;

There is obviously no error in the resulting mapping of values ​​for the top-level.

Mapping slopes (equation 3 and 4)

The slope equations are the only ones that have the hue as a parameter. The hue is the angle on the cylinder, which is divided into 6 sextants. The normalization of the angle ensures (0.0 le H _ { Delta} lt 1.0 ). Please note that the normalized angle never reaches one. Parameter (V ) and (S ) are already mapped to scale to 255. The scale of (H_ Delta ) is done in a similar way as (S ) and (V ), with the hue sextant mapping to fit into an 8-bit unsigned integer value, using (0 le H _ { Delta} lt 256 ). The mapping implies that 1.0 is mapped to 256 Because (H_ Delta ) never reaches one and therefore always fits in an 8-bit unsigned integer:

( begin {align}  slope_ {down} &=255 cdot Big [ {V over 255} cdot Big(1.0 – {S over 255} cdot {H_{Delta} over 256}Big) Big]  ={{V cdot Big (255 – {{ S cdot H _ { Delta}} over {256} } Big)} over {255}} \  & \  slope_ {up} &=255 cdot Big [ {V over 255} cdot Big(1.0 – {S over 255} cdot {{256 – H_{Delta}} over 256}Big) Big]  ={{V cdot Big (255 – {{ S cdot (256 – H _ { Delta})} over {256}} Big)} over {255}} end {align} )

The primary problem in these equations is the larger error when the division is changed to (1 / 256 ) instead of (1 / 255 ). This is because there are two multiplications and divisions in sequence, which cause a compounded calculation error. The resulting error would be in the /- 2 range and not in the /- 1 range as preferred if unchanged.

Rescaling of any of the parameters (V ), (S ) or (H _ { Delta} ) is not really an option because scaling to other values ​​than 0 … 255 would imply that an 8-bit type cannot hold the magnitude of the parameter. Therefore, a similar type of compensation as in the bottom level equation must be employed to balance the equations.

The first compensation is the same as in equation 1, where a self-similar correction is employed to compensate for the (255 / 256 ) factor. The second correction targets the compounded multiplication and lifts the numerator by a fraction of (V / 2 ). This second compensation balances the rounding error to be distributed evenly between under- and over-estimation by /- 1:

( begin {align}  slope_ {down} &={{V cdot Big (255 – {{S cdot H _ { Delta} {{S cdot H _ { Delta}} over {256}}} over {256}} Big) {V over 2}} over {256}} \  & \  slope_ {up} &={{V cdot Big (255 – {{S cdot (256 – H _ { Delta}) {{S cdot) 256 – H_ { Delta})} over {256}}} over {256}} Big) {V over 2}} over {256}} end {align} )

The computation (S cdot 256 – H _ { Delta }) ) can be shortcut to prevent an expensive multiplication. The problem is appears when (H _ { Delta}=0 ), where the subtraction results in a 9-bit value (0x 100). All other subtraction results are 8-bit in size. However, the case of (H _ { Delta}=0 ) is trivially detected and a 16 – bit multiplication can be prevented (S ⋅ 256=S

As before, the calculation (255 – X ) is replaced with the one’s complement. The calculation (256 – X ) is replaced with the two’s complement equivalent using 8-bit calculation, resulting in (- X ). Finally, the division by 256 is done by shifting. The resulting code is then:

(uint8_t)  v;uint8_ts;  (uint)  _th;  (uint)  _tnumerator;uint8_tslope_dn;uint8_tslope_up;  numerator=s * h; numerator =numerator>>8; numerator=v * (uint8_t) (~ (uint8_t) (numerator>>8)); numerator =v>>1; slope_dn=numerator>>8;  numerator=! h? (( (uint)  _t) suint8_t) (- h)); numerator =numerator>>8; numerator=v * (uint8_t) (~ (uint8_t) (numerator>>8)); numerator =v>>1; slope_up=numerator>>8;

An error analysis of the slopes reveals the following maps:

  

  

All values ​​have a maximum error of /- 1. The error of the slope down image is 6. 00% at the -1 level and 6. 16% at the 1 level. The error of the slope up image is 6. 36% at -1 and 6. 14% at 1. The map reveals that high value (V) will have the highest error concentrations. These errors are primarily in the more saturated color values ​​and are usually in the bright region of intensities. Therefore, it should not result in any significant visual artifacts.

The up- and down-slopes differ and seem to be moved /- 1 for the (H _ { Delta} ) parameter. The reason for this behavior is that the slopes are calculated from each end ( (H _ { Delta} | _ {0 rightarrow 1} ) versus (H _ { Delta} | _ {1 rightarrow 0} )). The resulting shift in calculation causes (H _ { Delta} ) to vary between 0 … 255 for the down-slope and 256 … 1 for the up-slope, which moves the equations slightly with respect to each other (a permanent fix is ​​to use thecalculation with larger resolution).

The artifact for (H _ { Delta}=0 ) in the down-slope, which is consistently at -1, can be compensated if necessary. However, fixing it requires an awkward conditional in the code that makes it run slower. The reason for the -1 error can be found in the fact that the compensation term (S cdot H _ { Delta} / 256 ) is zero for specifically this condition and the correct secondary compensation would then be (V ) instead of (V / 2 ). Therefore, the output is not compensated for this eventuality.

Hue sextants

The scaling of the hue to (6 cdot 256 ) allows for 8-bit selection of the sextant. Each 60 degree sextant is entirely located in the LSByte, whereas the sextant selection is entirely located in the MSByte. The RGB values ​​are assigned from the calculations in a specific order determined by the sextant (see (figure 2).

uint8_tsextant=h>>8;uint8_th_fraction=h & 0xff;  (if)  (sextantelse (if)  (sextant==1) {  }else{  } }else{   (if)  (sextant==4) {  }else(if)  (sextant==5) {  }else{  } }

The calculation’s hue value is based onh_fraction, which is the LSByte, and all comparisons are based in the MSByte.

Another advantage of the (6 cdot 256 ) hue scaling model is that the slope selection is rather easy. The least significant bit of the sextant indicates directly which slope calculation is relevant. All even sextants are slope-up and all odd sextants are slope-down:

uint8_tslope;  (if)  (sextant & 1) { slope=... slope_down ... }else{ slope=... slope_up ... }

A different construct is also possible. It is possible to bypass sextant comparisons and associated assignments and only use one fixed set of assignments. The sextant can be used to mangle the RGB return values. If the function uses pointers, then the sextant bits indicate how to exchange RGB pointers. The assignments are stable, but are mapped to the appropriate target color. The following table shows how up / down slopes (u, d) top level (v) and bottom level (c) are distributed and exchanged on the basis of the sextant bits:

(RGB) (R⇔B) (G⇔B) (R⇔G)

(0 0 1)

(uvc)

[ {256.0 – {H_{Delta}}} big]

(1 0 0)

()

(uvc)

Pointer swapping is done not at all (sextant 1), once (sextants 0, 2, 4) or twice (sextants 3, 5). Combining this idea with the previous code-snippet and completing an outline for the function could look like:

void  hsv2rgb ( (uint)  _th,uint8_t  (s,uint8_tV,uint8_t* r,uint8_t* g,uint8_t* b) {   (if)  (! s) { * r=* g=* b=v;return; }uint8_tsextant=h>>8;    (if)  (sextant>5) { sextant=5; }    (if)  (sextant & 2) { swapptr (r, b); }   (if)  (sextant & 4) { swapptr (g, b); }   (if)  (! (sextant & 6) {  (if)  (! (sextant & 1)) { swapptr (r, g); } }else{  (if)  (sextant & 1) { swapptr (r, g); } }   * g=v; * b=... lvl_bot ...;uint8_th_fraction=h & 0xff;  (if)  (! (sextant & 1)) { * r=... slope_up ...; }else{ * r=... slope_down ...; } }

An alternative to pointer swapping is assigning local variables to the calculated values ​​and then swap the values ​​accordingly. The advantage of such implementation is that no pointers need to be passed and a 32 – bit return-value could encapsulate the calculated RGB values. A reduced number of arguments is especially advantageous to reduce the call-setup time. A possible drawback is the extent for which intermediate values ​​must be held in registers. Anyhow, the caller must still handle the returned value (s) and act accordingly. It is very implementation specific which strategy ultimately will be better.

The setup again uses (0 le H _ { Delta} lt 256 ), (0 le S le 255 ) and (0 le V le 255 ). The equations are altered in a fashion to obtain balance in the rounding error as good as possible. That means that the order of the calculation is controlled more tightly to match the proper slope in both directions. The new slope calculations do compounded scaling divisions in one operation, thereby reducing the intermediate rounding error and keeping intermediate fractional results intact. This allows for a better result by effectively using a fixed-point calculation with higher resolution. Rewriting the equations for larger multiplication by eliminating the intermediate divisions:

( begin {align}  slope_ {down} &={{V cdot Big (255 – {{S cdot H _ { Delta}} over {256}} Big)} over {255}}  ={{V cdot Big (255 cdot 256) – {S cdot H _ { Delta}} Big)} over {255 cdot 256}} \  & \  slope_ {up} &={{V cdot Big (255 – {{S cdot (256 – H _ { Delta})} over {256}} Big)} over {255}}  ={{V cdot Big (255 cdot 256 – {S cdot (256 – H _ { Delta})} Big)} over {255 cdot 256}} end {align} )

Replacing the scaling factor (1 / 255 ) again with (1 / 256 ) for the final division means a similar correction as previously performed using a self-similar correction. The final division can then be performed using a 16 – bit shift operation on the nominator:

( begin {align} slope_ {down} &={{V cdot big (255 cdot 256) – {S cdot H_ { Delta}} big) {{V cdot big (255 cdot 256) – {S cdot H _ { Delta}} big)} over {256}} V} over {256 cdot 256}} \  & \  slope_ {up} &={{V cdot big (255 cdot 256) – {S cdot (256 – H _ { Delta})} big) {{ V cdot big ((255 cdot 256) – {S cdot (256 – H _ { Delta})} big)} over {256}} V} over {256 cdot 256}} end {align} )

The bottom level equation (equation 1) is unaltered with respect to the

Both slope equations use the same correction is before with the exception of a larger final addition before the division. The 3D slope of the equations is compensated by adding (V ) instead of (V / 2 ), to move the rounding error to be better balanced. It also assures a simpler calculation. The origin of this change can be found in the fix-point resolution, which is higher in this version and therefore slopes differently.

There seems to be an obvious simplification to the equations by taking the individual (V ) multiplications out in the numerator and calculate an intermediate result in parentheses before multiplication by (V ), like in: ((V cdot X) (V cdot X / 256) V=V cdot (X X / 256 1) ). However, this line of thought is flawed and fails because it effectively reduces the fix-point resolution by a factor of 2 (8) , thereby missing the target by a reasonable margin. Nor is there any significant calculation benefit in taking out (V ) in the multiplication, in terms of time, because the calculation is self-similar and only consists of additions and a shift operation.

The calculation can be written as:

uint8_ts;uint8_tv;uint8_th;   (uint) ******************************************************************************************************************************************************************************************************************************************** (_t)  d;uint8_tlvl_bot;uint8_tslope_up;uint8_tslope_dn;  d=(v * (uint8_t) (~ s)); d =d>>8; lvl_bot=d>>8;  d=s * h; D=(255>8; d =v; slope_dn=d>>16;  d=s * (256 - h); D=(255>8; d =v; slope_up=d>>16;

It should be noted that not all of the intermediate results are required to be written explicitly. Only the compensation term needs to be accessible for fast calculations. The compiler will use the CPU’s native size by default. This also means that, when reading above code,you mustknowhow and whenthe compiler converts its calculations to the CPU’s native types.

An error analysis of the slopes reveals the following maps:

  

  

The errors are all zero or 1. There are only very few errors in the map. The slope-down and slope-up have 0. 03% errors each. These are rather insignificant numbers.

HSV generic performance

The average performance of the implementations have been tested by calculating every possible conversion for the integer calculations (256 ⋅ 256 ⋅ 256 ⋅6, or ≈ 10 (8) ************************************************************************************************* (tests). The floating point implementation was reduced by about a factor of 25 to speedup the test (52 ⋅ 52 ⋅ 256 ⋅6 or ≈ 4⋅ 10 (6) tests). All tests were performed on an Arduino Uno R3, ATmega (P @) ******************************************************************************************************************************************************************************************************************************************************* (MHz.)

Test and call-overhead has been isolated by creating an empty function with the same prototype. The empty function only consists of a return. Therefore, all call-setup and test-related counters are encapsulated within this test. All performance-test values ​​have been subtracted this overhead to see the actual calculation timing for the algorithms as implemented. The non-integer clock-cycles count is due to conditional branching, which are not necessarily balanced to be an integer multiple when calculating averages.

The code is tested with both a hardware multiplier (HW mul) and a software implemented multiplication (SW mul). The C-compiler uses the default for the CPU, which is hardware multiplication.

Relative timing of the function body (actual calculation), sorted by timing:

It is clear that the floating point implementation is rather naive, ranging from 5.6 … 28 .8 times slower than any integer implementation. An interesting observation is that the 8-bit and 32 – bit implementations are relatively close (about a factor of two). This means that one probably could suffice with a pure C implementation most of the time and use generic portable code. Only timing constrained systems would need the assembly implementations. There is only about 12% difference between the 8-bit and 32 -bit assembly versions with hardware multiplication. Therefore, it is usually most logical to use the 32 – bit version because it is more accurate in the calculation. Only very tight systems, timing-wise, would require the 8-bit version.

The compiler (avr-gcc 4.9.3) is not very good at optimizing these specific calculations. The C-code, as it is implemented, underwent several iterations to improve on the compiler’s performance. The primary reason for the poorer performance is the C-language itself because the compiler does not know the range of values ​​it is handling and only knows about types. Therefore, shortcuts in the calculation (for example multiplying with values ​​≤0x 100) cannot be optimized optimally using a generic compiler.

HSV specific performance

Above tests show the performance of the generic algorithm averaged over the entire conversion space. This is not a real-world test. The most common, and therefore important conversions, are those where the saturation is at maximum (full color) or at minimum (gray-scale). A second test was performed to examine these conversions more closely.

The average calculations using software multiplication differ significantly between the generic case ( (0 le S le 255 )) and the specific cases ( (S=0 ) and (S=255 )). The differences can be found in the multiplication algorithm and the specific values. The multiplication uses an early-out algorithm and the number of ones in the multiplicands are significant for how long the algorithm loops. The case where S=0 prevents any multiplication (pure gray-scale) and (V ) is assigned directly to the RGB constituents. The case S=255 is reduced in complexity because thebottom leveluses (V cdot (255 – S) ), which means that one multiplicand becomes zero, thus performing a very fast multiplication.

Another interesting observation, for these specific conversions with S=255, is that the 8-bit assembly implementation with software multiplication is about as fast as the 32 – bit C version using hardware multiplications. The accuracy of both differ, but if speed is the motto, then it is a good option for the very small micro controllers.

There are timing differences in the calculations due to the special handling of the sextants. The assignments are swapped zero, once or twice, depending on the sextant (seeHue sextantsabove). However, the differences are rather small as indicated below:

The timing differences reflect the number of pointer-swaps that must be performed, except for sextants 2 … 6 in the 32 -bit assembly implementation. The branching in the 32 – bit assembly code is balanced between the pointer -swaps and the slope-up / -down calculation. This was no intentional coding, but simply turned out be be this way.  

HSV performance on ARM

No actual tests have been performed on ARM micro controllers. An inspection of the generated assembly from the C-code did reveal that the 32 -bit version performs better or as well as the 8-bit AVR assembly version. Estimating ARM performance is much more difficult because the CPU is pipelined and can execute more than one instruction at any given time, depending register scheduling and may apply out-of-order completion (Cortex-M7 or the big SOCs). Even the lower-end CPUs (M0, M1, M2, M3 and M4) have decoupled memory access and internal CPU operation by means of D- and I-cache. And then there is the C-compiler, which is actually very good at doing register scheduling (and most of the time better than humans).

That said, the most important factor for performance is the hardware multiplier, which can be a slow 32 – cycle version in M0 , M1 and M2, depending the chip-vendor’s implementation, or the fast 1-cycle version (also used in all M3 and better). You must check the vendor’s specification to see which version was implemented in the M [012] silicon.

Baseline performance of the 32 – bit HSV to RGB algorithm in C , based on ARM’s assembly inspection, is expected to follow the operating frequency of the core proportionally to the AVR’s performance above. A 32 MHz core should be able to execute the conversion in less than 2 µs and a MHz core should be able to do it in under 1 µs.

To put these numbers in perspective, a PC with an Intel Core-i5 at 2.4 GHz (64 – bit Linux OS (runs the 32 – bit conversion (in C) on average, over the entire conversion space, in 6. 39 ns per conversion including test- and call-overhead. That is more than three orders of magnitude faster than the AVR core. Rescaling the operating frequency to an Arduino at 16 MHz results in 958 ns per conversion, or still an order of magnitude faster than an AVR core. It is to be expected that an ARM is, silicon-wise, somewhere in between these extremes. Therefore, the above estimate can be considered valid, or maybe even slightly conservative.