image.cc

来自「用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵」· CC 代码 · 共 1,083 行 · 第 1/3 页

  void labNormalize (ImageStack& lab)
  {
    for(int x = 0; x < lab.size(1); x++)
    {
      for(int y = 0; y < lab.size(2); y++)
      {
        float L = lab(LAB_L,x,y);
        float a = lab(LAB_A,x,y);
        float b = lab(LAB_B,x,y);
        const float minab = -73;
        const float maxab = 95;
        const float range = maxab - minab;

        L = L / 100.0;
        a = (a - minab) / range;
        b = (b - minab) / range;
        L = Util::minmax(0.0f, L, 1.0f);
        a = Util::minmax(0.0f, a, 1.0f);
        b = Util::minmax(0.0f, b, 1.0f);
        lab(LAB_L,x,y) = L;
        lab(LAB_A,x,y) = a;
        lab(LAB_B,x,y) = b;
      }
    }
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // create a translated version of this image where 
  // the old image appears embedded in a new image of
  // size [newwidth x newheight].  undefined pixels 
  // are filled in with value fill.
  //
  void getTranslated (const Image& im, const int xoffset, const int yoffset,
                      const int newwidth, const int newheight, const float fill, 
                      Image& translated)
  {

      assert(newwidth >= 0);
      assert(newheight >= 0);
      
      int oldwidth = im.size(1);
      int oldheight = im.size(0);
      translated.resize(newwidth,newheight);
      translated.init(fill);

      for (int x = xoffset; x < oldwidth + xoffset; x++)
      {
        for (int y = yoffset; y < oldheight + yoffset; y++)
        {
          if ((y >= 0) && (y < newheight) && (x >= 0) && (x < newwidth))
          {
            translated(x,y) = im(x-xoffset,y-yoffset);
          }
        }
      }
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // utility method to compute an antialiasing filter of 
  // size 13 for the given downscaling factor
  //
  // used by getScaled
  //
  void createAntialiasFilter (const float scale, Image& filter) 
  {
      Util::Array1D<float> wind(13);
      Util::Array1D<float> b(6);
      float sum = 0;

      filter.resize(1,13);
      filter.init(0);

      for (int i = 0; i < 13; i++)
      {
          wind(i) = 0.54 - 0.46 * cos (2 * M_PI * i / 11);
      }

      for (int i = 0; i < 6; i++)
      {
          b(i) = sin (scale * M_PI * (i + 0.5)) / (M_PI * (i + 0.5));
      }

      for (int i = 0; i < 6; i++)
      {
          filter(1,i) = b(6 - i) * wind(i);
          filter(1,i + 6) = b(i) * wind(i + 6);
          sum = sum + filter(1,i) + filter(1,i + 6);
      }
      sum = fabs (sum);

      for (int i = 0; i < 12; i++)
      {
          filter(1,i) = filter(1,i) / sum;
      }
  }

  //
  // create a resized version of this image of size [newwidth x newheight]
  // if bilinear = true, use bilinear interpolation
  // otherwise use bicubic interpolation
  //
  void getScaled(const Image& im, const int newwidth, const int newheight, const bool bilinear, Image& scaled)
  {
      assert(newwidth >= 0);
      assert(newheight >= 0);

      //first compute the scale factor
      float oldheight = im.size(0);
      float oldwidth = im.size(1);
      float xscale = (float) newwidth / (float) oldwidth;
      float yscale = (float) newheight / (float) oldheight;

      //filter to prevent aliasing if necessary
      Image antialiased = im;
      
      if (xscale < 1)
      {
          Image filtx;
          Image filtered;
          createAntialiasFilter(xscale,filtx);
          getFiltered(antialiased,filtered,filtx);
          antialiased = filtered;
      }

      if (yscale < 1)
      {
          Image filty;
          Image filtered;
          createAntialiasFilter(yscale,filty);
          getFiltered(antialiased,filtered,filty);
          antialiased = filtered;
      }

      //build the affine matrix
      Matrix A(3,3);
      A.init(0);
      A(0,0) = xscale;
      A(1,1) = yscale;
      A(2,2) = 1.0;

      //transform the image
      getTransformed(antialiased,A,newwidth,newheight,0,0,bilinear,scaled);
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // create a rotated version of this image
  // using an appropriate affine transform.
  // if bilinear = true, use bilinear interpolation
  // otherwise use bicubic.
  //
  void getRotated(const Image& im, const float angle, const bool bilinear, Image& rotated) 
  {
      //put theta in 0 - 2pi
      float theta = angle;
      while (theta < 0)
          theta = (theta + 2 * M_PI);
      while (theta > 2 * M_PI)
          theta = (theta - 2 * M_PI);

      //build the affine matrix
      Matrix A(3,3);
      A(0,0) =  cos(theta); A(0,1) = sin(theta); A(0,2) = 0.0;
      A(1,0) = -sin(theta); A(1,1) = cos(theta); A(1,2) = 0.0;
      A(2,0) = 0.0;         A(2,1) = 0.0;        A(2,2) = 1.0;

      getTransformed(im, A, im.size(0), im.size(1), 0, 0, bilinear, rotated);
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // the cubic spline interpolation kernel
  //
  static float cubic_bspline (const float x)
  {
      float a, b, c, d;

      if ((x + 2.0) <= 0.0)
      {
          a = 0.0;
      }
      else
      {
          a = pow ((x + 2.0), 3.0);
      }

      if ((x + 1.0) <= 0.0)
      {
          b = 0.0;
      }
      else
      {
          b = pow ((x + 1.0), 3.0);
      }

      if (x <= 0)
      {
          c = 0.0;
      }
      else
      {
          c = pow (x, 3.0);
      }

      if ((x - 1.0) <= 0.0)
      {
          d = 0.0;
      }
      else
      {
          d = pow ((x - 1.0), 3.0);
      }

      return ((1.0 / 6.0) * (a - (4.0 * b) + (6.0 * c) - (4.0 * d)));
  }

  //
  // returns the inverse.  only works for 3x3 matricies
  //
  static void getInverse(const Matrix mat, Matrix& inv)
  {
      assert((mat.size(0) != 3) || (mat.size(1) != 3));
       
      inv.resize(3,3);

      float denom = mat(0, 0) * mat(1, 1) * mat(2, 2) -
                     mat(0, 0) * mat(1, 2) * mat(2, 1) -
                     mat(1, 0) * mat(0, 1) * mat(2, 2) +
                     mat(1, 0) * mat(0, 2) * mat(2, 1) +
                     mat(2, 0) * mat(0, 1) * mat(1, 2) -
                     mat(2, 0) * mat(0, 2) * mat(1, 1);
                                                                                                                         
      inv(0,0) = ( mat(1, 1) * mat(2, 2) - mat(1, 2) * mat(2, 1)) / denom;
      inv(0,1) = (-mat(0, 1) * mat(2, 2) + mat(0, 2) * mat(2, 1)) / denom;
      inv(0,2) = ( mat(0, 1) * mat(1, 2) - mat(0, 2) * mat(1, 1)) / denom;
      inv(1,0) = (-mat(1, 0) * mat(2, 2) + mat(1, 2) * mat(2, 0)) / denom;
      inv(1,1) = ( mat(0, 0) * mat(2, 2) - mat(0, 2) * mat(2, 0)) / denom;
      inv(1,2) = (-mat(0, 0) * mat(1, 2) + mat(0, 2) * mat(1, 0)) / denom;
      inv(2,0) = ( mat(1, 0) * mat(2, 1) - mat(1, 1) * mat(2, 0)) / denom;
      inv(2,1) = (-mat(0, 0) * mat(2, 1) + mat(0, 1) * mat(2, 0)) / denom;
      inv(2,2) = ( mat(0, 0) * mat(1, 1) - mat(0, 1) * mat(1, 0)) / denom;
  }

  //
  // returns a new image which is an affine
  // transformed version of this image.
  // newimage = A*image.  the new image
  // is of size (height, width) such that
  // the corners of the old image are transformed
  // to locations inside the new image
  //
  // if bilinear is TRUE then use bilinear interpolation
  // otherwise use bicubic B-spline interpolation
  //
  void getTransformed (const Image& im, const Matrix& A, const int width,
                         const int height, const int xoffset,
                         const int yoffset, const bool bilinear, 
                         Matrix& transformed) 
  {
      assert(width >= 0);
      assert(height >= 0);

      float oldwidth = (float)im.size(1);
      float oldheight = (float)im.size(0);

      Matrix B;  
      getInverse(A,B);

      //allocate the result
      transformed.resize(width,height);

      //transform the image
      for (float x = 0; x < width; x++)
      {
        for (float y = 0; y < height; y++)
        {
          //compute the coordinates in the original image plane
          float u = (x + xoffset) * B(0, 0) + (y + yoffset) * B(0, 1) + B(0, 2);
          float v = (x + xoffset) * B(1, 0) + (y + yoffset) * B(1, 1) + B(1, 2);

          //if it's outside the bounds of the 
          //source image, fill with zeros
          if ((u >= oldwidth) || (u < 0.0) || (v >= oldheight) || (v < 0.0))
          {
              transformed((int)x, (int)y) = 0.0;
          }
          else
          {
            //do bilinear or bicubic interpolation as required
            if (bilinear == true)
            {
              float u1 = floor (u);
              float u2 = u1 + 1;
              float v1 = floor (v);
              float v2 = v1 + 1;

              float du = u - u1;
              float dv = v - v1;

              u1 = Util::max (0.0f, u1);
              u2 = Util::min (oldwidth - 1, u2);
              v1 = Util::max (0.0f, v1);
              v2 = Util::min (oldheight - 1, v2);

              float val11 = im((int)u1, (int)v1);
              float val12 = im((int)u1, (int)v2);
              float val21 = im((int)u2, (int)v1);
              float val22 = im((int)u2, (int)v2);

              float val = (1 - dv) * (1 - du) * val11 + 
                            (1 - dv) * du * val12 + 
                             dv * (1 - du) * val21 + 
                              dv * du * val22;

              transformed((int)x, (int)y) = val;
            }
            else
            {
              float a = u - floor (u);
              float b = v - floor (v);
              float val = 0.0;
              for (int m = -1; m < 3; m++)
              {
                float r1 = cubic_bspline ((float) m - a);
                for (int n = -1; n < 3; n++)
                {
                  float r2 = cubic_bspline (-1.0 * ((float) n - b));
                  float u1 = floor (u) + m;
                  float v1 = floor (v) + n;
                  u1 = Util::min (oldwidth - 1, Util::max (0.0f, u1));
                  v1 = Util::min (oldheight - 1, Util::max (0.0f, v1));

                  val += im((int)u1,(int)v1) * r1 * r2;
                }
              }
              transformed((int)x,(int)y) = val;
            }
          }
        }
      }
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // filters the image via convolution with the given 
  // kernel and returns the resulting image.  kernel 
  // must have odd dimensions.
  //
  void getFiltered (const Image& im, const Image& kernel, Image& filtered) 
  {
      // image and kernel dimensions
      const int iwidth = im.size(0);
      const int iheight = im.size(1);
      const int kwidth = kernel.size(0);