From f33018188bf5bdf6f6b8af0d646e3f8c519d9d71 Mon Sep 17 00:00:00 2001 From: mwturvey Date: Thu, 23 Mar 2017 12:13:44 -0700 Subject: Added support for empty string in config.json & other cleanup. --- redist/json_helpers.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'redist') diff --git a/redist/json_helpers.c b/redist/json_helpers.c index 0267932..3b5cc0d 100644 --- a/redist/json_helpers.c +++ b/redist/json_helpers.c @@ -117,7 +117,8 @@ char* load_file_to_mem(const char* path) { fseek( f, 0, SEEK_END ); int len = ftell( f ); fseek( f, 0, SEEK_SET ); - char * JSON_STRING = malloc( len ); + char * JSON_STRING = malloc( len + 1); + memset(JSON_STRING,0,len+1); fread( JSON_STRING, len, 1, f ); fclose( f ); return JSON_STRING; -- cgit v1.3.1 From 5404526ae8da8c5fdff81b8ee8120ffe73647747 Mon Sep 17 00:00:00 2001 From: cnlohr Date: Fri, 24 Mar 2017 01:05:07 -0400 Subject: Dave's affine solve is getting close. --- dave/AffineSolve | Bin 40104 -> 40104 bytes redist/linmath.c | 2 +- src/poser_daveortho.c | 68 +++++++++++++++++++++++++++++++++++++------------- 3 files changed, 51 insertions(+), 19 deletions(-) (limited to 'redist') diff --git a/dave/AffineSolve b/dave/AffineSolve index 7d51b34..bd93cbd 100755 Binary files a/dave/AffineSolve and b/dave/AffineSolve differ diff --git a/redist/linmath.c b/redist/linmath.c index eefcd5f..1724a13 100644 --- a/redist/linmath.c +++ b/redist/linmath.c @@ -215,7 +215,7 @@ void quatfrommatrix( FLT * q, const FLT * matrix44 ) q[1] = (matrix44[9] - matrix44[6]) / S; q[2] = (matrix44[2] - matrix44[8]) / S; q[3] = (matrix44[4] - matrix44[1]) / S; - } else if ((matrix44[0] > matrix44[5])&(matrix44[0] > matrix44[10])) { + } else if ((matrix44[0] > matrix44[5])&&(matrix44[0] > matrix44[10])) { FLT S = FLT_SQRT(1.0 + matrix44[0] - matrix44[5] - matrix44[10]) * 2.; // S=4*qx q[0] = (matrix44[9] - matrix44[6]) / S; q[1] = 0.25f * S; diff --git a/src/poser_daveortho.c b/src/poser_daveortho.c index 80f65a9..906b21e 100644 --- a/src/poser_daveortho.c +++ b/src/poser_daveortho.c @@ -456,41 +456,73 @@ PRINT(ab,2,1); //------------------- // Orthogonalize the matrix //------------------- - - PRINT_MAT(T,4,4); -matrix44transpose(T2, T); -//matrix44copy(T2,T); - cross3d( &T2[0][0], &T2[1][0], &T2[2][0] ); - cross3d( &T2[2][0], &T2[0][0], &T2[1][0] ); - normalize3d( &T2[0][0], &T2[0][0] ); - normalize3d( &T2[1][0], &T2[1][0] ); - normalize3d( &T2[2][0], &T2[2][0] ); -//matrix44copy(T2,T); -matrix44transpose(T, T2); PRINT_MAT(T,4,4); +#if 1 +// matrix44transpose(T2, T); //Transpose so we are + matrix44copy(T2,T); + cross3d( &T2[1][0], &T2[0][0], &T2[2][0] ); + cross3d( &T2[2][0], &T2[1][0], &T2[0][0] ); //Replace axes in-place. + matrix44copy(T,T2); +// matrix44transpose(T, T2); + +#endif + + normalize3d( &T[0][0], &T[0][0] ); + normalize3d( &T[1][0], &T[1][0] ); + normalize3d( &T[2][0], &T[2][0] ); + //Change handedness + + T[1][0]*=-1; + T[1][1]*=-1; + T[1][2]*=-1; + +/* + //Check Orthogonality. Yep. It's orthogonal. + FLT tmp[3]; + cross3d( tmp, &T[0][0], &T[1][0] ); + printf( "M3: %f\n", magnitude3d( tmp ) ); + cross3d( tmp, &T[2][0], &T[1][0] ); + printf( "M3: %f\n", magnitude3d( tmp ) ); + cross3d( tmp, &T[2][0], &T[0][0] ); + printf( "M3: %f\n", magnitude3d( tmp ) ); +*/ + +// PRINT_MAT(T,4,4); +#if 1 - //memcpy(T2,T,16*sizeof(float)); -// matrix44copy(T2,T); + //matrix44copy(T2,T); matrix44transpose(T2,T); + quatfrommatrix( quat, &T2[0][0] ); + printf( "QM: %f\n", quatmagnitude( quat ) ); quatnormalize(quatNorm,quat); quattoeuler(euler,quatNorm); - quattomatrix( T2, quatNorm ); + quattomatrix( &T2[0][0], quatNorm ); + + PRINT_MAT(T2,4,4); printf("rot %f %f %f len %f\n", euler[0], euler[1], euler[2], quatmagnitude(quat)); - PRINT(T,4,4); +// PRINT(T,4,4); // matrix44copy(temp,T2); matrix44transpose(temp,T2); - memcpy(T2,temp,16*sizeof(float)); + +// matrix44transpose(T2, temp); +// memcpy(T2,temp,16*sizeof(float)); for (i=0; i<3; i++) { for (j=0; j<3; j++) { - T[i][j] = T2[j][i]; + T[i][j] = temp[i][j]; } } - PRINT(T2,4,4); + +/* PRINT(T2,4,4); */ +#endif + + T[1][0]*=-1; + T[1][1]*=-1; + T[1][2]*=-1; /* -- cgit v1.3.1 From 4dc1d72785c660c206f8def9d8c8aa32289c2709 Mon Sep 17 00:00:00 2001 From: mwturvey Date: Fri, 24 Mar 2017 15:19:59 -0700 Subject: More cleanup & finishing genericization of calibrator --- calibrate.c | 4 ++-- include/libsurvive/survive.h | 2 +- include/libsurvive/survive_types.h | 2 +- redist/CNFGWinDriver.c | 2 +- redist/json_helpers.c | 2 +- src/survive_cal.c | 30 +++++++++++++++++++----------- src/survive_cal.h | 6 ++++-- src/survive_data.c | 6 +++--- src/survive_process.c | 6 +++--- src/survive_vive.c | 6 ++---- 10 files changed, 37 insertions(+), 29 deletions(-) (limited to 'redist') diff --git a/calibrate.c b/calibrate.c index ee22abf..abf592a 100644 --- a/calibrate.c +++ b/calibrate.c @@ -99,9 +99,9 @@ void my_imu_process( struct SurviveObject * so, int mask, FLT * accelgyro, uint3 } -void my_angle_process( struct SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle ) +void my_angle_process( struct SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle, uint32_t lh) { - survive_default_angle_process( so, sensor_id, acode, timecode, length, angle ); + survive_default_angle_process( so, sensor_id, acode, timecode, length, angle, lh ); } char* sensor_name[32]; diff --git a/include/libsurvive/survive.h b/include/libsurvive/survive.h index 1165e9d..e13312d 100644 --- a/include/libsurvive/survive.h +++ b/include/libsurvive/survive.h @@ -132,7 +132,7 @@ void survive_cal_install( SurviveContext * ctx ); //XXX This will be removed if //Accept higher-level data. void survive_default_light_process( SurviveObject * so, int sensor_id, int acode, int timeinsweep, uint32_t timecode, uint32_t length , uint32_t lh); void survive_default_imu_process( SurviveObject * so, int mode, FLT * accelgyro, uint32_t timecode, int id ); -void survive_default_angle_process( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle ); +void survive_default_angle_process( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle, uint32_t lh ); ////////////////////// Survive Drivers //////////////////////////// diff --git a/include/libsurvive/survive_types.h b/include/libsurvive/survive_types.h index def30b8..bfd0b1d 100644 --- a/include/libsurvive/survive_types.h +++ b/include/libsurvive/survive_types.h @@ -30,7 +30,7 @@ typedef struct SurviveCalData SurviveCalData; //XXX Warning: This may be remov typedef void (*text_feedback_func)( SurviveContext * ctx, const char * fault ); typedef void (*light_process_func)( SurviveObject * so, int sensor_id, int acode, int timeinsweep, uint32_t timecode, uint32_t length, uint32_t lighthouse); typedef void (*imu_process_func)( SurviveObject * so, int mask, FLT * accelgyro, uint32_t timecode, int id ); -typedef void (*angle_process_func)( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle ); +typedef void (*angle_process_func)( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle, uint32_t lh); //Device drivers (prefix your drivers with "DriverReg") i.e. diff --git a/redist/CNFGWinDriver.c b/redist/CNFGWinDriver.c index c5da925..b1c1eb0 100644 --- a/redist/CNFGWinDriver.c +++ b/redist/CNFGWinDriver.c @@ -232,7 +232,7 @@ void CNFGHandleInput() case WM_MBUTTONUP: HandleButton( (msg.lParam & 0xFFFF), (msg.lParam>>16) & 0xFFFF, 3, 0 ); break; case WM_KEYDOWN: case WM_KEYUP: - HandleKey( tolower( msg.wParam ), (msg.message==WM_KEYDOWN) ); + HandleKey( tolower( (int)(msg.wParam) ), (msg.message==WM_KEYDOWN) ); break; default: DispatchMessage(&msg); diff --git a/redist/json_helpers.c b/redist/json_helpers.c index 3b5cc0d..29d48bd 100644 --- a/redist/json_helpers.c +++ b/redist/json_helpers.c @@ -174,7 +174,7 @@ void json_load_file(const char* path) { int16_t children = -1; - for (i=0; i<(int)items; i+=2) + for (i=0; i<(unsigned int)items; i+=2) { //increment i on each successful tag + values combination, not individual tokens jsmntok_t* tag_t = tokens+i; diff --git a/src/survive_cal.c b/src/survive_cal.c index 22a8eff..6c153b4 100755 --- a/src/survive_cal.c +++ b/src/survive_cal.c @@ -220,10 +220,13 @@ void survive_cal_light( struct SurviveObject * so, int sensor_id, int acode, int else if( acode < -4 ) break; int lh = (-acode) - 3; - if( strcmp( so->codename, "WM0" ) == 0 ) - sensor_id += 32; - if( strcmp( so->codename, "WM1" ) == 0 ) - sensor_id += 64; + for (int i=0; i < min(MAX_DEVICES_TO_CAL, cd->numPoseObjects); i++) + { + if( strcmp( so->codename, cd->poseobjects[i]->codename ) == 0 ) + { + sensor_id += i*32; + } + } cd->all_sync_times[sensor_id][lh][cd->all_sync_counts[sensor_id][lh]++] = length; break; @@ -233,7 +236,7 @@ void survive_cal_light( struct SurviveObject * so, int sensor_id, int acode, int } -void survive_cal_angle( struct SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle ) +void survive_cal_angle( struct SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle, uint32_t lh ) { struct SurviveContext * ctx = so->ctx; struct SurviveCalData * cd = ctx->calptr; @@ -241,14 +244,18 @@ void survive_cal_angle( struct SurviveObject * so, int sensor_id, int acode, uin if( !cd ) return; int sensid = sensor_id; - if( strcmp( so->codename, "WM0" ) == 0 ) - sensid += 32; - if( strcmp( so->codename, "WM1" ) == 0 ) - sensid += 64; + + for (int i=0; i < min(MAX_DEVICES_TO_CAL, cd->numPoseObjects); i++) + { + if( strcmp( so->codename, cd->poseobjects[i]->codename ) == 0 ) + { + sensid += i*32; + } + } if( sensid >= MAX_SENSORS_TO_CAL || sensid < 0 ) return; - int lighthouse = acode>>2; + int lighthouse = lh; int axis = acode & 1; switch( cd->stage ) @@ -292,7 +299,8 @@ void survive_cal_angle( struct SurviveObject * so, int sensor_id, int acode, uin int min_peaks = PTS_BEFORE_COMMON; int i, j, k; cd->found_common = 1; - for( i = 0; i < MAX_SENSORS_TO_CAL/SENSORS_PER_OBJECT; i++ ) + for( i = 0; i < cd->numPoseObjects; i++ ) + //for( i = 0; i < MAX_SENSORS_TO_CAL/SENSORS_PER_OBJECT; i++ ) for( j = 0; j < NUM_LIGHTHOUSES; j++ ) { int sensors_visible = 0; diff --git a/src/survive_cal.h b/src/survive_cal.h index 8f4e4de..ae644d1 100644 --- a/src/survive_cal.h +++ b/src/survive_cal.h @@ -30,9 +30,11 @@ int survive_cal_get_status( SurviveContext * ctx, char * description, int descri //Called from survive_default_light_process void survive_cal_light( SurviveObject * so, int sensor_id, int acode, int timeinsweep, uint32_t timecode, uint32_t length, uint32_t lighthouse); -void survive_cal_angle( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle ); +void survive_cal_angle( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle, uint32_t lh ); -#define MAX_SENSORS_TO_CAL 96 +#define MAX_SENSORS_PER_DEVICE 32 +#define MAX_DEVICES_TO_CAL 3 +#define MAX_SENSORS_TO_CAL (MAX_SENSORS_PER_DEVICE * MAX_DEVICES_TO_CAL) #define MIN_PTS_BEFORE_CAL 24 diff --git a/src/survive_data.c b/src/survive_data.c index 4e2479a..9447104 100644 --- a/src/survive_data.c +++ b/src/survive_data.c @@ -25,7 +25,7 @@ typedef struct typedef struct { - float acode_offset; + double acode_offset; } lightcap2_global_data; typedef struct @@ -40,11 +40,11 @@ typedef struct int handle_lightcap2_getAcodeFromSyncPulse(SurviveObject * so, int pulseLen) { - float oldOffset = ((lightcap2_data*)so->disambiguator_data)->global.acode_offset; + double oldOffset = ((lightcap2_data*)so->disambiguator_data)->global.acode_offset; int modifiedPulseLen = pulseLen - (int)oldOffset; - float newOffset = (((pulseLen) + 250) % 500) - 250; + double newOffset = (((pulseLen) + 250) % 500) - 250; ((lightcap2_data*)so->disambiguator_data)->global.acode_offset = oldOffset * 0.9 + newOffset * 0.1; diff --git a/src/survive_process.c b/src/survive_process.c index d4604d8..b58b344 100644 --- a/src/survive_process.c +++ b/src/survive_process.c @@ -37,16 +37,16 @@ void survive_default_light_process( SurviveObject * so, int sensor_id, int acode #endif FLT length_sec = length / (FLT)so->timebase_hz; - ctx->angleproc( so, sensor_id, acode, timecode, length_sec, angle ); + ctx->angleproc( so, sensor_id, acode, timecode, length_sec, angle, lh); } -void survive_default_angle_process( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle ) +void survive_default_angle_process( SurviveObject * so, int sensor_id, int acode, uint32_t timecode, FLT length, FLT angle, uint32_t lh) { SurviveContext * ctx = so->ctx; if( ctx->calptr ) { - survive_cal_angle( so, sensor_id, acode, timecode, length, angle ); + survive_cal_angle( so, sensor_id, acode, timecode, length, angle, lh ); } if( so->PoserFn ) { diff --git a/src/survive_vive.c b/src/survive_vive.c index 55e949a..a5c731d 100755 --- a/src/survive_vive.c +++ b/src/survive_vive.c @@ -834,7 +834,6 @@ static void handle_watchman( SurviveObject * w, uint8_t * readdata ) qty-=2; int propset = 0; int doimu = 0; - int i; if( (type & 0xf0) == 0xf0 ) { @@ -911,10 +910,9 @@ static void handle_watchman( SurviveObject * w, uint8_t * readdata ) *readdata = type; //Put 'type' back on stack. uint8_t * mptr = readdata + qty-3-1; //-3 for timecode, -1 to -//#define DEBUG_WATCHMAN #ifdef DEBUG_WATCHMAN printf( "_%s ", w->codename); - for( i = 0; i < qty; i++ ) + for(int i = 0; i < qty; i++ ) { printf( "%02x ", readdata[i] ); } @@ -1198,7 +1196,7 @@ void survive_data_cb( SurviveUSBInterface * si ) unsigned short *length = (unsigned short *)(&(readdata[2])); unsigned long *time = (unsigned long *)(&(readdata[4])); LightcapElement le; - le.sensor_id = POP2; + le.sensor_id = (uint8_t)POP2; le.length = POP2; le.timestamp = POP4; if( le.sensor_id == 0xff ) break; -- cgit v1.3.1 From 5ef6668c3d2bc0a4e7c364119a99fcd6c6e43017 Mon Sep 17 00:00:00 2001 From: Mike Turvey Date: Sun, 26 Mar 2017 19:08:44 -0700 Subject: Adding svd --- redist/svd.h | 450 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 450 insertions(+) create mode 100644 redist/svd.h (limited to 'redist') diff --git a/redist/svd.h b/redist/svd.h new file mode 100644 index 0000000..41708da --- /dev/null +++ b/redist/svd.h @@ -0,0 +1,450 @@ +/************************************************************************** +** +** svd3 +** +** Quick singular value decomposition as described by: +** A. McAdams, A. Selle, R. Tamstorf, J. Teran and E. Sifakis, +** "Computing the Singular Value Decomposition of 3x3 matrices +** with minimal branching and elementary floating point operations", +** University of Wisconsin - Madison technical report TR1690, May 2011 +** +** OPTIMIZED CPU VERSION +** Implementation by: Eric Jang +** +** 13 Apr 2014 +** +** This file originally retrieved from: +** https://github.com/ericjang/svd3/blob/master/svd3.h 3/26/2017 +** +** Original licesnse is MIT per: +** https://github.com/ericjang/svd3/blob/master/LICENSE.md +** +** Ported from C++ to C by Mike Turvey. All modifications also released +** under an MIT license. +**************************************************************************/ + + +#ifndef SVD3_H +#define SVD3_H + +#define _gamma 5.828427124 // FOUR_GAMMA_SQUARED = sqrt(8)+3; +#define _cstar 0.923879532 // cos(pi/8) +#define _sstar 0.3826834323 // sin(p/8) +#define EPSILON 1e-6 + +#include + +/* This is a novel and fast routine for the reciprocal square root of an +IEEE float (single precision). +http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf +http://playstation2-linux.com/download/p2lsd/fastrsqrt.pdf +http://www.beyond3d.com/content/articles/8/ +*/ +inline float rsqrt(float x) { + // int ihalf = *(int *)&x - 0x00800000; // Alternative to next line, + // float xhalf = *(float *)&ihalf; // for sufficiently large nos. + float xhalf = 0.5f*x; + int i = *(int *)&x; // View x as an int. + // i = 0x5f3759df - (i >> 1); // Initial guess (traditional). + i = 0x5f375a82 - (i >> 1); // Initial guess (slightly better). + x = *(float *)&i; // View i as float. + x = x*(1.5f - xhalf*x*x); // Newton step. + // x = x*(1.5008908 - xhalf*x*x); // Newton step for a balanced error. + return x; +} + +/* This is rsqrt with an additional step of the Newton iteration, for +increased accuracy. The constant 0x5f37599e makes the relative error +range from 0 to -0.00000463. +You can't balance the error by adjusting the constant. */ +inline float rsqrt1(float x) { + float xhalf = 0.5f*x; + int i = *(int *)&x; // View x as an int. + i = 0x5f37599e - (i >> 1); // Initial guess. + x = *(float *)&i; // View i as float. + x = x*(1.5f - xhalf*x*x); // Newton step. + x = x*(1.5f - xhalf*x*x); // Newton step again. + return x; +} + +inline float accurateSqrt(float x) +{ + return x * rsqrt1(x); +} + +inline void condSwap(bool c, float *X, float *Y) +{ + // used in step 2 + float Z = *X; + *X = c ? *Y : *X; + *Y = c ? Z : *Y; +} + +inline void condNegSwap(bool c, float *X, float *Y) +{ + // used in step 2 and 3 + float Z = -*X; + *X = c ? *Y : *X; + *Y = c ? Z : *Y; +} + +// matrix multiplication M = A * B +inline void multAB(float a11, float a12, float a13, + float a21, float a22, float a23, + float a31, float a32, float a33, + // + float b11, float b12, float b13, + float b21, float b22, float b23, + float b31, float b32, float b33, + // + float *m11, float *m12, float *m13, + float *m21, float *m22, float *m23, + float *m31, float *m32, float *m33) +{ + + *m11 = a11*b11 + a12*b21 + a13*b31; *m12 = a11*b12 + a12*b22 + a13*b32; *m13 = a11*b13 + a12*b23 + a13*b33; + *m21 = a21*b11 + a22*b21 + a23*b31; *m22 = a21*b12 + a22*b22 + a23*b32; *m23 = a21*b13 + a22*b23 + a23*b33; + *m31 = a31*b11 + a32*b21 + a33*b31; *m32 = a31*b12 + a32*b22 + a33*b32; *m33 = a31*b13 + a32*b23 + a33*b33; +} + +// matrix multiplication M = Transpose[A] * B +inline void multAtB(float a11, float a12, float a13, + float a21, float a22, float a23, + float a31, float a32, float a33, + // + float b11, float b12, float b13, + float b21, float b22, float b23, + float b31, float b32, float b33, + // + float *m11, float *m12, float *m13, + float *m21, float *m22, float *m23, + float *m31, float *m32, float *m33) +{ + *m11 = a11*b11 + a21*b21 + a31*b31; *m12 = a11*b12 + a21*b22 + a31*b32; *m13 = a11*b13 + a21*b23 + a31*b33; + *m21 = a12*b11 + a22*b21 + a32*b31; *m22 = a12*b12 + a22*b22 + a32*b32; *m23 = a12*b13 + a22*b23 + a32*b33; + *m31 = a13*b11 + a23*b21 + a33*b31; *m32 = a13*b12 + a23*b22 + a33*b32; *m33 = a13*b13 + a23*b23 + a33*b33; +} + +inline void quatToMat3(const float * qV, + float *m11, float *m12, float *m13, + float *m21, float *m22, float *m23, + float *m31, float *m32, float *m33 +) +{ + float w = qV[3]; + float x = qV[0]; + float y = qV[1]; + float z = qV[2]; + + float qxx = x*x; + float qyy = y*y; + float qzz = z*z; + float qxz = x*z; + float qxy = x*y; + float qyz = y*z; + float qwx = w*x; + float qwy = w*y; + float qwz = w*z; + + *m11 = 1 - 2 * (qyy + qzz); *m12 = 2 * (qxy - qwz); *m13 = 2 * (qxz + qwy); + *m21 = 2 * (qxy + qwz); *m22 = 1 - 2 * (qxx + qzz); *m23 = 2 * (qyz - qwx); + *m31 = 2 * (qxz - qwy); *m32 = 2 * (qyz + qwx); *m33 = 1 - 2 * (qxx + qyy); +} + +inline void approximateGivensQuaternion(float a11, float a12, float a22, float *ch, float *sh) +{ + /* + * Given givens angle computed by approximateGivensAngles, + * compute the corresponding rotation quaternion. + */ + *ch = 2 * (a11 - a22); + *sh = a12; + bool b = _gamma* (*sh)*(*sh) < (*ch)*(*ch); + // fast rsqrt function suffices + // rsqrt2 (https://code.google.com/p/lppython/source/browse/algorithm/HDcode/newCode/rsqrt.c?r=26) + // is even faster but results in too much error + float w = rsqrt((*ch)*(*ch) + (*sh)*(*sh)); + *ch = b ? w*(*ch) : (float)_cstar; + *sh = b ? w*(*sh) : (float)_sstar; +} + +inline void jacobiConjugation(const int x, const int y, const int z, + float *s11, + float *s21, float *s22, + float *s31, float *s32, float *s33, + float * qV) +{ + float ch, sh; + approximateGivensQuaternion(*s11, *s21, *s22, &ch, &sh); + + float scale = ch*ch + sh*sh; + float a = (ch*ch - sh*sh) / scale; + float b = (2 * sh*ch) / scale; + + // make temp copy of S + float _s11 = *s11; + float _s21 = *s21; float _s22 = *s22; + float _s31 = *s31; float _s32 = *s32; float _s33 = *s33; + + // perform conjugation S = Q'*S*Q + // Q already implicitly solved from a, b + *s11 = a*(a*_s11 + b*_s21) + b*(a*_s21 + b*_s22); + *s21 = a*(-b*_s11 + a*_s21) + b*(-b*_s21 + a*_s22); *s22 = -b*(-b*_s11 + a*_s21) + a*(-b*_s21 + a*_s22); + *s31 = a*_s31 + b*_s32; *s32 = -b*_s31 + a*_s32; *s33 = _s33; + + // update cumulative rotation qV + float tmp[3]; + tmp[0] = qV[0] * sh; + tmp[1] = qV[1] * sh; + tmp[2] = qV[2] * sh; + sh *= qV[3]; + + qV[0] *= ch; + qV[1] *= ch; + qV[2] *= ch; + qV[3] *= ch; + + // (x,y,z) corresponds to ((0,1,2),(1,2,0),(2,0,1)) + // for (p,q) = ((0,1),(1,2),(0,2)) + qV[z] += sh; + qV[3] -= tmp[z]; // w + qV[x] += tmp[y]; + qV[y] -= tmp[x]; + + // re-arrange matrix for next iteration + _s11 = *s22; + _s21 = *s32; _s22 = *s33; + _s31 = *s21; _s32 = *s31; _s33 = *s11; + *s11 = _s11; + *s21 = _s21; *s22 = _s22; + *s31 = _s31; *s32 = _s32; *s33 = _s33; + +} + +inline float dist2(float x, float y, float z) +{ + return x*x + y*y + z*z; +} + +// finds transformation that diagonalizes a symmetric matrix +inline void jacobiEigenanlysis( // symmetric matrix + float *s11, + float *s21, float *s22, + float *s31, float *s32, float *s33, + // quaternion representation of V + float * qV) +{ + qV[3] = 1; qV[0] = 0; qV[1] = 0; qV[2] = 0; // follow same indexing convention as GLM + for (int i = 0; i<4; i++) + { + // we wish to eliminate the maximum off-diagonal element + // on every iteration, but cycling over all 3 possible rotations + // in fixed order (p,q) = (1,2) , (2,3), (1,3) still retains + // asymptotic convergence + jacobiConjugation(0, 1, 2, s11, s21, s22, s31, s32, s33, qV); // p,q = 0,1 + jacobiConjugation(1, 2, 0, s11, s21, s22, s31, s32, s33, qV); // p,q = 1,2 + jacobiConjugation(2, 0, 1, s11, s21, s22, s31, s32, s33, qV); // p,q = 0,2 + } +} + + +inline void sortSingularValues(// matrix that we want to decompose + float *b11, float *b12, float *b13, + float *b21, float *b22, float *b23, + float *b31, float *b32, float *b33, + // sort V simultaneously + float *v11, float *v12, float *v13, + float *v21, float *v22, float *v23, + float *v31, float *v32, float *v33) +{ + float rho1 = dist2(*b11, *b21, *b31); + float rho2 = dist2(*b12, *b22, *b32); + float rho3 = dist2(*b13, *b23, *b33); + bool c; + c = rho1 < rho2; + condNegSwap(c, b11, b12); condNegSwap(c, v11, v12); + condNegSwap(c, b21, b22); condNegSwap(c, v21, v22); + condNegSwap(c, b31, b32); condNegSwap(c, v31, v32); + condSwap(c, &rho1, &rho2); + c = rho1 < rho3; + condNegSwap(c, b11, b13); condNegSwap(c, v11, v13); + condNegSwap(c, b21, b23); condNegSwap(c, v21, v23); + condNegSwap(c, b31, b33); condNegSwap(c, v31, v33); + condSwap(c, &rho1, &rho3); + c = rho2 < rho3; + condNegSwap(c, b12, b13); condNegSwap(c, v12, v13); + condNegSwap(c, b22, b23); condNegSwap(c, v22, v23); + condNegSwap(c, b32, b33); condNegSwap(c, v32, v33); +} + + +void QRGivensQuaternion(float a1, float a2, float *ch, float *sh) +{ + // a1 = pivot point on diagonal + // a2 = lower triangular entry we want to annihilate + float epsilon = (float)EPSILON; + float rho = accurateSqrt(a1*a1 + a2*a2); + + *sh = rho > epsilon ? a2 : 0; + *ch = fabsf(a1) + fmaxf(rho, epsilon); + bool b = a1 < 0; + condSwap(b, sh, ch); + float w = rsqrt((*ch)*(*ch) + (*sh)*(*sh)); + *ch *= w; + *sh *= w; +} + + +inline void QRDecomposition(// matrix that we want to decompose + float b11, float b12, float b13, + float b21, float b22, float b23, + float b31, float b32, float b33, + // output Q + float *q11, float *q12, float *q13, + float *q21, float *q22, float *q23, + float *q31, float *q32, float *q33, + // output R + float *r11, float *r12, float *r13, + float *r21, float *r22, float *r23, + float *r31, float *r32, float *r33) +{ + float ch1, sh1, ch2, sh2, ch3, sh3; + float a, b; + + // first givens rotation (ch,0,0,sh) + QRGivensQuaternion(b11, b21, &ch1, &sh1); + a = 1 - 2 * sh1*sh1; + b = 2 * ch1*sh1; + // apply B = Q' * B + *r11 = a*b11 + b*b21; *r12 = a*b12 + b*b22; *r13 = a*b13 + b*b23; + *r21 = -b*b11 + a*b21; *r22 = -b*b12 + a*b22; *r23 = -b*b13 + a*b23; + *r31 = b31; *r32 = b32; *r33 = b33; + + // second givens rotation (ch,0,-sh,0) + QRGivensQuaternion(*r11, *r31, &ch2, &sh2); + a = 1 - 2 * sh2*sh2; + b = 2 * ch2*sh2; + // apply B = Q' * B; + b11 = a*(*r11) + b*(*r31); b12 = a*(*r12) + b*(*r32); b13 = a*(*r13) + b*(*r33); + b21 = *r21; b22 = *r22; b23 = *r23; + b31 = -b*(*r11) + a*(*r31); b32 = -b*(*r12) + a*(*r32); b33 = -b*(*r13) + a*(*r33); + + // third givens rotation (ch,sh,0,0) + QRGivensQuaternion(b22, b32, &ch3, &sh3); + a = 1 - 2 * sh3*sh3; + b = 2 * ch3*sh3; + // R is now set to desired value + *r11 = b11; *r12 = b12; *r13 = b13; + *r21 = a*b21 + b*b31; *r22 = a*b22 + b*b32; *r23 = a*b23 + b*b33; + *r31 = -b*b21 + a*b31; *r32 = -b*b22 + a*b32; *r33 = -b*b23 + a*b33; + + // construct the cumulative rotation Q=Q1 * Q2 * Q3 + // the number of floating point operations for three quaternion multiplications + // is more or less comparable to the explicit form of the joined matrix. + // certainly more memory-efficient! + float sh12 = sh1*sh1; + float sh22 = sh2*sh2; + float sh32 = sh3*sh3; + + *q11 = (-1 + 2 * sh12)*(-1 + 2 * sh22); + *q12 = 4 * ch2*ch3*(-1 + 2 * sh12)*sh2*sh3 + 2 * ch1*sh1*(-1 + 2 * sh32); + *q13 = 4 * ch1*ch3*sh1*sh3 - 2 * ch2*(-1 + 2 * sh12)*sh2*(-1 + 2 * sh32); + + *q21 = 2 * ch1*sh1*(1 - 2 * sh22); + *q22 = -8 * ch1*ch2*ch3*sh1*sh2*sh3 + (-1 + 2 * sh12)*(-1 + 2 * sh32); + *q23 = -2 * ch3*sh3 + 4 * sh1*(ch3*sh1*sh3 + ch1*ch2*sh2*(-1 + 2 * sh32)); + + *q31 = 2 * ch2*sh2; + *q32 = 2 * ch3*(1 - 2 * sh22)*sh3; + *q33 = (-1 + 2 * sh22)*(-1 + 2 * sh32); +} + +void svd(// input A + float a11, float a12, float a13, + float a21, float a22, float a23, + float a31, float a32, float a33, + // output U + float *u11, float *u12, float *u13, + float *u21, float *u22, float *u23, + float *u31, float *u32, float *u33, + // output S + float *s11, float *s12, float *s13, + float *s21, float *s22, float *s23, + float *s31, float *s32, float *s33, + // output V + float *v11, float *v12, float *v13, + float *v21, float *v22, float *v23, + float *v31, float *v32, float *v33) +{ + // normal equations matrix + float ATA11, ATA12, ATA13; + float ATA21, ATA22, ATA23; + float ATA31, ATA32, ATA33; + + multAtB(a11, a12, a13, a21, a22, a23, a31, a32, a33, + a11, a12, a13, a21, a22, a23, a31, a32, a33, + &ATA11, &ATA12, &ATA13, &ATA21, &ATA22, &ATA23, &ATA31, &ATA32, &ATA33); + + // symmetric eigenalysis + float qV[4]; + jacobiEigenanlysis(&ATA11, &ATA21, &ATA22, &ATA31, &ATA32, &ATA33, qV); + quatToMat3(qV, v11, v12, v13, v21, v22, v23, v31, v32, v33); + + float b11, b12, b13; + float b21, b22, b23; + float b31, b32, b33; + multAB(a11, a12, a13, a21, a22, a23, a31, a32, a33, + *v11, *v12, *v13, *v21, *v22, *v23, *v31, *v32, *v33, + &b11, &b12, &b13, &b21, &b22, &b23, &b31, &b32, &b33); + + // sort singular values and find V + sortSingularValues(&b11, &b12, &b13, &b21, &b22, &b23, &b31, &b32, &b33, + v11, v12, v13, v21, v22, v23, v31, v32, v33); + + // QR decomposition + QRDecomposition(b11, b12, b13, b21, b22, b23, b31, b32, b33, + u11, u12, u13, u21, u22, u23, u31, u32, u33, + s11, s12, s13, s21, s22, s23, s31, s32, s33 + ); +} + +/// polar decomposition can be reconstructed trivially from SVD result +// A = UP +void pd(float a11, float a12, float a13, + float a21, float a22, float a23, + float a31, float a32, float a33, + // output U + float *u11, float *u12, float *u13, + float *u21, float *u22, float *u23, + float *u31, float *u32, float *u33, + // output P + float *p11, float *p12, float *p13, + float *p21, float *p22, float *p23, + float *p31, float *p32, float *p33) +{ + float w11, w12, w13, w21, w22, w23, w31, w32, w33; + float s11, s12, s13, s21, s22, s23, s31, s32, s33; + float v11, v12, v13, v21, v22, v23, v31, v32, v33; + + svd(a11, a12, a13, a21, a22, a23, a31, a32, a33, + &w11, &w12, &w13, &w21, &w22, &w23, &w31, &w32, &w33, + &s11, &s12, &s13, &s21, &s22, &s23, &s31, &s32, &s33, + &v11, &v12, &v13, &v21, &v22, &v23, &v31, &v32, &v33); + + // P = VSV' + float t11, t12, t13, t21, t22, t23, t31, t32, t33; + multAB(v11, v12, v13, v21, v22, v23, v31, v32, v33, + s11, s12, s13, s21, s22, s23, s31, s32, s33, + &t11, &t12, &t13, &t21, &t22, &t23, &t31, &t32, &t33); + + multAB(t11, t12, t13, t21, t22, t23, t31, t32, t33, + v11, v21, v31, v12, v22, v32, v13, v23, v33, + p11, p12, p13, p21, p22, p23, p31, p32, p33); + + // U = WV' + multAB(w11, w12, w13, w21, w22, w23, w31, w32, w33, + v11, v21, v31, v12, v22, v32, v13, v23, v33, + u11, u12, u13, u21, u22, u23, u31, u32, u33); +} + +#endif \ No newline at end of file -- cgit v1.3.1 From ad5afae7e264ffa72ed16d2ec2c6f6090e7aa7b5 Mon Sep 17 00:00:00 2001 From: mwturvey Date: Mon, 27 Mar 2017 13:02:23 -0700 Subject: Adding operation for rotating about an axis --- redist/linmath.c | 19 +++++++++++++++++++ redist/linmath.h | 1 + 2 files changed, 20 insertions(+) (limited to 'redist') diff --git a/redist/linmath.c b/redist/linmath.c index eefcd5f..0e06156 100644 --- a/redist/linmath.c +++ b/redist/linmath.c @@ -82,6 +82,25 @@ FLT anglebetween3d( FLT * a, FLT * b ) return FLT_ACOS(dot); } +// algorithm found here: http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/ +void rotatearoundaxis(FLT *outvec3, FLT *invec3, FLT *axis, FLT angle) +{ + FLT s = FLT_SIN(angle); + FLT c = FLT_COS(angle); + + FLT u=axis[0]; + FLT v=axis[1]; + FLT w=axis[2]; + + FLT x=invec3[0]; + FLT y=invec3[1]; + FLT z=invec3[2]; + + outvec3[0] = u*(u*x + v*y + w*z)*(1-c) + x*c + (-w*y + v*z)*s; + outvec3[1] = v*(u*x + v*y + w*z)*(1-c) + y*c + ( w*x - u*z)*s; + outvec3[2] = w*(u*x + v*y + w*z)*(1-c) + z*c + (-v*x + u*y)*s; +} + /////////////////////////////////////QUATERNIONS////////////////////////////////////////// //Originally from Mercury (Copyright (C) 2009 by Joshua Allen, Charles Lohr, Adam Lowman) //Under the mit/X11 license. diff --git a/redist/linmath.h b/redist/linmath.h index 4e0cb77..6f0bf60 100644 --- a/redist/linmath.h +++ b/redist/linmath.h @@ -70,6 +70,7 @@ FLT magnitude3d(const FLT * a ); FLT anglebetween3d( FLT * a, FLT * b ); +void rotatearoundaxis(FLT *outvec3, FLT *invec3, FLT *axis, FLT angle); //Quaternion things... void quatsetnone( FLT * q ); -- cgit v1.3.1 From e5c15af3af93356bb056624726e0b6068354690f Mon Sep 17 00:00:00 2001 From: mwturvey Date: Mon, 27 Mar 2017 14:20:29 -0700 Subject: Getting very close --- redist/linmath.c | 3 +++ src/poser_turveytori.c | 17 ++++++++++++++++- 2 files changed, 19 insertions(+), 1 deletion(-) (limited to 'redist') diff --git a/redist/linmath.c b/redist/linmath.c index 0e06156..69a70f6 100644 --- a/redist/linmath.c +++ b/redist/linmath.c @@ -85,6 +85,9 @@ FLT anglebetween3d( FLT * a, FLT * b ) // algorithm found here: http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/ void rotatearoundaxis(FLT *outvec3, FLT *invec3, FLT *axis, FLT angle) { + // TODO: this really should be external. + normalize3d(axis, axis); + FLT s = FLT_SIN(angle); FLT c = FLT_COS(angle); diff --git a/src/poser_turveytori.c b/src/poser_turveytori.c index 824fabb..d737723 100644 --- a/src/poser_turveytori.c +++ b/src/poser_turveytori.c @@ -779,8 +779,10 @@ static void RefineRotationEstimate(FLT *rotOut, Point lhPoint, FLT *initialEstim lastMatchFitness = newMatchFitness; quatcopy(rotOut, point4); //#ifdef TORI_DEBUG - printf("+ %8.8f, %f\n", newMatchFitness, point4[3]); + printf("+ %8.8f, (%8.8f, %8.8f, %8.8f) %f\n", newMatchFitness, point4[0], point4[1], point4[2], point4[3]); //#endif + g *= 1.03; + } else { @@ -796,6 +798,17 @@ static void RefineRotationEstimate(FLT *rotOut, Point lhPoint, FLT *initialEstim printf("\nRi=%d\n", i); } +static void WhereIsTheTrackedObject(FLT *rotation, Point lhPoint) +{ + FLT reverseRotation[4] = {rotation[0], rotation[1], rotation[2], -rotation[3]}; + FLT objPoint[3] = {lhPoint.x, lhPoint.y, lhPoint.z}; + + rotatearoundaxis(objPoint, objPoint, reverseRotation, reverseRotation[3]); + + printf("The tracked object is at location (%f, %f, %f)\n", objPoint[0], objPoint[1], objPoint[2]); +} + + void SolveForRotation(FLT rotOut[4], TrackedObject *obj, Point lh) { @@ -812,6 +825,8 @@ void SolveForRotation(FLT rotOut[4], TrackedObject *obj, Point lh) // Step 2, optimize the quaternion to match the data. RefineRotationEstimate(rotOut, lh, zAxis, obj); + WhereIsTheTrackedObject(rotOut, lh); + } -- cgit v1.3.1