<p>I have converted the algorithm for calculating the HealPix value from Lat/Lon values from C++ to GLSL. The result is saved in a BMP image and looks correct. Both algorithms (C++ and GLSL) use the same primitive types with the same program code. See figure 1.</p>
<p>I noticed that there are artifacts in GLSL, which you can see in Figure 2.
These artifacts only occur in the GLSL version. In the C++ version, there are no artifacts at all.
As you can see, I was initially concerned about accuracy and therefore added <code>highp</code> to all my floats. But nothing has changed.</p>
<p><a href="https://i.sstatic.net/TpBUvITJ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/TpBUvITJ.png" alt="enter image description here" /></a></p>
<p><em>Figure 1: Results of the HealPix to Lat/Lon algorithm in C++ and Vulkan(GLSL)</em></p>
<p><a href="https://i.sstatic.net/JFmL4A2C.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JFmL4A2C.png" alt="enter image description here" /></a></p>
<p><em>Figure 2: Enlarged view of the result. Here you can see the artifacts</em></p>
<p>Here is the GLSL code:</p>
<pre><code>#version 450
layout (local_size_x = 16, local_size_y = 8, local_size_z = 1) in;

const highp float HP_levelSetPow[14] = { 8192.0, 4096.0, 2048.0, 1024.0, 512.0, 256.0, 
    128.0, 64.0, 32.0, 16.0, 8.0, 4.0, 2.0, 1.0 };

const uint HP_levelSet[14] = { 26, 24, 22, 20, 18, 16, 14, 12, 10, 8, 6, 4, 2, 0};

const int HP_x2pix[128] = {0, 1, 4, 5, 16, 17, 20, 21, 64, 65, 68, 69, 80, 81, 84, 
    85, 256, 257, 260, 261, 272, 273, 276, 277, 320, 321, 324, 325, 336, 337, 340, 
    341, 1024, 1025, 1028, 1029, 1040, 1041, 1044, 1045, 1088, 1089, 1092, 1093, 
    1104, 1105, 1108, 1109, 1280, 1281, 1284, 1285, 1296, 1297, 1300, 1301, 1344, 
    1345, 1348, 1349, 1360, 1361, 1364, 1365, 4096, 4097, 4100, 4101, 4112, 4113, 
    4116, 4117, 4160, 4161, 4164, 4165, 4176, 4177, 4180, 4181, 4352, 4353, 4356, 
    4357, 4368, 4369, 4372, 4373, 4416, 4417, 4420, 4421, 4432, 4433, 4436, 4437, 
    5120, 5121, 5124, 5125, 5136, 5137, 5140, 5141, 5184, 5185, 5188, 5189, 5200, 
    5201, 5204, 5205, 5376, 5377, 5380, 5381, 5392, 5393, 5396, 5397, 5440, 5441, 
    5444, 5445, 5456, 5457, 5460, 5461};
const int HP_y2pix[128] = {0, 2, 8, 10, 32, 34, 40, 42, 128, 130, 136, 138, 160, 
    162, 168, 170, 512, 514, 520, 522, 544, 546, 552, 554, 640, 642, 648, 650, 672,
    674, 680, 682, 2048, 2050, 2056, 2058, 2080, 2082, 2088, 2090, 2176, 2178, 2184,
    2186, 2208, 2210, 2216, 2218, 2560, 2562, 2568, 2570, 2592, 2594, 2600, 2602, 
    2688, 2690, 2696, 2698, 2720, 2722, 2728, 2730, 8192, 8194, 8200, 8202, 8224, 
    8226, 8232, 8234, 8320, 8322, 8328, 8330, 8352, 8354, 8360, 8362, 8704, 8706, 
    8712, 8714, 8736, 8738, 8744, 8746, 8832, 8834, 8840, 8842, 8864, 8866, 8872, 
    8874, 10240, 10242, 10248, 10250, 10272, 10274, 10280, 10282, 10368, 10370, 
    10376, 10378, 10400, 10402, 10408, 10410, 10752, 10754, 10760, 10762, 10784,
    10786, 10792, 10794, 10880, 10882, 10888, 10890, 10912, 10914, 10920, 10922};

