Converting recursive function to iterative [Solved]

TheSnidr

Heavy metal viking dentist
GMC Elder
I've created a recursive script for casting a ray onto an octree and finding the nearest intersection. It performs recursion with the following format (pseudocode):
Code:
/// @description recursion(stuff)
/// @param stuff
if conditionA
{
    //End recursion
    return thing;
}
if conditionB
{
    something = recursion(stuff);
}
if conditionC
{
    something = recursion(stuff);
}
if conditionD
{
    something = recursion(stuff);
}
if conditionE
{
    something = recursion(stuff);
}
return something;
Here's the actual script in case you're interested:
Code:
/// @description smf_collision_cast_ray_recursive(collisionBuffer, lineStart, lineEnd, regionSize, regionPos)
/// @param collisionBuffer
/// @param lineStart
/// @param lineEnd
/// @param regionSize
/// @param regionPos
var i, j, k, r, t, n, a, b, colBuffer, rayStart, rayEnd, regionSize, halfRegionSize, regionPos, checkPos, intersect, newPos, region, tris, tri, ret;
colBuffer = argument0;
rayStart = argument1;
rayEnd = argument2;
regionSize = argument3;
halfRegionSize = regionSize / 2;
regionPos = argument4;

r = buffer_read(colBuffer, buffer_s32);

//Exit condition
if (!r)
{
    buffer_seek(colBuffer, buffer_seek_start, -r);
    n = buffer_read(colBuffer, buffer_u16);
    for (i = 0; i < n; i ++){tris[i] = buffer_read(colBuffer, buffer_u16);}
    for (i = 0; i < n; i ++)
    {
        //Find intersection with triangle plane
        tri = smf_get_triangle(colBuffer, tris[i]);
        t = dot_product_3d(tri[9], tri[10], tri[11], rayEnd[0] - rayStart[0], rayEnd[1] - rayStart[1], rayEnd[2] - rayStart[2]);
        if (t == 0){continue;}
        t = dot_product_3d(tri[9], tri[10], tri[11], tri[0] - rayStart[0], tri[1] - rayStart[1], tri[2] - rayStart[2]) / t;
        if ((t <= 0) or (t >= 1)){continue;}
        ret[0] = lerp(rayStart[0], rayEnd[0], t);
        ret[1] = lerp(rayStart[1], rayEnd[1], t);
        ret[2] = lerp(rayStart[2], rayEnd[2], t);
  
        //Check if the intersection is inside the triangle. If not, discard and continue.
        a[0] = ret[0] - tri[0]; a[1] = ret[1] - tri[1]; a[2] = ret[2] - tri[2];
        b[0] = tri[3] - tri[0]; b[1] = tri[4] - tri[1]; b[2] = tri[5] - tri[2];
        if (dot_product_3d(tri[9], tri[10], tri[11], a[2] * b[1] - a[1] * b[2], a[0] * b[2] - a[2] * b[0], a[1] * b[0] - a[0] * b[1]) < 0){continue;}
        a[0] = ret[0] - tri[3]; a[1] = ret[1] - tri[4]; a[2] = ret[2] - tri[5];
        b[0] = tri[6] - tri[3]; b[1] = tri[7] - tri[4]; b[2] = tri[8] - tri[5];
        if (dot_product_3d(tri[9], tri[10], tri[11], a[2] * b[1] - a[1] * b[2], a[0] * b[2] - a[2] * b[0], a[1] * b[0] - a[0] * b[1]) < 0){continue;}
        a[0] = ret[0] - tri[6]; a[1] = ret[1] - tri[7]; a[2] = ret[2] - tri[8];
        b[0] = tri[0] - tri[6]; b[1] = tri[1] - tri[7]; b[2] = tri[2] - tri[8];
        if (dot_product_3d(tri[9], tri[10], tri[11], a[2] * b[1] - a[1] * b[2], a[0] * b[2] - a[2] * b[0], a[1] * b[0] - a[0] * b[1]) < 0){continue;}
  
        //The line intersects the triangle. The ray is altered by a tiny amount to make sure the returned value is on the correct side of the triangle
        t += 1 / dot_product_3d(tri[9], tri[10], tri[11], ret[0] - rayStart[0], ret[1] - rayStart[1], ret[2] - rayStart[2]);
        rayEnd[0] = lerp(rayStart[0], rayEnd[0], t);
        rayEnd[1] = lerp(rayStart[1], rayEnd[1], t);
        rayEnd[2] = lerp(rayStart[2], rayEnd[2], t);
    }
    return rayEnd;
}

