params.cpp 25.2 KB
Newer Older
1 2 3
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions

   Copyright (C) 2014 CNRS
4
   Copyright (C) 2013, 2014, 2015, 2016 Inria
5 6 7 8 9 10 11

   This file is part of ALTA.

   This Source Code Form is subject to the terms of the Mozilla Public
   License, v. 2.0.  If a copy of the MPL was not distributed with this
   file, You can obtain one at http://mozilla.org/MPL/2.0/.  */

12
#include "params.h"
13
#include "common.h"
14
#include <cassert>
15

16 17
using namespace alta;

18 19
struct param_info
{
20
	param_info(std::string n, int d, std::string i) :
21
        name(n), dimension(d), info(i) { }
22 23 24 25 26 27

	std::string name;
	int dimension;
	std::string info;
};

28 29
std::map<params::input, const param_info> create_map()
{
30 31 32 33
	std::map<params::input, const param_info> _map;
	/* 1D Params */
	_map.insert(std::make_pair<params::input, const param_info>(params::COS_TH, param_info("COS_TH", 1, "Cosine of the Half angle")));
	_map.insert(std::make_pair<params::input, const param_info>(params::COS_TK, param_info("COS_TK", 1, "Cosine of the Back angle")));
34 35
	_map.insert(std::make_pair<params::input, const param_info>(params::COS_TLV, param_info("COS_TLV", 1, "Cosine of the Light and View directions")));
	_map.insert(std::make_pair<params::input, const param_info>(params::COS_TLR, param_info("COS_TLR", 1, "Cosine of the Light and Reflected directions")));
36

37 38
	/* 2D Params */
	_map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_TH_TD, param_info("RUSIN_TH_TD", 2, "Radialy symmetric Half angle parametrization")));
39
	_map.insert(std::make_pair<params::input, const param_info>(params::COS_TH_TD, param_info("COS_TH_TD", 2, "Cosines of the elevation angles of the Half angle parametrization")));
40
	_map.insert(std::make_pair<params::input, const param_info>(params::ISOTROPIC_TV_PROJ_DPHI, param_info("ISOTROPIC_TV_PROJ_DPHI", 2, "Isoptropic projected phi parametrization without a light direction.")));
41
	_map.insert(std::make_pair<params::input, const param_info>(params::STARK_2D, param_info("STARK_2D", 2, "Stark parametrization H, B but without third component.")));
42
	_map.insert(std::make_pair<params::input, const param_info>(params::NEUMANN_2D, param_info("NEUMANN_2D", 2, "Neumann parametrization H, B but without third component.")));
43

44 45 46
	/* 3D Params */
	_map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_TH_TD_PD, param_info("RUSIN_TH_TD_PD", 3, "Isotropic Half angle parametrization")));
	_map.insert(std::make_pair<params::input, const param_info>(params::ISOTROPIC_TV_TL_DPHI, param_info("ISOTROPIC_TV_TL_DPHI", 3, "Isotropic Light/View angle parametrization")));
47
	_map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_VH, param_info("RUSIN_VH", 3, "Vector representation of the Half angle only")));
48
	_map.insert(std::make_pair<params::input, const param_info>(params::SCHLICK_VK, param_info("SCHLICK_VH", 3, "Vector representation of the Back angle only")));
49
	_map.insert(std::make_pair<params::input, const param_info>(params::ISOTROPIC_TL_TV_PROJ_DPHI, param_info("ISOTROPIC_TL_TV_PROJ_DPHI", 3, "Isoptropic projected phi parametrization.")));
50
	_map.insert(std::make_pair<params::input, const param_info>(params::SCHLICK_TL_TK_PROJ_DPHI, param_info("SCHLICK_TL_TK_PROJ_DPHI", 3, "Isoptropic projected phi parametrization centered around the back vector.")));
51
	_map.insert(std::make_pair<params::input, const param_info>(params::RETRO_TL_TVL_PROJ_DPHI, param_info("RETRO_TL_TVL_PROJ_DPHI", 3, "Isoptropic retro projected phi parametrization.")));
52 53
	_map.insert(std::make_pair<params::input, const param_info>(params::STARK_3D, param_info("STARK_3D", 3, "Stark parametrization H, B.")));
	_map.insert(std::make_pair<params::input, const param_info>(params::NEUMANN_3D, param_info("NEUMANN_3D", 3, "Neumann parametrization H, B.")));
54