//theta: range [0, Pi)
//phi: range [0, 2Pi)
uint HP_latLonToHealPix(highp float theta, highp float phi) 
{       
    highp float z, za, z0, tt, tp, tmp;
    int    face_num, jp, jm;
    int   ifp, ifm;
    int    ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
    highp float piover2 = 0.5*3.1415926, twopi = 2.0*3.1415926;
    int    ns_max = 8192;

    z = cos(theta);
    za = abs(z);
    z0 = 2.0 / 3.0; 
    tt = phi / piover2; /* in [0,4[ */

    if (za &lt;= z0) 
    { /* equatorial region */

                    /* (the index of edge lines increase when the longitude=phi goes up) */
        jp = int(floor(8192*(0.5 + tt - z*0.75))); /* ascending edge line index */
        jm = int(floor(8192*(0.5 + tt + z*0.75))); /* descending edge line index */

                                                     /* finds the face */
        ifp = jp / 8192; /* in {0,4} */
        ifm = jm / 8192;

        if (ifp == ifm) face_num = int(mod(ifp, 4)) + 4; /* faces 4 to 7 */
        else if (ifp&lt;ifm) face_num = int(mod(ifp, 4)); /* (half-)faces 0 to 3 */
        else face_num = int(mod(ifm, 4)) + 8;           /* (half-)faces 8 to 11 */

        ix = int(mod(jm, 8192));
        iy = 8191 - int(mod(jp, 8192));
    }
    else 
    { /* polar region, za &gt; 2/3 */

        ntt = clamp(int(floor(tt)), 0, 3);
        tp = tt - ntt;
        tmp = sqrt(3.0*(1.0 - za)); /* in ]0,1] */

                                  /* (the index of edge lines increase when distance from the closest pole
                                  * goes up)
                                  */
                                  /* line going toward the pole as phi increases */
        jp = int(floor(8192 * tp * tmp));

        /* that one goes away of the closest pole */
        jm = int(floor(8192 * (1.0 - tp) * tmp));
        jp = clamp(jp, 0, 8191);
        jm = clamp(jm, 0, 8191);

        /* finds the face and pixel's (x,y) */
        if (z &gt;= 0) {
            face_num = ntt; /* in {0,3} */
            ix = 8191 - jm;
            iy = 8191 - jp;
        }
        else {
            face_num = ntt + 8; /* in {8,11} */
            ix = jp;
            iy = jm;
        }
    }

    ix_low = int(mod(ix, 128));
    ix_hi = ix / 128;
    iy_low = int(mod(iy, 128));
    iy_hi = iy / 128;

    ipf = (HP_x2pix[ix_hi] + HP_y2pix[iy_hi]) * (128 * 128) + (HP_x2pix[ix_low] + HP_y2pix[iy_low]);
    ipf = int((ipf / pow(ns_max / 8192, 2)));     /* in {0, nside**2 - 1} */
    return uint(ipf + face_num*pow(8192, 2)); /* in {0, 12*nside**2 - 1} */ 
}

uint HP_ignoreEverySecondBit(uint value, bool odd, uint countBits)
{
    uint result = 0u;
    uint mask = odd ? 1u : 2u;
    countBits = countBits / 2;
    for (uint i = 0; i &lt; countBits; ++i)
    {       
        result += clamp((value &amp; mask), 0u, 1u) * uint(pow(2u, i));     
        value = value &gt;&gt; 2;
    }
    return result;
}

highp float HP_getU(uint healPixNumber, uint level)
{
    return float(HP_ignoreEverySecondBit(healPixNumber, true, HP_levelSet[level])) / HP_levelSetPow[level];
}

highp float HP_getV(uint healPixNumber, uint level)
{
    return float(HP_ignoreEverySecondBit(healPixNumber, false, HP_levelSet[level])) / HP_levelSetPow[level];
}

struct healPixData
{
    uint healPixNumber;
    uint level;
    highp float u;
    highp float v;
};

