Anonymous avatar Anonymous committed 129dc03

lala

Comments (0)

Files changed (6)

project2/photomerge/source-1/feature.cpp

 #include "feature.h"
 #include "GaussianFilter.h"
 
-#include <ANN/ANN.h>
+#include "ANN/ANN.h"
 
 #include <algorithm>
 #include <vector>

project2/photomerge/source/GaussianFilter.cpp

     : m_Sigma(_sigma), m_Kernel(NULL), m_KernelSize(0)
   {
     int khs = 0;
-    while(gaussianWeight(khs, _sigma) >= 0.000001)
+    while(gaussianWeight(khs, _sigma) >= 0.001)
       khs = khs + 1;
     
     m_KernelSize = 2 * khs + 1;

project2/photomerge/source/Image.cpp

 namespace vfx
 {
   
-  // **************************************************************************
-  
   ByteImage::ByteImage()
-    : m_Data(NULL)
+  : m_Data(NULL)
   {
   }
   
   ByteImage::ByteImage(const ByteImage& obj)
-    : m_Data(obj.m_Data)
+  : m_Data(obj.m_Data)
   {
     if(m_Data)
-      ImageData_retain( static_cast< ImageData<byte>* >(m_Data) );
+      ImageData_retain( static_cast< ImageData<Byte>* >(m_Data) );
   }
   
   ByteImage& ByteImage::operator=(const ByteImage& obj)
   {
     if(m_Data)
-      ImageData_release( static_cast< ImageData<byte>* >(m_Data) );
+      ImageData_release( static_cast< ImageData<Byte>* >(m_Data) );
     
     m_Data = obj.m_Data;
     
     if(m_Data)
-      ImageData_retain( static_cast< ImageData<byte>* >(m_Data) );
+      ImageData_retain( static_cast< ImageData<Byte>* >(m_Data) );
     
     return *this;
   }
   ByteImage::~ByteImage()
   {
     if(m_Data)
-      ImageData_release( static_cast< ImageData<byte>* >(m_Data) );
+      ImageData_release( static_cast< ImageData<Byte>* >(m_Data) );
   }
   
   int ByteImage::width()const
   {
-    return static_cast< ImageData<byte>* >(m_Data)->width;
+    return static_cast< ImageData<Byte>* >(m_Data)->width;
   }
   
   int ByteImage::height()const
   {
-    return static_cast< ImageData<byte>* >(m_Data)->height;
+    return static_cast< ImageData<Byte>* >(m_Data)->height;
   }
   
   int ByteImage::row()const
   {
-    return static_cast< ImageData<byte>* >(m_Data)->row;
+    return static_cast< ImageData<Byte>* >(m_Data)->row;
   }
   
-  const byte* ByteImage::data()const
+  const Byte* ByteImage::data()const
   {
-    return ImageData_data(static_cast<const ImageData<byte>* >(m_Data));
+    return ImageData_data(static_cast<const ImageData<Byte>* >(m_Data));
   }
   
-  byte* ByteImage::data()
+  Byte* ByteImage::data()
   {
-    return ImageData_data(static_cast< ImageData<byte>* >(m_Data));
+    return ImageData_data(static_cast< ImageData<Byte>* >(m_Data));
   }
   
   ByteImage::ByteImage(void* data)
-    : m_Data(data)
+  : m_Data(data)
+  {
+  }
+  
+  // **************************************************************************
+  
+  ColorImage::ColorImage()
+  : m_Data(NULL)
+  {
+  }
+  
+  ColorImage::ColorImage(const ColorImage& obj)
+  : m_Data(obj.m_Data)
+  {
+    if(m_Data)
+      ImageData_retain( static_cast< ImageData<PixARGB>* >(m_Data) );
+  }
+  
+  ColorImage& ColorImage::operator=(const ColorImage& obj)
+  {
+    if(m_Data)
+      ImageData_release( static_cast< ImageData<PixARGB>* >(m_Data) );
+    
+    m_Data = obj.m_Data;
+    
+    if(m_Data)
+      ImageData_retain( static_cast< ImageData<PixARGB>* >(m_Data) );
+    
+    return *this;
+  }
+  
+  ColorImage::~ColorImage()
+  {
+    if(m_Data)
+      ImageData_release( static_cast< ImageData<PixARGB>* >(m_Data) );
+  }
+  
+  int ColorImage::width()const
+  {
+    return static_cast< ImageData<PixARGB>* >(m_Data)->width;
+  }
+  
+  int ColorImage::height()const
+  {
+    return static_cast< ImageData<PixARGB>* >(m_Data)->height;
+  }
+  
+  int ColorImage::row()const
+  {
+    return static_cast< ImageData<PixARGB>* >(m_Data)->row;
+  }
+  
+  const PixARGB* ColorImage::data()const
+  {
+    return ImageData_data(static_cast<const ImageData<PixARGB>* >(m_Data));
+  }
+  
+  PixARGB* ColorImage::data()
+  {
+    return ImageData_data(static_cast< ImageData<PixARGB>* >(m_Data));
+  }
+  
+  ColorImage::ColorImage(void* data)
+  : m_Data(data)
   {
   }
   
   
   ByteImage ImageCreator::createByteImage(int w, int h)
   {
-    return ByteImage( ImageData_create<byte>(w, h) );
+    return ByteImage( ImageData_create<Byte>(w, h) );
+  }
+  
+  ColorImage ImageCreator::createColorImage(int w, int h)
+  {
+    return ColorImage( ImageData_create<PixARGB>(w, h) );
   }
   
   FloatImage ImageCreator::createFloatImage(int w, int h)
     return FloatImage( ImageData_create<float>(w, h) );
   }
   
-  ByteImage ImageCreator::loadFromJpeg(const char* path)
+  ColorImage ImageCreator::loadFromJpeg(const char* path)
   {
-    return ByteImage(ImageData_loadJpegY<byte>(path));
+    return ColorImage(ImageData_loadJpeg<PixARGB>(path));
   }
   
-  bool ImageCreator::saveToJpeg(const char* path, const ByteImage& img)
+  bool ImageCreator::saveToJpeg(const char* path, const ColorImage& img)
   {
-    const ImageData<byte>* ptr = (const ImageData<byte>*)(img.m_Data);
-    return ImageData_saveJpegY(path, ptr);
+    const ImageData<PixARGB>* ptr = (const ImageData<PixARGB>*)(img.m_Data);
+    return ImageData_saveJpegRGB(path, ptr);
   }
   
-  bool ImageCreator::saveToJpeg(const char* path, 
-                                const ByteImage& imgR, 
-                                const ByteImage& imgG, 
-                                const ByteImage& imgB)
+  bool ImageCreator::saveToJpeg(const char* path, const FloatImage& img)
   {
-    const ImageData<byte>* pR = (const ImageData<byte>*)(imgR.m_Data);
-    const ImageData<byte>* pG = (const ImageData<byte>*)(imgG.m_Data);
-    const ImageData<byte>* pB = (const ImageData<byte>*)(imgB.m_Data);
-    return ImageData_saveJpegRGB(path, pR, pG, pB);
+    const ImageData<float>* ptr = (const ImageData<float>*)(img.m_Data);
+    return ImageData_saveJpegY(path, ptr);
   }
   
   namespace
   {
-    inline float _b2f(byte c)
+    inline float _argb2f(PixARGB c)
     {
-      return c / 255.0f;
+      return makePixel<float>( pixelR(c), pixelG(c), pixelB(c), 255 );
     }
     
-    inline byte _f2b(float c)
+    inline PixARGB _copypixel(PixARGB c)
     {
-      return c * 255;
+      return c;
     }
   }
   
-  ByteImage ImageCreator::toByteImage(const FloatImage& img)
+  FloatImage ImageCreator::toFloatImage(const ColorImage& img)
   {
-    const ImageData<float>* ptr = (const ImageData<float>*)(img.m_Data);
-    return ByteImage( ImageData_clone(ptr, byte(0), _f2b) );
+    const ImageData<PixARGB>* ptr = (const ImageData<PixARGB>*)(img.m_Data);
+    return FloatImage( ImageData_clone(ptr, float(0), _argb2f) );
   }
   
-  FloatImage ImageCreator::toFloatImage(const ByteImage& img)
+  ColorImage ImageCreator::clone(const ColorImage& img)
   {
-    const ImageData<byte>* ptr = (const ImageData<byte>*)(img.m_Data);
-    return FloatImage( ImageData_clone(ptr, float(0), _b2f) );
+    const ImageData<PixARGB>* ptr = (const ImageData<PixARGB>*)(img.m_Data);
+    return ColorImage( ImageData_clone(ptr, PixARGB(0), _copypixel) );
   }
   
   // **************************************************************************
     return FloatImage( ImageData_shrink(ptr) );
   }
   