55 56 57
	/* 4D Params */
	_map.insert(std::make_pair<params::input, const param_info>(params::RUSIN_TH_PH_TD_PD, param_info("RUSIN_TH_PH_TD_PD", 4, "Complete Half angle parametrization")));
	_map.insert(std::make_pair<params::input, const param_info>(params::SPHERICAL_TL_PL_TV_PV, param_info("SPHERICAL_TL_PL_TV_PV", 4, "Complete classical parametrization")));
58
	_map.insert(std::make_pair<params::input, const param_info>(params::STEREOGRAPHIC, param_info("STEREOGRAPHIC", 4, "Light/View vector in stereographic projection"))),
59 60 61 62 63 64 65 66

	/* 6D Param */
	_map.insert(std::make_pair<params::input, const param_info>(params::CARTESIAN, param_info("CARTESIAN", 6, "Complete vector parametrization")));

	return _map;
}
static const std::map<params::input, const param_info> input_map = create_map();

67

68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
static std::map<params::output, std::string> create_output_map()
{
		std::map<params::output, std::string> result;

#define STRINGIFY_(x) #x
#define STRINGIFY(x)  STRINGIFY_(x)
#define DEFINE_MAPPING(name)										\
		result[params::name] = STRINGIFY(name);

		DEFINE_MAPPING(INV_STERADIAN);
		DEFINE_MAPPING(INV_STERADIAN_COSINE_FACTOR);
		DEFINE_MAPPING(ENERGY);
		DEFINE_MAPPING(RGB_COLOR);
		DEFINE_MAPPING(XYZ_COLOR);

#undef DEFINE_MAPPING
#undef STRINGIFY
#undef STRINGIFY_

		return result;
}

static const std::map<params::output, std::string> output_map =
		create_output_map();