layout(binding = 0) uniform UniformBufferObject 
{
    uint level;
    uint scale;
} ubo;

layout(std430, binding = 2) buffer SSBO2
{
    healPixData valuesSSBO2[];
};

void main()
{
    uint width = 16 * ubo.scale;
    uint height = 8 * ubo.scale;
    uint arrayIndex = gl_GlobalInvocationID.x * height + height - gl_GlobalInvocationID.y - 1;
    
    vec2 currentIndex = vec2(gl_GlobalInvocationID.xy);
    vec2 latLon = currentIndex / vec2(width, height) * vec2(2 * 3.1415926, 3.1415926);
    
    healPixData data;
    data.healPixNumber = HP_latLonToHealPix(latLon.y, latLon.x); 
    data.u = HP_getU(data.healPixNumber, ubo.level);
    data.v = HP_getV(data.healPixNumber, ubo.level);
    
    valuesSSBO2[arrayIndex] = data;     
}
</code></pre>
<p>There are no race conditions that can cause artifacts. The artifacts do not change, there is no randomness. Any idea where and why this happened?</p>

CC BY-SA 4.0
2024-07-16T17:22:44.830Z
I don't know the cause, but can suggest some ideas for debugging.  You know the (x,y) of an artifact pixel.  Add shader [debug printf](https://github.com/KhronosGroup/Vulkan-ValidationLayers/blob/ac2e01fc1a67329da0c879ecd2296b76643c49c0/docs/debug_printf.md) statements that print out some of the intermediate values only for that (x,y) and compare them with good values from the C++ code.  In addition, I see some usages of `pow(x,2)` that could be changed to `(x * x)`.  It is possible that the shader compiler doesn't recognize the shortcut and uses logs to compute the function.
Karl Schultz
CC BY-SA 4.0
2024-07-17T07:09:07.913Z
If that can help: it seems that the artifacts are aligned horizontally, with equal spacing on a given row. Spacing may differ from row to row. It is not impossible that the rows are complete with some artifacts so little contrasted that one does not see them.
Yves Daoust
<glsl><vulkan><artifacts>
HealPix: Different results when comparing C++ and Vulkan
Thomas
# HealPix: Different results when comparing C++ and Vulkan

(Asked by: username_0 on 2024-07-14)

I have converted the algorithm for calculating the HealPix value from Lat/Lon values from C++ to GLSL. The result is saved in a BMP image and looks correct. Both algorithms (C++ and GLSL) use the same primitive types with the same program code. See figure 1.

I noticed that there are artifacts in GLSL, which you can see in Figure 2.
These artifacts only occur in the GLSL version. In the C++ version, there are no artifacts at all.
As you can see, I was initially concerned about accuracy and therefore added `highp` to all my floats. But nothing has changed.

