Distorted Android

//////////////////////////////////////////////////////////////////////////////////////////////

// Copyright 2016 Leszek Koltunski                                                          //

//                                                                                          //

// This file is part of Distorted.                                                          //

//                                                                                          //

// Distorted is free software: you can redistribute it and/or modify                        //

// it under the terms of the GNU General Public License as published by                     //

// the Free Software Foundation, either version 2 of the License, or                        //

// (at your option) any later version.                                                      //

//                                                                                          //

// Distorted is distributed in the hope that it will be useful,                             //

// but WITHOUT ANY WARRANTY; without even the implied warranty of                           //

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                            //

// GNU General Public License for more details.                                             //

//                                                                                          //

// You should have received a copy of the GNU General Public License                        // 

// along with Distorted.  If not, see <http://www.gnu.org/licenses/>.                       //

//////////////////////////////////////////////////////////////////////////////////////////////

uniform vec3 u_objD;                 // half of object width x half of object height X half the depth;

                                     // point (0,0,0) is the center of the object

uniform float u_Depth;               // max absolute value of v.z ; beyond that the vertex would be culled by the near or far planes.

                                     // I read OpenGL ES has a built-in uniform variable gl_DepthRange.near = n,

                                     // .far = f, .diff = f-n so maybe u_Depth is redundant

                                     // Update: this struct is only available in fragment shaders

uniform mat4 u_MVPMatrix;            // A constant representing the combined model/view/projection matrix.

uniform mat4 u_MVMatrix;             // A constant representing the combined model/view matrix.

attribute vec3 a_Position;           // Per-vertex position information we will pass in.

attribute vec3 a_Normal;             // Per-vertex normal information we will pass in.

attribute vec2 a_TexCoordinate;      // Per-vertex texture coordinate information we will pass in.

varying vec3 v_Position;             //

varying vec3 v_Normal;               //

varying vec2 v_TexCoordinate;        //

uniform int vNumEffects;             // total number of vertex effects

#if NUM_VERTEX>0

uniform int vType[NUM_VERTEX];       // their types.

uniform vec4 vUniforms[3*NUM_VERTEX];// i-th effect is 3 consecutive vec4's: [3*i], [3*i+1], [3*i+2].

                                     // The first vec4 is the Interpolated values,

                                     // next is half cache half Center, the third -  the Region.

#endif

#if NUM_VERTEX>0

//////////////////////////////////////////////////////////////////////////////////////////////

// HELPER FUNCTIONS

//////////////////////////////////////////////////////////////////////////////////////////////

// The trick below is the if-less version of the

//

// t = dx<0.0 ? (u_objD.x-v.x) / (u_objD.x-ux) : (u_objD.x+v.x) / (u_objD.x+ux);

// h = dy<0.0 ? (u_objD.y-v.y) / (u_objD.y-uy) : (u_objD.y+v.y) / (u_objD.y+uy);

// d = min(t,h);

//

// float d = min(-ps.x/(sign(ps.x)*u_objD.x+p.x),-ps.y/(sign(ps.y)*u_objD.y+p.y))+1.0;

//

// We still have to avoid division by 0 when p.x = +- u_objD.x or p.y = +- u_objD.y (i.e on the edge of the Object)

// We do that by first multiplying the above 'float d' with sign(denominator1*denominator2)^2.

//

//////////////////////////////////////////////////////////////////////////////////////////////

// return degree of the point as defined by the bitmap rectangle