//Find starting region of the line segment
region[0] = (rayStart[0] > regionPos[0] + halfRegionSize);
region[1] = (rayStart[1] > regionPos[1] + halfRegionSize);
region[2] = (rayStart[2] > regionPos[2] + halfRegionSize);
for (i = 0; i < 3; i ++)
{
    if (rayEnd[i] == rayStart[i]){continue;}
    if ((rayStart[i] < regionPos[i]) or (rayStart[i] > regionPos[i] + regionSize))
    {
        checkPos = regionPos;
        checkPos[i] += region[i] * regionSize;
        t = (checkPos[i] - rayStart[i]) / (rayEnd[i] - rayStart[i]);
        if ((t <= 0) or (t >= 1)){continue;}
        j = (i + 1) mod 3;
        k = (i + 2) mod 3;
        ret[j] = lerp(rayStart[j], rayEnd[j], t);
        ret[k] = lerp(rayStart[k], rayEnd[k], t);
        if  ((ret[j] < regionPos[j]) or (ret[j] > regionPos[j] + regionSize) or (ret[k] < regionPos[k]) or (ret[k] > regionPos[k] + regionSize)){continue;}
        region[j] = (ret[j] > regionPos[j] + halfRegionSize);
        region[k] = (ret[k] > regionPos[k] + halfRegionSize);
        break;
    }
}

//Check for collisions in the starting region of the line segment
buffer_seek(colBuffer, buffer_seek_start, r + region[0] * 4 + region[1] * 8 + region[2] * 16);
newPos[0] = regionPos[0] + halfRegionSize * region[0];
newPos[1] = regionPos[1] + halfRegionSize * region[1];
newPos[2] = regionPos[2] + halfRegionSize * region[2];
rayEnd = smf_collision_cast_ray_recursive(colBuffer, rayStart, rayEnd, halfRegionSize, newPos);

//Check for intersections along the middle axis in all three dimensions
for (i = 0; i < 3; i ++)
{
    if (rayEnd[i] == rayStart[i]){continue;}
    t = (regionPos[i] + halfRegionSize - rayStart[i]) / (rayEnd[i] - rayStart[i]);
    if ((t <= 0) or (t >= 1)){continue;}
    j = (i + 1) mod 3;
    k = (i + 2) mod 3;
    intersect[j] = lerp(rayStart[j], rayEnd[j], t) - regionPos[j];
    intersect[k] = lerp(rayStart[k], rayEnd[k], t) - regionPos[k];
    if ((intersect[j] < 0) or (intersect[j] > regionSize) or (intersect[k] < 0) or (intersect[k] > regionSize)){continue;}
    region[i] = (rayStart[i] < regionPos[i] + halfRegionSize);
    region[j] = (intersect[j] >= halfRegionSize);
    region[k] = (intersect[k] >= halfRegionSize);
    newPos[0] = regionPos[0] + halfRegionSize * region[0];
    newPos[1] = regionPos[1] + halfRegionSize * region[1];
    newPos[2] = regionPos[2] + halfRegionSize * region[2];
    buffer_seek(colBuffer, buffer_seek_start, r + region[0] * 4 + region[1] * 8 + region[2] * 16);
    rayEnd = smf_collision_cast_ray_recursive(colBuffer, rayStart, rayEnd, halfRegionSize, newPos);
}
return rayEnd;
I just can't see any ways of converting this to an iterative script. I've read multiple places that "any recursive function can be rewritten to an iterative one", but I'm kinda lost on this one.
And even if it's possible, is it feasible? The maximal recursive depth in this function is 10, it will never go deeper and will rarely even reach 10.

What I'm currently trying to get working is to store my progress in ds_stacks. Once I reach the exit condition, I pop the stack and contionue iterating. If I get it to work I'll post the solution, but I seem to be far away still...

