######################################################################
######################################################################
# 
# isosurface expression functions for houdini

# Matthew Lewis - http://www.accad.ohio-state.edu/~mlewis
# mlewis@accad.ohio-state.edu
#
# Most of the shapes were modified from those found at:
#     http://www.uib.no/People/nfytn/surfaces.htm  (Tore Nordstrand)
#
# version 01.06.05b
#

#################################################
lerp(a,b,t)
{
  return a+t*(b-a);
}

#################################################
smoothInts(v)
{
  fv = floor(v);
  return fv + smooth(v, fv+0.5, fv+1);
}

#################################################
# Perlin bias
bias(b,t)
{
    return t^(log(b)/log(0.5));
}


#################################################
# Perlin gain 
gain(g,t)
{
    return if(t<0.5,(.5*bias(1.0-(g),2.0*(t))),(1.0-.5*bias(1.0-(g),2.0-2.0*(t))));
}

#################################################
# Lawson Wade's reimplementation of Perlin's bias function
# makes it more "symmetric"
symBias(b,t) {
  v = 0;
  if(b<.5) { v = bias(b,t); }
  else { v = 1.0 - bias(1-b, 1-t); }
  return v;
}

#################################################
# Lawson Wade's reimplementation of Perlin's gain function
# makes it more "symmetric"
symGain(g,t) {
  v = 0;
  if(t<.5) { v = .5 * symBias(1-g, 2*t); }
  else { v = 1 - .5 * symBias(1-g, 2-2*t); }
  return v;
}

#################################################
iso_blend_old(f,g)
{
  v = min(f,g); # union
#  b = .001;  v = (1+b) - (b^(f)) - (b^(g)); # smooth blend
  return v;
}

#################################################
# 3sided topo thing
iso_chair(x,y,z)
{
  v = ((((x-.5)*5)^2)+((y-.5)*5^2)+((z-.5)*5^2)-.95*(2.5^2))^2-.8*(((z-.5)*5-2.5)^2-2*((x-.5)*5^2))*(((z-.5)*5+2.5)^2-2*((y-.5)*5^2));
  return v;
}

###########################
# cube
iso_cube(x,y,z)
{
  v = max(abs(x-.5)-.495,max(abs(y-.5)-.495,abs(z-.5)-.495));
  return v;
}

###########################
# radiating rods
iso_kummer(x,y,z)
{
  rods = (x-.5)*6^4+(y-.5)*6^4+(z-.5)*6^4-((x-.5)*6^2)-((y-.5)*6^2)-((z-.5)*6^2)-((x-.5)*6^2)*((y-.5)*6^2)-((x-.5)*6^2)*((z-.5)*6^2)-((y-.5)*6^2)*((z-.5)*6^2)+1;

  cube =  max(abs(x-.5)-.495,max(abs(y-.5)-.495,abs(z-.5)-.495));
  v = max(rods,cube);
  return v;
}

###########################
# hourglass
iso_lem(x,y,z)
{
  v = (y-.5)*2^4 - ((y-.5)*2^2) + ((x-.5)*2^2) + ((z-.5)*2^2);
  return v;
}

###########################
# octahedron
iso_oct(x,y,z)
{
  v = abs((x-.5))+abs((y-.5))+abs((z-.5))-.5;
  return v;
}

###########################
# pillow/tooth thing
iso_pillow(x,y,z)
{
  v = (x-.5)*2.35^4 + (y-.5)*2.35^4 + (z-.5)*2.35^4 - (((x-.5)*2.35^2) + ((y-.5)*2.35^2) + ((z-.5)*2.35^2));
  return v;
}

###########################
# top
iso_qcylinder(x,y,z)
{
  top = (((x-.5)*2^2) + ((z-.5)*2^2)) * ((y-.5)*2^2) + .01 * (((x-.5)*2^2) + ((z-.5)*2^2)) - .01 * 1;

  cube = iso_cube(x,y,z);
  v = max(top,cube);
  return v;
}

###########################
# sphere
iso_sphere(x,y,z)
{
  v = (x-.5)*(x-.5) + (y-.5)*(y-.5) + (z-.5)*(z-.5) - .25;
  return v;
}

###########################
# bulgy wireframe cube, kinda
iso_tangle(x,y,z)
{
  v = (x-.5)*4.6^4 - 5*((x-.5)*4.6^2) + (y-.5)*4.6^4 - 5*((y-.5)*4.6^2) + (z-.5)*4.6^4 - 5*((z-.5)*4.6^2) + 11.8;
  return v;
}

