 	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));
+	im_params.filt_scale = im_params.a_prime_square*im_params.ainv*im_params.ainv;
 
 	float value = 
-		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);
+		eval_cell(cpos, gpos, vec2(-1., -1.), params, im_params) +
+		eval_cell(cpos, gpos, vec2(-1.,  0.), params, im_params) +
+		eval_cell(cpos, gpos, vec2(-1.,  1.), params, im_params) +
+		eval_cell(cpos, gpos, vec2( 0., -1.), params, im_params) +
+		eval_cell(cpos, gpos, vec2( 0.,  0.), params, im_params) +
+		eval_cell(cpos, gpos, vec2( 0.,  1.), params, im_params) +
+		eval_cell(cpos, gpos, vec2( 1., -1.), params, im_params) +
+		eval_cell(cpos, gpos, vec2( 1.,  0.), params, im_params) +
+		eval_cell(cpos, gpos, vec2( 1.,  1.), params, im_params);
 	
 	
 	//ad hoc attempt to normalize
 	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;
+
 	
-	return value;
 }
  


struct gnoise_im_params{ //struct to pass intermediate values within gnoise
	mat2 filter, sigma_f_plus_g_inv;
	float ainv, a_prime_square;
};

ivec2 bound_grid(const ivec2 gpos){
	return gpos + IMOD*(1 - gpos/IMOD);
}

//hash based on Blum, Blum & Shub 1986
//and Sharpe http://briansharpe.wordpress.com/2011/10/01/gpu-texture-free-noise/
const float bbsm = 137023.;//magic product of two primes chosen to have high period without float precision issues
vec2 bbsmod( const vec2 a ) {
	return a - floor( a * ( 1.0 / bbsm ) ) * bbsm;
}
vec2 bbs(const vec2 a) {
	return bbsmod(a*a);
}
float bbsopt( const float a ) {
	return fract( a * a * ( 1.0 / bbsm ) ) * bbsm;
}
float seed(const ivec2 p){
	vec2 h = bbs(vec2(p.xy));
	return bbsopt(h.x+bbsopt(h.y))*(1./bbsm);
}

//permutation polynomial
//based on Gustavson/McEwan https://github.com/ashima/webgl-noise/
//and Sharpe http://briansharpe.wordpress.com/2011/10/01/gpu-texture-free-noise/
const float pp_epsilon = .01;
float nextRand(inout float u){//rng
	u = fract(((u*34.0*289.)+1.0)*u+pp_epsilon);
	return fract(7.*u);
}

//approximate poisson distribution uniform u

	//compute positions for this fragment
	vec2 temp = pos*params.a; 
	vec2 cpos = fract(temp);
	vec2 gpos = floor(temp);
		
	mat2 jacob = params.jacob;
	mat2 jacob_t = mat2(jacob[0][0], jacob[1][0], jacob[0][1], jacob[1][1]);