EDIT:
OMG I managed to turn it into an iterative function! :D
The stack method I mentioned works like a charm. Instead of recursively performing the function, I made a "while true"-loop and go back and forth in this with the help of a ds_stack. It's ridiculously much faster than the recursive function, and can be performed with a single script call.
In essence, this is what it does:
Code:
stack = ds_stack_create();
progress = 0;
conditionA = false;
stuff = something;
while true
{
    if conditionA
    {
        //If we're at the lowest level, perform collision checks
        //This is the same as the exit condition for the recursive function
        perform_collision_check(stuff);
    
        //Go back up to the previous level
        progress = ds_stack_pop(stack);
        stuff = ds_stack_pop(stack);
    }
    if progress == 0 and functionA(stuff)
    {
        ds_stack_push(stack, progress + 1);
        ds_stack_push(stack, stuff);
        if functionB(stuff){conditionA = true;} //If certain conditions are met, perform the recursive exit condition in the next loop
        progress = 0;
        continue; //Continue makes the loop start over again, effectively jumping back to the start - but with the new variables
    }
    if progress == 1 and functionC(stuff)
    {
        ds_stack_push(stack, progress + 1);
        ds_stack_push(stack, stuff);
        if functionD(stuff){conditionA = true;} //If certain conditions are met, perform the recursive exit condition in the next loop
        progress = 0;
        continue; //Continue makes the loop start over again, effectively jumping back to the start - but with the new variables
    }
    if progress == 2 and functionE(stuff)
    {
        ds_stack_push(stack, progress + 1);
        ds_stack_push(stack, stuff);
        if functionF(stuff){conditionA = true;} //If certain conditions are met, perform the recursive exit condition in the next loop
        progress = 0;
        continue; //Continue makes the loop start over again, effectively jumping back to the start - but with the new variables
    }
    if progress == 3 and functionG(stuff)
    {
        ds_stack_push(stack, progress + 1);
        ds_stack_push(stack, stuff);
        if functionH(stuff){conditionA = true;} //If certain conditions are met, perform the recursive exit condition in the next loop
        progress = 0;
        continue; //Continue makes the loop start over again, effectively jumping back to the start - but with the new variables
    }
 
    //If the stack is empty, we're done iterating
    if ds_stack_size(stack) == 0{break;}
 
    //If we're done in this level, go back up one level
    progress = ds_stack_pop(stack);
    stuff = ds_stack_pop(stack);
}
ds_stack_destroy(stack);
And here's the iterative script in its entirety:
Code:
/// @description smf_collision_cast_ray(collisionBuffer, x1, y1, z1, x2, y2, z2)
/// @param collisionBuffer
/// @param x1
/// @param y1
/// @param z1
/// @param x2
/// @param y2
/// @param z2
/*
Casts a ray from one point to another and returns the position of the first collision with geometry

Script made by TheSnidr
www.TheSnidr.com
*/
var colBuffer = argument0, regionSize = 65535, progress = 0, stack = ds_stack_create();
var a, b, n, i, j, k, l, t, tri, tris, modelX, modelY, modelZ, modelSize, transformScale, lineStart, lineEnd, regionPos, bufferPos, halfRegionSize, region, intersect;
buffer_seek(colBuffer, buffer_seek_start, 0);
bufferPos = buffer_read(colBuffer, buffer_s32);
modelX = buffer_read(colBuffer, buffer_f32);
modelY = buffer_read(colBuffer, buffer_f32);
modelZ = buffer_read(colBuffer, buffer_f32);
modelSize = buffer_read(colBuffer, buffer_f32);
transformScale = 65535 / modelSize;
lineStart[0] = transformScale * (argument1 - modelX);
lineStart[1] = transformScale * (argument2 - modelY);
lineStart[2] = transformScale * (argument3 - modelZ);
lineEnd[0] = transformScale * (argument4 - modelX);
lineEnd[1] = transformScale * (argument5 - modelY);
lineEnd[2] = transformScale * (argument6 - modelZ);
regionPos[2] = 0;
returnNormal = -1;
returnNormal[2] = 1;
while true
{
    if (bufferPos > 0)
    {   //Iterate through the octree
        halfRegionSize = regionSize / 2;
        for (l = progress; l < 4; l ++)
        {
            if l == 0
            {   //Check either the starting region of the ray or the first region it intersects
                region[0] = (lineStart[0] > regionPos[0] + halfRegionSize);
                region[1] = (lineStart[1] > regionPos[1] + halfRegionSize);
                region[2] = (lineStart[2] > regionPos[2] + halfRegionSize);
                for (i = 0; i < 3; i ++)
                {
                    if (lineEnd[i] == lineStart[i]) continue;
                    if (lineStart[i] >= regionPos[i] and lineStart[i] <= regionPos[i] + regionSize) continue;
                    t = (regionPos[i] + region[i] * regionSize - lineStart[i]) / (lineEnd[i] - lineStart[i]);
                    if (t < 0 or t > 1) continue;
                    j = (i + 1) mod 3; k = (i + 2) mod 3;
                    intersect[j] = lerp(lineStart[j], lineEnd[j], t) - regionPos[j];
                    intersect[k] = lerp(lineStart[k], lineEnd[k], t) - regionPos[k];
                    if (intersect[j] < 0 or intersect[j] > regionSize or intersect[k] < 0 or intersect[k] > regionSize) continue;
                    region[j] = (intersect[j] > halfRegionSize);
                    region[k] = (intersect[k] > halfRegionSize);
                    break;
                }
            }
            else
            {   //Check for intersections with the middle plane of each dimension
                i = l - 1;
                if (lineEnd[i] == lineStart[i]) continue;
                t = (regionPos[i] + halfRegionSize - lineStart[i]) / (lineEnd[i] - lineStart[i]);
                if (t < 0 or t > 1) continue;
                j = (i + 1) mod 3; k = (i + 2) mod 3;
                intersect[j] = lerp(lineStart[j], lineEnd[j], t) - regionPos[j];
                intersect[k] = lerp(lineStart[k], lineEnd[k], t) - regionPos[k];
                if (intersect[j] < 0 or intersect[j] > regionSize or intersect[k] < 0 or intersect[k] > regionSize) continue;
                region[i] = lineStart[i] < regionPos[i] + halfRegionSize;
                region[j] = intersect[j] >= halfRegionSize;
                region[k] = intersect[k] >= halfRegionSize;
            }
       
            //Push this region to stack
            ds_stack_push(stack, bufferPos);
            ds_stack_push(stack, regionPos[0]);
            ds_stack_push(stack, regionPos[1]);
            ds_stack_push(stack, regionPos[2]);
            ds_stack_push(stack, l);
       
            //Go to intersected child region
            if (region[0]){regionPos[0] += halfRegionSize; bufferPos += 4;}
            if (region[1]){regionPos[1] += halfRegionSize; bufferPos += 8;}
            if (region[2]){regionPos[2] += halfRegionSize; bufferPos += 16;}
            buffer_seek(colBuffer, buffer_seek_start, bufferPos);
            bufferPos = buffer_read(colBuffer, buffer_s32);
            regionSize /= 2;
            progress = 0;
            break;
        }
        if (l < 4) continue; //If we ended the for-loop prematurely, we should also restart the while-loop
    }
    else
    {   //If this is a leaf region, check for intersections with the triangles in this leaf
        l = 0;
        buffer_seek(colBuffer, buffer_seek_start, -bufferPos);
        repeat buffer_read(colBuffer, buffer_u16) 
            tris[l++] = buffer_read(colBuffer, buffer_u16);
        repeat l
        {   //Find intersection with triangle plane
            tri = smf_get_triangle(colBuffer, tris[--l]);
            t = dot_product_3d(tri[9], tri[10], tri[11], lineEnd[0] - lineStart[0], lineEnd[1] - lineStart[1], lineEnd[2] - lineStart[2]);
            if (t == 0) continue;
            t = dot_product_3d(tri[9], tri[10], tri[11], tri[0] - lineStart[0], tri[1] - lineStart[1], tri[2] - lineStart[2]) / t;
            if (t <= 0 or t >= 1) continue;
            intersect[0] = lerp(lineStart[0], lineEnd[0], t);
            intersect[1] = lerp(lineStart[1], lineEnd[1], t);
            intersect[2] = lerp(lineStart[2], lineEnd[2], t);
           
            //Check if the intersection is inside the triangle. If not, discard and continue.
            a[0] = intersect[0] - tri[0]; a[1] = intersect[1] - tri[1]; a[2] = intersect[2] - tri[2];
            b[0] = tri[3] - tri[0]; b[1] = tri[4] - tri[1]; b[2] = tri[5] - tri[2];
            if (dot_product_3d(tri[9], tri[10], tri[11], a[2] * b[1] - a[1] * b[2], a[0] * b[2] - a[2] * b[0], a[1] * b[0] - a[0] * b[1]) < 0) continue;
            a[0] = intersect[0] - tri[3]; a[1] = intersect[1] - tri[4]; a[2] = intersect[2] - tri[5];
            b[0] = tri[6] - tri[3]; b[1] = tri[7] - tri[4]; b[2] = tri[8] - tri[5];
            if (dot_product_3d(tri[9], tri[10], tri[11], a[2] * b[1] - a[1] * b[2], a[0] * b[2] - a[2] * b[0], a[1] * b[0] - a[0] * b[1]) < 0) continue;
            a[0] = intersect[0] - tri[6]; a[1] = intersect[1] - tri[7]; a[2] = intersect[2] - tri[8];
            b[0] = tri[0] - tri[6]; b[1] = tri[1] - tri[7]; b[2] = tri[2] - tri[8];
            if (dot_product_3d(tri[9], tri[10], tri[11], a[2] * b[1] - a[1] * b[2], a[0] * b[2] - a[2] * b[0], a[1] * b[0] - a[0] * b[1]) < 0) continue;
       
            //The line intersects the triangle. Save the triangle normal and intersection.
            returnNormal[0] = tri[9];
            returnNormal[1] = tri[10];
            returnNormal[2] = tri[11];
            lineEnd = intersect;
            intersect = -1;
        }
    }
    if !ds_stack_size(stack) break; //If the stack is empty, break the loop
   
    //Pop the previous region from stack
    progress = ds_stack_pop(stack) + 1;
    regionPos[2] = ds_stack_pop(stack);
    regionPos[1] = ds_stack_pop(stack);
    regionPos[0] = ds_stack_pop(stack);
    bufferPos = ds_stack_pop(stack);
    regionSize *= 2;
}
ds_stack_destroy(stack);
returnX = modelX + lineEnd[0] / transformScale;
returnY = modelY + lineEnd[1] / transformScale;
returnZ = modelZ + lineEnd[2] / transformScale;
It is in essence still a recursive script, but instead of jumping deeper and deeper into recursion depth, it jumps back and forth within the while-loop itself in a controlled manner. It performs the same number of checks in the exact same order as it did while it was recursive, but it's a lot faster since it doesn't have to pass variables back and forth.
 
