Search Unity

  1. Good news ✨ We have more Unite Now videos available for you to watch on-demand! Come check them out and ask our experts any questions!
    Dismiss Notice
  2. Ever participated in one our Game Jams? Want pointers on your project? Our Evangelists will be available on Friday to give feedback. Come share your games with us!
    Dismiss Notice

Infinite parallax hair volume using wrapping grid tracing (attempt)

Discussion in 'Shaders' started by neoshaman, Dec 5, 2018.

  1. neoshaman


    Feb 11, 2011
    There is this thing I need to get done, it's an extension and end goal of experiment with relief function shader. Basically I want to trace in fragment shader hair that had parallax property. Since I'm having problem with simple math due to a past burn out, my brain blanking out on thing I should know, I figure out that chronicling the attempt publicly and verbalizing it to other people might help me get thing fall into place, and maybe some people might be interested and drop some clue.

    Hair rendering is notoriously hard and costly, due to the high number of occluded elements, which may or not contribute to the final image. I propose a technique that use a simplified model of hair, as parallel line in an extruded grid, and a mathematical representation that would allow to find the first visible intersection to a line without having to pay the cost for the occluded many, and potentially infer how many hair has been traversed with the same operation at the exit of the volume.

    The main idea is to realize that the simplified hair is "coherent", ie all are in the same position in the grid (idea to shuffle it will be explored later) and therefore all cells are indistinguishable from another, they can be superposed. This mean that a ray entering at one cell boundary will have a predictable trajectory wrapping in the superposed space. The challenge will be to figure out the math of intersection in that wrapping space, which would be the equivalent of determining occlusion in one hit. To simplify further, we will move the hair position to the origin corner of the cell so it's 0,0. We will adapt the DDA algorithm to get to that result.
    Further analysis show we can simplified the problem more by simply projecting to 2d position as it's extruded, therefore we can align the grid with one axis, then project to 1d, where it become a problem similar to the " how many timed clock hand overlap in a day" or planetary alignment puzzle

    The end goal is to have a version of the shader that can also represent curly hair (intersection with helix line).

    Given I found all the mathematical foundation it should have been trivial from there ... I'll add image later to illustrate the technique.
    AlanMattano likes this.
  2. neoshaman


    Feb 11, 2011
    Oups, end of years + a loss in the family + a new occupation didn't make me live up to those statement. Slowly coming back to it...

    a. the basic
    I was trying to model afro textured hair cut, but wasn't satisfied with the result. Afro textured is inherently volumetric and I wanted to reflect that fact.
    There is multiple way to implement that idea but my first instinct was to try a parallax shader, to model the superposition of strand, basically representing multiple layer displaced virtually by the view vector.

    Going through the math I realized there was opportunity to expand on that. But let's review the math of layer parallax (I hope there is no mistake, I'm not the most competent mathematician here).

    Let's consider an UV surface, with layer at regular interval d, and a view normal vector v. Let's simplify the representation of a 2d cross section.

    hair parallax.png
    Due to thales theorem, we can see than the offset is a ratio k of distance to the plane. If we set d=1 then k<1 because v is a unit vector as such k is equivalent to the depth when v is looking from the top.

    Once we see the result we understand that the sample are regularly space in a predictable way, we might be able to derive mathematically an intersection with a fixed position, and skip sampling multiple time an occlusion assuming it happen at the same regular interval.

    If we assume hair strands are regularly space line in uv space, we can approximate the hair volume occlusion mathematically. While strands have thickness, we will assume they are mathematical line and try to find the basis of the math from there.

    The insight is that multiple plane at the same distance is equivalent to a wrapping plane, also UV wrap around the 0-1 range. SO if we can port the math to wrapping space ... it would be equivalent to tracing through an infinite regular grid. Also since we have wrapping, we can realize it's equivalent to looping, stepping at regular interval is therefore equivalent to the clock needle overlap problem.

    The other trick is to realize that the support (the triangle face) can bias the visual such has the straight line become curved while the math stay the same (operating on UV space not in visual space). Of course it has limitation (such has curving beyond the horizon wouldn't work, as the ray would need to intersect twice the same strand.

    It's not really a new idea, some people have experimented with modulo space, though I need to find one that give interesting result for hair strands occlusion.

    I'll explore that next time.
  3. neoshaman


    Feb 11, 2011
    Wrapping the space

    First let’s assume the hair strands are all situated at position 1 in x in all layers, and the view rays start at position 0 in x on the first. We will trace rays through all layers such as they intercept the strand position at each layer.

    If we trace back to the first layer, every intersections with a layer before hitting a strand position, we can observe a pattern. Because the ray slope is constant, and each layer distance is the same, the projected intersection to the first layer are equidistant. We can see there is a correlation between the depth of the layer, where the ray hit the strand is, and the number of projected intersection. Therefore the intersections number give us the depth. We then know we can just divide the distance, between the ray x origin and the strand x position, by the number of intersection to get the depth.

    The corollary is that we can guess which depth a strand is hit, by looking at the number of step it takes, which is a function of the slope. Figuring the slope at the first layer give you the delta of the step, you then just to divide the distance.
    test hair drawing.png
    Wrapping horizontally is similar, the difference is that the number of step happen vertically. But really we just need to observe, that the distance to the targeted strand horizontally, is a multiple of the distance to the reference strand on the first later.

    Those simplified observation are here to demonstrate how we can infer positions of virtual strands (ie strand beyond the first layer and horizontal wrapping), based on the rays slopes vs the reference strand position, in a kind of wrapping trigonometry. And also how we can basicaly simplify the math to a linear convergence using division or modulo.

    However to generalize further, we still need to demonstrate how it works when we start changing the slopes and the offset of the ray, instead of assuming all ray point to a strand from the zero position. And we need to verify when the math allow convergence to a virtual position of a strands and when it does not. Doing so would allow us to guess the first hit in a virtual volume.

    Once we have done that, we have effectively found a way to trace an infinite volume of primitive, however hair has thickness, we will need to find a way to intercept more than an infinitely small position. Ideally we would also generalized to different kind of wrapping shape, such has helix, which would allow us to emulate curly hair.
  4. neoshaman


    Feb 11, 2011
    Up until we cheated by locking the ray slopes and position to the a single place into the receiving surface and only to slopes guaranteed to hit a target point. It was done to demonstrate how it would translate to a 1d problem of interval overlapping. But now we need to introduce two new parameter:
    1. arbitrary slopes
    2. arbitrary starting position of the ray.

    1. variation of slopes
    As we have shown earlier, the slope of the ray translate into discrete steps, with same size, on the "surface", and there is collision when they are basically a multiple of the wrapping size, and that multiple give you the depth or distance of the collision within the infinite space. But what happen when we have an arbitrary slope, that give us an arbitrary step, and we want to find the multiple to the wrapping size, does that multiple even exist?

    Basically we can reformulate the problem as beat synchronicity, let say you have two beat, one that goes on twos and the other on three, when will the beat pulse together? In this simple set up it's easy, they will be synchrone every 6 beats, that's the "collision depth". But what does that mean in mathematical form, how can we generalize it? seems simple though isn't it?

    2. variation of position
    Up until we cheated (again), by fixing the position at the start we made it easier, since the space is wrapping the starting and the target are basically the same position, and since we normalized the wrapping distance, it also mean we had neat number to play with. But what happen when we move the ray to anywhere within the wrap space, as it should happen in a real case? Well let's see the beat analogy again, let's say we have two beat on twos, but one start half way after the first beat, they will never be in sync again, ie similar beat with an offset, they will never pulse, can beat we translate that into mathematical form, to check that two beat will never converge to the same time?

    I use beat analogy, but we can basically use other one, since the space is wrapping, we can also see the surface "segment" as a circle, since we have stuff moving at different rate on that circle is analogous to a clock needles, which mean we can relate to the clock needle superposition problem (itself analogous to the planet alignement problem), with the added complication that we are moving by "discrete steps" rather than a "continuous rate".

    So basically it's this, given two elements moving at discrete step, can they overlap and when? The minimal distance will gave us the depth of the collision, and probably we can define the periodicity of the collision after that.

    So here is where my competence end, I need to translate those insight into mathematical form to resolve.
  5. neoshaman


    Feb 11, 2011
    Up until now I showed that the problem could express as a 1d line with interval overlapping. But since the space is wrapping on itself, that is the first line and teh last line is basically the same, instead of an infinite line with repeated "beat", we can represent it as a circle with interval stop around it:
    In this example the black tick overshoot the target after one wrap, the red tick is the resulting position of overshooting. We need to calculate how many wrap before it can actually overlap the target.

    Using this representation allow us to find equivalence, it looks like a clock, it also look like circular planet orbital. Turns out we can probably get inspire by that to find mathematical equivalence. If you think a bit about it, it turns out that's a similar problem than planets alignment, or the clock needle superposition riddle. The hair target

    Let's use planetary alignment.

    We find this googling on the net:

    which gave us this:
    If it's correct we have this:
    let say we have two point moving at two rate,
    - the position of each point = rate * t.
    - Therefore position1 = position2 it mean that rate1 * t = rate2 * t,
    - that is (rate1-rate2) * t .
    What's missing is the bigger period in which rate1 and rate2 are inscribed,
    - we can get it by having rate1 * rate2.
    - So rate1 and rate2 are sync for t= k(rate1*rate2) / (rate1 -rate2) with k being the number of period.

    However we have a specific case where rate2 = the normalized wraping size, rate2 = 1, for k = 1 we have t = rate1 / (rate1-1).

    But we only assume no rate as an offset,
    - so if we check with offset we have position1+ offset = position2,
    - which translate (rate1 *t) +offset = rate2 *t
    - which lead to t * (rate1 - rate2) + offset
    and we can probably write the final result as t = k(rate1*rate2) / (rate1-rate2+offset)
    which finalized as
    t = k(rate1) / (rate1+offset-1) which has no solution when rate1+offset = 1

    Probably? Any mistake?

    Extruding 2d intersections points
    More math (help given by someone else in another forum), this would help us define the line in 3d, since hair is axis align with the support geometry UV space (the deformation/orientation of that support geometry being the approximation of the hair direction) once we have the virtual depth, we can compute the virtual position as a function of the origin and the slope.

    Jittering the depth alignment
    Anyway, we can have some more insight, we have assumed the wrapping itself was fixed, if the wrapping windows itself move, it's equivalent to find 3 planets alignment ... What's the use? well if we stay in the two beat system, we get all hair perfectly aligned in regular interval, if the wrapping move too, we can jitter the depth predictably by multiplying by clever number.

    However this isn't that easy, if we go back to the 2d representation, let say we have an angle where the beat is bigger than the wrapping size, we would wrap horizontally first, so the jittering position will not match a ray sampling the nearby cells position, because the order of traversal would matter, so we need to define the wrapping displacement as a 2d problem. Which I haven't figure out yet.

    More works is needed
    Up until then we have assume mathematical object, ie point that have no size. In order to properly compute hair, we need thickness, we need to find a way to detect when our beat is within a distance of the target to get a hit, not just right on the target. Once we can do that, we can probably extend to a 3d volume cross section (the circle shape of the hair), the problem still being finding a single point, so the extrusion should remain the same.

    We also need to find a way to intersect helicoidal shape, they are too repeating, so a wrapping function is fine, the challenge is that there isn't cheap formula to intersect sin wave, which is needed to get the helix.

    The end game of this exercice, infinite field of infinite helix shape!

    If we were to just do infinite straight hair, probably a beat base solution is overkill, on shadertoy, there is already wrapping tracing shader, we could just extend them to get hair approximation.
    Last edited: Jul 12, 2019
    Invertex likes this.
  6. neoshaman


    Feb 11, 2011
  7. neoshaman


    Feb 11, 2011
    Doing recursive interception on wrapping helix is hard, the equation of the helix is x = r cos t, y =r sin t, z=a t ... all those sin and cos create problem, there is no simple equation to intercept them. And it's a curve that is not convex, it's hard to find simplification.

    I try simplifying by boxing the shape abusing the inherent symmetry, and then approximate it with a x² curves that we can reliably intercept but I wasn't able to see anything that will put me closer to the goal.

    Up until I decided to step back and instead of trying to intercept the coil, intercept the support cylinder, it can be reduce to a 2D problem of line circle intersection, which doesn't have any sin or cos, but then that lead to a realization that the cylinder unroll into a 2D wrapping space where the coil is basically just the diagonal line ...

    Now the idea is to try and translate the cylinder intersection equation into an equation of all potential wrapping hit points, given that all points are in the form of (x+n, y+m, z+o), with n,m,o being the wrapping period. Convert that into cylinder space, hope that no sincos appear and trying to intercept the diagonal for hit on the helix...
    Last edited: Oct 1, 2019
  8. neoshaman


    Feb 11, 2011
    I spend some time with math equation ... that lead nowhere :oops: Trying to express thing into varying variable, couldn't do it, maybe I wasn't using the proper visualization:


    Oh ...
    You can turn the whole 2D cylinder projection into a single dimension where the line move at integer wrapping step ... where I have seen that before?

    That also open the possibility of doing 3D sampling of pseudo noise integral, ie doing cloud without raymarching a noise?

    Anyway, now I can cleanly project the point back to the cylinder with a cos and have a definitive way to measure hit and miss of rays. Also it solve the vertical on the unwrapping of the cylinder, but now we need to find the delta ie where on that vertical the infinite set of point fall and how to find that intersection with the helix ... We have the horizontal delta but we need the vertical delta, if possible as a predictable pattern. My guess is that it will be nasty SIN, my hope is that it is quadric instead ...

    Anybody has any idea?
    Last edited: Oct 1, 2019
  9. neoshaman


    Feb 11, 2011
  10. neoshaman


    Feb 11, 2011
    Code (CSharp):
    1. Hairtech flat only
    2. [code]
    3. #define DARKEN_OVER_DISTANCE 1  // This makes it easier to see the different layers in a static image, not to hide a max distance.
    4. #define SHOW_2D_SHAPE        1  // This makes the 2d shape be shown in upper left corner
    5. const float c_camearaDistance    = 6.0;
    6. const float c_cameraViewWidth    = 24.0;
    7. //============================================================
    8. //============================================================
    9. bool BooleanFunction_Square (vec2 current)
    10. {
    11.     return (current.x < 0.25 || current.x > 0.75 || current.x < 0.25 || current.x > 0.75);
    12. }
    13. //============================================================
    14. float NumberStepsFunction_Square (vec2 current, vec2 stepValue)
    15. {
    16.     // if it's already in the shape, no steps to take
    17.     if (BooleanFunction_Square(current))
    18.         return 0.0;
    20.     float stepsX;
    21.     if (stepValue.x < 0.0)
    22.         stepsX = (current.x - 0.25) / -stepValue.x;
    23.     else
    24.         stepsX = (0.75 - current.x) / stepValue.x;
    26.     return ceil(stepsX);
    27. }
    28. //============================================================
    29. void mainImage( out vec4 fragColor, in vec2 fragCoord )
    30. {
    31.     // set up the camera
    32.     vec3 cameraPos;
    33.     vec3 rayDir;
    34.     {
    35.         vec2 percent = (fragCoord / iResolution.xy) - vec2(0.5,0.5);
    36.         vec3 offset = vec3(0.5, 0.5, iTime+0.01);
    37.         float angleX = 0.0;
    38.         float angleY = 0.0;
    39.         if (iMouse.z > 0.0) {
    40.             vec2 mouse = iMouse.xy / iResolution.xy;
    41.             angleX = 3.14 + 6.28 * mouse.x;
    42.             angleY = (mouse.y - 0.5) * 3.14;//(mouse.y * 3.90) - 0.4;
    43.         }
    44.         vec3 cameraFwd    = (vec3(sin(angleX)*cos(angleY), sin(angleY), cos(angleX)*cos(angleY)));      
    45.         vec3 cameraRight = normalize(cross(vec3(0.0,1.0,0.0),cameraFwd));
    46.         vec3 cameraUp = normalize(cross(cameraFwd, cameraRight));
    47.         cameraPos = vec3(0.0, 0.0, -1.0) + offset;
    48.         vec3 cameraTarget = vec3(0.0, 0.0, 0.0) + offset;
    49.         float cameraViewHeight    = c_cameraViewWidth * iResolution.y / iResolution.x;
    50.         vec3 rayTarget = cameraPos +  cameraFwd * c_cameraDistance + cameraRight * c_cameraViewWidth * percent.x + cameraUp * cameraViewHeight * percent.y;
    51.         rayDir = normalize(rayTarget - cameraPos);
    52.     }
    53.     // keep the camera in a unit cube
    54.     cameraPos = fract(cameraPos);
    56.     // If ray facing negative on z axis, just flip direction and invert where we are in the cube on the z axis.
    57.     // Now we only have to deal with positive z directions.
    58.     if (rayDir.z < 0.0) {
    59.         rayDir *= -1.0;
    60.         cameraPos.z = 1.0 - cameraPos.z;
    61.     }
    63.     // calculate the 3d position of the next two plane intersections
    64.     float intersection1Distance = (1.0 - cameraPos.z) / rayDir.z;
    65.     float intersection2Distance = (2.0 - cameraPos.z) / rayDir.z;
    66.     vec3 intersection1 = fract(cameraPos + rayDir * intersection1Distance);
    67.     vec3 intersection2 = fract(cameraPos + rayDir * intersection2Distance);
    69.     // Calculate how much the uv changes from intersection1 to intersection2.
    70.     // Convert it from [0,1] to [-0.5, 0.5].
    71.     // We need to know this to know if the uvs are going positive or negative and by how much, on each axis.
    72.     vec2 uvStep = intersection2.xy - intersection1.xy;
    73.     uvStep = fract(uvStep + 0.5) - 0.5;
    75.     // calculate how many steps it takes to hit something on the X and Y axis and take whichever hits first.
    76.     float steps = 0.0;
    77.     steps = NumberStepsFunction_Square(intersection1.xy, uvStep);
    78.     // calculate how far it is to the intersection we found
    79.     float dist = (1.0 - cameraPos.z) / rayDir.z + steps / rayDir.z;
    81.     #if DARKEN_OVER_DISTANCE
    82.     float tint = clamp(1.0 - dist / 5.0, 0.0, 1.0);
    83.     #else
    84.     float tint = 1.0;
    85.     #endif
    87.     // calculate the hit point
    88.     vec3 hitPoint = cameraPos + rayDir * dist;
    89.     vec2 uv = hitPoint.xy;
    91.     // sample the texture
    92.     fragColor = vec4(texture(iChannel0, uv).rgb * tint, 1.0);
    93. }


    Code (CSharp):
    1. //============================================================
    2. float binarySign (float v)
    3. {
    4.     return step(0.0, v) * 2.0 - 1.0;
    5. }
    7. //============================================================
    8. // returns t
    9. // circle xy = position, z = radius
    10. // Adapted from "real time collision detection" IntersectRaySphere()
    11. float RayIntersectCircle (in vec2 rayPos, in vec2 rayDir, in vec3 circle)
    12. {
    13.     // rayDir isn't normalized, so normalize it but remember it's length
    14.     float rayLen = length(rayDir);
    15.     rayDir = normalize(rayDir);
    17.     vec2 m = rayPos - circle.xy;
    18.     float b = dot(m, rayDir);
    19.     float c = dot(m, m) - circle.z*circle.z;
    21.     // Exit if the ray is outside the circle and pointing away from the circle
    22.     if (c > 0.0 && b > 0.0)
    23.         return -1.0;
    25.     float discr = b*b - c;
    27.     // A negative discriminant means it missed the sphere
    28.     if (discr < 0.0)
    29.         return -1.0;
    31.     float t = -b - sqrt(discr);
    32.     if (t < 0.0)
    33.         t = -b + sqrt(discr);
    35.     return t / rayLen;
    36. }
    Code (CSharp):
    2. //============================================================
    3. float NumberStepsFunction_L_OneAxis (float current, float stepValue)
    4. {
    5.     float steps;
    6.     if (stepValue < 0.0)
    7.         steps = (current - 0.5) / -stepValue;
    8.     else
    9.         steps = (1.0 - current) / stepValue;
    10.     return steps;
    11. }
    13. //============================================================
    14. bool BooleanFunction_L (vec2 current)
    15. {
    16.     return (current.x < 0.5 || current.y < 0.5);
    17. }
    19. //============================================================
    20. float NumberStepsFunction_L (vec2 current, vec2 stepValue)
    21. {
    22.     if (BooleanFunction_L(current))
    23.         return 0.0;
    25.     float stepsX = NumberStepsFunction_L_OneAxis(current.x, stepValue.x);
    26.     float stepsY = NumberStepsFunction_L_OneAxis(current.y, stepValue.y);
    27.     return ceil(min(stepsX, stepsY));
    28. }
    Code (CSharp):
    2. //============================================================
    3. bool BooleanFunction_Checker (vec2 current)
    4. {
    5.     if (current.x >= 0.5)
    6.     {
    7.         return current.y < 0.5;
    8.     }
    9.     else
    10.     {
    11.         return current.y >= 0.5;
    12.     }
    13. }
    15. //============================================================
    16. float NumberStepsFunction_Checker (vec2 current, vec2 stepValue)
    17. {
    18.     // if it's already in the shape, no steps to take
    19.     if (BooleanFunction_Checker(current))
    20.         return 0.0;
    22.     // there are different values to reach based on where we are in the pattern.
    23.     float lower = step(0.5, current.x) * 0.5;
    24.     float upper = lower + 0.5;
    26.     // see how long to escape the box on each axis, take sooner event
    27.     float stepsX;
    28.     if (stepValue.x < 0.0)
    29.         stepsX = (current.x - lower) / -stepValue.x;
    30.     else
    31.         stepsX = (upper - current.x) / stepValue.x;
    33.     float stepsY;
    34.     if (stepValue.y < 0.0)
    35.         stepsY = (current.y - lower) / -stepValue.y;
    36.     else
    37.         stepsY = (upper - current.y) / stepValue.y;    
    39.     return ceil(min(stepsX, stepsY));  
    40. }
    Code (CSharp):
    2. //============================================================
    3. bool BooleanFunction_Square (vec2 current)
    4. {
    5.     return (current.x < 0.25 || current.x > 0.75 || current.y < 0.25 || current.y > 0.75);
    6. }
    8. //============================================================
    9. float NumberStepsFunction_Square (vec2 current, vec2 stepValue)
    10. {
    11.     // if it's already in the shape, no steps to take
    12.     if (BooleanFunction_Square(current))
    13.         return 0.0;
    15.     float stepsX;
    16.     if (stepValue.x < 0.0)
    17.         stepsX = (current.x - 0.25) / -stepValue.x;
    18.     else
    19.         stepsX = (0.75 - current.x) / stepValue.x;
    21.     float stepsY;
    22.     if (stepValue.y < 0.0)
    23.         stepsY = (current.y - 0.25) / -stepValue.y;
    24.     else
    25.         stepsY = (0.75 - current.y) / stepValue.y;  
    27.     return ceil(min(stepsX, stepsY));
    28. }
    Code (CSharp):
    2. //============================================================
    3. bool BooleanFunction_Circle (vec2 current)
    4. {
    5.     const float radius = 0.25;
    6.     const float radiusSq = radius * radius;
    7.     vec2 rel = current - vec2(0.5, 0.5);
    8.     return rel.x*rel.x + rel.y*rel.y > radiusSq;
    9. }
    11. //============================================================
    12. float NumberStepsFunction_Circle (vec2 current, vec2 stepValue)
    13. {
    14.     // if it's already in the shape, no steps to take
    15.     if (BooleanFunction_Circle(current))
    16.         return 0.0;
    18.     float steps = RayIntersectCircle(current, stepValue, vec3(0.5, 0.5, 0.25));
    19.     return ceil(steps);
    20. }
    Code (CSharp):
    2. //============================================================
    3. bool BooleanFunction_ThinSquare (vec2 current)
    4. {
    5.     return (current.x < 0.05 || current.x > 0.95 || current.y < 0.05 || current.y > 0.95);
    6. }
    8. //============================================================
    9. float NumberStepsFunction_ThinSquare (vec2 current, vec2 stepValue)
    10. {
    11.     // if it's already in the shape, no steps to take
    12.     if (BooleanFunction_ThinSquare(current))
    13.         return 0.0;
    15.     float stepsX;
    16.     if (stepValue.x < 0.0)
    17.         stepsX = (current.x - 0.05) / -stepValue.x;
    18.     else
    19.         stepsX = (0.95 - current.x) / stepValue.x;
    21.     float stepsY;
    22.     if (stepValue.y < 0.0)
    23.         stepsY = (current.y - 0.05) / -stepValue.y;
    24.     else
    25.         stepsY = (0.95 - current.y) / stepValue.y;  
    27.     return ceil(min(stepsX, stepsY));
    28. }
    Last edited: May 22, 2020 at 5:01 AM
    Invertex likes this.