 	mat2 jacob_t = mat2(jacob[0][0], jacob[1][0], jacob[0][1], jacob[1][1]);
 	mat2 sigma_f_inv = (4.*PI*PI*params.filterSigma*params.filterSigma)*(jacob*jacob_t);
 	mat2 sigma_f = inv2x2(sigma_f_inv);
-	mat2 sigma_g_inv = (2.*PI*params.ainv*params.ainv)* id2x2();
+	mat2 sigma_g_inv = (2.*PI*im_params.ainv*im_params.ainv)* id2x2();
 	mat2 sigma_g = inv2x2(sigma_g_inv);
 	mat2 sigma_fg_inv = sigma_f_inv + sigma_g_inv;
 	mat2 sigma_fg = inv2x2(sigma_fg_inv);
 	
 	//filter params
-	params.filter = sigma_fg * sigma_g_inv;	
-	params.sigma_f_plus_g_inv = inv2x2(sigma_f + sigma_g);
-	params.a_prime_square = 2.*PI*sqrt(det2x2(sigma_fg));
-	
-	params.octaves = params.sector.y;//log2(params.sector.y/params.sector.x);
+	im_params.filter = sigma_fg * sigma_g_inv;	
+	im_params.sigma_f_plus_g_inv = inv2x2(sigma_f + sigma_g);
+	im_params.a_prime_square = 2.*PI*sqrt(det2x2(sigma_fg));
 
 	float value = 
-		eval_cell(cpos, gpos, ivec2(-1, -1), params) +
-		eval_cell(cpos, gpos, ivec2(-1,  0), params) +
-		eval_cell(cpos, gpos, ivec2(-1,  1), params) +
-		eval_cell(cpos, gpos, ivec2( 0, -1), params) +
-		eval_cell(cpos, gpos, ivec2( 0,  0), params) +
-		eval_cell(cpos, gpos, ivec2( 0,  1), params) +
-		eval_cell(cpos, gpos, ivec2( 1, -1), params) +
-		eval_cell(cpos, gpos, ivec2( 1,  0), params) +
-		eval_cell(cpos, gpos, ivec2( 1,  1), params);
+		eval_cell(cpos, gpos, ivec2(-1, -1), params, im_params) +
+		eval_cell(cpos, gpos, ivec2(-1,  0), params, im_params) +
+		eval_cell(cpos, gpos, ivec2(-1,  1), params, im_params) +
+		eval_cell(cpos, gpos, ivec2( 0, -1), params, im_params) +
+		eval_cell(cpos, gpos, ivec2( 0,  0), params, im_params) +
+		eval_cell(cpos, gpos, ivec2( 0,  1), params, im_params) +
+		eval_cell(cpos, gpos, ivec2( 1, -1), params, im_params) +
+		eval_cell(cpos, gpos, ivec2( 1,  0), params, im_params) +
+		eval_cell(cpos, gpos, ivec2( 1,  1), params, im_params);
 	
-	//float lambda = params.a*params.a*params.density;
 	
 	//ad hoc attempt to normalize
-	//value/=log2(2.+2.*params.density);
-	value*=(1./PI)*pow(params.density+1., -.5);
-	//value*=(1./3.)/log2(params.density+2.);
-	float octexp = pow(2., params.octaves);
-	value*= (1.+params.octaves)*octexp/(2.*octexp-1.);
+	value*=.5*pow(params.density+1., -.5);
+	float octexp = pow(2., params.sector.y);
+	value*= (1.+params.sector.y)*octexp/(2.*octexp-1.);
 	
 	return value;
 }

			//position of impulse in cell space - uniform distribution
			vec2 ipos = vec2(nextRand(u), nextRand(u));
			//displacement to fragment
			vec2 delta = (fpos - ipos)*params.ainv;
			//impulse frequency, orientation - uniform distribution on input ranges
			float mfreq = pow(2., nextRand(u)*params.octaves);
			float ifreq = h.x*mfreq; 
			float iorientation = mix(h.z-.5*h.w, h.z+.5*h.w, nextRand(u));
			//evaluate kernel, accumulate fragment value
			vec2 mu = ifreq*vec2(cos(iorientation), sin(iorientation));
			float phi = nextRand(u); //phase - uniform dist [0, 1]
			float filt_exp = -.5*dot(mu, params.sigma_f_plus_g_inv*mu);
			acc+= filt_scale/mfreq*exp(-PI*aps*dot(delta,delta)+filt_exp)*cos(2.*PI*(dot(delta, params.filter*mu)+phi));
		}else {break;}
	}
	return acc;
 }
 
 float det2x2(mat2 m){
	return (m[0][0]*m[1][1] - m[0][1]*m[1][0]);
 }
 mat2 inv2x2(mat2 m){
	return (1./det2x2(m))*mat2(m[1][1], -m[0][1], -m[1][0], m[0][0]);
 }
 mat2 id2x2(){

//approximate poisson distribution uniform u
//from Galerne, Lagae, Lefebvre, Drettakis
const float poisson_epsilon = .001;
int poisson(inout float u, const float m){
	float u1 = nextRand(u);
	float u2 = nextRand(u);
	float x = sqrt(-2.*log(u1+poisson_epsilon))*cos(2.*PI*u2);