93 94 95
void params::to_cartesian(const double* invec, params::input intype,
		double* outvec)
{
96 97 98 99
	#ifdef DEBUG_PARAM
	std::cout << " ENTERING to_cartesian with inttype = " << get_name(intype) << std::endl;
	#endif

100 101 102 103 104
	switch(intype)
	{
		// 1D Parametrizations
		case params::COS_TH:
			half_to_cartesian(acos(invec[0]), 0.0, 0.0, 0.0, outvec);
105 106 107 108 109 110 111 112
			break;
		case params::COS_TK:
			outvec[0] = invec[0]*invec[0]-1.0;
			outvec[1] = 0.0;
			outvec[2] = invec[0];
			outvec[3] = 1.0-invec[0]*invec[0];
			outvec[4] = 0.0;
			outvec[5] = invec[0];
113
			break;
114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
		case params::COS_TLV:
			outvec[0] = sqrt(1.0 - invec[0]*invec[0]);
			outvec[1] = 0.0;
			outvec[2] = invec[0];
			outvec[3] = 0.0;
			outvec[4] = 0.0;
			outvec[5] = 1.0;
			break;
		case params::COS_TLR:
			outvec[0] = - sqrt(1.0 - invec[0]*invec[0]);
			outvec[1] = 0.0;
			outvec[2] = invec[0];
			outvec[3] = 0.0;
			outvec[4] = 0.0;
			outvec[5] = 1.0;
			break;
130 131 132 133 134 135

			// 2D Parametrizations
		case params::COS_TH_TD:
			half_to_cartesian(acos(invec[0]), 0.0, acos(invec[1]), 0.0, outvec);
			break;
		case params::RUSIN_TH_TD:
136 137 138 139
         half_to_cartesian(invec[0], 0.0, invec[1], 0.5*M_PI, outvec);
			break;
		case ISOTROPIC_TV_PROJ_DPHI:
		{
Laurent Belcour's avatar
Laurent Belcour committed
140
			const double theta = sqrt(invec[0]*invec[0] + invec[1]*invec[1]);
141 142
			if(theta > 0.0)
			{
143 144
				outvec[0] = (invec[0]/theta)*sin(theta);
				outvec[1] = (invec[1]/theta)*sin(theta);
145 146 147
			}
			else
			{
148 149
				outvec[0] = 0.0;
				outvec[1] = 0.0;
150
			}
151 152 153 154
			outvec[2] = cos(theta);
			outvec[3] = 0.0;
			outvec[4] = 0.0;
			outvec[5] = 1.0;
155 156
		}
			break;
157 158 159
		// invec[0] = ||Hp|| Norm of the projected normalized Half vector
		// H=(L+V)/2 invec[1] = ||B||  Norm of the unormalized Back vector
		// B=(L-V)/2
160 161 162 163
		case STARK_2D:
		{
			const double Hx = invec[0];
			const double Hy = 0;
164 165 166
         const double sum = Hx*Hx + invec[1]*invec[1];
			const double Hz = (sum <= 1.0) ? sqrt(1.0 - std::min(sum, 1.0)) : -1.0;

167
			// Ensuring that <H,B> = 0.
168 169 170
			const double Bx = 0.0;
			const double By = invec[1];
			const double Bz = 0.0;
171

172
         //assert(Hz >= 0 && Hz <= 1);
173 174 175
			outvec[0] = Hx-Bx;
			outvec[1] = Hy-By;
			outvec[2] = Hz-Bz;
176
			outvec[3] = Hx+Bx;
177 178
			outvec[4] = Hy+By;
			outvec[5] = Hz+Bz;
179
		}
180
		break;
181 182 183 184 185 186 187 188
		case NEUMANN_2D:
		{
			outvec[0] =   invec[0];
			outvec[1] =   invec[1];
			outvec[2] =   sqrt(1.0 - outvec[0]*outvec[0] - outvec[1]*outvec[1]);
			outvec[3] =   invec[0];
			outvec[4] = - invec[1];
			outvec[5] =   sqrt(1.0 - outvec[3]*outvec[3] - outvec[4]*outvec[4]);
189
		}
190 191
		break;

192 193 194 195 196 197 198

			// 3D Parametrization
		case params::RUSIN_TH_PH_TD:
			half_to_cartesian(invec[0], invec[1], invec[2], 0.0, outvec);
			break;
		case params::RUSIN_TH_TD_PD:
			half_to_cartesian(invec[0], 0.0, invec[1], invec[2], outvec);
Laurent Belcour's avatar
Laurent Belcour committed
199 200 201
#ifdef DEBUG
			std::cout << outvec[2] << std::endl;
#endif
202 203 204 205
			break;
		case params::ISOTROPIC_TV_TL_DPHI:
			classical_to_cartesian(invec[0], 0.0, invec[1], invec[2], outvec);
			break;
206
		case params::RUSIN_VH:
207
            half_to_cartesian(acos(invec[2]), atan2(invec[1], invec[0]), 0.0, 0.0, outvec);
208
			break;
209
        // \todo I should handle the phi_k in the conversion to CARTESIAN
210
      case params::SCHLICK_VK:
211 212 213 214 215 216 217
            outvec[0] = invec[2]*invec[2]-1.0;
            outvec[1] = 0.0;
            outvec[2] = invec[2];
            outvec[3] = 1.0-invec[2]*invec[2];
            outvec[4] = 0.0;
            outvec[5] = invec[2];
            break;
218 219
		case ISOTROPIC_TL_TV_PROJ_DPHI:
		{
220
			const double theta = sqrt(invec[1]*invec[1] + invec[2]*invec[2]);
Laurent Belcour's avatar
Laurent Belcour committed
221 222 223 224 225
			// theta can equals zero. In that case, you should not
			// divide by its value, and the relative orientation does
			// not matter (normal direction).
			if(theta > 0.0)
			{
226 227
				outvec[0] = (invec[1]/theta)*sin(theta);
				outvec[1] = (invec[2]/theta)*sin(theta);
Laurent Belcour's avatar
Laurent Belcour committed
228 229 230 231 232 233
			}
			else
			{
				outvec[0] = 0.0;
				outvec[1] = 0.0;
			}
234 235 236 237
			outvec[2] = cos(theta);
			outvec[3] = sin(invec[0]);
			outvec[4] = 0.0;
			outvec[5] = cos(invec[0]);
238 239
		}
			break;
240
		case SCHLICK_TL_TK_PROJ_DPHI:
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
		{
			// Set the light direction
			outvec[3] = sin(invec[0]);
			outvec[4] = 0.0;
			outvec[5] = cos(invec[0]);

			// The view direction is the symmetric of the reflected direction
			// with respect to the back direction:
			//     v = 2 |r.k| k - r
			//     r = 2 |l.n| n - l = {-lx, -ly, lz }
			const double theta = sqrt(invec[1]*invec[1] + invec[2]*invec[2]);
			const double Kx = (theta > 0.0) ? (invec[1]/theta)*sin(theta) : 0.0;
			const double Ky = (theta > 0.0) ? (invec[2]/theta)*sin(theta) : 0.0;
			const double Kz = cos(theta);
			const double dotKR = outvec[5]*Kz - outvec[3]*Kx;
			outvec[0] = 2.0*dotKR*Kx + outvec[3];
			outvec[1] = 2.0*dotKR*Ky;
			outvec[2] = 2.0*dotKR*Kz - outvec[5];
		}
260
			break;
261
		case RETRO_TL_TVL_PROJ_DPHI:
262 263 264 265 266 267
		{
			const double theta_vl = sqrt(invec[1]*invec[1] + invec[2]*invec[2]);
			const double cos_vl   = cos(theta_vl);
			const double sin_vl   = sin(theta_vl);
			const double cos_l    = cos(invec[0]);
			const double sin_l    = sin(invec[0]);
268 269
			const double cos_phi  = (theta_vl > 0.0) ? invec[1] / theta_vl : 0.0;
			const double sin_phi  = (theta_vl > 0.0) ? invec[2] / theta_vl : 0.0;
270 271 272 273

			// Compute the cosine of the outgoing vector using the
			// spherical law of cosines
			const double cos_v = cos_l*cos_vl + cos_phi*sin_l*sin_vl;
274
			const double sin_v = sqrt(1.0 - std::min(cos_v*cos_v, 1.0));
275 276 277 278 279

			if(sin_v != 0.0)
			{
				const double sin_dphi = (sin_vl*sin_phi) / sin_v;

280
				outvec[0] = sin_v * sqrt(1.0 - std::min(sin_dphi*sin_dphi, 1.0));
281 282 283 284 285 286 287 288 289 290 291 292
				outvec[1] = sin_v * sin_dphi;
			}
			else
			{
				outvec[0] = 0.0;
				outvec[1] = 0.0;
			}
			outvec[2] = cos_v;
			outvec[3] = sin(invec[0]);
			outvec[4] = 0.0;
			outvec[5] = cos(invec[0]);
		}
293
			break;
294 295
		case STARK_3D:
		{
296
			// Constructing the Half vector
297 298
			const double Hx = invec[0];
			const double Hy = 0;
299
			const double Hz = sqrt(1.0 - invec[0]*invec[0] - invec[1]*invec[1]);
300

301
			// Constructing the Back vector using azimuth and elevation angles
302 303
			const double cosPhi = cos(invec[2]);
			const double sinPhi = sin(invec[2]);
304

305
#define STARK_BUILDING_THETA
306

307 308 309 310 311 312 313 314 315
#ifdef STARK_BUILDING_THETA
			const double Theta  = atan2(-Hz, invec[0]*cosPhi);
			const double cosThe = cos(Theta);
			const double sinThe = sin(Theta);
#else
			// Pascal Barla's way to do it
			const double cosThe = -cosPhi*invec[0] / sqrt(Hz*Hz + cosPhi*cosPhi*invec[0]*invec[0]);
			const double sinThe = Hz / (sqrt(Hz*Hz + cosPhi*cosPhi * invec[0]*invec[0]));
#endif
316 317 318
			const double Bx = invec[1]*sinThe*cosPhi;
			const double By = invec[1]*sinThe*sinPhi;
			const double Bz = invec[1]*cosThe;
319

320 321 322
			outvec[0] = Hx-Bx;
			outvec[1] = Hy-By;
			outvec[2] = Hz-Bz;
323
			outvec[3] = Hx+Bx;
324 325 326
			outvec[4] = Hy+By;
			outvec[5] = Hz+Bz;
		}
327
			break;
328 329 330 331 332
		case NEUMANN_3D:
		{
			const double cosPhi = cos(invec[2]);
			const double sinPhi = sin(invec[2]);

333
			// Build the projected Half and Back vectors
334 335
			outvec[0] =   invec[0] + cosPhi*invec[1];
			outvec[1] =   sinPhi*invec[1];
336
			outvec[3] =   invec[0] - cosPhi*invec[1];
337
			outvec[4] = - sinPhi*invec[1];
338

339 340 341 342 343 344 345 346 347 348 349 350 351 352 353
			// Safeguard, if the vectors are not under unit length return an
			// invalid configuration
			if(outvec[0]*outvec[0]+outvec[1]*outvec[1] > 1.0 ||
				outvec[3]*outvec[3]+outvec[4]*outvec[4] > 1.0) {
				outvec[0] =  0.0;
				outvec[1] =  0.0;
				outvec[2] = -1.0;
				outvec[3] =  0.0;
				outvec[4] =  0.0;
				outvec[5] = -1.0;
				break;
			}

			// Project the vectorx on the hemisphere.
			outvec[2] =   sqrt(1.0 - outvec[0]*outvec[0] - outvec[1]*outvec[1]);
354 355
			outvec[5] =   sqrt(1.0 - outvec[3]*outvec[3] - outvec[4]*outvec[4]);
		}
356
			break;
357 358 359 360 361 362 363

			// 4D Parametrization
		case params::RUSIN_TH_PH_TD_PD:
			half_to_cartesian(invec[0], invec[1], invec[2], invec[3], outvec);
			break;

		case params::SPHERICAL_TL_PL_TV_PV:
364 365 366 367 368 369
			outvec[0] = cos(invec[3])*sin(invec[2]);
			outvec[1] = sin(invec[3])*sin(invec[2]);
			outvec[2] = cos(invec[2]);
			outvec[3] = cos(invec[1])*sin(invec[0]);
			outvec[4] = sin(invec[1])*sin(invec[0]);
			outvec[5] = cos(invec[0]);
370 371
			break;

372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390
		case params::STEREOGRAPHIC:
      {
			// Project the 2D direction on the surface invec[0,1] to the View
         // vector outvec[0,1,2]
         const double normL = invec[0]*invec[0] + invec[1]*invec[1];
         const double alphL = (1.0 - normL) / (1.0 + normL);
         outvec[0] = invec[0] + alphL;
         outvec[1] = invec[1] + alphL;
         outvec[2] = alphL;

         // Project the 2D direction on the surface invec[2,3] to the Light
         // vector outvec[3,4,5]
         const double normV = invec[2]*invec[2] + invec[3]*invec[3];
         const double alphV = (1.0 - normV) / (1.0 + normV);
         outvec[3] = invec[2] + alphV;
         outvec[4] = invec[3] + alphV;
         outvec[5] = alphV;
		}
         break;
391

392 393
		// 6D Parametrization
      case params::CARTESIAN:
394 395 396 397
			memcpy(outvec, invec, 6*sizeof(double));
			break;

		default:
398 399
			std::cerr << "<<ERROR>> Transformation not implemented, " << get_name(intype) << " " << __FILE__ << ":" << __LINE__ << std::endl;
			throw;
400 401 402 403 404
			break;
	}

}

