PageRenderTime 109ms CodeModel.GetById 10ms app.highlight 90ms RepoModel.GetById 1ms app.codeStats 1ms

/src/core/Lut3DOp.cpp

http://github.com/imageworks/OpenColorIO
C++ | 982 lines | 708 code | 146 blank | 128 comment | 61 complexity | 07ed781d1b994b1f5987ab679f606137 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 "HashUtils.h"
 32#include "Lut3DOp.h"
 33#include "MathUtils.h"
 34
 35#include <cmath>
 36#include <limits>
 37#include <sstream>
 38
 39OCIO_NAMESPACE_ENTER
 40{
 41    Lut3D::Lut3D()
 42    {
 43        for(int i=0; i<3; ++i)
 44        {
 45            from_min[i] = 0.0f;
 46            from_max[i] = 1.0f;
 47            size[i] = 0;
 48        }
 49    }
 50    
 51    Lut3DRcPtr Lut3D::Create()
 52    {
 53        return Lut3DRcPtr(new Lut3D());
 54    }
 55
 56    std::string Lut3D::getCacheID() const
 57    {
 58        AutoMutex lock(m_cacheidMutex);
 59        
 60        if(lut.empty())
 61            throw Exception("Cannot compute cacheID of invalid Lut3D");
 62        
 63        if(!m_cacheID.empty())
 64            return m_cacheID;
 65        
 66        md5_state_t state;
 67        md5_byte_t digest[16];
 68        
 69        md5_init(&state);
 70        
 71        md5_append(&state, (const md5_byte_t *)from_min, 3*sizeof(float));
 72        md5_append(&state, (const md5_byte_t *)from_max, 3*sizeof(float));
 73        md5_append(&state, (const md5_byte_t *)size, 3*sizeof(int));
 74        md5_append(&state, (const md5_byte_t *)&lut[0], (int) (lut.size()*sizeof(float)));
 75        
 76        md5_finish(&state, digest);
 77        
 78        m_cacheID = GetPrintableHash(digest);
 79        
 80        return m_cacheID;
 81    }
 82    
 83    
 84    namespace
 85    {
 86        // Linear
 87        inline float lerp(float a, float b, float z)
 88            { return (b - a) * z + a; }
 89        
 90        inline void lerp_rgb(float* out, float* a, float* b, float* z)
 91        {
 92            out[0] = (b[0] - a[0]) * z[0] + a[0];
 93            out[1] = (b[1] - a[1]) * z[1] + a[1];
 94            out[2] = (b[2] - a[2]) * z[2] + a[2];
 95        }
 96        
 97        // Bilinear
 98        inline float lerp(float a, float b, float c, float d, float y, float z)
 99            { return lerp(lerp(a, b, z), lerp(c, d, z), y); }
100        
101        inline void lerp_rgb(float* out, float* a, float* b, float* c,
102                             float* d, float* y, float* z)
103        {
104            float v1[3];
105            float v2[3];
106            
107            lerp_rgb(v1, a, b, z);
108            lerp_rgb(v2, c, d, z);
109            lerp_rgb(out, v1, v2, y);
110        }
111        
112        // Trilinear
113        inline float lerp(float a, float b, float c, float d,
114                          float e, float f, float g, float h,
115                          float x, float y, float z)
116            { return lerp(lerp(a,b,c,d,y,z), lerp(e,f,g,h,y,z), x); }
117        
118        inline void lerp_rgb(float* out, float* a, float* b, float* c, float* d,
119                             float* e, float* f, float* g, float* h,
120                             float* x, float* y, float* z)
121        {
122            float v1[3];
123            float v2[3];
124            
125            lerp_rgb(v1, a,b,c,d,y,z);
126            lerp_rgb(v2, e,f,g,h,y,z);
127            lerp_rgb(out, v1, v2, x);
128        }
129        
130        inline float lookupNearest_3D(int rIndex, int gIndex, int bIndex,
131                                      int size_red, int size_green, int size_blue,
132                                      const float* simple_rgb_lut, int channelIndex)
133        {
134            return simple_rgb_lut[ GetLut3DIndex_B(rIndex, gIndex, bIndex,
135                size_red, size_green, size_blue) + channelIndex];
136        }
137        
138        inline void lookupNearest_3D_rgb(float* rgb,
139                                         int rIndex, int gIndex, int bIndex,
140                                         int size_red, int size_green, int size_blue,
141                                         const float* simple_rgb_lut)
142        {
143            int offset = GetLut3DIndex_B(rIndex, gIndex, bIndex, size_red, size_green, size_blue);
144            rgb[0] = simple_rgb_lut[offset];
145            rgb[1] = simple_rgb_lut[offset+1];
146            rgb[2] = simple_rgb_lut[offset+2];
147        }
148        
149        // Note: This function assumes that minVal is less than maxVal
150        inline int clamp(float k, float minVal, float maxVal)
151        {
152            return static_cast<int>(roundf(std::max(std::min(k, maxVal), minVal)));
153        }
154        
155        ///////////////////////////////////////////////////////////////////////
156        // Nearest Forward
157        
158        void Lut3D_Nearest(float* rgbaBuffer, long numPixels, const Lut3D & lut)
159        {
160            float maxIndex[3];
161            float mInv[3];
162            float b[3];
163            float mInv_x_maxIndex[3];
164            int lutSize[3];
165            const float* startPos = &(lut.lut[0]);
166            
167            for(int i=0; i<3; ++i)
168            {
169                maxIndex[i] = (float) (lut.size[i] - 1);
170                mInv[i] = 1.0f / (lut.from_max[i] - lut.from_min[i]);
171                b[i] = lut.from_min[i];
172                mInv_x_maxIndex[i] = (float) (mInv[i] * maxIndex[i]);
173                
174                lutSize[i] = lut.size[i];
175            }
176            
177            int localIndex[3];
178            
179            for(long pixelIndex=0; pixelIndex<numPixels; ++pixelIndex)
180            {
181                if(isnan(rgbaBuffer[0]) || isnan(rgbaBuffer[1]) || isnan(rgbaBuffer[2]))
182                {
183                    rgbaBuffer[0] = std::numeric_limits<float>::quiet_NaN();
184                    rgbaBuffer[1] = std::numeric_limits<float>::quiet_NaN();
185                    rgbaBuffer[2] = std::numeric_limits<float>::quiet_NaN();
186                }
187                else
188                {
189                    localIndex[0] = clamp(mInv_x_maxIndex[0] * (rgbaBuffer[0] - b[0]), 0.0f, maxIndex[0]);
190                    localIndex[1] = clamp(mInv_x_maxIndex[1] * (rgbaBuffer[1] - b[1]), 0.0f, maxIndex[1]);
191                    localIndex[2] = clamp(mInv_x_maxIndex[2] * (rgbaBuffer[2] - b[2]), 0.0f, maxIndex[2]);
192                    
193                    lookupNearest_3D_rgb(rgbaBuffer, localIndex[0], localIndex[1], localIndex[2],
194                                                     lutSize[0], lutSize[1], lutSize[2], startPos);
195                }
196                
197                rgbaBuffer += 4;
198            }
199        }
200        
201        ///////////////////////////////////////////////////////////////////////
202        // Linear Forward
203        
204        void Lut3D_Linear(float* rgbaBuffer, long numPixels, const Lut3D & lut)
205        {
206            float maxIndex[3];
207            float mInv[3];
208            float b[3];
209            float mInv_x_maxIndex[3];
210            int lutSize[3];
211            const float* startPos = &(lut.lut[0]);
212            
213            for(int i=0; i<3; ++i)
214            {
215                maxIndex[i] = (float) (lut.size[i] - 1);
216                mInv[i] = 1.0f / (lut.from_max[i] - lut.from_min[i]);
217                b[i] = lut.from_min[i];
218                mInv_x_maxIndex[i] = (float) (mInv[i] * maxIndex[i]);
219                
220                lutSize[i] = lut.size[i];
221            }
222            
223            for(long pixelIndex=0; pixelIndex<numPixels; ++pixelIndex)
224            {
225                
226                if(isnan(rgbaBuffer[0]) || isnan(rgbaBuffer[1]) || isnan(rgbaBuffer[2]))
227                {
228                    rgbaBuffer[0] = std::numeric_limits<float>::quiet_NaN();
229                    rgbaBuffer[1] = std::numeric_limits<float>::quiet_NaN();
230                    rgbaBuffer[2] = std::numeric_limits<float>::quiet_NaN();
231                }
232                else
233                {
234                    float localIndex[3];
235                    int indexLow[3];
236                    int indexHigh[3];
237                    float delta[3];
238                    float a[3];
239                    float b_[3];
240                    float c[3];
241                    float d[3];
242                    float e[3];
243                    float f[3];
244                    float g[3];
245                    float h[4];
246                    float x[4];
247                    float y[4];
248                    float z[4];
249                    
250                    localIndex[0] = std::max(std::min(mInv_x_maxIndex[0] * (rgbaBuffer[0] - b[0]), maxIndex[0]), 0.0f);
251                    localIndex[1] = std::max(std::min(mInv_x_maxIndex[1] * (rgbaBuffer[1] - b[1]), maxIndex[1]), 0.0f);
252                    localIndex[2] = std::max(std::min(mInv_x_maxIndex[2] * (rgbaBuffer[2] - b[2]), maxIndex[2]), 0.0f);
253                    
254                    indexLow[0] =  static_cast<int>(std::floor(localIndex[0]));
255                    indexLow[1] =  static_cast<int>(std::floor(localIndex[1]));
256                    indexLow[2] =  static_cast<int>(std::floor(localIndex[2]));
257                    
258                    indexHigh[0] =  static_cast<int>(std::ceil(localIndex[0]));
259                    indexHigh[1] =  static_cast<int>(std::ceil(localIndex[1]));
260                    indexHigh[2] =  static_cast<int>(std::ceil(localIndex[2]));
261                    
262                    delta[0] = localIndex[0] - static_cast<float>(indexLow[0]);
263                    delta[1] = localIndex[1] - static_cast<float>(indexLow[1]);
264                    delta[2] = localIndex[2] - static_cast<float>(indexLow[2]);
265                    
266                    // Lookup 8 corners of cube
267                    lookupNearest_3D_rgb(a, indexLow[0],  indexLow[1],  indexLow[2],  lutSize[0], lutSize[1], lutSize[2], startPos);
268                    lookupNearest_3D_rgb(b_, indexLow[0],  indexLow[1],  indexHigh[2], lutSize[0], lutSize[1], lutSize[2], startPos);
269                    lookupNearest_3D_rgb(c, indexLow[0],  indexHigh[1], indexLow[2],  lutSize[0], lutSize[1], lutSize[2], startPos);
270                    lookupNearest_3D_rgb(d, indexLow[0],  indexHigh[1], indexHigh[2], lutSize[0], lutSize[1], lutSize[2], startPos);
271                    lookupNearest_3D_rgb(e, indexHigh[0], indexLow[1],  indexLow[2],  lutSize[0], lutSize[1], lutSize[2], startPos);
272                    lookupNearest_3D_rgb(f, indexHigh[0], indexLow[1],  indexHigh[2], lutSize[0], lutSize[1], lutSize[2], startPos);
273                    lookupNearest_3D_rgb(g, indexHigh[0], indexHigh[1], indexLow[2],  lutSize[0], lutSize[1], lutSize[2], startPos);
274                    lookupNearest_3D_rgb(h, indexHigh[0], indexHigh[1], indexHigh[2], lutSize[0], lutSize[1], lutSize[2], startPos);
275                    
276                    // Also store the 3d interpolation coordinates
277                    x[0] = delta[0]; x[1] = delta[0]; x[2] = delta[0];
278                    y[0] = delta[1]; y[1] = delta[1]; y[2] = delta[1];
279                    z[0] = delta[2]; z[1] = delta[2]; z[2] = delta[2];
280                    
281                    // Do a trilinear interpolation of the 8 corners
282                    // 4726.8 scanlines/sec
283                    
284                    lerp_rgb(rgbaBuffer, a, b_, c, d, e, f, g, h,
285                                         x, y, z);
286                }
287                
288                rgbaBuffer += 4;
289            }
290        }
291    }
292    
293    
294    void Lut3D_Tetrahedral(float* rgbaBuffer, long numPixels, const Lut3D & lut)
295    {
296        // Tetrahedral interoplation, as described by:
297        // http://www.filmlight.ltd.uk/pdf/whitepapers/FL-TL-TN-0057-SoftwareLib.pdf
298        // http://blogs.mathworks.com/steve/2006/11/24/tetrahedral-interpolation-for-colorspace-conversion/
299        // http://www.hpl.hp.com/techreports/98/HPL-98-95.html
300
301        float maxIndex[3];
302        float mInv[3];
303        float b[3];
304        float mInv_x_maxIndex[3];
305        int lutSize[3];
306        const float* startPos = &(lut.lut[0]);
307
308        for(int i=0; i<3; ++i)
309        {
310            maxIndex[i] = (float) (lut.size[i] - 1);
311            mInv[i] = 1.0f / (lut.from_max[i] - lut.from_min[i]);
312            b[i] = lut.from_min[i];
313            mInv_x_maxIndex[i] = (float) (mInv[i] * maxIndex[i]);
314
315            lutSize[i] = lut.size[i];
316        }
317
318        for(long pixelIndex=0; pixelIndex<numPixels; ++pixelIndex)
319        {
320
321            if(isnan(rgbaBuffer[0]) || isnan(rgbaBuffer[1]) || isnan(rgbaBuffer[2]))
322            {
323                rgbaBuffer[0] = std::numeric_limits<float>::quiet_NaN();
324                rgbaBuffer[1] = std::numeric_limits<float>::quiet_NaN();
325                rgbaBuffer[2] = std::numeric_limits<float>::quiet_NaN();
326            }
327            else
328            {
329                float localIndex[3];
330                int indexLow[3];
331                int indexHigh[3];
332                float delta[3];
333
334                // Same index/delta calculation as linear interpolation
335                localIndex[0] = std::max(std::min(mInv_x_maxIndex[0] * (rgbaBuffer[0] - b[0]), maxIndex[0]), 0.0f);
336                localIndex[1] = std::max(std::min(mInv_x_maxIndex[1] * (rgbaBuffer[1] - b[1]), maxIndex[1]), 0.0f);
337                localIndex[2] = std::max(std::min(mInv_x_maxIndex[2] * (rgbaBuffer[2] - b[2]), maxIndex[2]), 0.0f);
338
339                indexLow[0] =  static_cast<int>(std::floor(localIndex[0]));
340                indexLow[1] =  static_cast<int>(std::floor(localIndex[1]));
341                indexLow[2] =  static_cast<int>(std::floor(localIndex[2]));
342
343                indexHigh[0] =  static_cast<int>(std::ceil(localIndex[0]));
344                indexHigh[1] =  static_cast<int>(std::ceil(localIndex[1]));
345                indexHigh[2] =  static_cast<int>(std::ceil(localIndex[2]));
346
347                delta[0] = localIndex[0] - static_cast<float>(indexLow[0]);
348                delta[1] = localIndex[1] - static_cast<float>(indexLow[1]);
349                delta[2] = localIndex[2] - static_cast<float>(indexLow[2]);
350
351                // Rebind for consistency with Truelight paper
352                float fx = delta[0];
353                float fy = delta[1];
354                float fz = delta[2];
355
356                // Compute index into LUT for surrounding corners
357                const int n000 = GetLut3DIndex_B(indexLow[0], indexLow[1], indexLow[2],
358                                                 lutSize[0], lutSize[1], lutSize[2]);
359                const int n100 = GetLut3DIndex_B(indexHigh[0], indexLow[1], indexLow[2],
360                                                 lutSize[0], lutSize[1], lutSize[2]);
361                const int n010 = GetLut3DIndex_B(indexLow[0], indexHigh[1], indexLow[2],
362                                                 lutSize[0], lutSize[1], lutSize[2]);
363                const int n001 = GetLut3DIndex_B(indexLow[0], indexLow[1], indexHigh[2],
364                                                 lutSize[0], lutSize[1], lutSize[2]);
365                const int n110 = GetLut3DIndex_B(indexHigh[0], indexHigh[1], indexLow[2],
366                                                 lutSize[0], lutSize[1], lutSize[2]);
367                const int n101 = GetLut3DIndex_B(indexHigh[0], indexLow[1], indexHigh[2],
368                                                 lutSize[0], lutSize[1], lutSize[2]);
369                const int n011 = GetLut3DIndex_B(indexLow[0], indexHigh[1], indexHigh[2],
370                                                 lutSize[0], lutSize[1], lutSize[2]);
371                const int n111 = GetLut3DIndex_B(indexHigh[0], indexHigh[1], indexHigh[2],
372                                                 lutSize[0], lutSize[1], lutSize[2]);
373
374                if (fx > fy) {
375                    if (fy > fz) {
376                       rgbaBuffer[0] =
377                           (1-fx)  * startPos[n000] +
378                           (fx-fy) * startPos[n100] +
379                           (fy-fz) * startPos[n110] +
380                           (fz)    * startPos[n111];
381
382                       rgbaBuffer[1] =
383                           (1-fx)  * startPos[n000+1] +
384                           (fx-fy) * startPos[n100+1] +
385                           (fy-fz) * startPos[n110+1] +
386                           (fz)    * startPos[n111+1];
387
388                       rgbaBuffer[2] =
389                           (1-fx)  * startPos[n000+2] +
390                           (fx-fy) * startPos[n100+2] +
391                           (fy-fz) * startPos[n110+2] +
392                           (fz)    * startPos[n111+2];
393                    }
394                    else if (fx > fz)
395                    {
396                        rgbaBuffer[0] =
397                            (1-fx)  * startPos[n000] +
398                            (fx-fz) * startPos[n100] +
399                            (fz-fy) * startPos[n101] +
400                            (fy)    * startPos[n111];
401
402                        rgbaBuffer[1] =
403                            (1-fx)  * startPos[n000+1] +
404                            (fx-fz) * startPos[n100+1] +
405                            (fz-fy) * startPos[n101+1] +
406                            (fy)    * startPos[n111+1];
407
408                        rgbaBuffer[2] =
409                            (1-fx)  * startPos[n000+2] +
410                            (fx-fz) * startPos[n100+2] +
411                            (fz-fy) * startPos[n101+2] +
412                            (fy)    * startPos[n111+2];
413                    }
414                    else
415                    {
416                        rgbaBuffer[0] =
417                            (1-fz)  * startPos[n000] +
418                            (fz-fx) * startPos[n001] +
419                            (fx-fy) * startPos[n101] +
420                            (fy)    * startPos[n111];
421
422                        rgbaBuffer[1] =
423                            (1-fz)  * startPos[n000+1] +
424                            (fz-fx) * startPos[n001+1] +
425                            (fx-fy) * startPos[n101+1] +
426                            (fy)    * startPos[n111+1];
427
428                        rgbaBuffer[2] =
429                            (1-fz)  * startPos[n000+2] +
430                            (fz-fx) * startPos[n001+2] +
431                            (fx-fy) * startPos[n101+2] +
432                            (fy)    * startPos[n111+2];
433                    }
434                }
435                else
436                {
437                    if (fz > fy)
438                    {
439                        rgbaBuffer[0] =
440                            (1-fz)  * startPos[n000] +
441                            (fz-fy) * startPos[n001] +
442                            (fy-fx) * startPos[n011] +
443                            (fx)    * startPos[n111];
444
445                        rgbaBuffer[1] =
446                            (1-fz)  * startPos[n000+1] +
447                            (fz-fy) * startPos[n001+1] +
448                            (fy-fx) * startPos[n011+1] +
449                            (fx)    * startPos[n111+1];
450
451                        rgbaBuffer[2] =
452                            (1-fz)  * startPos[n000+2] +
453                            (fz-fy) * startPos[n001+2] +
454                            (fy-fx) * startPos[n011+2] +
455                            (fx)    * startPos[n111+2];
456                    }
457                    else if (fz > fx)
458                    {
459                        rgbaBuffer[0] =
460                            (1-fy)  * startPos[n000] +
461                            (fy-fz) * startPos[n010] +
462                            (fz-fx) * startPos[n011] +
463                            (fx)    * startPos[n111];
464
465                        rgbaBuffer[1] =
466                            (1-fy)  * startPos[n000+1] +
467                            (fy-fz) * startPos[n010+1] +
468                            (fz-fx) * startPos[n011+1] +
469                            (fx)    * startPos[n111+1];
470
471                        rgbaBuffer[2] =
472                            (1-fy)  * startPos[n000+2] +
473                            (fy-fz) * startPos[n010+2] +
474                            (fz-fx) * startPos[n011+2] +
475                            (fx)    * startPos[n111+2];
476                    }
477                    else
478                    {
479                        rgbaBuffer[0] =
480                            (1-fy)  * startPos[n000] +
481                            (fy-fx) * startPos[n010] +
482                            (fx-fz) * startPos[n110] +
483                            (fz)    * startPos[n111];
484
485                        rgbaBuffer[1] =
486                            (1-fy)  * startPos[n000+1] +
487                            (fy-fx) * startPos[n010+1] +
488                            (fx-fz) * startPos[n110+1] +
489                            (fz)    * startPos[n111+1];
490
491                        rgbaBuffer[2] =
492                            (1-fy)  * startPos[n000+2] +
493                            (fy-fx) * startPos[n010+2] +
494                            (fx-fz) * startPos[n110+2] +
495                            (fz)    * startPos[n111+2];
496                    }
497                }
498            } // !isnan
499            
500            rgbaBuffer += 4;
501        }
502    }
503
504
505    void GenerateIdentityLut3D(float* img, int edgeLen, int numChannels, Lut3DOrder lut3DOrder)
506    {
507        if(!img) return;
508        if(numChannels < 3)
509        {
510            throw Exception("Cannot generate idenitity 3d lut with less than 3 channels.");
511        }
512        
513        float c = 1.0f / ((float)edgeLen - 1.0f);
514        
515        if(lut3DOrder == LUT3DORDER_FAST_RED)
516        {
517            for(int i=0; i<edgeLen*edgeLen*edgeLen; i++)
518            {
519                img[numChannels*i+0] = (float)(i%edgeLen) * c;
520                img[numChannels*i+1] = (float)((i/edgeLen)%edgeLen) * c;
521                img[numChannels*i+2] = (float)((i/edgeLen/edgeLen)%edgeLen) * c;
522            }
523        }
524        else if(lut3DOrder == LUT3DORDER_FAST_BLUE)
525        {
526            for(int i=0; i<edgeLen*edgeLen*edgeLen; i++)
527            {
528                img[numChannels*i+0] = (float)((i/edgeLen/edgeLen)%edgeLen) * c;
529                img[numChannels*i+1] = (float)((i/edgeLen)%edgeLen) * c;
530                img[numChannels*i+2] = (float)(i%edgeLen) * c;
531            }
532        }
533        else
534        {
535            throw Exception("Unknown Lut3DOrder.");
536        }
537    }
538    
539    
540    int Get3DLutEdgeLenFromNumPixels(int numPixels)
541    {
542        int dim = static_cast<int>(roundf(powf((float) numPixels, 1.0f/3.0f)));
543        
544        if(dim*dim*dim != numPixels)
545        {
546            std::ostringstream os;
547            os << "Cannot infer 3D Lut size. ";
548            os << numPixels << " element(s) does not correspond to a ";
549            os << "unform cube edge length. (nearest edge length is ";
550            os << dim << ").";
551            throw Exception(os.str().c_str());
552        }
553        
554        return dim;
555    }
556    
557    namespace
558    {
559        class Lut3DOp : public Op
560        {
561        public:
562            Lut3DOp(Lut3DRcPtr lut,
563                    Interpolation interpolation,
564                    TransformDirection direction);
565            virtual ~Lut3DOp();
566            
567            virtual OpRcPtr clone() const;
568            
569            virtual std::string getInfo() const;
570            virtual std::string getCacheID() const;
571            
572            virtual bool isNoOp() const;
573            virtual bool isSameType(const OpRcPtr & op) const;
574            virtual bool isInverse(const OpRcPtr & op) const;
575            virtual bool hasChannelCrosstalk() const;
576            virtual void finalize();
577            virtual void apply(float* rgbaBuffer, long numPixels) const;
578            
579            virtual bool supportsGpuShader() const;
580            virtual void writeGpuShader(std::ostream & shader,
581                                        const std::string & pixelName,
582                                        const GpuShaderDesc & shaderDesc) const;
583            
584        private:
585            Lut3DRcPtr m_lut;
586            Interpolation m_interpolation;
587            TransformDirection m_direction;
588            
589            // Set in finalize
590            std::string m_cacheID;
591        };
592        
593        typedef OCIO_SHARED_PTR<Lut3DOp> Lut3DOpRcPtr;
594        
595        
596        Lut3DOp::Lut3DOp(Lut3DRcPtr lut,
597                         Interpolation interpolation,
598                         TransformDirection direction):
599                            Op(),
600                            m_lut(lut),
601                            m_interpolation(interpolation),
602                            m_direction(direction)
603        {
604        }
605        
606        OpRcPtr Lut3DOp::clone() const
607        {
608            OpRcPtr op = OpRcPtr(new Lut3DOp(m_lut, m_interpolation, m_direction));
609            return op;
610        }
611        
612        Lut3DOp::~Lut3DOp()
613        { }
614        
615        std::string Lut3DOp::getInfo() const
616        {
617            return "<Lut3DOp>";
618        }
619        
620        std::string Lut3DOp::getCacheID() const
621        {
622            return m_cacheID;
623        }
624        
625        // TODO: compute real value for isNoOp
626        bool Lut3DOp::isNoOp() const
627        {
628            return false;
629        }
630        
631        bool Lut3DOp::isSameType(const OpRcPtr & op) const
632        {
633            Lut3DOpRcPtr typedRcPtr = DynamicPtrCast<Lut3DOp>(op);
634            if(!typedRcPtr) return false;
635            return true;
636        }
637        
638        bool Lut3DOp::isInverse(const OpRcPtr & op) const
639        {
640            Lut3DOpRcPtr typedRcPtr = DynamicPtrCast<Lut3DOp>(op);
641            if(!typedRcPtr) return false;
642            
643            if(GetInverseTransformDirection(m_direction) != typedRcPtr->m_direction)
644                return false;
645            
646            return (m_lut->getCacheID() == typedRcPtr->m_lut->getCacheID());
647        }
648        
649        // TODO: compute real value for hasChannelCrosstalk
650        bool Lut3DOp::hasChannelCrosstalk() const
651        {
652            return true;
653        }
654        
655        void Lut3DOp::finalize()
656        {
657            if(m_direction != TRANSFORM_DIR_FORWARD)
658            {
659                std::ostringstream os;
660                os << "3D Luts can only be applied in the forward direction. ";
661                os << "(" << TransformDirectionToString(m_direction) << ")";
662                os << " specified.";
663                throw Exception(os.str().c_str());
664            }
665            
666            // Validate the requested interpolation type
667            switch(m_interpolation)
668            {
669                 // These are the allowed values.
670                case INTERP_NEAREST:
671                case INTERP_LINEAR:
672                case INTERP_TETRAHEDRAL:
673                    break;
674                case INTERP_BEST:
675                    m_interpolation = INTERP_LINEAR;
676                    break;
677                case INTERP_UNKNOWN:
678                    throw Exception("Cannot apply Lut3DOp, unspecified interpolation.");
679                    break;
680                default:
681                    throw Exception("Cannot apply Lut3DOp, invalid interpolation specified.");
682            }
683            
684            for(int i=0; i<3; ++i)
685            {
686                if(m_lut->size[i] == 0)
687                {
688                    throw Exception("Cannot apply Lut3DOp, lut object is empty.");
689                }
690                // TODO if from_min[i] == from_max[i]
691            }
692            
693            if(m_lut->size[0]*m_lut->size[1]*m_lut->size[2] * 3 != (int)m_lut->lut.size())
694            {
695                throw Exception("Cannot apply Lut3DOp, specified size does not match data.");
696            }
697            
698            // Create the cacheID
699            std::ostringstream cacheIDStream;
700            cacheIDStream << "<Lut3DOp ";
701            cacheIDStream << m_lut->getCacheID() << " ";
702            cacheIDStream << InterpolationToString(m_interpolation) << " ";
703            cacheIDStream << TransformDirectionToString(m_direction) << " ";
704            cacheIDStream << ">";
705            m_cacheID = cacheIDStream.str();
706        }
707        
708        void Lut3DOp::apply(float* rgbaBuffer, long numPixels) const
709        {
710            if(m_interpolation == INTERP_NEAREST)
711            {
712                Lut3D_Nearest(rgbaBuffer, numPixels, *m_lut);
713            }
714            else if(m_interpolation == INTERP_LINEAR)
715            {
716                Lut3D_Linear(rgbaBuffer, numPixels, *m_lut);
717            }
718            else if(m_interpolation == INTERP_TETRAHEDRAL)
719            {
720                Lut3D_Tetrahedral(rgbaBuffer, numPixels, *m_lut);
721            }
722        }
723        
724        bool Lut3DOp::supportsGpuShader() const
725        {
726            return false;
727        }
728        
729        void Lut3DOp::writeGpuShader(std::ostream & /*shader*/,
730                                     const std::string & /*pixelName*/,
731                                     const GpuShaderDesc & /*shaderDesc*/) const
732        {
733            throw Exception("Lut3DOp does not support analytical shader generation.");
734        }
735    }
736    
737    void CreateLut3DOp(OpRcPtrVec & ops,
738                       Lut3DRcPtr lut,
739                       Interpolation interpolation,
740                       TransformDirection direction)
741    {
742        ops.push_back( Lut3DOpRcPtr(new Lut3DOp(lut, interpolation, direction)) );
743    }
744}
745OCIO_NAMESPACE_EXIT
746
747///////////////////////////////////////////////////////////////////////////////
748
749#ifdef OCIO_UNIT_TEST
750
751#include <cstring>
752#include <cstdlib>
753#include <sys/time.h>
754
755namespace OCIO = OCIO_NAMESPACE;
756#include "UnitTest.h"
757
758OIIO_ADD_TEST(Lut3DOp, NanInfValueCheck)
759{
760    OCIO::Lut3DRcPtr lut = OCIO::Lut3D::Create();
761    
762    lut->from_min[0] = 0.0f;
763    lut->from_min[1] = 0.0f;
764    lut->from_min[2] = 0.0f;
765    
766    lut->from_max[0] = 1.0f;
767    lut->from_max[1] = 1.0f;
768    lut->from_max[2] = 1.0f;
769    
770    lut->size[0] = 3;
771    lut->size[1] = 3;
772    lut->size[2] = 3;
773    
774    lut->lut.resize(lut->size[0]*lut->size[1]*lut->size[2]*3);
775    
776    GenerateIdentityLut3D(&lut->lut[0], lut->size[0], 3, OCIO::LUT3DORDER_FAST_RED);
777    for(unsigned int i=0; i<lut->lut.size(); ++i)
778    {
779        lut->lut[i] = powf(lut->lut[i], 2.0f);
780    }
781    
782    const float reference[4] = {  std::numeric_limits<float>::signaling_NaN(),
783                                  std::numeric_limits<float>::quiet_NaN(),
784                                  std::numeric_limits<float>::infinity(),
785                                  -std::numeric_limits<float>::infinity() };
786    float color[4];
787    
788    memcpy(color, reference, 4*sizeof(float));
789    OCIO::Lut3D_Nearest(color, 1, *lut);
790    
791    memcpy(color, reference, 4*sizeof(float));
792    OCIO::Lut3D_Linear(color, 1, *lut);
793}
794
795
796OIIO_ADD_TEST(Lut3DOp, ValueCheck)
797{
798    OCIO::Lut3DRcPtr lut = OCIO::Lut3D::Create();
799    
800    lut->from_min[0] = 0.0f;
801    lut->from_min[1] = 0.0f;
802    lut->from_min[2] = 0.0f;
803    
804    lut->from_max[0] = 1.0f;
805    lut->from_max[1] = 1.0f;
806    lut->from_max[2] = 1.0f;
807    
808    lut->size[0] = 32;
809    lut->size[1] = 32;
810    lut->size[2] = 32;
811    
812    lut->lut.resize(lut->size[0]*lut->size[1]*lut->size[2]*3);
813    GenerateIdentityLut3D(&lut->lut[0], lut->size[0], 3, OCIO::LUT3DORDER_FAST_RED);
814    for(unsigned int i=0; i<lut->lut.size(); ++i)
815    {
816        lut->lut[i] = powf(lut->lut[i], 2.0f);
817    }
818    
819    const float reference[] = {  0.0f, 0.2f, 0.3f, 1.0f,
820                                 0.1234f, 0.4567f, 0.9876f, 1.0f,
821                                 11.0f, -0.5f, 0.5010f, 1.0f
822                                };
823    const float nearest[] = { 0.0f, 0.03746097535f, 0.0842871964f, 1.0f,
824                              0.01664932258f, 0.2039542049f, 1.0f, 1.0f,
825                              1.0f, 0.0f, 0.2663891613f, 1.0f
826                            };
827    const float linear[] = { 0.0f, 0.04016649351f, 0.09021852165f, 1.0f,
828                              0.01537752338f, 0.2087130845f, 0.9756000042f, 1.0f,
829                              1.0f, 0.0f, 0.2512601018f, 1.0f
830                            };
831    float color[12];
832    
833    // Check nearest
834    memcpy(color, reference, 12*sizeof(float));
835    OCIO::Lut3D_Nearest(color, 3, *lut);
836    for(unsigned int i=0; i<12; ++i)
837    {
838        OIIO_CHECK_CLOSE(color[i], nearest[i], 1e-8);
839    }
840    
841    // Check linear
842    memcpy(color, reference, 12*sizeof(float));
843    OCIO::Lut3D_Linear(color, 3, *lut);
844    for(unsigned int i=0; i<12; ++i)
845    {
846        OIIO_CHECK_CLOSE(color[i], linear[i], 1e-8);
847    }
848
849    // Check tetrahedral
850    memcpy(color, reference, 12*sizeof(float));
851    OCIO::Lut3D_Tetrahedral(color, 3, *lut);
852    for(unsigned int i=0; i<12; ++i)
853    {
854        OIIO_CHECK_CLOSE(color[i], linear[i], 1e-7); // Note, max delta lowered from 1e-8
855    }
856}
857
858
859
860OIIO_ADD_TEST(Lut3DOp, InverseComparisonCheck)
861{
862    OCIO::Lut3DRcPtr lut_a = OCIO::Lut3D::Create();
863    lut_a->from_min[0] = 0.0f;
864    lut_a->from_min[1] = 0.0f;
865    lut_a->from_min[2] = 0.0f;
866    lut_a->from_max[0] = 1.0f;
867    lut_a->from_max[1] = 1.0f;
868    lut_a->from_max[2] = 1.0f;
869    lut_a->size[0] = 32;
870    lut_a->size[1] = 32;
871    lut_a->size[2] = 32;
872    lut_a->lut.resize(lut_a->size[0]*lut_a->size[1]*lut_a->size[2]*3);
873    GenerateIdentityLut3D(&lut_a->lut[0], lut_a->size[0], 3, OCIO::LUT3DORDER_FAST_RED);
874    
875    OCIO::Lut3DRcPtr lut_b = OCIO::Lut3D::Create();
876    lut_b->from_min[0] = 0.5f;
877    lut_b->from_min[1] = 0.5f;
878    lut_b->from_min[2] = 0.5f;
879    lut_b->from_max[0] = 1.0f;
880    lut_b->from_max[1] = 1.0f;
881    lut_b->from_max[2] = 1.0f;
882    lut_b->size[0] = 32;
883    lut_b->size[1] = 32;
884    lut_b->size[2] = 32;
885    lut_b->lut.resize(lut_b->size[0]*lut_b->size[1]*lut_b->size[2]*3);
886    GenerateIdentityLut3D(&lut_b->lut[0], lut_b->size[0], 3, OCIO::LUT3DORDER_FAST_RED);
887    
888    OCIO::OpRcPtrVec ops;
889    CreateLut3DOp(ops, lut_a, OCIO::INTERP_NEAREST, OCIO::TRANSFORM_DIR_FORWARD);
890    CreateLut3DOp(ops, lut_a, OCIO::INTERP_LINEAR, OCIO::TRANSFORM_DIR_INVERSE);
891    CreateLut3DOp(ops, lut_b, OCIO::INTERP_LINEAR, OCIO::TRANSFORM_DIR_FORWARD);
892    CreateLut3DOp(ops, lut_b, OCIO::INTERP_LINEAR, OCIO::TRANSFORM_DIR_INVERSE);
893    
894    OIIO_CHECK_EQUAL(ops.size(), 4);
895    
896    OIIO_CHECK_ASSERT(ops[0]->isSameType(ops[1]));
897    OIIO_CHECK_ASSERT(ops[0]->isSameType(ops[2]));
898    OIIO_CHECK_ASSERT(ops[0]->isSameType(ops[3]->clone()));
899    
900    OIIO_CHECK_EQUAL( ops[0]->isInverse(ops[1]), true);
901    OIIO_CHECK_EQUAL( ops[0]->isInverse(ops[2]), false);
902    OIIO_CHECK_EQUAL( ops[0]->isInverse(ops[2]), false);
903    OIIO_CHECK_EQUAL( ops[0]->isInverse(ops[3]), false);
904    OIIO_CHECK_EQUAL( ops[2]->isInverse(ops[3]), true);
905}
906
907
908OIIO_ADD_TEST(Lut3DOp, PerformanceCheck)
909{
910    /*
911    OCIO::Lut3D lut;
912    
913    lut.from_min[0] = 0.0f;
914    lut.from_min[1] = 0.0f;
915    lut.from_min[2] = 0.0f;
916    
917    lut.from_max[0] = 1.0f;
918    lut.from_max[1] = 1.0f;
919    lut.from_max[2] = 1.0f;
920    
921    lut.size[0] = 32;
922    lut.size[1] = 32;
923    lut.size[2] = 32;
924    
925    lut.lut.resize(lut.size[0]*lut.size[1]*lut.size[2]*3);
926    GenerateIdentityLut3D(&lut.lut[0], lut.size[0], 3, OCIO::LUT3DORDER_FAST_RED);
927    
928    std::vector<float> img;
929    int xres = 2048;
930    int yres = 1;
931    int channels = 4;
932    img.resize(xres*yres*channels);
933    
934    srand48(0);
935    
936    // create random values from -0.05 to 1.05
937    // (To simulate clipping performance)
938    
939    for(unsigned int i=0; i<img.size(); ++i)
940    {
941        float uniform = (float)drand48();
942        img[i] = uniform*1.1f - 0.05f;
943    }
944    
945    timeval t;
946    gettimeofday(&t, 0);
947    double starttime = (double) t.tv_sec + (double) t.tv_usec / 1000000.0;
948    
949    int numloops = 1024;
950    for(int i=0; i<numloops; ++i)
951    {
952        //OCIO::Lut3D_Nearest(&img[0], xres*yres, lut);
953        OCIO::Lut3D_Linear(&img[0], xres*yres, lut);
954    }
955    
956    gettimeofday(&t, 0);
957    double endtime = (double) t.tv_sec + (double) t.tv_usec / 1000000.0;
958    double totaltime_a = (endtime-starttime)/numloops;
959    
960    printf("Linear: %0.1f ms  - %0.1f fps\n", totaltime_a*1000.0, 1.0/totaltime_a);
961
962
963    // Tetrahedral
964    gettimeofday(&t, 0);
965    starttime = (double) t.tv_sec + (double) t.tv_usec / 1000000.0;
966
967    for(int i=0; i<numloops; ++i)
968    {
969        OCIO::Lut3D_Tetrahedral(&img[0], xres*yres, lut);
970    }
971
972    gettimeofday(&t, 0);
973    endtime = (double) t.tv_sec + (double) t.tv_usec / 1000000.0;
974    double totaltime_b = (endtime-starttime)/numloops;
975
976    printf("Tetra: %0.1f ms  - %0.1f fps\n", totaltime_b*1000.0, 1.0/totaltime_b);
977
978    double speed_diff = totaltime_a/totaltime_b;
979    printf("Tetra is %.04f speed of Linear\n", speed_diff);
980    */
981}
982#endif // OCIO_UNIT_TEST