Last edited:
M

MishMash

Guest
Just wanted to spice this topic up a bit. First of all, it is not possible for every recursive function to be represented iteratively (it is possible for any iterative function to be represented recursively however). What I mean here isn't so much to do with "recursion" in terms of function calls, but recursion in terms of using a stack. (When I say recursion, I'm talking about either recursion using the system stack with function calls, and also iteration over some form of stack structure).

Fundamentally, recursion is essential if you NEED to maintain your old stack frames, that is, if the computation in older frames needs to be re-visited and re-used later on. Any function where this information can be thrown away can be made iterative. In some cases, you can also cheat the system with Memoisation in the case of dynamic programming optimisation, which essentially stores previous results in an alternative datastructure like a grid.

In the case of navigating an octree, it is a little harder given the multiple possible end-states you can end up in, however given that you have full information of the structure, and once you have reached a particular layer, you no longer care about the previous layer, it is reasonable to assume an iterative solution that is simpler than emulating recursion exists.

A more efficient iterative solution would require knowledge of how you implemented your octree, as there are many different ways of doing it:
- Whether you have a fixed depth, and at each depth, you store a list of intersecting primitives:
- Even for this one, each node in the octree could either store a cascading subset of primitives, or only the final layer could store primitives, with it being possible for a triangle to exist in the list belonging to multiple nodes.
(In this case, empty nodes would just end, however any node with a triangle in it would reach maximum spacial depth).
- Whether you have an infinite depth, and each branch ends on one single primitive.
- Whether you have a minimal spacial size, and you're quantizing your model down to a voxel style structure where each node simply represents a "yes-no" collision, and you don't actually test on the original geometry, just whether there is or isn't something in a given cube at each level.