405
// Return in HALF the half vector for the two vectors in INVEC.
406 407
template<class Vector, typename Real>
static void half_vector(Real const *invec, Vector &half)
408 409 410 411
{
    half[0] = invec[0] + invec[3];
    half[1] = invec[1] + invec[4];
    half[2] = invec[2] + invec[5];
412
    Real sqnorm = half[0]*half[0] + half[1]*half[1] + half[2]*half[2];
413 414 415 416 417 418 419 420 421 422 423 424 425 426

    if (sqnorm <= 0.) {
        half[0] = 0;
        half[1] = 0;
        half[2] = 1;
    }
    else {
        half /= sqrt(sqnorm);
    }

    assert(half[2] <= 1.);
    assert(half[2] >= 0.);
}

427 428 429
// Return true if X is essentially zero.  We have to resort to this hack
// before using functions that have a high derivative around zero, such as
// 'atan2'.
430 431
template<typename Real>
static bool is_zero(Real x)
432
{
433
    return close_to(x, 0., 1e-14);
434 435
}

436 437 438 439
// Return true if VEC is alongside z⃗.
template<class Vector>
static bool alongside_z(const Vector& vec)
{
440 441
    return (is_zero(vec[0]) && is_zero(vec[1]))
        || close_to(vec[2], 1., 1e-14);
442 443 444
}