###########################
# torus
iso_torus(x,y,z)
{
  v = ((x-.5)^2+(y-.5)^2+(z-.5)^2+0.325^2-0.175^2)^2 - 4*(0.325^2)*((x-0.5)^2 + (y-0.5)^2);
  return v;
}

#################################################
# blend t of the way between two shapes s1 and s2
# use bias and gain to control transition
iso_blend(s1,s2,t,b,g)
{
  v = lerp(s1, s2, symBias(b, symGain(g, t)));
  return v;
}

######################################################################
# t: interpolation shape, [0..1]
iso_shape(x,y,z,t)
{
  v = 0;
  b = (t*10)%1;   # map tenths into [0..1]

  if(t<.1) {
    v = iso_blend(iso_kummer(x,y,z),iso_tangle(x,y,z),b,.85, .5);
  }
  if((t>=.1)&&(t<.2)) {
    v = iso_blend(iso_tangle(x,y,z),iso_pillow(x,y,z),b,.9, .3);
  }
  if((t>=.2)&&(t<.3)) {
    v = iso_blend(iso_pillow(x,y,z),iso_cube(x,y,z),b,.75, .5);
  }
  if((t>=.3)&&(t<.4)) {
    v = iso_blend(iso_cube(x,y,z),iso_sphere(x,y,z),b,.5, .5);
  }
  if((t>=.4)&&(t<.5)) {
    v = iso_blend(iso_sphere(x,y,z),iso_oct(x,y,z),b,.5, .5);
  }
  if((t>=.5)&&(t<.6)) {
    v = iso_blend(iso_oct(x,y,z),iso_qcylinder(x,y,z),b,.9, .5);
  }
  if((t>=.6)&&(t<.7)) {
    v = iso_blend(iso_qcylinder(x,y,z),iso_lem(x,y,z),b,.1, .5);
  }
  if((t>=.7)&&(t<.8)) {
    v = iso_blend(iso_lem(x,y,z),iso_torus(x,y,z),b,.975, .15);
  }
  if((t>=.8)&&(t<.9)) {
    v = iso_blend(iso_torus(x,y,z),iso_chair(x,y,z),b,.01, .3);
  }
  if(t>=.9) {
    v = iso_blend(iso_chair(x,y,z),iso_kummer(x,y,z),b,.6, .8);
  }

  return v;
}

######################################################################
# translated from rmannotes
######################################################################
repeat(x,freq)
{
  return (x * freq) % 1.0;
}

######################################################################
odd(x)           
{
  return ((x%2) == 1);
}

######################################################################
even(x)          
{
  return ((x%2) == 0);
}


######################################################################
whichtile(x,freq)
{
  return  floor((x) * (freq));
}


######################################################################
nnoise(x)
{
  return noise(x,.5,.5) * .5 + .5;
}


######################################################################
nnoise2(x,y)
{
  return noise(x,y,.5) * .5 + .5;
}

######################################################################
nnoise3(x,y,z)
{
  return noise(x,y,z) * .5 + .5;
}

######################################################################
udn(x,lo,hi)
{
  return  smooth(nnoise(x), .25, .75) * (hi - lo) + lo;
}

######################################################################
udn2(x,y,lo,hi)
{
  return smooth(nnoise2(x,y), .25, .75) * (hi-lo)+lo;
}

######################################################################
udn3(x,y,z,lo,hi)
{
  return smooth(nnoise3(x,y,z), .25, .75) * (hi-lo)+lo;
}

######################################################################
iso_pattern(x,y,z,s,fx,fy,fz)
{
  xx = repeat(x,fx);
  yy = repeat(y,fy);
  zz = repeat(z,fz);

  v = iso_shape(xx,yy,zz,s);
  return v;
}

######################################################################
# s: base shape
# fx,fy,fz: pattern freq
# su: uniform scale bombing [0..1]
iso_bomb(x,y,z,s,fx,fy,fz,su)
{
  xx = repeat(x,fx);
  yy = repeat(y,fy);
  zz = repeat(z,fz);

  xt = whichtile(x,fx);
  yt = whichtile(y,fx);
  zt = whichtile(z,fx);
  sc = udn3(xt+23.5,yt+345.5,zt+892.5,0,1);
  sc = (sc-0.5)*su+0.5;
  if(sc<0.5) { sc = 1-sc; }
  else { sc *= 2; }

  xx = (xx-0.5)*sc + 0.5;
  yy = (yy-0.5)*sc + 0.5;
  zz = (zz-0.5)*sc + 0.5;

  v = iso_shape(xx,yy,zz,s);

  cube = iso_cube(x,y,z,1);
  v = max(v,cube);

  return v;
}

######################################################################