-  void ImageTool::fill(ByteImage& img, byte c)
+  void ImageTool::fill(ByteImage& img, Byte c)
   {
-    ImageData<byte>* ptr = (ImageData<byte>*)(img.m_Data);
+    ImageData<Byte>* ptr = (ImageData<Byte>*)(img.m_Data);
+    ImageData_fill(ptr, c);
+  }
+  
+  void ImageTool::fill(ColorImage& img, PixARGB c)
+  {
+    ImageData<PixARGB>* ptr = (ImageData<PixARGB>*)(img.m_Data);
     ImageData_fill(ptr, c);
   }
   
   }
 
   void ImageTool::line
-  (ByteImage& img, int x1, int y1, int x2, int y2, byte c)
+  (ColorImage& img, int x1, int y1, int x2, int y2, PixARGB c)
   {
-    ImageData<byte>* ptr = (ImageData<byte>*)(img.m_Data);
+    ImageData<PixARGB>* ptr = (ImageData<PixARGB>*)(img.m_Data);
     ImageData_line(ptr, x1, y1, x2, y2, c);
   }
   
     return dst;
   }
   
+  namespace
+  {
+    template<bool x_min, bool x_max, bool y_min, bool y_max>
+    inline float _convolute3x3(const float* s, int rs, 
+                               float k00, float k01, float k02,
+                               float k10, float k11, float k12,
+                               float k20, float k21, float k22)
+    {
+      float sum = 0.0f;
+      
+      if(!y_min)
+      {
+        sum += k00 * ( !x_min ? s[-rs-1] : s[-rs] );
+        sum += k01 * ( s[-rs] );
+        sum += k02 * ( !x_max ? s[-rs+1] : s[-rs] );
+      }
+      else
+      {
+        sum += k00 * ( !x_min ? s[-1] : s[0] );
+        sum += k01 * ( s[0] );
+        sum += k02 * ( !x_max ? s[1] : s[0] );
+      }
+      
+      {
+        sum += k10 * ( !x_min ? s[-1] : s[0] );
+        sum += k11 * ( s[0] );
+        sum += k12 * ( !x_max ? s[1] : s[0] );
+      }
+      
+      if(!y_max)
+      {
+        sum += k20 * ( !x_min ? s[rs-1] : s[rs] );
+        sum += k21 * ( s[rs] );
+        sum += k22 * ( !x_max ? s[rs+1] : s[rs] );
+      }
+      else
+      {
+        sum += k20 * ( !x_min ? s[-1] : s[0] );
+        sum += k21 * ( s[0] );
+        sum += k22 * ( !x_max ? s[1] : s[0] );
+      }
+      
+      return sum;
+    }
+  }
+  
   FloatImage ImageTool::convolute3x3(const FloatImage& src, 
                                      float k00, float k01, float k02,
                                      float k10, float k11, float k12,
       // x
       for(int x = 1; x < w-1; ++x)
       {
-        float sum = 0.0f;
-        sum += k00 * s[x-rs-1] + k01 * s[x-rs] + k02 * s[x-rs+1];
-        sum += k10 * s[x   -1] + k11 * s[x   ] + k12 * s[x   +1];
-        sum += k20 * s[x+rs-1] + k21 * s[x+rs] + k22 * s[x+rs+1];
+        //float sum = 0.0f;
+        //sum += k00 * s[x-rs-1] + k01 * s[x-rs] + k02 * s[x-rs+1];
+        //sum += k10 * s[x   -1] + k11 * s[x   ] + k12 * s[x   +1];
+        //sum += k20 * s[x+rs-1] + k21 * s[x+rs] + k22 * s[x+rs+1];
         
-        d[x] = sum;
+        float c;
+        c = _convolute3x3<false, false, false, false>(s+x, rs,
+                                                      k00, k01, k02,
+                                                      k10, k11, k12,
+                                                      k20, k21, k22);
+        d[x] = c;
       }
       
       // x = w-1

project2/photomerge/source/Image.h

 #pragma once
 
+#include "ImagePixel.h"
+
 #include <algorithm>
 
 #include <cstddef>
 
 namespace vfx
 {
-  typedef unsigned char byte;
-  
-  // *************************************************************************
   
   class ByteImage
   {
     int height()const;
     int row()const;
     
-    const byte* data()const;
-    byte* data();
+    const Byte* data()const;
+    Byte* data();
     
     bool valid()const { return m_Data != NULL; }
     
     ByteImage(void* data);
   };
   
+  class ColorImage
+  {
+    friend class ImageCreator;
+    friend class ImageTool;
+    
+  public:
+    ColorImage();
+    ColorImage(const ColorImage& obj);
+    ColorImage& operator=(const ColorImage& obj);
+    ~ColorImage();
+    
+    int width()const;
+    int height()const;
+    int row()const;
+    
+    const PixARGB* data()const;
+    PixARGB* data();
+    
+    bool valid()const { return m_Data != NULL; }
+    
+  private:
+    void* m_Data;
+    ColorImage(void* data);
+  };
+  
   class FloatImage
   {
     friend class ImageCreator;
   {
   public:
     static ByteImage createByteImage(int w, int h);
+    static ColorImage createColorImage(int w, int h);
     static FloatImage createFloatImage(int w, int h);
     
-    static ByteImage loadFromJpeg(const char* path);
-    static bool saveToJpeg(const char* path, const ByteImage& img);
-    static bool saveToJpeg(const char* path, 
-                           const ByteImage& imgR, 
-                           const ByteImage& imgG, 
-                           const ByteImage& imgB);
+    static ColorImage loadFromJpeg(const char* path);
+    static bool saveToJpeg(const char* path, const ColorImage& img);
+    static bool saveToJpeg(const char* path, const FloatImage& img);
     
-    static ByteImage toByteImage(const FloatImage& img);
-    static FloatImage toFloatImage(const ByteImage& img);
+    static FloatImage toFloatImage(const ColorImage& img);
+    static ColorImage clone(const ColorImage& img);
     
   };
   
   class ImageTool
   {
   public:
+    
     static FloatImage shrink(const FloatImage& src);
     
-    static void fill(ByteImage& img, byte c);
+    static void fill(ByteImage& img, Byte c);
+    static void fill(ColorImage& img, PixARGB c);
     static void fill(FloatImage& img, float c);
 
     static void line
-    (ByteImage& img, int x1, int y1, int x2, int y2, byte c);
+    (ColorImage& img, int x1, int y1, int x2, int y2, PixARGB c);
     
     static void line
     (FloatImage& img, int x1, int y1, int x2, int y2, float c);

project2/photomerge/source/feature.cpp

 #include "feature.h"
+#include "GaussianFilter.h"
 
-#include <cassert>
+#include "ANN/ANN.h"
+
+#include <algorithm>
+#include <vector>
 #include <cstdio>
 #include <cstdlib>
 #include <cstring>
-#include <cmath>
-
+#include <ctime>
 #include <list>
 
 namespace vfx
 {
   
+  enum{ 
+    FEATURE_DESC_WIDTH = 8,
+    FEATURE_DESC_SIZE = 64,
+    FEATURE_DESC_WINDOW_SIZE = 40
+  };
+  
+  const int FEATURE_PYRAMID_LEVELS = 4;
+  const float FEATURE_PYRAMID_SIGMA = 1.0f;
+  const float FEATURE_SIGMA_D = 1.0f;
+  const float FEATURE_SIGMA_I = 1.5f;
+  const float FEATURE_SIGMA_O = 4.5f;
+  const float FEATURE_RESPONSE_THRESHOLD = 10.0f;
+  const int FEATURE_NONMAXIMAL_RADIUS_DENOM = 50;
+  const int FEATURE_NUMBER_LOWERBOUND = 1000;
+  
+  const float FEATURE_MATCH_THRESHOLD = 0.6f;
+  
+  const int FEATURE_RANSAC_SAMPLES = 2;
+  const float FEATURE_RANSAC_THRESHOLD = 9.0f;
+  const float FEATURE_RANSAC_SUCESSS_PROBABILITY = 0.99;
+  const float FEATURE_RANSAC_INLIER_PROBABILITY = 0.6;
+  
+  struct Features
+  {
+    int levels, count;
+    std::vector<int> level;
+    std::vector<float> orient;
+    std::vector<float> x;
+    std::vector<float> y;
+    std::vector<float> desc[FEATURE_DESC_SIZE];
+  };
+  
+  static 
   FloatImage computeGradientX(const FloatImage& source)
   {
     return ImageTool::convolute3x3(source,
                                    );
   }
   
+  static 
   FloatImage computeGradientY(const FloatImage& source)
   {
     return ImageTool::convolute3x3(source,
                                    );
   }
   
+  static 
   FloatImage harrisResponse(const FloatImage& source, 
                             const GaussianFilter& filterD, 
                             const GaussianFilter& filterI)
     return response;
   }
   
+  static 
   FloatImage harrisOrient(const FloatImage& source, 
                           const GaussianFilter& filter)
   {
         float* d = pDst + y * rd;
         for(int x = 0; x < w; ++x)
           d[x] = std::atan2( sx[x], sy[x] ); // NOTE: It is strange.
-          //d[x] = std::atan2( sy[x], sx[x] );
+        //d[x] = std::atan2( sy[x], sx[x] );
       }
     }
     
     return orient;
   }
   
+  static 
   ByteImage findCorners(const FloatImage& response, float threshold)
   {
     const int w = response.width(), h = response.height();
     ByteImage mask = ImageCreator::createByteImage(w, h);
     
     const float* pSrc = response.data();
-    byte* pDst = mask.data();
+    Byte* pDst = mask.data();
     const int rs = response.row(), rd = mask.row();
     
     ImageTool::fill(mask, 0);
     for(int y = 1; y < h-1; ++y)
     {
       const float* s = pSrc + y * rs;
-      byte* d = pDst + y * rd;
+      Byte* d = pDst + y * rd;
       
       for(int x = 1; x < w-1; ++x)
       {
     return mask;
   }
   
+  static 
   ByteImage processSubPixelAccuracy(const ByteImage& featmask, 
                                     const FloatImage& response)
   {
     ImageTool::fill(newmask, 0);
     
     const int iFeatSrcRow = featmask.row();
-    const byte* pFeatSrc = featmask.data();
+    const Byte* pFeatSrc = featmask.data();
     const int rs = response.row();
     const float* pSrc = response.data();
     const int rd = newmask.row();
-    byte* pDst = newmask.data();
+    Byte* pDst = newmask.data();
     
     for(int y = 1; y < h-1; ++y)
     {
-      const byte* fs = pFeatSrc + y * iFeatSrcRow;
+      const Byte* fs = pFeatSrc + y * iFeatSrcRow;
       const float* s = pSrc + y * rs;
       
       for(int x = 1; x < w-1; ++x)
     
   }
   
+  static 
   void detectFeatures(const FloatImage& inputImg, Features& feats)
   {
     const bool debug = false;
     
     // Step 1: Build Pyramid
-    const int levels = 4;
-    const float sigma = 1.0f;
+    const int levels = FEATURE_PYRAMID_LEVELS;
+    const float sigma = FEATURE_PYRAMID_SIGMA;
     std::vector<FloatImage> pyramid(levels, inputImg);
     {
       GaussianFilter filter(sigma);
     }
     
     // Step 2: Calculate Harris Corner Response
-    const float sigmaD = 1.0f, sigmaI = 1.5f;
+    const float sigmaD = FEATURE_SIGMA_D, sigmaI = FEATURE_SIGMA_I;
     std::vector<FloatImage> response(levels);
     {
       GaussianFilter filterD(sigmaD), filterI(sigmaI);
     }
     
     // Step 3: Harris Feature Orientation
-    const float sigmaO = 4.5;
+    const float sigmaO = FEATURE_SIGMA_O;
     std::vector<FloatImage> orient(levels);
     {
       GaussianFilter filter(sigmaO);
     }
     
     // Step 4: Local Maxima and Thresholding
-    const float responseThreshold = 10.0;
+    const float responseThreshold = FEATURE_RESPONSE_THRESHOLD;
     std::vector<ByteImage> featmask(levels);
     {
       for(int lv = 0; lv < levels; ++lv)
         
         for(int y = 0; y < h; ++y)
         {
-          byte* pMask = featmask[lv].data() + y * featmask[lv].row();
+          Byte* pMask = featmask[lv].data() + y * featmask[lv].row();
           const float* pRespon = response[lv].data() + y * response[lv].row();
           const float* pOrient = orient[lv].data() + y * orient[lv].row();
-          byte* pFeatLv = featLevel.data() + sc * y * featLevel.row();
+          Byte* pFeatLv = featLevel.data() + sc * y * featLevel.row();
           float* pFeatOr = featOrient.data() + sc * y * featOrient.row();
           float* pFeatRp = featRespon.data() + sc * y * featRespon.row();
           
     }
     
     // Step 7: Convert Feature Map to Feature List
+    int totalFeatures = 0;
     {
-      // TODO:
+      const int w = inputImg.width(), h = inputImg.height();
+      
+      for(int y = 0; y < h; ++y)
+      {
+        const Byte* pLv = featLevel.data() + y * featLevel.row();
+        
+        for(int x = 0; x < w; ++x)
+        {
+          if(pLv[x] == 255)
+            continue;
+          totalFeatures++;
+        }
+      }
     }
     
     // Step 8: Non-Maximal Suppression
     //         (Adaptive Non-Maximal Suppression
     // TODO: set a variable to constraint the number of key points
+    const int nonmaximal_radius_ratio = FEATURE_NONMAXIMAL_RADIUS_DENOM;
+    const int lowerbound_nFeatures = FEATURE_NUMBER_LOWERBOUND;
     {
       const int w = inputImg.width(), h = inputImg.height();
       
       // TODO: tunable factor 
-      int maxRadius = std::min(w, h) / 100;
+      int maxRadius = std::min(w, h) / nonmaximal_radius_ratio;
       
-      for(int radius = 8; radius < maxRadius; ++radius)
+      for(int radius = 8; 
+          radius < maxRadius && totalFeatures > lowerbound_nFeatures; 
+          ++radius)
       {
         int r2 = radius * radius;
         
         int iRpRow = featRespon.row();
         ByteImage newLevel = ImageCreator::createByteImage(w, h);
+        ImageTool::fill(newLevel, 255);
         
         for(int y = 0; y < h; ++y)
         {
-          const byte* pLv = featLevel.data() + y * featLevel.row();
+          const Byte* pLv = featLevel.data() + y * featLevel.row();
           const float* pRp = featRespon.data() + y * featRespon.row();
-          byte* pDst = newLevel.data() + y * newLevel.row();
+          Byte* pDst = newLevel.data() + y * newLevel.row();
           
           for(int x = 0; x < w; ++x)
           {
             if(pLv[x] == 255)
-            {
-              pDst[x] = 255;
               continue;
-            }
             
             bool valid = true;
             
                   valid = false;
               }
             
-            if(!valid)
-              pDst[x] = 255;
+            if(!valid && totalFeatures > lowerbound_nFeatures)
+              pDst[x] = 255, --totalFeatures ;
             else
               pDst[x] = pLv[x];
           }
       
       for(int y = 0; y < h; ++y)
       {
-        const byte* pLv = featLevel.data() + y * featLevel.row();
+        const Byte* pLv = featLevel.data() + y * featLevel.row();
         const float* pOr = featOrient.data() + y * featOrient.row();
         
         for(int x = 0; x < w; ++x)
             }
           
           // Normalize Descriptor
-          if(false){
+          if(true)
+          {
             float mean = 0.0f, stdv = 0.0f;
             
             for(int i = 0; i < FEATURE_DESC_SIZE; ++i)
         FloatImage tmp;
         
         std::sprintf(buf, "pyramid-%d.jpg", i);
-        ImageCreator::saveToJpeg(buf, ImageCreator::toByteImage(pyramid[i]));
+        ImageCreator::saveToJpeg(buf, pyramid[i]);
         
         std::sprintf(buf, "response-%d.jpg", i);
-        ImageCreator::saveToJpeg(buf, ImageCreator::toByteImage(response[i]));
+        ImageCreator::saveToJpeg(buf, response[i]);
         
         std::sprintf(buf, "orient-%d.jpg", i);
-        ImageCreator::saveToJpeg(buf, ImageCreator::toByteImage(orient[i]));
-        
-        std::sprintf(buf, "featmask-%d.jpg", i);
-        ImageCreator::saveToJpeg(buf, (featmask[i]));
+        ImageCreator::saveToJpeg(buf, orient[i]);
         
         std::sprintf(buf, "grad-x-%d.jpg", i);
         tmp = computeGradientX(pyramid[i]);
-        ImageCreator::saveToJpeg(buf, ImageCreator::toByteImage(tmp));
+        ImageCreator::saveToJpeg(buf, tmp);
         
         std::sprintf(buf, "grad-y-%d.jpg", i);
         tmp = computeGradientY(pyramid[i]);
-        ImageCreator::saveToJpeg(buf, ImageCreator::toByteImage(tmp));
+        ImageCreator::saveToJpeg(buf, tmp);
         
       }
     }
   }
   
   //
+  
+  static 
   void outputFeatureMap(const char* path, 
-                        const FloatImage& inputImg, 
+                        const ColorImage& src, 
                         const Features& feats)
   {
-    ByteImage finalR = ImageCreator::toByteImage(inputImg);
-    ByteImage finalGB = ImageCreator::toByteImage(inputImg);
+    ColorImage final = ImageCreator::clone(src);
     
     std::vector<float>::const_iterator itX = feats.x.begin();
     std::vector<float>::const_iterator itY = feats.y.begin();
       
       int x = *itX, y = *itY, sc = feats.levels - lv - 1;
       
+      int siz = FEATURE_DESC_WINDOW_SIZE / 2;
       float cor[4][2] = {
-        { -20<<sc, -20<<sc }, { 19<<sc, -20<<sc }, 
-        { 19<<sc, 19<<sc }, { -20<<sc, 19<<sc }
+        { -siz<<sc, -siz<<sc }, { siz<<sc, -siz<<sc }, 
+        { siz<<sc, siz<<sc }, { -siz<<sc, siz<<sc }
       };
       
       for(int i = 0; i < 4; ++i)
       {
         int x1 = cor[i][0], y1 = cor[i][1];
         int x2 = cor[(i+1)%4][0], y2 = cor[(i+1)%4][1];
-        ImageTool::line(finalR, x1, y1, x2, y2, 255);
-        ImageTool::line(finalGB, x1, y1, x2, y2, 0);
+        ImageTool::line(final, x1, y1, x2, y2, 0xff0000);
         //printf("%d %d %d %d\n", x1, y1, x2, y2);
       }
       
       float tx, ty;
       _rotate(rot, 19<<sc, 0, tx, ty);
       
-      ImageTool::line(finalR, x, y, x+tx, y+ty, 255);
-      ImageTool::line(finalGB, x, y, x+tx, y+ty, 0);
+      ImageTool::line(final, x, y, x+tx, y+ty, 0xff0000);
+    }
+    
+    ImageCreator::saveToJpeg(path, final);
+  }
+  
+  // ********************************************************************
+  
+  static 
+  std::vector<ImageMatchingPair> matchFeatures
+  (const std::vector<Features> feats, int nearestNeighborNum)
+  {
+    const float threshold = FEATURE_MATCH_THRESHOLD;
+    
+    int numImgs = feats.size();
+    const int descsize = FEATURE_DESC_SIZE;
+    
+    // Init ANN Resources
+    
+    std::vector< ANNpointArray > descs(numImgs);
+    std::vector< ANNkd_tree* > kdtrees(numImgs);
+    
+    for(int i = 0; i < numImgs; ++i)
+    {
+      descs[i] = annAllocPts( feats[i].count, descsize );
+      for(int n = 0; n < feats[i].count; ++n)
+        for(int k = 0; k < descsize; ++k)
+          descs[i][n][k] = feats[i].desc[k][n];
+    }
+    
+    for(int i = 0; i < numImgs; ++i)
+      kdtrees[i] = new ANNkd_tree(descs[i], feats[i].count, descsize);
+    
+    ANNpoint testTuple = annAllocPt(descsize);
+    ANNidxArray nnIdxs = new ANNidx[nearestNeighborNum];
+    ANNdistArray dists = new ANNdist[nearestNeighborNum];
+    
+    // Create Match List
+    std::vector<ImageMatchingPair> matchList;
+    matchList.reserve( numImgs * (numImgs-1) / 2 );
+    
+    // Do Dual-Way Matching
+    for(int second = 0; second < numImgs; ++second)
+      for(int first = 0; first < second; ++first)
+      {
+        int cnt = feats[first].count;
+        
+        std::list< std::pair<int, int> > feature_pairs;
+        
+        for(int n = 0; n < cnt; ++n)
+        {
+          for(int d = 0; d < descsize; ++d)
+            testTuple[d] = feats[first].desc[d][n];
+          
+          kdtrees[second]->annkSearch(testTuple, nearestNeighborNum, 
+                                      nnIdxs, dists, 0);
+          
+          if(dists[0] < threshold * dists[1])
+          {
+            int f1 = n, f2 = nnIdxs[0];
+            
+            for(int d = 0; d < descsize; ++d)
+              testTuple[d] = feats[second].desc[d][f2];
+            
+            kdtrees[first]->annkSearch(testTuple, nearestNeighborNum, 
+                                       nnIdxs, dists, 0);
+            
+            if(dists[0] < threshold * dists[1] && nnIdxs[0] == f1)
+              feature_pairs.push_back( std::make_pair(f1, f2) );
+          }
+        }
+        
+        ImageMatchingPair data;
+        data.first = first;
+        data.second = second;
+        data.x1.reserve( feature_pairs.size() );
+        data.y1.reserve( feature_pairs.size() );
+        data.x2.reserve( feature_pairs.size() );
+        data.y2.reserve( feature_pairs.size() );
+        
+        std::list< std::pair<int, int> >::iterator it;
+        for(it = feature_pairs.begin(); it != feature_pairs.end(); ++it)
+        {
+          int f1 = it->first, f2 = it->second;
+          data.x1.push_back( feats[first].x[f1] );
+          data.y1.push_back( feats[first].y[f1] );
+          data.x2.push_back( feats[second].x[f2] );
+          data.y2.push_back( feats[second].y[f2] );
+        }
+        
+        data.count = data.x1.size();
+        
+        matchList.push_back(data);
+      }
+    
+    
+    // Release ANN Resources
+    
+    annDeallocPt(testTuple);
+    delete[] nnIdxs;
+    delete[] dists;
+    
+    for(int i = 0; i < numImgs; ++i)
+      delete kdtrees[i];
+    for(int i = 0; i < numImgs; ++i)
+      annDeallocPts( descs[i] );
+    
+    // Return
+    return matchList;
+  }
+  
+  // ********************************************************************
+  
+  static 
+  void ransac(ImageMatchingPair& mat)
+  {
+    const int nsamples = FEATURE_RANSAC_SAMPLES;
+    const float ransac_threshold = FEATURE_RANSAC_THRESHOLD;
+    const float ransac_success_prob = FEATURE_RANSAC_SUCESSS_PROBABILITY;
+    const float ransac_inlier_prob = FEATURE_RANSAC_INLIER_PROBABILITY;
+    const int nfeatpairs = mat.count;
+    
+    std::srand(time(NULL));
+    int nRANSAC = (int)(log(1.0 - ransac_success_prob) / 
+                        log(1.0 - pow(ransac_inlier_prob, nsamples)) );
+    
+    int best_nInliers = 0;
+    float best_dx = 0.0f, best_dy = 0.0f;
+    std::vector<bool> best_selected(nfeatpairs, false);
+    
+    std::vector<bool> isSelected(nfeatpairs);
+    
+    for(int k = 0; k < nRANSAC; ++k)
+    {
+      isSelected.assign(nfeatpairs, false);
+      
+      float dx = 0.0, dy = 0.0;
+      
+      for(int i = 0; i < nsamples; ++i)
+      {
+        int pidx;
+        do {
+          pidx = std::rand() % nfeatpairs;
+        } while(isSelected[pidx]);
+        isSelected[pidx] = true;
+        
+        float x1 = mat.x1[pidx];
+        float y1 = mat.y1[pidx];
+        float x2 = mat.x2[pidx];
+        float y2 = mat.y2[pidx];
+        
+        dx += x2 - x1, dy += y2 - y1;
+      }
+      
+      int nInliers = nsamples;
+      dx /= nsamples, dy /= nsamples;
+      
+      for(int i = 0; i < nfeatpairs; ++i)
+        if(!isSelected[i])
+        {
+          float x1 = mat.x1[i];
+          float y1 = mat.y1[i];
+          float x2 = mat.x2[i];
+          float y2 = mat.y2[i];
+          
+          float dist = (x2-x1-dx)*(x2-x1-dx)+(y2-y1-dy)*(y2-y1-dy);
+          if(dist < ransac_threshold)
+            nInliers++, isSelected[i] = true;
+        }
+      
+      //
+      if(nInliers > best_nInliers)
+      {
+        best_nInliers = nInliers;
+        best_dx = dx, best_dy = dy;
+        best_selected.assign(isSelected.begin(), isSelected.end());
+      }
+    }
+    
+    std::vector<float> temp;
+    temp.reserve(nfeatpairs);
+    
+    for(int i = 0; i < nfeatpairs; ++i)
+    {
+      if(best_selected[i])
+        temp.push_back( mat.x1[i] );
+    }
+    std::swap( mat.x1, temp );
+    temp.clear();
+    
+    for(int i = 0; i < nfeatpairs; ++i)
+    {
+      if(best_selected[i])
+        temp.push_back( mat.y1[i] );
+    }
+    std::swap( mat.y1, temp );
+    temp.clear();
+    
+    for(int i = 0; i < nfeatpairs; ++i)
+    {
+      if(best_selected[i])
+        temp.push_back( mat.x2[i] );
+    }
+    std::swap( mat.x2, temp );
+    temp.clear();
+    
+    for(int i = 0; i < nfeatpairs; ++i)
+    {
+      if(best_selected[i])
+        temp.push_back( mat.y2[i] );
+    }
+    std::swap( mat.y2, temp );
+    temp.clear();
+    
+    // update
+    mat.count = best_nInliers;
+    
+    // done
+  }
+  
+  // ********************************************************************
+  
+  namespace
+  {
+    inline bool _comparator(const ImageMatchingPair& a, 
+                            const ImageMatchingPair& b)
+    {
+      return a.count > b.count;
+    }
+  }
+  
+  std::vector<ImageMatchingPair> 
+  featureDetectionAndMatching(int numImages, const ColorImage* images,
+                              bool dump)
+  {
+    std::vector<Features> features(numImages);
+    
+    for(int i = 0; i < numImages; ++i)
+    {
+      FloatImage img = ImageCreator::toFloatImage(images[i]);
+      detectFeatures(img, features[i]);
+      std::printf("Image %d: %d features.\n", i, (int)features[i].count);
+    }
+    
+    if(dump)
+      for(int i = 0; i < numImages; ++i)
+      {
+        char buf[128];
+        std::sprintf(buf, "feature-%d.jpg", i);
+        outputFeatureMap(buf, images[i], features[i]);
+      }
+    
+    // ANN matching
+    std::vector<ImageMatchingPair> pairs = matchFeatures(features, 2);
+    
+    // RANSAC
+    for(std::vector<ImageMatchingPair>::iterator it = pairs.begin();
+        it != pairs.end(); ++it )
+    {
+      ransac(*it);
+    }
+    
+    // Sort Matching List
+    std::sort(pairs.begin(), pairs.end(), _comparator);
+    
+    // Remove Invalid Matching
+    for(int i = pairs.size()-1; i >= 0; --i)
+    {
+      if(pairs[i].count <= FEATURE_RANSAC_SAMPLES)
+        pairs.pop_back();
+    }
+    
+    if(dump)
+    {
+      int n = pairs.size();
+      for(int i = 0; i < n; ++i)
+      {
+        const int id0 = pairs[i].first, id1 = pairs[i].second;
+        const int w0 = images[id0].width(), h0 = images[id0].height();
+        const int w1 = images[id1].width(), h1 = images[id1].height();
+        const int tw = std::max(w0, w1);
+        const int th = h0 + h1;
+        
+        ColorImage tmp = ImageCreator::createColorImage(tw, th);
+        
+        for(int y = 0; y < h0; ++y)
+        {
+          const PixARGB* pSrc0 = images[id0].data() + y * images[id0].row();
+          PixARGB* pDst = tmp.data() + y * tmp.row();
+          
+          for(int x = 0; x < w0; ++x)
+            pDst[x] = pSrc0[x];
+        }
+        
+        for(int y = 0; y < h1; ++y)
+        {
+          const PixARGB* pSrc1 = images[id1].data() + y * images[id1].row();
+          PixARGB* pDst = tmp.data() + (y+h0) * tmp.row();
+          
+          for(int x = 0; x < w0; ++x)
+            pDst[x] = pSrc1[x];
+        }
+        
+        for(int k = 0; k < pairs[i].count; ++k)
+        {
+          int x0 = pairs[i].x1[k];
+          int y0 = pairs[i].y1[k];
+          int x1 = pairs[i].x2[k];
+          int y1 = h0 + pairs[i].y2[k];
+          
+          ImageTool::line(tmp, x0, y0, x1, y1, 0xff0000);
+        }
+        
+        char buf[128];
+        std::sprintf(buf, "matching-%d-%d.jpg", id0, id1);
+        ImageCreator::saveToJpeg(buf, tmp);
+      }
     }
     
-    ImageCreator::saveToJpeg(path, finalR, finalGB, finalGB);
+    return pairs;
   }
   
 }

project2/photomerge/source/feature.h

 #pragma once
 
 #include "Image.h"
-#include "GaussianFilter.h"
 
 #include <vector>
 
 namespace vfx
 {
   
-  enum{ 
-    FEATURE_DESC_WIDTH = 20,
-    FEATURE_DESC_SIZE = 400,
-    FEATURE_DESC_WINDOW_SIZE = 40
-  };
-  
-  struct Features
+  struct ImageMatchingPair
   {
-    int levels, count;
-    std::vector<int> level;
-    std::vector<float> orient;
-    std::vector<float> x;
-    std::vector<float> y;
-    std::vector<float> desc[FEATURE_DESC_SIZE];
+    int first; // First Image
+    int second; // Second Image
+    
+    //
+    int count; // Number of matched key points
+    
+    // Matched Feature Points
+    std::vector<float> x1; // On First Image
+    std::vector<float> y1; // On First Image
+    std::vector<float> x2; // On Second Image
+    std::vector<float> y2; // On Second Image
   };
   
-  void detectFeatures(const FloatImage& inputImg, Features& feats);
-  
-  void outputFeatureMap(const char* path, 
-                        const FloatImage& inputImg, 
-                        const Features& feats);
+  std::vector<ImageMatchingPair> 
+  featureDetectionAndMatching(int numImages, const ColorImage* images, 
+                              bool dump);
   
 }
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.