445 446 447
void params::from_cartesian(const double* invec, params::input outtype,
		double* outvec)
{
448 449 450 451
	#ifdef DEBUG_PARAM
	std::cout << " ENTERING from_cartesian with outtype = " << get_name(outtype) << std::endl;
	std::cout << " invec = " << invec[0] <<  " " << invec[1] << " " << invec[2] << " " << invec[3]
						<< "  " << invec[4] << " " << invec[5] << std::endl;
452
	#endif
453

454
	// Compute the half vector.
455
  Eigen::Vector3d half;
456
  half_vector(invec, half);
457

458
	// Difference vector.
459
  Eigen::Vector3d diff;
460 461 462 463 464 465 466

	switch(outtype)
	{
		// 1D Parametrizations
		case params::COS_TH:
			outvec[0] = half[2];
			break;
467
		case params::COS_TK:
468 469 470 471 472 473
		{
			const double Kx = invec[0]-invec[3];
			const double Ky = invec[1]-invec[4];
			const double Kz = invec[2]+invec[5];
			outvec[0] = (invec[2] + invec[5]) / sqrt(Kx*Kx + Ky*Ky + Kz*Kz);
		}
474
			break;
475 476 477 478 479 480
		case params::COS_TLV:
			outvec[0] = invec[0]*invec[3] + invec[1]*invec[4] + invec[2]*invec[5];
			break;
		case params::COS_TLR:
			outvec[0] = invec[0]*invec[3] - (invec[1]*invec[4] + invec[2]*invec[5]);
			break;
481 482 483 484

			// 2D Parametrizations
		case params::COS_TH_TD:
			outvec[0] = half[2];
485
			outvec[1] = half[0]*invec[0] + half[1]*invec[1] + half[2]*invec[2];
486 487 488
			break;
		case params::RUSIN_TH_TD:
			outvec[0] = acos(half[2]);
489
			outvec[1] = acos(clamp(half[0]*invec[0] + half[1]*invec[1] + half[2]*invec[2], -1.0, 1.0));
490
			break;
491 492
		case ISOTROPIC_TV_PROJ_DPHI:
		{
493
			const double theta = acos(invec[2]);
494 495 496
			const double dphi  = atan2(invec[1], invec[0]) - atan2(invec[4], invec[3]);
			outvec[0] = theta * cos(dphi);
			outvec[1] = theta * sin(dphi);
497 498
		}
			break;
499
		// outvec[0] = ||Hp|| Norm of the projected unormalized Half vector
500
		//             H = (V+L)/2
501
		// outvec[1] = ||B|| Norm of the unormalized Back vector B = (L-V)/2
502 503 504 505 506 507 508 509 510 511
		case STARK_2D:
		{
			double Hx = 0.5*(invec[0]+invec[3]);
			double Hy = 0.5*(invec[1]+invec[4]);
			outvec[0] = sqrt(Hx*Hx + Hy*Hy);

			double Bx = 0.5*(invec[3]-invec[0]);
			double By = 0.5*(invec[4]-invec[1]);
			double Bz = 0.5*(invec[5]-invec[2]);
			outvec[1] = sqrt(Bx*Bx + By*By + Bz*Bz);
512 513 514 515
		}
			break;
		case NEUMANN_2D:
		{
516 517 518
				double Hx = 0.5*(invec[0]+invec[3]);
				double Hy = 0.5*(invec[1]+invec[4]);
				outvec[0] = sqrt(Hx*Hx + Hy*Hy);
519

520 521 522
				double Bx = 0.5*(invec[3]-invec[0]);
				double By = 0.5*(invec[4]-invec[1]);
				outvec[1] = sqrt(Bx*Bx + By*By);
523 524
		}
			break;
525

526

527 528 529
			// 3D Parametrization
		case params::RUSIN_TH_PH_TD:
			outvec[0] = acos(half[2]);
530
			outvec[1] = atan2(half[1], half[0]);
531
			outvec[2] = acos(half[0]*invec[0] + half[1]*invec[1] + half[2]*invec[2]);
532
			break;
533 534
		case params::RUSIN_TH_TD_PD:
			outvec[0] = acos(half[2]);
535

536 537 538 539
			// Compute the diff vector
			diff[0] = invec[0];
			diff[1] = invec[1];
			diff[2] = invec[2];
540

541
      rotate(diff, -atan2(half[1], half[0]), -outvec[0]);
542

543 544 545 546 547 548 549 550 551 552
      // Note: the following approach is numerically more stable (up to
      // 1e-12) but slower.
      // {
      //     vec view(3);
      //     view[0] = invec[0];
      //     view[1] = invec[1];
      //     view[2] = invec[2];
      //     outvec[1] = acos(clamp(half.dot(view),0.,1.));
      // }
      outvec[1] = acos(diff[2]);
553 554

      // By convention, when DIFF is alongside z⃗, return φd = 0.
555
      if (alongside_z(diff))
556 557 558
          outvec[2] = 0;
      else
          outvec[2] = atan2(diff[1], diff[0]);
559

560
			break;
561 562
		case params::ISOTROPIC_TV_TL_DPHI:
			outvec[0] = acos(invec[2]);
563
			outvec[1] = acos(invec[5]);
564
			outvec[2] = atan2(invec[4], invec[3]) - atan2(invec[1], invec[0]);
565
			break;
566
		case params::RUSIN_VH:
567 568 569
			outvec[0] = half[0];
			outvec[1] = half[1];
			outvec[2] = half[2];
570
			break;
571 572
      case params::SCHLICK_VK:
      {
573 574 575
			const double Kx = invec[3]-invec[0];
         const double Ky = invec[4]-invec[1];
         const double Kz = invec[5]+invec[2];
576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591

         const double norm =  sqrt(Kx*Kx + Ky*Ky + Kz*Kz);
			if(norm > 1.0E-10)
			{
				outvec[0] = Kx / norm;
		      outvec[1] = Ky / norm;
			   outvec[2] = Kz / norm;
			}
			else
			{
				outvec[0] = 0.0;
		      outvec[1] = 0.0;
			   outvec[2] = 1.0;
			}
      }
			break;
592 593
		case ISOTROPIC_TL_TV_PROJ_DPHI:
		{
594 595 596
			const double theta_l = acos(invec[5]);
			const double theta_v = acos(invec[2]);
			const double dphi    = atan2(invec[4], invec[3]) - atan2(invec[1], invec[0]);
597 598 599
			outvec[0] = theta_l;
			outvec[1] = theta_v * cos(dphi);
			outvec[2] = theta_v * sin(dphi);
600 601
		}
			break;
602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629
		case RETRO_TL_TVL_PROJ_DPHI:
		{
			const double cos_v  = invec[2];
			const double cos_l  = invec[5];
			const double cos_vl = invec[0]*invec[3] + invec[1]*invec[4] + invec[2]*invec[5];
			const double sin_l  = sqrt(1.0 - cos_l * cos_l);
			const double sin_v  = sqrt(1.0 - cos_v * cos_v);
			const double sin_vl = sqrt(1.0 - cos_vl* cos_vl);

			// Sine of the \Delta_{\phi} angle
			const double sin_dp = sin(atan2(invec[4], invec[3]) - atan2(invec[1], invec[0]));

			if(sin_vl != 0.0 && sin_l != 0.0)
			{
				// Cosine and sine of the phi_{L,V} angle: the spherical
				// angle between the arclength {L,N} and {V,N}
				const double cos_phi = (cos_v - cos_l*cos_vl) / (sin_l * sin_vl);
				const double sin_phi = (sin_dp * sin_v) / sin_vl;

				const double theta = acos(cos_vl);

				outvec[0] = acos(cos_l);
				outvec[1] = theta * cos_phi;
				outvec[2] = theta * sin_phi;
			}
			else
			{
				outvec[0] = acos(cos_l);
630
				outvec[1] = acos(cos_v) - acos(cos_l);
631 632 633 634
				outvec[2] = 0.0;
			}
		}
			break;
635 636 637 638 639 640 641 642 643 644 645
		case SCHLICK_TL_TK_PROJ_DPHI:
		{
			const double vkx     = invec[0]-invec[3];
			const double vky     = invec[1]-invec[4];
			const double vkz     = invec[2]+invec[5];
			const double norm    = sqrt(vkx*vkx + vky*vky + vkz*vkz);
			const double theta_k = acos(vkz / norm);
			const double dphi    = atan2(invec[4], invec[3]) - atan2(vky/norm, vkx/norm);
			outvec[0] = acos(invec[5]);
			outvec[1] = theta_k * cos(dphi);
			outvec[2] = theta_k * sin(dphi);
646 647 648 649 650 651 652 653 654 655 656 657
		}
			break;
		case STARK_3D:
		{
			double Hx = 0.5*(invec[0]+invec[3]);
			double Hy = 0.5*(invec[1]+invec[4]);
			outvec[0] = sqrt(Hx*Hx + Hy*Hy);

			double Bx = 0.5*(invec[3]-invec[0]);
			double By = 0.5*(invec[4]-invec[1]);
			double Bz = 0.5*(invec[5]-invec[2]);
			outvec[1] = sqrt(Bx*Bx + By*By + Bz*Bz);
658
			outvec[2] = atan2(By, Bx) -  atan2(Hy, Hx);
659 660 661
         if(outvec[2] < 0.0) {
            outvec[2] = 2.0*M_PI + outvec[2];
         }
662 663 664 665 666 667 668 669 670 671 672
		}
			break;
		case NEUMANN_3D:
		{
				double Hx = 0.5*(invec[0]+invec[3]);
				double Hy = 0.5*(invec[1]+invec[4]);
				outvec[0] = sqrt(Hx*Hx + Hy*Hy);

				double Bx = 0.5*(invec[3]-invec[0]);
				double By = 0.5*(invec[4]-invec[1]);
				outvec[1] = sqrt(Bx*Bx + By*By);
673
				outvec[2] = atan2(By, Bx) - atan2(Hy, Hx);
674 675
		}
			break;
676 677 678 679

			// 4D Parametrization
		case params::RUSIN_TH_PH_TD_PD:
			outvec[0] = acos(half[2]);
680
			outvec[1] = atan2(half[1], half[0]);
681 682 683 684 685

			// Compute the diff vector
			diff[0] = invec[0];
			diff[1] = invec[1];
			diff[2] = invec[2];
686
      rotate(diff, -outvec[1], -outvec[0]);
687

688
      // By convention, when DIFF is alongside z⃗, return φd = θd = 0.
689
      if (alongside_z(diff))
690 691 692 693
      {
          outvec[2] = 0.;
          outvec[3] = 0.;
      }
694
      else
695 696
      {
          outvec[2] = acos(diff[2]);
697
          outvec[3] = atan2(diff[1], diff[0]);
698
      }
699 700 701
			break;

		case params::SPHERICAL_TL_PL_TV_PV:
702 703 704 705
			outvec[0] = acos(invec[5]);
			outvec[1] = atan2(invec[4], invec[3]);
			outvec[2] = acos(invec[2]);
			outvec[3] = atan2(invec[1], invec[0]);
Laurent Belcour's avatar
Laurent Belcour committed
706 707 708
#ifdef DEBUG
			std::cout << invec[2] << " - acos -> " << outvec[0] << std::endl;
#endif
709 710
			break;

711 712 713 714 715 716 717
      case params::STEREOGRAPHIC:
      {
			// Project the View vector invec[0,1,2] on a 2D direction on the
         // surface outvec[0,1]
         const double dotVN = invec[2];
         outvec[0] = invec[0] / (1.0+dotVN);
         outvec[1] = invec[1] / (1.0+dotVN);
718

719 720 721 722 723
         // Project the Light vector invec[0,1,2] on a 2D direction on the
         // surface outvec[2,3]
         const double dotLN = invec[5];
         outvec[2] = invec[3] / (1.0+dotLN);
         outvec[3] = invec[4] / (1.0+dotLN);
724

725 726
		}
         break;
727

728 729 730 731 732 733
			// 6D Parametrization
		case params::CARTESIAN:
			memcpy(outvec, invec, 6*sizeof(double));
			break;

		default:
734
			std::cerr << "<<ERROR>> Transformation not implemented, n°" << outtype << ", " << __FILE__ << ":" << __LINE__ << std::endl;
735 736 737 738
			assert(false);
			break;
	}
}
739 740
params::input params::parse_input(const std::string& txt)
{
741

742 743 744 745
	for(std::map<params::input, const param_info>::const_iterator it=input_map.begin(); it != input_map.end(); ++it)
	{
		if(txt.compare(it->second.name) == 0)
		{
746
			std::cout << "<<INFO>> parsed input parametrization " << it->second.name << " from name \"" << txt << "\"" << std::endl;
747 748 749 750
			return it->first;
		}
	}

751
	std::cout << "<<INFO>> the input parametrization is UNKNOWN_INPUT" << std::endl;
752 753
	return params::UNKNOWN_INPUT;
}
754

