Voxel <> IFS

[David Makin]

To some it will be obvious to others maybe not, but if you take 8 simple affine transforms all scaling by 0.5 in all 3D with translations to the corners of a cube centred on the origin then the IFS tree from that can be adapted to generate the equivalent of a voxel object using distance estimation - in fact because the transforms are simple scale and translate any voxel object done this way could be rendered using KIFS i.e. just 8 tests per depth through the IFS with only one path to the next depth.

On this basis I decided to investigate how much memory would be required for how much detail and, as a bonus, check that my understanding of box counting was correct.

So here is a little code - just a test global section and one user parameter as a UF formula (standard ufm not class-based).

At the moment the surftest function simply checks for surface points of a sphere by checking each "voxel" to see if there are corner points on both sides of the spherical surface - if so it's a surface cube, if not then not - this could simply be adapted to a similar routine for any object e.g. a Mandelbulb just testing to see if there are corner points both "inside" and "outside" (with appropriate scaling instead of the current -1...+1 range). Of course this method assumes the object concerned is fully connected.

Code:

voxelIFStest {
global:
    $define debug
    func surftest(float &x,float &y,float &z,float &s)
      if sqr(x-s)+sqr(y-s)+sqr(z-s)<=1
        if sqr(x+s)+sqr(y-s)+sqr(z-s)>1
          return
        endif
        if sqr(x-s)+sqr(y+s)+sqr(z-s)>1
          return
        endif
        if sqr(x+s)+sqr(y+s)+sqr(z-s)>1
          return
        endif
        if sqr(x-s)+sqr(y-s)+sqr(z+s)>1
          return
        endif
        if sqr(x+s)+sqr(y-s)+sqr(z+s)>1
          return
        endif
        if sqr(x-s)+sqr(y+s)+sqr(z+s)>1
          return
        endif
        if sqr(x+s)+sqr(y+s)+sqr(z+s)>1
          return
        endif
      else
        if sqr(x+s)+sqr(y-s)+sqr(z-s)<=1
          return
        endif
        if sqr(x-s)+sqr(y+s)+sqr(z-s)<=1
          return
        endif
        if sqr(x+s)+sqr(y+s)+sqr(z-s)<=1
          return
        endif
        if sqr(x-s)+sqr(y-s)+sqr(z+s)<=1
          return
        endif
        if sqr(x+s)+sqr(y-s)+sqr(z+s)<=1
          return
        endif
        if sqr(x-s)+sqr(y+s)+sqr(z+s)<=1
          return
        endif
        if sqr(x+s)+sqr(y+s)+sqr(z+s)<=1
          return
        endif
      endif
      s = 0
      return
    endfunc
    int a[@depth]
    float b[@depth,8,3]
    float n[@depth]
    b[0,0,0]=-0.5
    b[0,0,1]=-0.5
    b[0,0,2]=-0.5
    b[0,1,0]=0.5
    b[0,1,1]=-0.5
    b[0,1,2]=-0.5
    b[0,2,0]=-0.5
    b[0,2,1]=0.5
    b[0,2,2]=-0.5
    b[0,3,0]=0.5
    b[0,3,1]=0.5
    b[0,3,2]=-0.5
    b[0,4,0]=-0.5
    b[0,4,1]=-0.5
    b[0,4,2]=0.5
    b[0,5,0]=0.5
    b[0,5,1]=-0.5
    b[0,5,2]=0.5
    b[0,6,0]=-0.5
    b[0,6,1]=0.5
    b[0,6,2]=0.5
    b[0,7,0]=0.5
    b[0,7,1]=0.5
    b[0,7,2]=0.5
    float u = 0
    float v = 0
    float w = 0
    int i = 0
    int j = 0
    int k = 0
    float c = 0
    repeat
      a[i] = 0
      n[i] = 0
    until (i=i+1)>=@depth
    i = 0
    j = 0
    float t = 0.5
    float s
    repeat
      s = t
      surftest(b[i,j,0],b[i,j,1],b[i,j,2],s)
      if s!=0
        c = c + 1
        n[i] = n[i] + 1
        a[i]=j
        if (i=i+1)<@depth
          t = 0.5*t
          u = b[i-1,j,0]
          v = b[i-1,j,1]
          w = b[i-1,j,2]
          b[i,0,0]=u-t
          b[i,0,1]=v-t
          b[i,0,2]=w-t
          b[i,1,0]=u+t
          b[i,1,1]=v-t
          b[i,1,2]=w-t
          b[i,2,0]=u-t
          b[i,2,1]=v+t
          b[i,2,2]=w-t
          b[i,3,0]=u+t
          b[i,3,1]=v+t
          b[i,3,2]=w-t
          b[i,4,0]=u-t
          b[i,4,1]=v-t
          b[i,4,2]=w+t
          b[i,5,0]=u+t
          b[i,5,1]=v-t
          b[i,5,2]=w+t
          b[i,6,0]=u-t
          b[i,6,1]=v+t
          b[i,6,2]=w+t
          b[i,7,0]=u+t
          b[i,7,1]=v+t
          b[i,7,2]=w+t
          j = -1
        else
          i = @depth-1
        endif
      endif
      if (j=j+1)>=8
        repeat
          i = i-1
          if i>=0
            j = a[i] + 1
            t = 2.0*t
          endif
        until i<0 || j<8
        if i==0 || i==1
          print("Depth:",i," Box:",j)
        endif
      endif
    until i<0
    print("Total nodes: ",c)
    print("Box dimensions")
    i = 0
    t = 0.5
    repeat
      print(log(n[i])/log(1.0/t))
      t = 0.5*t
    until (i=i+1)>=@depth
default:
  param depth
    default = 7
  endparam
}

And here's the important output for depth=16 for those not wanting to wait:

Total nodes: 26986060952
Box dimensions
(3, 0)
(2.9036774610288021, 0)
(2.6958209470834465, 0)
(2.5449772725037336, 0)
(2.4433491716390612, 0)
(2.3718701967851976, 0)
(2.3191954842509092, 0)
(2.2794082915946075, 0)
(2.2484719798223884, 0)
(2.223631789897514, 0)
(2.2033102274863545, 0)
(2.1863709202239249, 0)
(2.1720347292799661, 0)
(2.1597468979876545, 0)
(2.1490971855356314, 0)
(2.1397786385977124, 0)

So the definition of a sphere in this form to this accuracy requires around 25GB of data (26986060952 bytes) using 8 bits per node (i.e. occupied/not) - the accuracy would give say earth to every 100m or so but only if earth were a perfect sphere, more memory would be required for the real thing since although as can be gleaned from the above the dimension of a sphere's surface is 2, that of real earth is somewhat larger

In perspective that 25GB data could be stored in a disk rendered image output from UF as RGBA at a resolution of 81920*81920 or so.

[David Makin]

Here's the data from the same code but testing the cubic volume around a main bulb of the degree 8 Mandelbulb - of course we're hoping for a dimension close to 3.....


Total nodes: 1106504
Box dimensions
3
2.5849625007211562
2.6741226043428182
2.7464604842508351
2.7850367038868417
2.8161736407773771
2.8403596681867006


I know so far only size 256*256*256 but encouraging on the dimension front - 2.84 and still rising - especially because if anything the calculated value using this method will always be lower than the real value.

[Patryk Kizny]

If you can save the voxel in particle data (just a set of points in space with coordinates and any extra properties, can be ASCI/CSV, whatever) I'd be happy to play with it. I can easily visualize these kinds of datasets up to 1 000 000 000 points in reasonable time.