[![enter image description here](https://i.sstatic.net/TpBUvITJ.png)](https://i.sstatic.net/TpBUvITJ.png)

*Figure 1: Results of the HealPix to Lat/Lon algorithm in C++ and Vulkan(GLSL)*

[![enter image description here](https://i.sstatic.net/JFmL4A2C.png)](https://i.sstatic.net/JFmL4A2C.png)

*Figure 2: Enlarged view of the result. Here you can see the artifacts*

Here is the GLSL code:

```
#version 450
layout (local_size_x = 16, local_size_y = 8, local_size_z = 1) in;

const highp float HP_levelSetPow[14] = { 8192.0, 4096.0, 2048.0, 1024.0, 512.0, 256.0, 
    128.0, 64.0, 32.0, 16.0, 8.0, 4.0, 2.0, 1.0 };

const uint HP_levelSet[14] = { 26, 24, 22, 20, 18, 16, 14, 12, 10, 8, 6, 4, 2, 0};

const int HP_x2pix[128] = {0, 1, 4, 5, 16, 17, 20, 21, 64, 65, 68, 69, 80, 81, 84, 
    85, 256, 257, 260, 261, 272, 273, 276, 277, 320, 321, 324, 325, 336, 337, 340, 
    341, 1024, 1025, 1028, 1029, 1040, 1041, 1044, 1045, 1088, 1089, 1092, 1093, 
    1104, 1105, 1108, 1109, 1280, 1281, 1284, 1285, 1296, 1297, 1300, 1301, 1344, 
    1345, 1348, 1349, 1360, 1361, 1364, 1365, 4096, 4097, 4100, 4101, 4112, 4113, 
    4116, 4117, 4160, 4161, 4164, 4165, 4176, 4177, 4180, 4181, 4352, 4353, 4356, 
    4357, 4368, 4369, 4372, 4373, 4416, 4417, 4420, 4421, 4432, 4433, 4436, 4437, 
    5120, 5121, 5124, 5125, 5136, 5137, 5140, 5141, 5184, 5185, 5188, 5189, 5200, 
    5201, 5204, 5205, 5376, 5377, 5380, 5381, 5392, 5393, 5396, 5397, 5440, 5441, 
    5444, 5445, 5456, 5457, 5460, 5461};
const int HP_y2pix[128] = {0, 2, 8, 10, 32, 34, 40, 42, 128, 130, 136, 138, 160, 
    162, 168, 170, 512, 514, 520, 522, 544, 546, 552, 554, 640, 642, 648, 650, 672,
    674, 680, 682, 2048, 2050, 2056, 2058, 2080, 2082, 2088, 2090, 2176, 2178, 2184,
    2186, 2208, 2210, 2216, 2218, 2560, 2562, 2568, 2570, 2592, 2594, 2600, 2602, 
    2688, 2690, 2696, 2698, 2720, 2722, 2728, 2730, 8192, 8194, 8200, 8202, 8224, 
    8226, 8232, 8234, 8320, 8322, 8328, 8330, 8352, 8354, 8360, 8362, 8704, 8706, 
    8712, 8714, 8736, 8738, 8744, 8746, 8832, 8834, 8840, 8842, 8864, 8866, 8872, 
    8874, 10240, 10242, 10248, 10250, 10272, 10274, 10280, 10282, 10368, 10370, 
    10376, 10378, 10400, 10402, 10408, 10410, 10752, 10754, 10760, 10762, 10784,
    10786, 10792, 10794, 10880, 10882, 10888, 10890, 10912, 10914, 10920, 10922};

//theta: range [0, Pi)
//phi: range [0, 2Pi)
uint HP_latLonToHealPix(highp float theta, highp float phi) 
{       
    highp float z, za, z0, tt, tp, tmp;
    int    face_num, jp, jm;
    int   ifp, ifm;
    int    ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
    highp float piover2 = 0.5*3.1415926, twopi = 2.0*3.1415926;
    int    ns_max = 8192;

    z = cos(theta);
    za = abs(z);
    z0 = 2.0 / 3.0; 
    tt = phi / piover2; /* in [0,4[ */

    if (za <= z0) 
    { /* equatorial region */

                    /* (the index of edge lines increase when the longitude=phi goes up) */
        jp = int(floor(8192*(0.5 + tt - z*0.75))); /* ascending edge line index */
        jm = int(floor(8192*(0.5 + tt + z*0.75))); /* descending edge line index */

                                                     /* finds the face */
        ifp = jp / 8192; /* in {0,4} */
        ifm = jm / 8192;

        if (ifp == ifm) face_num = int(mod(ifp, 4)) + 4; /* faces 4 to 7 */
        else if (ifp<ifm) face_num = int(mod(ifp, 4)); /* (half-)faces 0 to 3 */
        else face_num = int(mod(ifm, 4)) + 8;           /* (half-)faces 8 to 11 */

        ix = int(mod(jm, 8192));
        iy = 8191 - int(mod(jp, 8192));
    }
    else 
    { /* polar region, za > 2/3 */

        ntt = clamp(int(floor(tt)), 0, 3);
        tp = tt - ntt;
        tmp = sqrt(3.0*(1.0 - za)); /* in ]0,1] */

                                  /* (the index of edge lines increase when distance from the closest pole
                                  * goes up)
                                  */
                                  /* line going toward the pole as phi increases */
        jp = int(floor(8192 * tp * tmp));

        /* that one goes away of the closest pole */
        jm = int(floor(8192 * (1.0 - tp) * tmp));
        jp = clamp(jp, 0, 8191);
        jm = clamp(jm, 0, 8191);

        /* finds the face and pixel's (x,y) */
        if (z >= 0) {
            face_num = ntt; /* in {0,3} */
            ix = 8191 - jm;
            iy = 8191 - jp;
        }
        else {
            face_num = ntt + 8; /* in {8,11} */
            ix = jp;
            iy = jm;
        }
    }

    ix_low = int(mod(ix, 128));
    ix_hi = ix / 128;
    iy_low = int(mod(iy, 128));
    iy_hi = iy / 128;

    ipf = (HP_x2pix[ix_hi] + HP_y2pix[iy_hi]) * (128 * 128) + (HP_x2pix[ix_low] + HP_y2pix[iy_low]);
    ipf = int((ipf / pow(ns_max / 8192, 2)));     /* in {0, nside**2 - 1} */
    return uint(ipf + face_num*pow(8192, 2)); /* in {0, 12*nside**2 - 1} */ 
}

uint HP_ignoreEverySecondBit(uint value, bool odd, uint countBits)
{
    uint result = 0u;
    uint mask = odd ? 1u : 2u;
    countBits = countBits / 2;
    for (uint i = 0; i < countBits; ++i)
    {       
        result += clamp((value & mask), 0u, 1u) * uint(pow(2u, i));     
        value = value >> 2;
    }
    return result;
}

highp float HP_getU(uint healPixNumber, uint level)
{
    return float(HP_ignoreEverySecondBit(healPixNumber, true, HP_levelSet[level])) / HP_levelSetPow[level];
}

highp float HP_getV(uint healPixNumber, uint level)
{
    return float(HP_ignoreEverySecondBit(healPixNumber, false, HP_levelSet[level])) / HP_levelSetPow[level];
}

struct healPixData
{
    uint healPixNumber;
    uint level;
    highp float u;
    highp float v;
};

layout(binding = 0) uniform UniformBufferObject 
{
    uint level;
    uint scale;
} ubo;

layout(std430, binding = 2) buffer SSBO2
{
    healPixData valuesSSBO2[];
};

void main()
{
    uint width = 16 * ubo.scale;
    uint height = 8 * ubo.scale;
    uint arrayIndex = gl_GlobalInvocationID.x * height + height - gl_GlobalInvocationID.y - 1;
    
    vec2 currentIndex = vec2(gl_GlobalInvocationID.xy);
    vec2 latLon = currentIndex / vec2(width, height) * vec2(2 * 3.1415926, 3.1415926);
    
    healPixData data;
    data.healPixNumber = HP_latLonToHealPix(latLon.y, latLon.x); 
    data.u = HP_getU(data.healPixNumber, ubo.level);
    data.v = HP_getV(data.healPixNumber, ubo.level);
    
    valuesSSBO2[arrayIndex] = data;     
}

```

There are no race conditions that can cause artifacts. The artifacts do not change, there is no randomness. Any idea where and why this happened?

**Tags:** glsl, vulkan, artifacts

> **Comments:**
> **username_1**: I don't know the cause, but can suggest some ideas for debugging.  You know the (x,y) of an artifact pixel.  Add shader [debug printf](https://github.com/KhronosGroup/Vulkan-ValidationLayers/blob/ac2e01fc1a67329da0c879ecd2296b76643c49c0/docs/debug_printf.md) statements that print out some of the intermediate values only for that (x,y) and compare them with good values from the C++ code.  In addition, I see some usages of `pow(x,2)` that could be changed to `(x * x)`.  It is possible that the shader compiler doesn't recognize the shortcut and uses logs to compute the function.  
> **username_2**: If that can help: it seems that the artifacts are aligned horizontally, with equal spacing on a given row. Spacing may differ from row to row. It is not impossible that the rows are complete with some artifacts so little contrasted that one does not see them.  

---

## 0 answers

