PageRenderTime 108ms CodeModel.GetById 49ms app.highlight 54ms RepoModel.GetById 1ms app.codeStats 0ms

/src/core/MathUtils.cpp

http://github.com/imageworks/OpenColorIO
C++ | 622 lines | 470 code | 110 blank | 42 comment | 41 complexity | e8f424aa4f1d31cc19edd210490bced5 MD5 | raw file
  1/*
  2Copyright (c) 2003-2010 Sony Pictures Imageworks Inc., et al.
  3All Rights Reserved.
  4
  5Redistribution and use in source and binary forms, with or without
  6modification, are permitted provided that the following conditions are
  7met:
  8* Redistributions of source code must retain the above copyright
  9  notice, this list of conditions and the following disclaimer.
 10* Redistributions in binary form must reproduce the above copyright
 11  notice, this list of conditions and the following disclaimer in the
 12  documentation and/or other materials provided with the distribution.
 13* Neither the name of Sony Pictures Imageworks nor the names of its
 14  contributors may be used to endorse or promote products derived from
 15  this software without specific prior written permission.
 16THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 17"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 18LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 19A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 20OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 21SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 22LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 23DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 24THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 25(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 26OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 27*/
 28
 29#include <OpenColorIO/OpenColorIO.h>
 30
 31#include <cstring>
 32
 33#include "MathUtils.h"
 34
 35OCIO_NAMESPACE_ENTER
 36{
 37    namespace
 38    {
 39        const float FLTMIN = std::numeric_limits<float>::min();
 40    }
 41    
 42    bool IsScalarEqualToZero(float v)
 43    {
 44        return equalWithAbsError(v, 0.0f, FLTMIN);
 45    }
 46
 47    bool IsScalarEqualToZeroFlt(double v)
 48    {
 49        return equalWithAbsError(float(v), 0.0f, FLTMIN);
 50    }
 51
 52    bool IsScalarEqualToOne(float v)
 53    {
 54        return equalWithAbsError(v, 1.0f, FLTMIN);
 55    }
 56    
 57    bool IsScalarEqualToOneFlt(double v)
 58    {
 59        return equalWithAbsError(float(v), 1.0f, FLTMIN);
 60    }
 61
 62    float GetSafeScalarInverse(float v, float defaultValue)
 63    {
 64        if(IsScalarEqualToZero(v)) return defaultValue;
 65        return 1.0f / v;
 66    }
 67    
 68    bool IsVecEqualToZero(const float* v, int size)
 69    {
 70        for(int i=0; i<size; ++i)
 71        {
 72            if(!IsScalarEqualToZero(v[i])) return false;
 73        }
 74        return true;
 75    }
 76    
 77    bool IsVecEqualToOne(const float* v, int size)
 78    {
 79        for(int i=0; i<size; ++i)
 80        {
 81            if(!IsScalarEqualToOne(v[i])) return false;
 82        }
 83        return true;
 84    }
 85
 86    bool IsVecEqualToOneFlt(const double* v, int size)
 87    {
 88        for(int i=0; i<size; ++i)
 89        {
 90            if(!IsScalarEqualToOneFlt(v[i])) return false;
 91        }
 92        return true;
 93    }
 94
 95    bool VecContainsZero(const float* v, int size)
 96    {
 97        for(int i=0; i<size; ++i)
 98        {
 99            if(IsScalarEqualToZero(v[i])) return true;
100        }
101        return false;
102    }
103    
104    bool VecContainsOne(const float* v, int size)
105    {
106        for(int i=0; i<size; ++i)
107        {
108            if(IsScalarEqualToOne(v[i])) return true;
109        }
110        return false;
111    }
112    
113    bool VecsEqualWithRelError(const float* v1, int size1,
114                               const float* v2, int size2,
115                               float e)
116    {
117        if(size1 != size2) return false;
118        for(int i=0; i<size1; ++i)
119        {
120            if(!equalWithRelError(v1[i], v2[i], e)) return false;
121        }
122        
123        return true;
124    }
125    
126    double ClampToNormHalf(double val)
127    {
128        if(val < -GetHalfMax())
129        {
130            return -GetHalfMax();
131        }
132
133        if(val > -GetHalfNormMin() && val<GetHalfNormMin())
134        {
135            return 0.0;
136        }
137
138        if(val > GetHalfMax())
139        {
140            return GetHalfMax();
141        }
142
143        return val;
144    }
145
146    bool IsM44Identity(const float* m44)
147    {
148        int index=0;
149
150        for(unsigned int j=0; j<4; ++j)
151        {
152            for(unsigned int i=0; i<4; ++i)
153            {
154                index = 4*j+i;
155
156                if(i==j)
157                {
158                    if(!IsScalarEqualToOne(m44[index])) return false;
159                }
160                else
161                {
162                    if(!IsScalarEqualToZero(m44[index])) return false;
163                }
164            }
165        }
166        
167        return true;
168    }
169    
170    bool IsM44Diagonal(const float* m44)
171    {
172        for(int i=0; i<16; ++i)
173        {
174            if((i%5)==0) continue; // If we're on the diagonal, skip it
175            if(!IsScalarEqualToZero(m44[i])) return false;
176        }
177        
178        return true;
179    }
180    
181    void GetM44Diagonal(float* out4, const float* m44)
182    {
183        for(int i=0; i<4; ++i)
184        {
185            out4[i] = m44[i*5];
186        }
187    }
188    
189    // We use an intermediate double representation to make sure
190    // there is minimal float precision error on the determinant's computation
191    // (We have seen IsScalarEqualToZero sensitivities here on 32-bit
192    // virtual machines)
193    
194    bool GetM44Inverse(float* inverse_out, const float* m_)
195    {
196        double m[16];
197        for(unsigned int i=0; i<16; ++i) m[i] = (double)m_[i];
198        
199        double d10_21 = m[4]*m[9] - m[5]*m[8];
200        double d10_22 = m[4]*m[10] - m[6]*m[8];
201        double d10_23 = m[4]*m[11] - m[7]*m[8];
202        double d11_22 = m[5]*m[10] - m[6]*m[9];
203        double d11_23 = m[5]*m[11] - m[7]*m[9];
204        double d12_23 = m[6]*m[11] - m[7]*m[10];
205        
206        double a00 = m[13]*d12_23 - m[14]*d11_23 + m[15]*d11_22;
207        double a10 = m[14]*d10_23 - m[15]*d10_22 - m[12]*d12_23;
208        double a20 = m[12]*d11_23 - m[13]*d10_23 + m[15]*d10_21;
209        double a30 = m[13]*d10_22 - m[14]*d10_21 - m[12]*d11_22;
210        
211        double det = a00*m[0] + a10*m[1] + a20*m[2] + a30*m[3];
212        
213        if(IsScalarEqualToZero((float)det)) return false;
214        
215        det = 1.0/det;
216        
217        double d00_31 = m[0]*m[13] - m[1]*m[12];
218        double d00_32 = m[0]*m[14] - m[2]*m[12];
219        double d00_33 = m[0]*m[15] - m[3]*m[12];
220        double d01_32 = m[1]*m[14] - m[2]*m[13];
221        double d01_33 = m[1]*m[15] - m[3]*m[13];
222        double d02_33 = m[2]*m[15] - m[3]*m[14];
223        
224        double a01 = m[9]*d02_33 - m[10]*d01_33 + m[11]*d01_32;
225        double a11 = m[10]*d00_33 - m[11]*d00_32 - m[8]*d02_33;
226        double a21 = m[8]*d01_33 - m[9]*d00_33 + m[11]*d00_31;
227        double a31 = m[9]*d00_32 - m[10]*d00_31 - m[8]*d01_32;
228        
229        double a02 = m[6]*d01_33 - m[7]*d01_32 - m[5]*d02_33;
230        double a12 = m[4]*d02_33 - m[6]*d00_33 + m[7]*d00_32;
231        double a22 = m[5]*d00_33 - m[7]*d00_31 - m[4]*d01_33;
232        double a32 = m[4]*d01_32 - m[5]*d00_32 + m[6]*d00_31;
233        
234        double a03 = m[2]*d11_23 - m[3]*d11_22 - m[1]*d12_23;
235        double a13 = m[0]*d12_23 - m[2]*d10_23 + m[3]*d10_22;
236        double a23 = m[1]*d10_23 - m[3]*d10_21 - m[0]*d11_23;
237        double a33 = m[0]*d11_22 - m[1]*d10_22 + m[2]*d10_21;
238        
239        inverse_out[0] = (float) (a00*det);
240        inverse_out[1] = (float) (a01*det);
241        inverse_out[2] = (float) (a02*det);
242        inverse_out[3] = (float) (a03*det);
243        inverse_out[4] = (float) (a10*det);
244        inverse_out[5] = (float) (a11*det);
245        inverse_out[6] = (float) (a12*det);
246        inverse_out[7] = (float) (a13*det);
247        inverse_out[8] = (float) (a20*det);
248        inverse_out[9] = (float) (a21*det);
249        inverse_out[10] = (float) (a22*det);
250        inverse_out[11] = (float) (a23*det);
251        inverse_out[12] = (float) (a30*det);
252        inverse_out[13] = (float) (a31*det);
253        inverse_out[14] = (float) (a32*det);
254        inverse_out[15] = (float) (a33*det);
255        
256        return true;
257    }
258    
259    void GetM44M44Product(float* mout, const float* m1_, const float* m2_)
260    {
261        float m1[16];
262        float m2[16];
263        memcpy(m1, m1_, 16*sizeof(float));
264        memcpy(m2, m2_, 16*sizeof(float));
265        
266        mout[ 0] = m1[ 0]*m2[0] + m1[ 1]*m2[4] + m1[ 2]*m2[ 8] + m1[ 3]*m2[12];
267        mout[ 1] = m1[ 0]*m2[1] + m1[ 1]*m2[5] + m1[ 2]*m2[ 9] + m1[ 3]*m2[13];
268        mout[ 2] = m1[ 0]*m2[2] + m1[ 1]*m2[6] + m1[ 2]*m2[10] + m1[ 3]*m2[14];
269        mout[ 3] = m1[ 0]*m2[3] + m1[ 1]*m2[7] + m1[ 2]*m2[11] + m1[ 3]*m2[15];
270        mout[ 4] = m1[ 4]*m2[0] + m1[ 5]*m2[4] + m1[ 6]*m2[ 8] + m1[ 7]*m2[12];
271        mout[ 5] = m1[ 4]*m2[1] + m1[ 5]*m2[5] + m1[ 6]*m2[ 9] + m1[ 7]*m2[13];
272        mout[ 6] = m1[ 4]*m2[2] + m1[ 5]*m2[6] + m1[ 6]*m2[10] + m1[ 7]*m2[14];
273        mout[ 7] = m1[ 4]*m2[3] + m1[ 5]*m2[7] + m1[ 6]*m2[11] + m1[ 7]*m2[15];   
274        mout[ 8] = m1[ 8]*m2[0] + m1[ 9]*m2[4] + m1[10]*m2[ 8] + m1[11]*m2[12];
275        mout[ 9] = m1[ 8]*m2[1] + m1[ 9]*m2[5] + m1[10]*m2[ 9] + m1[11]*m2[13];
276        mout[10] = m1[ 8]*m2[2] + m1[ 9]*m2[6] + m1[10]*m2[10] + m1[11]*m2[14];
277        mout[11] = m1[ 8]*m2[3] + m1[ 9]*m2[7] + m1[10]*m2[11] + m1[11]*m2[15];
278        mout[12] = m1[12]*m2[0] + m1[13]*m2[4] + m1[14]*m2[ 8] + m1[15]*m2[12];
279        mout[13] = m1[12]*m2[1] + m1[13]*m2[5] + m1[14]*m2[ 9] + m1[15]*m2[13];
280        mout[14] = m1[12]*m2[2] + m1[13]*m2[6] + m1[14]*m2[10] + m1[15]*m2[14];
281        mout[15] = m1[12]*m2[3] + m1[13]*m2[7] + m1[14]*m2[11] + m1[15]*m2[15];
282    }
283    
284    namespace
285    {
286    
287    void GetM44V4Product(float* vout, const float* m, const float* v_)
288    {
289        float v[4];
290        memcpy(v, v_, 4*sizeof(float));
291        
292        vout[0] = m[ 0]*v[0] + m[ 1]*v[1] + m[ 2]*v[2] + m[ 3]*v[3];
293        vout[1] = m[ 4]*v[0] + m[ 5]*v[1] + m[ 6]*v[2] + m[ 7]*v[3];
294        vout[2] = m[ 8]*v[0] + m[ 9]*v[1] + m[10]*v[2] + m[11]*v[3];
295        vout[3] = m[12]*v[0] + m[13]*v[1] + m[14]*v[2] + m[15]*v[3];
296    }
297    
298    void GetV4Sum(float* vout, const float* v1, const float* v2)
299    {
300        for(int i=0; i<4; ++i)
301        {
302            vout[i] = v1[i] + v2[i];
303        }
304    }
305    
306    } // anon namespace
307    
308    // All m(s) are 4x4.  All v(s) are size 4 vectors.
309    // Return mout, vout, where mout*x+vout == m2*(m1*x+v1)+v2
310    // mout = m2*m1
311    // vout = m2*v1 + v2
312    void GetMxbCombine(float* mout, float* vout,
313                       const float* m1_, const float* v1_,
314                       const float* m2_, const float* v2_)
315    {
316        float m1[16];
317        float v1[4];
318        float m2[16];
319        float v2[4];
320        memcpy(m1, m1_, 16*sizeof(float));
321        memcpy(v1, v1_, 4*sizeof(float));
322        memcpy(m2, m2_, 16*sizeof(float));
323        memcpy(v2, v2_, 4*sizeof(float));
324        
325        GetM44M44Product(mout, m2, m1);
326        GetM44V4Product(vout, m2, v1);
327        GetV4Sum(vout, vout, v2);
328    }
329    
330    namespace
331    {
332    
333    void GetMxbResult(float* vout, float* m, float* x, float* v)
334    {
335        GetM44V4Product(vout, m, x);
336        GetV4Sum(vout, vout, v);
337    }
338    
339    } // anon namespace
340    
341    bool GetMxbInverse(float* mout, float* vout,
342                       const float* m_, const float* v_)
343    {
344        float m[16];
345        float v[4];
346        memcpy(m, m_, 16*sizeof(float));
347        memcpy(v, v_, 4*sizeof(float));
348
349        if(!GetM44Inverse(mout, m)) return false;
350        
351        for(int i=0; i<4; ++i)
352        {
353            v[i] = -v[i];
354        }
355        GetM44V4Product(vout, mout, v);
356        
357        return true;
358    }
359    
360}
361
362OCIO_NAMESPACE_EXIT
363
364
365
366
367///////////////////////////////////////////////////////////////////////////////
368
369#ifdef OCIO_UNIT_TEST
370
371OCIO_NAMESPACE_USING
372
373#include "UnitTest.h"
374
375OIIO_ADD_TEST(MathUtils, M44_is_diagonal)
376{
377    {
378        float m44[] = { 1.0f, 0.0f, 0.0f, 0.0f,
379                        0.0f, 1.0f, 0.0f, 0.0f,
380                        0.0f, 0.0f, 1.0f, 0.0f,
381                        0.0f, 0.0f, 0.0f, 1.0f };
382        bool isdiag = IsM44Diagonal(m44);
383        OIIO_CHECK_EQUAL(isdiag, true);
384
385        m44[1] += 1e-8f;
386        isdiag = IsM44Diagonal(m44);
387        OIIO_CHECK_EQUAL(isdiag, false);
388    }
389}
390
391
392OIIO_ADD_TEST(MathUtils, IsScalarEqualToZero)
393{
394        OIIO_CHECK_EQUAL(IsScalarEqualToZero(0.0f), true);
395        OIIO_CHECK_EQUAL(IsScalarEqualToZero(-0.0f), true);
396        
397        OIIO_CHECK_EQUAL(IsScalarEqualToZero(-1.072883670794056e-09f), false);
398        OIIO_CHECK_EQUAL(IsScalarEqualToZero(1.072883670794056e-09f), false);
399        
400        OIIO_CHECK_EQUAL(IsScalarEqualToZero(-1.072883670794056e-03f), false);
401        OIIO_CHECK_EQUAL(IsScalarEqualToZero(1.072883670794056e-03f), false);
402        
403        OIIO_CHECK_EQUAL(IsScalarEqualToZero(-1.072883670794056e-01f), false);
404        OIIO_CHECK_EQUAL(IsScalarEqualToZero(1.072883670794056e-01f), false);
405}
406
407OIIO_ADD_TEST(MathUtils, GetM44Inverse)
408{
409        // This is a degenerate matrix, and shouldnt be invertible.
410        float m[] = { 0.3f, 0.3f, 0.3f, 0.0f,
411                      0.3f, 0.3f, 0.3f, 0.0f,
412                      0.3f, 0.3f, 0.3f, 0.0f,
413                      0.0f, 0.0f, 0.0f, 1.0f };
414        
415        float mout[16];
416        bool invertsuccess = GetM44Inverse(mout, m);
417        OIIO_CHECK_EQUAL(invertsuccess, false);
418}
419
420
421OIIO_ADD_TEST(MathUtils, M44_M44_product)
422{
423    {
424        float mout[16];
425        float m1[] = { 1.0f, 2.0f, 0.0f, 0.0f,
426                       0.0f, 1.0f, 1.0f, 0.0f,
427                       1.0f, 0.0f, 1.0f, 0.0f,
428                       0.0f, 1.0f, 3.0f, 1.0f };
429        float m2[] = { 1.0f, 1.0f, 0.0f, 0.0f,
430                       0.0f, 1.0f, 0.0f, 0.0f,
431                       0.0f, 0.0f, 1.0f, 0.0f,
432                       2.0f, 0.0f, 0.0f, 1.0f };
433        GetM44M44Product(mout, m1, m2);
434        
435        float mcorrect[] = { 1.0f, 3.0f, 0.0f, 0.0f,
436                       0.0f, 1.0f, 1.0f, 0.0f,
437                       1.0f, 1.0f, 1.0f, 0.0f,
438                       2.0f, 1.0f, 3.0f, 1.0f };
439        
440        for(int i=0; i<16; ++i)
441        {
442            OIIO_CHECK_EQUAL(mout[i], mcorrect[i]);
443        }
444    }
445}
446
447OIIO_ADD_TEST(MathUtils, M44_V4_product)
448{
449    {
450        float vout[4];
451        float m[] = { 1.0f, 2.0f, 0.0f, 0.0f,
452                      0.0f, 1.0f, 1.0f, 0.0f,
453                      1.0f, 0.0f, 1.0f, 0.0f,
454                      0.0f, 1.0f, 3.0f, 1.0f };
455        float v[] = { 1.0f, 2.0f, 3.0f, 4.0f };
456        GetM44V4Product(vout, m, v);
457        
458        float vcorrect[] = { 5.0f, 5.0f, 4.0f, 15.0f };
459        
460        for(int i=0; i<4; ++i)
461        {
462            OIIO_CHECK_EQUAL(vout[i], vcorrect[i]);
463        }
464    }
465}
466
467OIIO_ADD_TEST(MathUtils, V4_add)
468{
469    {
470        float vout[4];
471        float v1[] = { 1.0f, 2.0f, 3.0f, 4.0f };
472        float v2[] = { 3.0f, 1.0f, 4.0f, 1.0f };
473        GetV4Sum(vout, v1, v2);
474        
475        float vcorrect[] = { 4.0f, 3.0f, 7.0f, 5.0f };
476        
477        for(int i=0; i<4; ++i)
478        {
479            OIIO_CHECK_EQUAL(vout[i], vcorrect[i]);
480        }
481    }
482}
483
484OIIO_ADD_TEST(MathUtils, mxb_eval)
485{
486    {
487        float vout[4];
488        float m[] = { 1.0f, 2.0f, 0.0f, 0.0f,
489                      0.0f, 1.0f, 1.0f, 0.0f,
490                      1.0f, 0.0f, 1.0f, 0.0f,
491                      0.0f, 1.0f, 3.0f, 1.0f };
492        float x[] = { 1.0f, 1.0f, 1.0f, 1.0f };
493        float v[] = { 1.0f, 2.0f, 3.0f, 4.0f };
494        GetMxbResult(vout, m, x, v);
495        
496        float vcorrect[] = { 4.0f, 4.0f, 5.0f, 9.0f };
497        
498        for(int i=0; i<4; ++i)
499        {
500            OIIO_CHECK_EQUAL(vout[i], vcorrect[i]);
501        }
502    }
503}
504
505OIIO_ADD_TEST(MathUtils, Combine_two_mxb)
506{
507    float m1[] = { 1.0f, 0.0f, 2.0f, 0.0f,
508                   2.0f, 1.0f, 0.0f, 1.0f,
509                   0.0f, 1.0f, 2.0f, 0.0f,
510                   1.0f, 0.0f, 0.0f, 1.0f };
511    float v1[] = { 1.0f, 2.0f, 3.0f, 4.0f };
512    float m2[] = { 2.0f, 1.0f, 0.0f, 0.0f,
513                   0.0f, 1.0f, 0.0f, 0.0f,
514                   1.0f, 0.0f, 3.0f, 0.0f,
515                   1.0f,1.0f, 1.0f, 1.0f };
516    float v2[] = { 0.0f, 2.0f, 1.0f, 0.0f };
517    float tolerance = 1e-9f;
518
519    {
520        float x[] = { 1.0f, 1.0f, 1.0f, 1.0f };
521        float vout[4];
522
523        // Combine two mx+b operations, and apply to test point
524        float mout[16];
525        float vcombined[4];      
526        GetMxbCombine(mout, vout, m1, v1, m2, v2);
527        GetMxbResult(vcombined, mout, x, vout);
528        
529        // Sequentially apply the two mx+b operations.
530        GetMxbResult(vout, m1, x, v1);
531        GetMxbResult(vout, m2, vout, v2);
532        
533        // Compare outputs
534        for(int i=0; i<4; ++i)
535        {
536            OIIO_CHECK_CLOSE(vcombined[i], vout[i], tolerance);
537        }
538    }
539    
540    {
541        float x[] = { 6.0f, 0.5f, -2.0f, -0.1f };
542        float vout[4];
543
544        float mout[16];
545        float vcombined[4];
546        GetMxbCombine(mout, vout, m1, v1, m2, v2);
547        GetMxbResult(vcombined, mout, x, vout);
548        
549        GetMxbResult(vout, m1, x, v1);
550        GetMxbResult(vout, m2, vout, v2);
551        
552        for(int i=0; i<4; ++i)
553        {
554            OIIO_CHECK_CLOSE(vcombined[i], vout[i], tolerance);
555        }
556    }
557    
558    {
559        float x[] = { 26.0f, -0.5f, 0.005f, 12.1f };
560        float vout[4];
561
562        float mout[16];
563        float vcombined[4];
564        GetMxbCombine(mout, vout, m1, v1, m2, v2);
565        GetMxbResult(vcombined, mout, x, vout);
566        
567        GetMxbResult(vout, m1, x, v1);
568        GetMxbResult(vout, m2, vout, v2);
569        
570        // We pick a not so small tolerance, as we're dealing with
571        // large numbers, and the error for CHECK_CLOSE is absolute.
572        for(int i=0; i<4; ++i)
573        {
574            OIIO_CHECK_CLOSE(vcombined[i], vout[i], 1e-3);
575        }
576    }
577}
578
579OIIO_ADD_TEST(MathUtils, mxb_invert)
580{
581    {
582        float m[] = { 1.0f, 2.0f, 0.0f, 0.0f,
583                      0.0f, 1.0f, 1.0f, 0.0f,
584                      1.0f, 0.0f, 1.0f, 0.0f,
585                      0.0f, 1.0f, 3.0f, 1.0f };
586        float x[] = { 1.0f, 0.5f, -1.0f, 60.0f };
587        float v[] = { 1.0f, 2.0f, 3.0f, 4.0f };
588        
589        float vresult[4];
590        float mout[16];
591        float vout[4];
592        
593        GetMxbResult(vresult, m, x, v);
594        bool invertsuccess = GetMxbInverse(mout, vout, m, v);
595        OIIO_CHECK_EQUAL(invertsuccess, true);
596        
597        GetMxbResult(vresult, mout, vresult, vout);
598        
599        float tolerance = 1e-9f;
600        for(int i=0; i<4; ++i)
601        {
602            OIIO_CHECK_CLOSE(vresult[i], x[i], tolerance);
603        }
604    }
605    
606    {
607        float m[] = { 0.3f, 0.3f, 0.3f, 0.0f,
608                      0.3f, 0.3f, 0.3f, 0.0f,
609                      0.3f, 0.3f, 0.3f, 0.0f,
610                      0.0f, 0.0f, 0.0f, 1.0f };
611        float v[] = { 0.0f, 0.0f, 0.0f, 0.0f };
612        
613        float mout[16];
614        float vout[4];
615        
616        bool invertsuccess = GetMxbInverse(mout, vout, m, v);
617        OIIO_CHECK_EQUAL(invertsuccess, false);
618    }
619}
620
621#endif
622