As a note, neither of these would be massively far off a stack, implementations I have done in the past have worked off of building up a list of visited nodes in the tree, and then maintaining a list of final triangle IDs to test against. However, they can be faster than a stack-based method as you have both control over the search pattern and the type of data you are storing.

As a rough example algorithm, assuming each node in your octree is accessible via an index, and that index subsequently points to 8 other indices, you can simply maintain a list of nodes to visit, then at each node, test whether any of the 8 sub-sections are valid (i.e. there is an intersection with the cube), if they are, add them to the "to-visit" list. Loop through until the to-visit list is empty whilst simultaneously building up a list of triangle IDs to test against later on (presumably found at the lowest level nodes).

Then once you have accumulated all of the triangles in all of the boxes that a ray intersects with through the first pass, you can then perform your nearest intersection on that small subset of triangles extracted from the model. Or alternatively, you dont nede to make a list of triangles, but can instead do the test in the previous step, maintaining a global reference to the currently closest intersection.

Using this is more optimal for a few reasons:
- We can minimise iterations, and therefore jumps by removal of the "continue" steps.
- We can control whether the structure is searched depth-first or breadth-first based on where we insert parameters in the "to-visit" list. (This may have implications depending on the situation or re-use of algorithm. Depth first for example is faster if you just want to check if the ray collides with anything at all, whereas breadth first might be a more pragmatic of finding the specific nearest intersection).
- Minimise the amount of data that needs to be stored in the "stack", as the list only needs to store IDs. This can save a fair bit of time from function calls.

