mitsuba-af602c6fd98a/src/bsdfs/irawan.cpp
/** parameters:
	 *	u	 for staple, we compute v(u)
	 *	v	 to be compared to v(u) in texturing
	 *	om_i  incident direction
	 *	om_r  exitant direction
	 *  fiber properties
	 *	alpha uniform scattering
	 *	beta  forward scattering
	 *  yarn geometry
	 *	psi   fiber twist angle
	 *	umax  maximum inclination angle
	 *	kappa spine curvature parameter
	 *  weave pattern
	 *	w	 width of segment rectangle
	 *	l	 length of segment rectangle
	 */
u        for staple, we compute v(u)
v        to be compared to v(u) in texturing
om_i  incident direction
om_r  exitant direction
fiber properties
alpha uniform scattering
beta  forward scattering
yarn geometry
psi   fiber twist angle
umax  maximum inclination angle
kappa spine curvature parameter
weave pattern
w        width of segment rectangle
l        length of segment rectangle
u for staple we compute v ( u ) v to be compared to v ( u ) in texturing om_i incident direction om_r exitant direction fiber properties alpha uniform scattering beta forward scattering yarn geometry psi fiber twist angle umax maximum inclination angle kappa spine curvature parameter weave pattern w width of segment rectangle l length of segment rectangle
Float evalStapleIntegrand(Float u, Float v, const Vector &om_i,
			const Vector &om_r, Float alpha, Float beta, Float psi,
			Float umax, Float kappa, Float w, Float l) const {
		if (w * std::sin(umax) >= l)
			return 0.0f;
		if (kappa < -1.0f)
			return 0.0f;
		Vector h = normalize(om_i + om_r);
		Float D = (h.y*std::cos(u) - h.z*std::sin(u))
			/ (std::sqrt(h.x * h.x + std::pow(h.y * std::sin(u) + h.z * std::cos(u), (Float) 2.0f)) * std::tan(psi));
		Float v_of_u = std::atan2(-h.y * std::sin(u) - h.z * std::cos(u), h.x) + math::safe_acos(D);
		if (std::abs(D) < 1.0f && std::abs(v_of_u) < M_PI / 2.0f) {
			Normal n = normalize(Normal(std::sin(v_of_u), std::sin(u)
					* std::cos(v_of_u), std::cos(u) * std::cos(v_of_u)));
			Float R = radiusOfCurvature(std::abs(u), umax, kappa, w, l);
			Float a = 0.5f * w;
			Vector om_i_plus_om_r(om_i + om_r);
			Float Gv = a * (R + a * std::cos(v_of_u))
				/ (om_i_plus_om_r.length() * dot(n, h) * std::abs(std::sin(psi)));
			Float fc = alpha + vonMises(-dot(om_i, om_r), beta);
			Float A = seeliger(dot(n, om_i), dot(n, om_r), 0, 1);
			Float fs = Gv * fc * A;
			fs = fs * 2.0f * w * umax;
			Float delta_x = w * m_pattern.hWidth;
			Float x_of_u = v_of_u * w / M_PI;
			if (x_of_u > 0.5f * (w - delta_x))
				x_of_u = 0.5f * (w - delta_x);
			else if (x_of_u < 0.5f * (delta_x - w))
				x_of_u = 0.5f * (delta_x - w);
			if (std::abs(x_of_u - v * w / M_PI) < 0.5f * delta_x)
				return fs / delta_x;
		}
		return 0.0f;
	}
Float evalStapleIntegrand ( Float u , Float v , const Vector & om_i , const Vector & om_r , Float alpha , Float beta , Float psi , Float umax , Float kappa , Float w , Float l ) const { if ( w * std :: sin ( umax ) >= l ) return 0.0f ; if ( kappa < -1.0f ) return 0.0f ; Vector h = normalize ( om_i + om_r ) ; Float D = ( h . y * std :: cos ( u ) - h . z * std :: sin ( u ) ) / ( std :: sqrt ( h . x * h . x + std :: pow ( h . y * std :: sin ( u ) + h . z * std :: cos ( u ) , ( Float ) 2.0f ) ) * std :: tan ( psi ) ) ; Float v_of_u = std :: atan2 ( - h . y * std :: sin ( u ) - h . z * std :: cos ( u ) , h . x ) + math :: safe_acos ( D ) ; if ( std :: abs ( D ) < 1.0f && std :: abs ( v_of_u ) < M_PI / 2.0f ) { Normal n = normalize ( Normal ( std :: sin ( v_of_u ) , std :: sin ( u ) * std :: cos ( v_of_u ) , std :: cos ( u ) * std :: cos ( v_of_u ) ) ) ; Float R = radiusOfCurvature ( std :: abs ( u ) , umax , kappa , w , l ) ; Float a = 0.5f * w ; Vector om_i_plus_om_r ( om_i + om_r ) ; Float Gv = a * ( R + a * std :: cos ( v_of_u ) ) / ( om_i_plus_om_r . length ( ) * dot ( n , h ) * std :: abs ( std :: sin ( psi ) ) ) ; Float fc = alpha + vonMises ( - dot ( om_i , om_r ) , beta ) ; Float A = seeliger ( dot ( n , om_i ) , dot ( n , om_r ) , 0 , 1 ) ; Float fs = Gv * fc * A ; fs = fs * 2.0f * w * umax ; Float delta_x = w * m_pattern . hWidth ; Float x_of_u = v_of_u * w / M_PI ; if ( x_of_u > 0.5f * ( w - delta_x ) ) x_of_u = 0.5f * ( w - delta_x ) ; else if ( x_of_u < 0.5f * ( delta_x - w ) ) x_of_u = 0.5f * ( delta_x - w ) ; if ( std :: abs ( x_of_u - v * w / M_PI ) < 0.5f * delta_x ) return fs / delta_x ; } return 0.0f ; }
parameters:
      u        for staple, we compute v(u)
      v        to be compared to v(u) in texturing
      om_i  incident direction
      om_r  exitant direction
  fiber properties
      alpha uniform scattering
      beta  forward scattering
  yarn geometry
      psi   fiber twist angle
      umax  maximum inclination angle
      kappa spine curvature parameter
  weave pattern
      w        width of segment rectangle
      l        length of segment rectangle
parameters : u for staple we compute v ( u ) v to be compared to v ( u ) in texturing om_i incident direction om_r exitant direction fiber properties alpha uniform scattering beta forward scattering yarn geometry psi fiber twist angle umax maximum inclination angle kappa spine curvature parameter weave pattern w width of segment rectangle l length of segment rectangle
// w * sin(umax) < l // -1 < kappa < inf // h is the half vector // v_of_u is location of specular reflection. // Check if v_of_u within the range of valid v values // n is normal to the yarn surface. // t is tangent of the fibers. /*Vector t = normalize(Vector(-std::cos(v_of_u) * std::sin(psi),
					std::cos(u) * std::cos(psi) + std::sin(u) * std::sin(v_of_u) * std::sin(psi),
					-std::sin(u) * std::cos(psi) + std::cos(u) * std::sin(v_of_u) * std::sin(psi))); */ // R is radius of curvature. // G is geometry factor. // fc is phase function. // A is attenuation function without smoothing. // fs is scattering function. // Domain transform. // Highlight has constant width delta_x on screen. // Clamp x_of_u between (w - delta_x)/2 and -(w - delta_x)/2. // Check if |x(v(u)) - x(v)| < delta_x/2.