image.cc

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

      const int kheight = kernel.size(1);

      //output is same size as input
      filtered.resize(iwidth,iheight);
      filtered.init(0);

      // the kernel must not be larger than the image
      assert (kwidth <= iwidth);
      assert (kheight <= iheight);

      // the kernel must be odd in each dimension so it can be centered
      // over a pixel
      assert ((kwidth % 2) == 1);
      assert ((kheight % 2) == 1);

      // radius of kernel in each dimension; also the coordinates of the
      // kernel's center pixel
      const int xr = kwidth / 2;
      const int yr = kheight / 2;

      //flip the kernel left-right and up-down
      Util::Array2D<float> kern(kwidth,kheight);
      kern.init(0);
      for (int x = 0; x < kwidth; x++)
      {
        const int xx = kwidth - 1 - x;
        for (int y = 0; y < kheight; y++)
        {
          const int yy = kheight - 1 - y;
          kern(xx,yy) = kernel(x,y);
        }
      }

      // padded image dimensions
      const int pwidth = iwidth + 2 * xr;
      const int pheight = iheight + 2 * yr;

      //create image with reflective padding
      Util::Array2D<float> pim(pwidth,pheight);

      // top left
      for (int y = 0; y < yr; y++)
      {
          const int py = yr - 1 - y;
          for (int x = 0; x < xr; x++)
            {
                const int px = xr - 1 - x;
                pim(px,py) = im(x,y);
            }
      }
      // top right
      for (int y = 0; y < yr; y++)
      {
          const int py = yr - 1 - y;
          for (int x = 0; x < xr; x++)
            {
                const int xs = iwidth - 1 - x;
                const int xd = xr + iwidth + x;
                pim(xd,py) = im(xs,y);
            }
      }
      // bottom left
      for (int y = 0; y < yr; y++)
      {
          const int ys = iheight - 1 - y;
          const int yd = yr + iheight + y;
          for (int x = 0; x < xr; x++)
            {
                const int px = xr - 1 - x;
                pim(px,yd) = im(x,ys);
            }
      }
      // bottom right
      for (int y = 0; y < yr; y++)
      {
          const int ys = iheight - 1 - y;
          const int yd = yr + iheight + y;
          for (int x = 0; x < xr; x++)
            {
                const int xs = iwidth - 1 - x;
                const int xd = xr + iwidth + x;
                pim(xd,yd) = im(xs,ys);
            }
      }
      // top
      for (int y = 0; y < yr; y++)
      {
          const int py = yr - 1 - y;
          for (int x = 0; x < iwidth; x++)
            {
                const int px = x + xr;
                pim(px,py) = im(x,y);
            }
      }
      // bottom
      for (int y = 0; y < yr; y++)
      {
          const int ys = iheight - 1 - y;
          const int yd = yr + iheight + y;
          for (int x = 0; x < iwidth; x++)
            {
                const int px = x + xr;
                pim(px,yd) = im(x,ys);
            }
      }
      // left 
      for (int y = 0; y < iheight; y++)
      {
          const int py = yr + y;
          for (int x = 0; x < xr; x++)
            {
                const int px = xr - 1 - x;
                pim(px,py) = im(x,y);
            }
      }
      // right
      for (int y = 0; y < iheight; y++)
      {
          const int py = yr + y;
          for (int x = 0; x < xr; x++)
            {
                const int xs = iwidth - 1 - x;
                const int xd = xr + iwidth + x;
                pim(xd,py) = im(xs,y);
            }
      }
      // center
      for (int y = 0; y < iheight; y++)
      {
          const int py = yr + y;
          for (int x = 0; x < iwidth; x++)
          {
              const int px = xr + x;
              pim(px,py) = im(x,y);
          }
      }


      // use direct access to underlying arrays for speed
      float *p_pim = pim.data();
      float *p_kern = kern.data();
      float *p_filtered = filtered.data();

      // do the convolution
      // interchange y and ky loops, and unroll ky loop
      // gets 371 MFLOPS (including all overhead) on 700MHz PIII (53% of peak)
      for (int x = 0; x < iwidth; x++)
      {
        for (int kx = 0; kx < kwidth; kx++)
        {
          const int pcol = (x + kx) * pheight;
          const int kcol = kx * kheight;
          int ky = 0;
          while (ky < (kheight & ~0x7))
          {
            const float k0 = p_kern[kcol + ky + 0];
            const float k1 = p_kern[kcol + ky + 1];
            const float k2 = p_kern[kcol + ky + 2];
            const float k3 = p_kern[kcol + ky + 3];
            const float k4 = p_kern[kcol + ky + 4];
            const float k5 = p_kern[kcol + ky + 5];
            const float k6 = p_kern[kcol + ky + 6];
            const float k7 = p_kern[kcol + ky + 7];
            float in0 = p_pim[pcol + 0 + ky + 0];
            float in1 = p_pim[pcol + 0 + ky + 1];
            float in2 = p_pim[pcol + 0 + ky + 2];
            float in3 = p_pim[pcol + 0 + ky + 3];
            float in4 = p_pim[pcol + 0 + ky + 4];
            float in5 = p_pim[pcol + 0 + ky + 5];
            float in6 = p_pim[pcol + 0 + ky + 6];
            for (int y = 0; y < iheight; y++)
            {
                const float in7 = p_pim[pcol + y + ky + 7];
                p_filtered[x*iheight + y] += k0*in0 + k1*in1 + k2*in2 + k3*in3 +
                                             k4*in4 + k5*in5 + k6*in6 + k7*in7;
                in0 = in1;
                in1 = in2;
                in2 = in3;
                in3 = in4;
                in4 = in5;
                in5 = in6;
                in6 = in7;
            }
            ky += 8;
          }
          while (ky < (kheight & ~0x3))
          {
            const float k0 = p_kern[kcol + ky + 0];
            const float k1 = p_kern[kcol + ky + 1];
            const float k2 = p_kern[kcol + ky + 2];
            const float k3 = p_kern[kcol + ky + 3];
            float in0 = p_pim[pcol + 0 + ky + 0];
            float in1 = p_pim[pcol + 0 + ky + 1];
            float in2 = p_pim[pcol + 0 + ky + 2];
            for (int y = 0; y < iheight; y++)
            {
                const float in3 = p_pim[pcol + y + ky + 3];
                p_filtered[x*iheight + y] += k0*in0 + k1*in1 + k2*in2 + k3*in3;
                in0 = in1;
                in1 = in2;
                in2 = in3;
            }
            ky += 4;
          }
          while (ky < (kheight & ~0x1))
          {
            const float k0 = p_kern[kcol + ky + 0];
            const float k1 = p_kern[kcol + ky + 1];
            float in0 = p_pim[pcol + 0 + ky + 0];
            for (int y = 0; y < iheight; y++)
            {
                const float in1 = p_pim[pcol + y + ky + 1];
                p_filtered[x*iheight + y] += k0*in0 + k1*in1;
                in0 = in1;
            }
            ky += 2;
          }
          while (ky < kheight)
          {
            const float k0 = p_kern[kcol + ky];
            for (int y = 0; y < iheight; y++)
            {
                p_filtered[x*iheight + y] += k0*p_pim[pcol + y + ky];
            }
            ky += 1;
          }
          assert (ky == kheight);
        }
      }
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // filter at the given radius.
  // the resulting value at a pixel is the n'th-order value
  // in a circular window cetnered at the pixel
  // 0 <= order <= 1, 0 is smallest value, 1 is largest.
  // Call with 0.5 to get median filtering.
  //
  void getPctFiltered(const Image& im, const float radius, const float order, Image& filtered) 
  {
      assert (order >= 0 && order <= 1);

      if (radius < 1)
      {
        filtered = im;
        return;
      }

      //allocate window and filtered image
      const int windowRadius = (int) ceil (radius);
      const int windowDiam = 2 * windowRadius + 1;
      const int windowPixels = windowDiam * windowDiam;
      Util::Array1D<float> values(windowPixels);

      const int iwidth = im.size(0);
      const int iheight = im.size(1);
      filtered.resize(iwidth,iheight);

      //loop over the image
      for (int x = 0; x < iwidth; x++)
      {
        for (int y = 0; y < iheight; y++)
        {
          //copy values out of the window
          int count = 0;
          for (int u = -windowRadius; u <= windowRadius; u++)
          {
            for (int v = -windowRadius; v <= windowRadius; v++)
            {
              if ((u * u + v * v) > radius * radius)
              {
                continue;
              }
              int yi = y + u;
              int xi = x + v;
              if (yi < 0 || yi >= iheight)
              {
                continue;
              }
              if (xi < 0 || xi >= iwidth)
              {
                continue;
              }
              assert (count < windowPixels);
              values(count++) = im(xi,yi);
            }
          }
          //sort the values in ascending order
          assert (count > 0);
          Util::sort(values.data(), count);
          assert(values(0) <= values(count - 1));
          //pick out percentile value
          int index = (int) rint (order * (count - 1));
          assert (index >= 0 && index < count);
          float pctVal = values(index);
          assert (pctVal >= values(0));
          assert (pctVal <= values(count - 1));
          filtered(x,y) = pctVal;
        }
      }
  }

  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // filter at the given radius.
  // the resulting value at a pixel is the maximum value
  // in a circular window cetnered at the pixel
  //
  void getMaxFiltered (const Image& im, const float radius, Image& filtered)
  {
    if (radius < 1)
    {
      filtered = im;
      return;
    }

    const int iwidth = im.size(0);
    const int iheight = im.size(1);
    const int windowRadius = (int)ceil(radius);
    filtered.resize(iwidth,iheight);

    //loop over the image
    for (int x = 0; x < iwidth; x++)
    {
      for (int y = 0; y < iheight; y++)
      {
        //extract max from window
        float maxVal = im(x,y);
        for (int u = -windowRadius; u <= windowRadius; u++)
        {
          for (int v = -windowRadius; v <= windowRadius; v++)
          {
            if ((u * u + v * v) > radius * radius)
            {
              continue;
            }
            int xi = x + u;
            int yi = y + v;
            if (yi < 0 || yi >= iheight)
            {
              continue;
            }
            if (xi < 0 || xi >= iwidth)
            {
              continue;
            }
            const float val = im(xi,yi);
            maxVal = Util::max(val,maxVal);
          }
        }
        // save max
        filtered(x,y) = maxVal;
      }
    }
  }

} //namespace Util