755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781
params::output params::parse_output(const std::string& txt)
{
	if(txt == std::string("ENERGY"))
	{
		return params::ENERGY;
	}
	else if(txt == std::string("INV_STERADIAN"))
	{
		return params::INV_STERADIAN;
	}
	else if(txt == std::string("INV_STERADIAN_COSINE_FACTOR"))
	{
		return params::INV_STERADIAN_COSINE_FACTOR;
	}
	else if(txt == std::string("RGB_COLOR"))
	{
		return params::RGB_COLOR;
	}
	else if(txt == std::string("XYZ_COLOR"))
	{
		return params::XYZ_COLOR;
	}
	else
	{
		return params::UNKNOWN_OUTPUT;
	}
}
782 783

const std::string& params::get_name(const params::input param)
784 785 786 787 788 789 790
{
	std::map<params::input, const param_info>::const_iterator it = input_map.find(param);
	if(it != input_map.end())
	{
		return it->second.name;
	}

791 792 793
#ifdef DEBUG
	std::cerr << "<<WARNING>> Unknown parametrization, n°" << param << ", "<< __FILE__ << ":" << __LINE__ << std::endl;
#endif
794 795 796

	static const std::string unknown = "UNKNOWN_INPUT";
	return unknown;
797 798
}

799
const std::string& params::get_name(const params::output param)
800 801 802 803 804 805 806 807 808
{
		std::map<params::output, std::string>::const_iterator it = output_map.find(param);
		if (it != output_map.end())
		{
				return it->second;
		}

		static const std::string unknown = "UNKNOWN_OUTPUT";
		return unknown;
809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825
}

int  params::dimension(params::input t)
{
	std::map<params::input, const param_info>::const_iterator it = input_map.find(t);
	if(it != input_map.end())
	{
		return it->second.dimension;
	}
	else
	{
		return -1;
	}
}

void params::print_input_params()
{
826
	std::cout << "List of available input parametrizations in the ALTA library:" << std::endl;
827 828
	for(std::map<params::input, const param_info>::const_iterator it=input_map.begin(); it != input_map.end(); ++it)
	{
829 830 831
		std::cout << "  - " << it->second.name;
		for(int i=it->second.name.size(); i<26; ++i) { std::cout << " "; }
		std::cout << it->second.info << std::endl;
832 833
	}
}