float degree_bitmap(in vec2 S, in vec2 PS)

  {

  vec2 A = sign(PS)*u_objD.xy + S;

  vec2 signA = sign(A);                           //

  vec2 signA_SQ = signA*signA;                    // div = PS/A if A!=0, 0 otherwise.

  vec2 div = signA_SQ*PS/(A-(vec2(1,1)-signA_SQ));//

  return 1.0-max(div.x,div.y);

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// Return degree of the point as defined by the Region. Currently only supports circular regions.

//

// Let us first introduce some notation.

// Let 'PS' be the vector from point P (the current vertex) to point S (the center of the effect).

// Let region.xy be the vector from point S to point O (the center point of the region circle)

// Let region.z be the radius of the region circle.

// (This all should work regardless if S is inside or outside of the circle).

//

// Then, the degree of a point with respect to a given (circular!) Region is defined by:

//

// If P is outside the circle, return 0.

// Otherwise, let X be the point where the halfline SP meets the region circle - then return |PX|/||SX|,

// aka the 'degree' of point P.

//

// We solve the triangle OPX.

// We know the lengths |PO|, |OX| and the angle OPX, because cos(OPX) = cos(180-OPS) = -cos(OPS) = -PS*PO/(|PS|*|PO|)

// then from the law of cosines PX^2 + PO^2 - 2*PX*PO*cos(OPX) = OX^2 so PX = -a + sqrt(a^2 + OX^2 - PO^2)

// where a = PS*PO/|PS| but we are really looking for d = |PX|/(|PX|+|PS|) = 1/(1+ (|PS|/|PX|) ) and

// |PX|/|PS| = -b + sqrt(b^2 + (OX^2-PO^2)/PS^2) where b=PS*PO/|PS|^2 which can be computed with only one sqrt.

float degree_region(in vec4 region, in vec2 PS)

  {

  vec2 PO  = PS + region.xy;

  float D = region.z*region.z-dot(PO,PO);      // D = |OX|^2 - |PO|^2

  if( D<=0.0 ) return 0.0;

  float ps_sq = dot(PS,PS);

  float one_over_ps_sq = 1.0/(ps_sq-(sign(ps_sq)-1.0));  // return 1.0 if ps_sq = 0.0

                                                         // Important: if we want to write

                                                         // b = 1/a if a!=0, b=1 otherwise

                                                         // we need to write that as

                                                         // b = 1 / ( a-(sign(a)-1) )

                                                         // [ and NOT 1 / ( a + 1 - sign(a) ) ]

                                                         // because the latter, if 0<a<2^-24,

                                                         // will suffer from round-off error and in this case

                                                         // a + 1.0 = 1.0 !! so 1 / (a+1-sign(a)) = 1/0 !

  float DOT  = dot(PS,PO)*one_over_ps_sq;

  return 1.0 / (1.0 + 1.0/(sqrt(DOT*DOT+D*one_over_ps_sq)-DOT));

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// return min(degree_bitmap,degree_region). Just like degree_region, currently only supports circles.

float degree(in vec4 region, in vec2 S, in vec2 PS)

  {

  vec2 PO  = PS + region.xy;

  float D = region.z*region.z-dot(PO,PO);      // D = |OX|^2 - |PO|^2

  if( D<=0.0 ) return 0.0;

  vec2 A = sign(PS)*u_objD.xy + S;

  vec2 signA = sign(A);

  vec2 signA_SQ = signA*signA;

  vec2 div = signA_SQ*PS/(A-(vec2(1,1)-signA_SQ));

  float E = 1.0-max(div.x,div.y);

  float ps_sq = dot(PS,PS);

  float one_over_ps_sq = 1.0/(ps_sq-(sign(ps_sq)-1.0));  // return 1.0 if ps_sq = 0.0

  float DOT  = dot(PS,PO)*one_over_ps_sq;

  return min(1.0/(1.0 + 1.0/(sqrt(DOT*DOT+D*one_over_ps_sq)-DOT)),E);

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// Clamp v.z to (-u_Depth,u_Depth) with the following function:

// define h to be, say, 0.7; let H=u_Depth

//      if v.z < -hH then v.z = (-(1-h)^2 * H^2)/(v.z+(2h-1)H) -H   (function satisfying f(-hH)=-hH, f'(-hH)=1, lim f(x) = -H)

// else if v.z >  hH then v.z = (-(1-h)^2 * H^2)/(v.z-(2h-1)H) +H   (function satisfying f(+hH)=+hH, f'(+hH)=1, lim f(x) = +H)

// else v.z = v.z

void restrictZ(inout float v)

  {

  const float h = 0.7;

  float signV = 2.0*max(0.0,sign(v))-1.0;

  float c = ((1.0-h)*(h-1.0)*u_Depth*u_Depth)/(v-signV*(2.0*h-1.0)*u_Depth) +signV*u_Depth;

  float b = max(0.0,sign(abs(v)-h*u_Depth));

  v = b*c+(1.0-b)*v; // Avoid branching: if abs(v)>h*u_Depth, then v=c; otherwise v=v.

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// DEFORM EFFECT

//

// Deform the whole shape of the Object by force V

// 

// If the point of application (Sx,Sy) is on the upper edge of the Object, then:

// a) ignore Vz

// b) change shape of the whole Object in the following way:

//    Suppose the upper-left corner of the Object rectangle is point L, upper-right - R, force vector V

//    is applied to point M on the upper edge, width of the Object = w, height = h, |LM| = Wl, |MR| = Wr,

//    force vector V=(Vx,Vy). Also let H = h/(h+Vy)

//

//    Let now L' and R' be points such that vec(LL') = Wr/w * vec(V) and vec(RR') = Wl/w * vec(V)

//    now let Vl be a point on the line segment L --> M+vec(V) such that Vl(y) = L'(y)

//    and let Vr be a point on the line segment R --> M+vec(V) such that Vr(y) = R'(y)

//    

//    Now define points Fl and Fr, the points L and R will be moved to under force V, with Fl(y)=L'(y)

//    and Fr(y)=R'(y) and |VrFr|/|VrR'| = |VlFl|/|VlL'| = H

//    Now notice that |VrR'| = |VlL'| = Wl*Wr / w   ( a little geometric puzzle! )

//

//    Then points L,R under force V move by vectors vec(Fl), vec(Fr) where

//    vec(Fl) = (Wr/w) * [ (Vx+Wl)-Wl*H, Vy ] = (Wr/w) * [ Wl*Vy / (h+Vy) + Vx, Vy ]

//    vec(Fr) = (Wl/w) * [ (Vx-Wr)+Wr*H, Vy ] = (Wl/w) * [-Wr*Vy / (h+Vy) + Vx, Vy ]

//

//    Lets now denote M+vec(V) = M'. The line segment LMR gets distorted to the curve Fl-M'-Fr. Let's

//    now arbitrarilly decide that:

//    a) at point Fl the curve has to be parallel to line LM'

//    b) at point M' - to line LR

//    c) at point Fr - to line M'R

//

//    Now if Fl=(flx,fly) , M'=(mx,my) , Fr=(frx,fry); direction vector at Fl is (vx,vy) and at M'

//    is (+c,0) where +c is some positive constant, then  the parametric equations of the Fl--->M'

//    section of the curve (which has to satisfy (X(0),Y(0)) = Fl, (X(1),Y(1))=M',

//    (X'(0),Y'(0)) = (vx,vy), (X'(1),Y'(1)) = (+c,0) ) is

//

//    X(t) = ( (mx-flx)-vx )t^2 + vx*t + flx                                  (*)

//    Y(t) = ( vy - 2(my-fly) )t^3 + ( 3(my-fly) -2vy )t^2 + vy*t + fly

//

//    Here we have to have X'(1) = 2(mx-flx)-vx which is positive <==> vx<2(mx-flx). We also have to

//    have vy<2(my-fly) so that Y'(t)>0 (this is a must otherwise we have local loops!)

//    Similarly for the Fr--->M' part of the curve we have the same equation except for the fact that

//    this time we have to have X'(1)<0 so now we have to have vx>2(mx-frx).

//

//    If we are stretching the left or right edge of the bitmap then the only difference is that we

//    have to have (X'(1),Y'(1)) = (0,+-c) with + or - c depending on which part of the curve

//    we are tracing. Then the parametric equation is

//

//    X(t) = ( vx - 2(mx-flx) )t^3 + ( 3(mx-flx) -2vx )t^2 + vx*t + flx

//    Y(t) = ( (my-fly)-vy )t^2 + vy*t + fly

//

//    If we are dragging the top edge:    

//

//    Then point (x,h/2) on the top edge will move by vector (X(t),Y(t)) where those functions are

//    given by (*) and t =  x < dSx ? (w/2+x)/(w/2+dSx) : (w/2-x)/(w/2-dSx)    (-w/2 < x < +w/2 !)

//    (this is 'vec2 time' below in the code).

//    Any point (x,y) will move by vector (a*X(t),a*Y(t)) where a is (y+h/2)/h

void deform(in int effect, inout vec4 v)

  {

  const float A = 0.5;

  const float B = 0.3;

  vec2 center = vUniforms[effect+1].yz;

  vec2 force  = vUniforms[effect].xy;

  vec2 dist   = center-v.xy;

  vec2 aDist  = abs(dist.xy);

  vec2 maxdist= u_objD.xy + abs(center);

  vec2 aForce = abs(force);

  vec2 Aw = A*maxdist;

  vec2 quot = dist / maxdist;

  float mvXvert =-((B*dist.x*aForce.y)/(aForce.y + Aw.x))*(1.0-quot.y*quot.y);

  float mvYvert = force.y*( 1.0 - quot.x*quot.x*(Aw.x/(aForce.y+Aw.x)) ) * aForce.y/(aForce.y+aDist.y);

  float mvXhorz =-force.x*( 1.0 - quot.y*quot.y*(Aw.y/(aForce.x+Aw.y)) ) * aForce.x/(aForce.x+aDist.x);

  float mvYhorz =-((B*dist.y*aForce.x)/(aForce.x + Aw.y))*(1.0-quot.x*quot.x);

  v.x -= (mvXvert+mvXhorz);

  v.y -= (mvYvert+mvYhorz);

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// DISTORT EFFECT

//

// Point (Px,Py) gets moved by vector (Wx,Wy,Wz) where Wx/Wy = Vx/Vy i.e. Wx=aVx and Wy=aVy where 

// a=Py/Sy (N --> when (Px,Py) is above (Sx,Sy)) or a=Px/Sx (W) or a=(w-Px)/(w-Sx) (E) or a=(h-Py)/(h-Sy) (S) 

// It remains to be computed which of the N,W,E or S case we have: answer: a = min[ Px/Sx , Py/Sy , (w-Px)/(w-Sx) , (h-Py)/(h-Sy) ]

// Computations above are valid for screen (0,0)x(w,h) but here we have (-w/2,-h/2)x(w/2,h/2)

//  

// the vertical part

// Let |(v.x,v.y),(ux,uy)| = |PS|, ux-v.x=dx,uy-v.y=dy, f(x) (0<=x<=|SX|) be the shape of the side of the bubble.

// H(v.x,v.y) = |PS|>|SX| ? 0 : f(|PX|)

// N(v.x,v.y) = |PS|>|SX| ? (0,0,1) : ( -(dx/|PS|)sin(beta), -(dy/|PS|)sin(beta), cos(beta) ) where tan(beta) is f'(|PX|) 

// ( i.e. normalize( dx, dy, -|PS|/f'(|PX|))         

//

// Now we also have to take into account the effect horizontal move by V=(u_dVx[i],u_dVy[i]) will have on the normal vector.

// Solution: 

// 1. Decompose the V into two subcomponents, one parallel to SX and another perpendicular.

// 2. Convince yourself (draw!) that the perpendicular component has no effect on normals.

// 3. The parallel component changes the length of |SX| by the factor of a=(|SX|-|Vpar|)/|SX| (where the length

//    can be negative depending on the direction)

// 4. that in turn leaves the x and y parts of the normal unchanged and multiplies the z component by a!

//

// |Vpar| = (u_dVx[i]*dx - u_dVy[i]*dy) / sqrt(ps_sq) = (Vx*dx-Vy*dy)/ sqrt(ps_sq)  (-Vy because y is inverted)

// a =  (|SX| - |Vpar|)/|SX| = 1 - |Vpar|/((sqrt(ps_sq)/(1-d)) = 1 - (1-d)*|Vpar|/sqrt(ps_sq) = 1-(1-d)*(Vx*dx-Vy*dy)/ps_sq 

//

// Side of the bubble

// 

// choose from one of the three bubble shapes: the cone, the thin bubble and the thick bubble          

// Case 1: 

// f(t) = t, i.e. f(x) = uz * x/|SX|   (a cone)

// -|PS|/f'(|PX|) = -|PS|*|SX|/uz but since ps_sq=|PS|^2 and d=|PX|/|SX| then |PS|*|SX| = ps_sq/(1-d)

// so finally -|PS|/f'(|PX|) = -ps_sq/(uz*(1-d))

//                    

// Case 2: 

// f(t) = 3t^2 - 2t^3 --> f(0)=0, f'(0)=0, f'(1)=0, f(1)=1 (the bell curve)

// here we have t = x/|SX| which makes f'(|PX|) = 6*uz*|PS|*|PX|/|SX|^3.

// so -|PS|/f'(|PX|) = (-|SX|^3)/(6uz|PX|) =  (-|SX|^2) / (6*uz*d) but

// d = |PX|/|SX| and ps_sq = |PS|^2 so |SX|^2 = ps_sq/(1-d)^2

// so finally -|PS|/f'(|PX|) = -ps_sq/ (6uz*d*(1-d)^2)

//                  

// Case 3:

// f(t) = 3t^4-8t^3+6t^2 would be better as this satisfies f(0)=0, f'(0)=0, f'(1)=0, f(1)=1,

// f(0.5)=0.7 and f'(t)= t(t-1)^2 >=0 for t>=0 so this produces a fuller, thicker bubble!

// then -|PS|/f'(|PX|) = (-|PS|*|SX)) / (12uz*d*(d-1)^2) but |PS|*|SX| = ps_sq/(1-d) (see above!) 

// so finally -|PS|/f'(|PX|) = -ps_sq/ (12uz*d*(1-d)^3)  

//

// Now, new requirement: we have to be able to add up normal vectors, i.e. distort already distorted surfaces.

// If a surface is given by z = f(x,y), then the normal vector at (x0,y0) is given by (-df/dx (x0,y0), -df/dy (x0,y0), 1 ).

// so if we have two surfaces defined by f1(x,y) and f2(x,y) with their normals expressed as (f1x,f1y,1) and (f2x,f2y,1) 

// then the normal to g = f1+f2 is simply given by (f1x+f2x,f1y+f2y,1), i.e. if the third components are equal, then we

// can simply add up the first and second components.

//

// Thus we actually want to compute N(v.x,v.y) = a*(-(dx/|PS|)*f'(|PX|), -(dy/|PS|)*f'(|PX|), 1) and keep adding

// the first two components. (a is the horizontal part)

void distort(in int effect, inout vec4 v, inout vec4 n)

  {

  vec2 center = vUniforms[effect+1].yz;

  vec2 ps = center-v.xy;

  vec3 force = vUniforms[effect].xyz;

  float d = degree(vUniforms[effect+2],center,ps);

  float denom = dot(ps+(1.0-d)*force.xy,ps);

  float one_over_denom = 1.0/(denom-0.001*(sign(denom)-1.0));          // = denom==0 ? 1000:1/denom;

  //v.z += force.z*d;                                                  // cone

  //b = -(force.z*(1.0-d))*one_over_denom;                             //

  //v.z += force.z*d*d*(3.0-2.0*d);                                    // thin bubble

  //b = -(6.0*force.z*d*(1.0-d)*(1.0-d))*one_over_denom;               //

  v.z += force.z*d*d*(3.0*d*d -8.0*d +6.0);                            // thick bubble

  float b = -(12.0*force.z*d*(1.0-d)*(1.0-d)*(1.0-d))*one_over_denom;  //

  v.xy += d*force.xy;

  n.xy += n.z*b*ps;

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// SINK EFFECT

//

// Pull P=(v.x,v.y) towards center of the effect with P' = P + (1-h)*dist(S-P)

// when h>1 we are pushing points away from S: P' = P + (1/h-1)*dist(S-P)

void sink(in int effect,inout vec4 v)

  {

  vec2 center = vUniforms[effect+1].yz;

  vec2 ps = center-v.xy;

  float h = vUniforms[effect].x;

  float t = degree(vUniforms[effect+2],center,ps) * (1.0-h)/max(1.0,h);

  v.xy += t*ps;           

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// PINCH EFFECT

//

// Pull P=(v.x,v.y) towards the line that

// a) passes through the center of the effect

// b) forms angle defined in the 2nd interpolated value with the X-axis

// with P' = P + (1-h)*dist(line to P)

// when h>1 we are pushing points away from S: P' = P + (1/h-1)*dist(line to P)

void pinch(in int effect,inout vec4 v)

  {

  vec2 center = vUniforms[effect+1].yz;

  vec2 ps = center-v.xy;

  float h = vUniforms[effect].x;

  float t = degree(vUniforms[effect+2],center,ps) * (1.0-h)/max(1.0,h);

  float angle = vUniforms[effect].y;

  vec2 dir = vec2(sin(angle),-cos(angle));

  v.xy += t*dot(ps,dir)*dir;

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// SWIRL EFFECT

//

// Let d be the degree of the current vertex V with respect to center of the effect S and Region vRegion.

// This effect rotates the current vertex V by vInterpolated.x radians clockwise around the circle dilated 

// by (1-d) around the center of the effect S.

void swirl(in int effect, inout vec4 v)

  {

  vec2 center  = vUniforms[effect+1].yz;

  vec2 PS = center-v.xy;

  vec4 SO = vUniforms[effect+2];

  float d1_circle = degree_region(SO,PS);

  float d1_bitmap = degree_bitmap(center,PS);

  float alpha = vUniforms[effect].x;

  float sinA = sin(alpha);

  float cosA = cos(alpha);

  vec2 PS2 = vec2( PS.x*cosA+PS.y*sinA,-PS.x*sinA+PS.y*cosA ); // vector PS rotated by A radians clockwise around center.

  vec4 SG = (1.0-d1_circle)*SO;                                // coordinates of the dilated circle P is going to get rotated around

  float d2 = max(0.0,degree(SG,center,PS2));                   // make it a max(0,deg) because otherwise when center=left edge of the

                                                               // bitmap some points end up with d2<0 and they disappear off view.

  v.xy += min(d1_circle,d1_bitmap)*(PS - PS2/(1.0-d2));        // if d2=1 (i.e P=center) we should have P unchanged. How to do it?

  }

//////////////////////////////////////////////////////////////////////////////////////////////

// WAVE EFFECT

//

// Directional sinusoidal wave effect.

//

// This is an effect from a (hopefully!) generic family of effects of the form (vec3 V: |V|=1 , f(x,y) )  (*)

// i.e. effects defined by a unit vector and an arbitrary function. Those effects are defined to move each

// point (x,y,0) of the XY plane to the point (x,y,0) + V*f(x,y).

//

// In this case V is defined by angles A and B (sines and cosines of which are precomputed in

// EffectQueueVertex and passed in the uniforms).

// Let's move V to start at the origin O, let point C be the endpoint of V, and let C' be C's projection

// to the XY plane. Then A is defined to be the angle C0C' and angle B is the angle C'O(axisY).

//

// Also, in this case f(x,y) = amplitude*sin(x/length), with those 2 parameters passed in uniforms.

//

//////////////////////////////////////////////////////////////////////////////////////////////

// How to compute any generic effect of type (*)

//////////////////////////////////////////////////////////////////////////////////////////////

//

// By definition, the vertices move by f(x,y)*V.

//

// Normals are much more complicated.

// Let angle X be the angle (0,Vy,Vz)(0,Vy,0)(Vx,Vy,Vz).

// Let angle Y be the angle (Vx,0,Vz)(Vx,0,0)(Vx,Vy,Vz).

//

// Then it can be shown that the resulting surface, at point to which point (x0,y0,0) got moved to,

// has 2 tangent vectors given by

//

// SX = (1.0+cosX*fx , cosY*sinX*fx , |sinY|*sinX*fx);  (**)

// SY = (cosX*sinY*fy , 1.0+cosY*fy , |sinX|*sinY*fy);  (***)

//

// and then obviously the normal N is given by N= SX x SY .

//

// We still need to remember the note from the distort function about adding up normals:

// we first need to 'normalize' the normals to make their third components equal, and then we

// simply add up the first and the second component while leaving the third unchanged.

//

// How to see facts (**) and (***) ? Briefly:

// a) compute the 2D analogon and conclude that in this case the tangent SX is given by

//    SX = ( cosA*f'(x) +1, sinA*f'(x) )    (where A is the angle vector V makes with X axis )

// b) cut the resulting surface with plane P which

//    - includes vector V

//    - crosses plane XY along line parallel to X axis

// c) apply the 2D analogon and notice that the tangent vector to the curve that is the common part of P

//    and our surface (I am talking about the tangent vector which belongs to P) is given by

//    (1+cosX*fx,0,sinX*fx) rotated by angle (90-|Y|) (where angles X,Y are defined above) along vector (1,0,0).

//

//    Matrix of rotation:

//

//    |sinY|  cosY

//    -cosY  |sinY|

//

// d) compute the above and see that this is equal precisely to SX from (**).

// e) repeat points b,c,d in direction Y and come up with (***).

//

//////////////////////////////////////////////////////////////////////////////////////////////

// Note: we should avoid passing certain combinations of parameters to this function. One such known

// combination is ( A: small but positive, B: any, amplitude >= length ).

// In this case, certain 'unlucky' points have their normals almost horizontal (they got moved by (almost!)

// amplitude, and other point length (i.e. <=amplitude) away got moved by 0, so the slope in this point is

// very steep). Visual effect is: vast majority of surface pretty much unchanged, but random 'unlucky'

// points very dark)

//

// Generally speaking I'd keep to amplitude < length, as the opposite case has some other problems as well.

void wave(in int effect, inout vec4 v, inout vec4 n)

  {

  vec2 center     = vUniforms[effect+1].yz;

  float amplitude = vUniforms[effect  ].x;

  float length    = vUniforms[effect  ].y;

  vec2 ps = center - v.xy;

  float deg = amplitude*degree_region(vUniforms[effect+2],ps);

  if( deg != 0.0 && length != 0.0 )

    {

    float phase = vUniforms[effect  ].z;

    float alpha = vUniforms[effect  ].w;

    float beta  = vUniforms[effect+1].x;

    float sinA = sin(alpha);

    float cosA = cos(alpha);

    float sinB = sin(beta);

    float cosB = cos(beta);

    float angle= 1.578*(ps.x*cosB-ps.y*sinB) / length + phase;

    vec3 dir= vec3(sinB*cosA,cosB*cosA,sinA);

    v.xyz += sin(angle)*deg*dir;

    if( n.z != 0.0 )

      {

      float sqrtX = sqrt(dir.y*dir.y + dir.z*dir.z);

      float sqrtY = sqrt(dir.x*dir.x + dir.z*dir.z);

      float sinX = ( sqrtY==0.0 ? 0.0 : dir.z / sqrtY);

      float cosX = ( sqrtY==0.0 ? 1.0 : dir.x / sqrtY);

      float sinY = ( sqrtX==0.0 ? 0.0 : dir.z / sqrtX);

      float cosY = ( sqrtX==0.0 ? 1.0 : dir.y / sqrtX);

      float abs_z = dir.z <0.0 ? -(sinX*sinY) : (sinX*sinY);

      float tmp = 1.578*cos(angle)*deg/length;

      float fx =-cosB*tmp;

      float fy = sinB*tmp;

      vec3 sx = vec3 (1.0+cosX*fx,cosY*sinX*fx,abs_z*fx);

      vec3 sy = vec3 (cosX*sinY*fy,1.0+cosY*fy,abs_z*fy);

      vec3 normal = cross(sx,sy);

      if( normal.z<=0.0 )                   // Why this bizarre shit rather than the straightforward

        {                                   //

        normal.x= 0.0;                      // if( normal.z>0.0 )

        normal.y= 0.0;                      //   {

        normal.z= 1.0;                      //   n.x = (n.x*normal.z + n.z*normal.x);

        }                                   //   n.y = (n.y*normal.z + n.z*normal.y);

                                            //   n.z = (n.z*normal.z);

                                            //   }

      n.x = (n.x*normal.z + n.z*normal.x);  //

      n.y = (n.y*normal.z + n.z*normal.y);  // ? Because if we do the above, my shitty Nexus4 crashes

      n.z = (n.z*normal.z);                 // during shader compilation!

      }

    }

  }

#endif

//////////////////////////////////////////////////////////////////////////////////////////////

void main()                                                 	

  {              

  vec4 v = vec4( 2.0*u_objD*a_Position,1.0 );

  vec4 n = vec4(a_Normal,0.0);

#if NUM_VERTEX>0

  for(int i=0; i<vNumEffects; i++)

    {

         if( vType[i]==DISTORT) distort(3*i,v,n);

    else if( vType[i]==DEFORM ) deform (3*i,v);

    else if( vType[i]==SINK   ) sink   (3*i,v);

    else if( vType[i]==PINCH  ) pinch  (3*i,v);

    else if( vType[i]==SWIRL  ) swirl  (3*i,v);

    else if( vType[i]==WAVE   ) wave   (3*i,v,n);

    }

  restrictZ(v.z);

#endif

  v_Position      = v.xyz;

  v_TexCoordinate = a_TexCoordinate;

  v_Normal        = normalize(vec3(u_MVMatrix*n));

  gl_Position     = u_MVPMatrix*v;      

  }