As a note, if you happen to be working with OpenGL for newer (2001+) graphics cards, and you are using the GL_ARB_vertex_program extension, you can:
1. create simple bulge effects (like the cast modifier in blender), for things eating spheres
2. create a 'push cast', which looks like something struggling inside (think pushing against a cloth, only in this case, an arbitrary predator represented by a polygonal mesh)

I'm working on skinning with it at the moment, which is rather difficult but highly necessary; since the iPhone can use Matrix Palette, it's worthwhile to build an analog in the C++ world I have, especially since hardware accelerated armatures + nifty bulging effects means the CPU won't die!

As a side note, this is one of my favorite things, 'Vore Mathematics'. As you all probably know and have tried, it's conceptually simple, but practically quite difficult. It applies a lot in rendering, but is also extremely kick-ass in realtime games. It's a lot easier now thanks to newer graphics cards and things like shaders and GLSL though.

Anyways, the proof is in the pudding:

#Bulge distortion (can be furthur optimized to use less TEMPS)###################
TEMP wNorm, wPos, blPr, blDm, blPf;
PARAM blZ = {0,0,0,0}; #Zero vector
PARAM blB = program.local[4]; #Bulge position in local space, w is radius. Pass this through your program
PARAM blS = {0.5,0.5,0.5,0.5}; #Bulge scaling factor (could also pass this as well, but 0.5 looks best

#Set initial position + normal (iPos is the vertex.position value, initially, iNormal is the same thing)
MOV wPos, iPos;
MOV wNorm, iNormal;

#Bulge distortion from a center, radius, and factor
SUB blPr, wPos, blB; #D = P - B,
DP3 blDm.w, blPr, blPr; #D dot D = |D|^2
RSQ blDm.w, blDm.w; #Dm = 1/sqrt(mag^2)
MUL blPr.xyz, blDm.w, blPr; #D * 1/sqrt(mag^2) => D = normal = Pr2
RCP blDm.w, blDm.w; #Dm = sqrt(mag^2)...
SUB blDm.w, blB.w, blDm.w;#Get invert distance (to move by), R - Dmag (R to -inf)
SGE blPf.w, blDm.w, blZ.x; #(R - sqrt(mag)) >= 0 ? in sphere?
MUL blPr, blPr, blPf.wwww; #Dm *= [1,0...] zero out current normal, or keep it

#Should fix normal here, via adding dist.normal * factor to current normal, and normalizing... but this is optional
#MOV temp, blPr; #Grab input normal
#MUL temp, temp, blS; #Scale input normal
#ADD wNorm, wNorm, temp; #Add normal vectors
#DP3 temp, wNorm, wNorm; #D dot D = |D|^2
#RSQ temp, temp.x; #Dm = 1/sqrt(mag^2)
#MUL wNorm, wNorm, temp; #Renormalize normal

MUL blPr, blPr, blS; #Dm *= [scale] Apply factor to normal vector (smoothing)
MUL blPf, blPr, blDm.wwww; #Dm * Prn Apply distance factor to normal

ADD wPos.xyz, blPf, wPos; #pos += Pf Add distance vec to position

#...continue with normal vertex program here, end of bulge#############

And my personal favorite, the Push Cast bulge (imagine this working along a spinal curve, it's pretty glorious in practice)

############################
#Push-Cast deformation requires 2 points and 2 params, 5 temps:
PARAM B = program.local[4]; #B.x,B.y,B.z, Radius;
PARAM C = program.local[5]; #C.x,C.y,C.z, Factor;
PARAM K = {1, -1, 0, 1.656854249492380195206754896837}; #1,-1,0,K for bezier curve

TEMP FVB, DV, RV, RN, FVT; #FVT: x = RS, y = DT, z = DM, w = t; & temp1, temp2

MOV FVB, K.zzzz; #Zero out flags
MUL FVT.x, B.wwww, B.wwww; #Calculate R^2
SUB DV, B, wPos; #Vector from point to bulge center
SUB RV, C, B; #RV = C - B Vector from point int direction of bulge
DP3 temp1, RV, RV; #RN = RV.Normalize(); Normalize RV
RSQ temp1, temp1.x;
MUL RN, RV, temp1;
DP3 temp1, RV, DV; #RV dot DV
DP3 temp2, RV, RV; #|RV|^2
RCP temp2, temp2.x; #1/|RV|^2
MUL temp2, temp2, K.yyyy; #-1/|RV|^2
MUL FVT.y, temp1, temp2; #DT = (-1/|RV|^2) * (RV dot DV)
SGE FVB, FVT.yyyy, K.zzzz; #if( DT > 0 )
MUL FVT.y, FVT.yyyy, FVB; #float DT = (DT > 1) ? 1 : ((DT < 0) ? 0 : DT) Time clamp
MIN FVT.y, FVT.yyyy, K.xxxx;
MUL temp1, FVT.yyyy, RV; #DV += DT * RV;
ADD DV, DV, temp1;
DP3 FVT.z, DV, DV; #float DM = DV.MagnitudeSquared();
SLT temp1, FVT.zzzz, FVT.xxxx; #if( DM < RadiusSq ){
MUL FVB, FVB, temp1.xxxx;
SLT temp2, FVT.yyyy, K.xxxx; #temp2 if( DT < 1 ){ TLT = temp2
SUB temp1, K.xxxx, FVT.yyyy; #float temp = (Factor*(1 - DT));
MUL temp1, temp1, C.wwww;
MUL temp1, temp1, RV; #wPos += temp*RV; //Push to end of cylinder (flat area) #QUACK, dangerous?
MUL temp1, temp1, temp2;
MUL temp1, temp1, FVB.xxxx;
ADD wPos.xyz, wPos, temp1; #}else{
DP3 temp1, DV, RN; #float temp = -DV.DotProduct( RN );
MUL temp1, temp1, temp1; #temp *= temp;
SUB temp2, K.xxxx, temp2;
MUL temp1, temp1, temp2;
SUB FVT.z, FVT.zzzz, temp1; #DM -= temp; } }
MOV temp1, FVT.xxxx; #float n = DM / RadiusSq;
RCP temp1, temp1.x;
MUL FVT.z, FVT.zzzz, temp1; #n = DM
SUB FVT.w, K.xxxx, FVT.zzzz; #float t = 1 - n;
MUL temp1, FVT.wwww, vOne.zzzz; #float F = t*( n*( t*3 + n*K ) + t*t );
MUL temp2, FVT.zzzz, K.w;
ADD temp1, temp1, temp2;
MUL temp1, temp1, FVT.zzzz;
MUL temp2, FVT.wwww, FVT.wwww;
ADD temp1, temp1, temp2;
MUL temp2, temp1, FVT.wwww;
MUL temp2, temp2, C.wwww; #F *= Factor * Radius;
MUL temp2, temp2, B.wwww;
MUL temp2.x, temp2.xxxx, FVB.xxxx;
MUL FVB, RN, temp2.xxxx; #wPos += F * RN; Add to final deformation

#Fix normal vector here?

ADD wPos.xyz, wPos, FVB; #Apply
############################

If you care, or are interested, I have all these and more implemented in C++ code so you can actually read them easier.

Note these vertex programs are NOT complete, because they are missing the header and OUTPUT values required. refer to any reference online to see what is missing, or message me for more info.

Wish me luck!

-Z

Hey Z, whatcha cooking up in there man? Are you working on a vore game as well? I tried looking through your blog but I can't seem to find exactly where you describe it all. What's the specs on it and how's it coming along?