There's probably no need for you to actually go away and implement this, your current solution is perfectly reasonable anyway! Just wanted to share some thoughts on the matter!
 

GMWolf

aka fel666
First of all, it is not possible for every recursive function to be represented iteratively
Pretty sure you can! After all, How does the compiler do it?
you have to use a stack most of the time though.
There is a theorem about it. dont remember which one.
 

FrostyCat

Redemption Seeker
First of all, it is not possible for every recursive function to be represented iteratively (it is possible for any iterative function to be represented recursively however).
Pretty sure you can! After all, How does the compiler do it?
you have to use a stack most of the time though.
There is a theorem about it. dont remember which one.
All recursive functions can be represented iteratively, it is a direct consequence of the Church-Turing Thesis. However, the choice of language can make it more difficult to convert some of them. With imperative languages, functions that are already tail-recursive generally have the easiest time converting, and the opening post demonstrates the general form that the conversion uses. Those that aren't are typically rewritten in tail-recursive form and then converted from there, which can either be manual or compiler-mediated.
 
M

MishMash

Guest
Pretty sure you can! After all, How does the compiler do it?
you have to use a stack most of the time though.
There is a theorem about it. dont remember which one.
In the very next sentence: "What I mean here isn't so much to do with "recursion" in terms of function calls, but recursion in terms of using a stack. (When I say recursion, I'm talking about either recursion using the system stack with function calls, and also iteration over some form of stack structure)." I clarify that I was talking about creating an iterative representation that does not use a stack. Though it depends on how you look at it, a stack could be considered any structure that requires space complexity to store progressive information that the algorithm uses during computation. For example, you can traverse tree structures non-recursively and without using a specific stack datastructure, but each node of the tree might then need to be extended to store some progressive value. This is not true for all algorithms, and I may be wrong, but I think the Ackerman function is one definitive example of a problem that needs some form of structure to store progressive values.

So to conclude, what I mean is that the concept of what "recursion" is, is merely a form of describing an algorithm. Iterative forms are one type of function that can be represented in this form. So I guess a better discussion would be on the space complexity of an algorithm, rather than this. (As some algorithms require more space to store intermediate values than the size of the structure. Not necessarily in terms of bytes, but in terms of scalability)


Edit: So yeah, I realise my definition of recursion is incorrect, I was talking about what most programmers refer to as recursion and its implications, but the point I meant to make was with regard to the space complexity of the algorithm. Where say an algorithm with a running counter has a space complexity of O(1), however there are many algorithm's which demand a stack, which could be O(n^2) space complexity. (As my thinking is, if you are just implementing it with a stack anyway, that's functionally the same as what the hardware is doing with nested function calls and creating stack frames of the local variables, which is what you are trying to avoid by converting how the algorithm is represented.)
 
Top