filterbank.cc

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

  float
  FilterBank::getSigmaY (int index) const
  {
    assert (index >= 0 && index < m_nfilters);
    return m_sigmasY(index);
  }

  float
  FilterBank::getOrient (int index) const
  {
    assert (index >= 0 && index < m_nfilters);
    return m_orients(index);
  }

  int FilterBank::getNscales () const
  {
    return m_nscales;
  }
                                                                                                       
  int FilterBank::getNorient () const
  {
    return m_norient;
  }
                                                                                                       
  int FilterBank::getIndexElong (int scale, int orient, int evenodd) const
  {
    return calcIndexElong(m_nscales, m_norient, scale, orient, evenodd);
  }
                                                                                                       
  int FilterBank::getIndexCenter (int scale) const
  {
    return calcIndexCenter(m_nscales, m_norient, scale);
  }


  int FilterBank::calcNumFilters (const int nscales, const int norient)
  {
    return nscales * (2 * norient + 1);
  }

  int FilterBank::calcIndexElong (const int nscales, const int norient,
                              const int scale, const int orient,
                              const int evenodd)
  {
    assert (scale >= 0 && scale < nscales);
    assert (orient >= 0 && orient < norient);
    assert (evenodd == 0 || evenodd == 1);
    return scale * (2 * norient + 1) + 2 * orient + evenodd;
  }

  int FilterBank::calcIndexCenter (const int nscales, const int norient,
                               const int scale)
  {
    assert (scale >= 0 && scale < nscales);
    return scale * (2 * norient + 1) + 2 * norient;
  }

  ////////////////////////////////////////////////////////////////////////////////////////////////

  //
  //run the filterbank on an image and return
  //an array of images containing numFilters elements.
  //
  void FilterBank::filter (const Util::Image& im, Util::ImageStack& filtered) const
  {
    filtered.resize(m_nfilters,im.size(0),im.size(1));
    filtered.init(0);
    Util::Image f;
    Util::Message::startBlock(m_nfilters,"filtering",1);
    for (int i = 0; i < m_nfilters; i++)
    {
      Util::Message::stepBlock(1);
      getFiltered(im,getFilter(i),f);
      assert(f.size(0) == filtered.size(1));
      assert(f.size(1) == filtered.size(2));
      for (int x = 0; x < f.size(0); x++)
      {
        for (int y = 0; y < f.size(1); y++)
        {
          filtered(i,x,y) = f(x,y);
        } 
      }
    }
    Util::Message::endBlock(1);
  }

  // given a set of quadrature pair images, returns 
  // half as many new images which give the orientation
  // energy at each point
  void FilterBank::orientationEnergy (const Util::ImageStack& filtered, Util::ImageStack& energy) const
  {
    int width = filtered.size(1);
    int height = filtered.size(2);
    energy.resize(m_nscales*m_norient,width,height);
    for (int scale = 0; scale < m_nscales; scale++)
    {
      for (int orient = 0; orient < m_norient; orient++)
      {
        for (int x = 0; x < width; x++)
        {
          for (int y = 0; y < height; y++)
          {

            const float e1 = filtered(getIndexElong (scale, orient,0), x, y);
            const float e2 = filtered(getIndexElong (scale, orient,1), x, y);
            energy(scale*m_norient + orient,x,y) = sqrt(e1*e1 + e2*e2);
          }
        }
      }
    }
  }

  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // create a 2D gaussian kernel with sigmaX and sigmaY where
  //  orient - the orientation in radians.
  //  deriv - 2 or 0 depending on whether we want an edge filter
  //          or just a gaussian bump
  //  hilbert - whether to take the hilbert transform in the Y direction
  //  
  void FilterBank::createGaussianKernel (const float sigmaX, const float sigmaY,
                                    const float support, const float orient,
                                    const int deriv, const bool hilbert, 
                                    Util::Image& kernelImage)
  {
      assert (deriv >= 0 && deriv <= 2);

      //make an odd size kernel
      int kernelSize = (int) rint (2.0 * support * Util::max (sigmaX, sigmaY));
      kernelSize = kernelSize + (1 - (kernelSize % 2));
      assert ((kernelSize % 2) == 1);
      int halfSize = (kernelSize - 1) / 2;

      //create a grayscale image
      kernelImage.resize(kernelSize, kernelSize);

      //compute the kernel value along the x and y directions at a fine
      //resolution
      int maxSize = 2 * kernelSize;
      int zoom = 64;
      int nsamples = zoom * maxSize + 1;
      int halfnsamples = zoom * maxSize / 2;
      Util::Array1D<float> xval(nsamples);
      Util::Array1D<float> yval(nsamples);
      for (int y = -halfnsamples; y <= halfnsamples; y++)
      {
        const float sigmaYsqr = (sigmaY * sigmaY);
        const float sigmaXsqr = (sigmaX * sigmaX);
        const float u = (((float) y) / ((float) zoom));
        xval(y + halfnsamples) = exp (-u * u / (2 * sigmaXsqr));
        if (deriv == 0)
        {
            yval(y + halfnsamples) = exp (-u * u / (2 * sigmaYsqr));
        }
        else if (deriv == 1)
        {
            yval(y + halfnsamples) = exp (-u * u / (2 * sigmaYsqr)) * (-u / sigmaYsqr);
        }
        else if (deriv == 2)
        {
            yval(y + halfnsamples) = exp (-u * u / (2 * sigmaYsqr)) * ((u * u / sigmaYsqr) - 1);
        }
      }

      //take the hilbert transform in the Y direction if necessary
      if (hilbert == true)
      {
        computeHilbert (yval);
      }

      //now create the filter by averaging values out of the sampled
      //kernels.
      for (int x = -halfSize; x <= halfSize; x++)
      {
        for (int y = -halfSize; y <= halfSize; y++)
        {
          float u = x * cos(orient) + y * sin(orient);
          float v = x * sin(orient) - y * cos(orient);

          // compute the x value
          int u1 = (int) rint (((u - 0.5) * zoom) + halfnsamples);
          int u2 = (int) rint (((u + 0.5) * zoom) + halfnsamples);
          u1 = Util::max (u1, 0);
          u2 = Util::min (u2, nsamples - 1);
          assert (u2 >= u1);
          float xvalue = 0;
          for (int i = u1; i <= u2; i++)
          {
              xvalue += xval(i);
          }
          xvalue /= (u2 - u1 + 1);

          // compute the y value
          int v1 = (int) rint (((v - 0.5) * zoom) + halfnsamples);
          int v2 = (int) rint (((v + 0.5) * zoom) + halfnsamples);
          v1 = Util::max (v1, 0);
          v2 = Util::min (v2, nsamples - 1);
          assert (v2 >= v1);
          float yvalue = 0;
          for (int i = v1; i <= v2; i++)
          {
              yvalue += yval(i);
          }
          yvalue /= (v2 - v1 + 1);

          float val = xvalue * yvalue;
          kernelImage(x + halfSize, y + halfSize) = val;
        }
      }

      //make it zero mean if it's a derivative of a gaussian
      if (deriv != 0)
      {
        makeZeroMean(kernelImage);
      }

      //make sure the filter has unit L1 norm 
      makeL1Normalized(kernelImage);
  }


  /////////////////////////////////////////////////////////////////////////////////////////////////////

  //
  // paste kernel into a bigger image
  //
  static void drawFilter(const Util::Image& filt, int cx, int cy, Util::Image& im)
  {
    int iwidth = im.size(0);
    int iheight = im.size(1);
    int w = filt.size(0);
    int h = filt.size(1);
    // (i,j) are coorindates local to the filter
    for (int i = 0; i < w; i++)
    {
      for (int j = 0; j < h; j++)
      {
        // (x,y) are image coordinates for (i,j)
        int x = cx + i - w / 2;
        int y = cy + j - h / 2;
        assert (x >= 0 && x < iwidth);
        assert (y >= 0 && y < iheight);
        im(x,y) = filt(i,j);
      }
    }
  }

  //
  // create a nice image of the filterbank for presentation or debugging
  //
  void FilterBank::getFilterBankImage (Util::Image& fbim) const
  {
      // figure out the max filter dimensions
      int maxd = 0;
      for (int i = 0; i < m_nfilters; i++)
      {
        maxd = Util::max (maxd,getFilter(i).size(0));
        maxd = Util::max (maxd,getFilter(i).size(1));
      }
      maxd++;

      // allocate an image large enough to hold images of all the filters
      int nscales = getNscales ();
      int norient = getNorient ();
      int iheight = maxd * nscales;
      int iwidth = maxd * (2 * norient + 1);
      fbim.resize(iwidth,iheight);
      fbim.init(0);

      // write the filters into the big image
      // (cx,cy) give the image coordinates for the center of
      // the this filter
      for (int scale = 0; scale < nscales; scale++)
      {
        // lay out filters for this scale like this: eoeoeoeoeoeoc
        for (int orient = 0; orient < norient; orient++)
        {
          for (int evenodd = 0; evenodd < 2; evenodd++)
          {
            int cx = (2 * orient + evenodd) * maxd + maxd / 2;
            int cy = scale * maxd + maxd / 2;
            const Util::Image filt = getFilter(getIndexElong (scale, orient, evenodd));
            drawFilter (filt, cx, cy, fbim);
          }
        }
        int cx = 2 * norient * maxd + maxd / 2;
        int cy = scale * maxd + maxd / 2;
        const Util::Image filt = getFilter(getIndexCenter (scale));
        drawFilter (filt, cx, cy, fbim);
      }
  }

} //namespace Group