nonmax.cc

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

      r = Util::max(r,1.5f); 

      //compute the best paraboloid fit at each point (x,y)
      for (int x = 0; x < width; x++)
      {
        for (int y = 0; y < height; y++)
        {
              float eta = theta(x,y); 
              float maxe = mag(x,y);
              assert (eta >= 0 && eta < M_PI);

              // save the maximal orientation and energy at all pixels
              parabs.maxo(x,y) = eta;
              parabs.maxe(x,y) = maxe;

              // these are the least squares paraboloid coeff. estimates
              // for z = a*d^2 + b*d + c, where d = -x*sint + y*cost
              // error is the normalized fit error
              float a = 0, b = 0, c = 0;
              float error = 0;

              fitParabolaLeastSqr ( x, y, eta, r, mag, a, b, c, error);

              // compute various parameters of fitted parabola
              // curvature is normalized to the window radius
              const float curvature = a * r * r;
              const float dist = -b / (2*a);
              const float energy = a*dist*dist + b*dist + c;
              const float smooth = Util::max(c,0.0f);
              const float ycoord = dist*cos(eta);
              const float xcoord = -dist*sin(eta);

              // save fit info for all putative needles
              parabs.X(x,y) = xcoord;
              parabs.Y(x,y) = ycoord;
              parabs.energy(x,y) = energy;
              parabs.orientation(x,y) = eta;
              parabs.curvature(x,y) = curvature;
              parabs.fiterror(x,y) = error;
              parabs.dist(x,y) = dist;
              parabs.smooth(x,y) = smooth;

              // if any of the values are infinite or NaN, then skip
              if (!finite(xcoord)) { continue; }
              if (!finite(ycoord)) { continue; }
              if (!finite(energy)) { continue; }
              if (!finite(eta)) { continue; }
              if (!finite(curvature)) { continue; }
              if (!finite(error)) { continue; }

              // if curvature is positive, then this is not a max; skip
              if (curvature > -1e-10) { continue; }

              // if the center is out of this pixel's neighborhood, then skip
              if (fabs(dist) > 1) { continue; }

              // only mark valid those needles that are maxima in the
              // pixel's immediate neighborhood
              parabs.edgemap(x,y) = 1;
          }
      }
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // interpolate the orientation matrix and energies,
  // fit clyndrical parabolas to each and return the resulting subpixel info
  // r = radius of square to which we fit the cylindrical parabola
  //
  void fitSubpixelParabolas (const Util::ImageStack& mag, float r, ParabolaSet& parabs)
  {
    r = Util::max(r,1.5f); 
    const int norient = mag.size(0);
    const int width = mag.size(1);
    const int height = mag.size(2);

    parabs.resize(width,height);
    parabs.init();

    //compute the best paraboloid fit at each point (x,y)
    for (int x = 0; x < width; x++)
    {
      for (int y = 0; y < height; y++)
      {
        // compute the max energy orientation at this point
        int maxOrient = 0;
        float maxEnergy = mag(0,x,y);
        for (int orient = 0; orient < norient; orient++) {
            const float e = mag(orient,x,y);
            if (e > maxEnergy) {
                maxOrient = orient;
                maxEnergy = e;
            }
        }
        const float oStep = M_PI / norient;
        float theta = maxOrient * oStep;
        float maxe = mag(maxOrient,x,y);
        assert (theta >= 0 && theta < M_PI);

        if (norient > 2) {
          // interpolate around the max energy point to get a more
          // accurate theta; we will use this orientation for the
          // parabola fit
          const float x1 = theta - oStep;
          const float x2 = theta;
          const float x3 = theta + oStep;
          const float y1 = mag(Util::wrap(maxOrient-1,norient),x,y);
          const float y2 = mag(maxOrient,x,y);
          const float y3 = mag(Util::wrap(maxOrient+1,norient),x,y);
          float maxo,maxeFit;
          interpolateMaxLagrange (x1, y1, x2, y2, x3, y3, maxo, maxeFit);
          assert(finite(maxo));
          assert(finite(maxeFit));
          if (fabs(theta-maxo) < oStep) 
          {
            theta = maxo;
            maxe = maxeFit;
            if (theta < 0) { theta += M_PI; }
            if (theta >= M_PI) { theta -= M_PI; }
            if (theta < 0) { theta = 0; }
            if (theta >= M_PI) { theta = 0; }
            assert (theta >= 0 && theta < M_PI);
          } // else bad numerical failure!
        }
        assert (theta >= 0 && theta < M_PI);

        // save the interpolated maximal orientation and energy 
        // at all pixels
        parabs.maxo(x,y) = theta;
        parabs.maxe(x,y) = maxe;

        // these are the least squares paraboloid coeff. estimates
        // for z = a*d^2 + b*d + c, where d = -x*sint + y*cost
        // error is the normalized fit error
        float a = 0, b = 0, c = 0;
        float error = 0;

        fitParabolaLeastSqr ( x, y, theta, r, *mag.slice(maxOrient), a, b, c, error);

        // compute various parameters of fitted parabola
        // curvature is normalized to the window radius
        const float curvature = a * r * r;
        const float dist = -b / (2*a);
        const float energy = a*dist*dist + b*dist + c;
        const float smooth = Util::max(c,0.0f);
        const float ycoord = dist*cos(theta);
        const float xcoord = -dist*sin(theta);

        // save fit info for all putative needles
        parabs.X(x,y) = xcoord;
        parabs.Y(x,y) = ycoord;
        parabs.energy(x,y) = energy;
        parabs.orientation(x,y) = theta;
        parabs.curvature(x,y) = curvature;
        parabs.fiterror(x,y) = error;
        parabs.dist(x,y) = dist;
        parabs.smooth(x,y) = smooth;

        // if any of the values are infinite or NaN, then skip
        if (!finite(xcoord)) { continue; }
        if (!finite(ycoord)) { continue; }
        if (!finite(energy)) { continue; }
        if (!finite(theta)) { continue; }
        if (!finite(curvature)) { continue; }
        if (!finite(error)) { continue; }

        // if curvature is positive, then this is not a max; skip
        if (curvature > -1e-10) { continue; }

        // if the center is out of this pixel's neighborhood, then skip
        if (fabs(dist) > 1) { continue; }

        // only mark valid those needles that are maxima in the
        // pixel's immediate neighborhood
        parabs.edgemap(x,y) = 1;
      }
    }
  }


  //////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // project an edgel with energy pb onto the dual lattice by tracing the
  // breshenham line that lies underneath it.
  // 
  void traceBoundary(const float pb, const float ex1, const float ey1, 
                     const float ex2, const float ey2, DualLattice& L)
  {
      const int width = L.H.size(0);
      const int height = L.V.size(1);

      int dx = (int)rintf(ex2) - (int)rintf(ex1);
      int dy = (int)rintf(ey2) - (int)rintf(ey1);
      float x1 = 0;
      float y1 = 0;
      float x2 = 0;
      float y2 = 0;
      if (dx > 0)
      {
        x1 = ex1; 
        y1 = ey1; 
        x2 = ex2; 
        y2 = ey2; 
      }
      else
      {
        x1 = ex2;
        y1 = ey2;
        x2 = ex1;
        y2 = ey1;
        dx = -dx;
        dy = -dy;
      }

      const int steps = Util::max(abs(dx),abs(dy));
      if (steps == 0) { return; }
                                                                                                                 
      const float xincr = (x2-x1) / (float) steps;
      const float yincr = (y2-y1) / (float) steps;
                                                                                                                 
      float x = x1;
      float y = y1;
      int oldx = (int)rint(x);
      int oldy = (int)rint(y);
      for (int k = 0; k < steps; k++)
      {
        float mx = x + 0.5*xincr; 
        float lmx = mx - floor(mx);
        float my = y + 0.5*yincr; 
        float lmy = my - floor(my);

        x += xincr;
        y += yincr;
        const int xi = (int) rintf(x);
        const int yi = (int) rintf(y);

        //don't try to accumulate out of bounds
        if ((oldx < 0) || (oldx >= width)) {continue;}
        if ((xi < 0) || (xi >= width)) {continue;}
        if ((oldy < 0) || (oldy >= height)) {continue;}
        if ((yi < 0) || (yi >= height)) {continue;}

        if ((oldx != xi) && (oldy != yi))
        {
          //diagonal move 
          if (dy > 0)
          {
            if ( (lmx-lmy) > 0 )
            {
              L.H(oldx,oldy) = pb;
              L.V(xi,oldy) = pb;
            }
            else
            {
              L.V(oldx,oldy) = pb;
              L.H(oldx,yi) = pb;
            }
          }
          else
          {
            if ( (lmx+lmy) > 1 )
            {
              L.H(oldx,oldy) = pb;
              L.V(xi,yi) = pb;
            }
            else
            {
              L.V(oldx,yi) = pb;
              L.H(oldx,yi) = pb;
            }
          }
        }
        else if (oldy != yi)
        {
          //vertical move
          if (dy > 0)
          {
            L.V(oldx,oldy) = pb;    
          }
          else
          {
            L.V(oldx,yi) = pb;    
          }
        }
        else if (oldx != xi)
        {
          //horizontal move
          L.H(oldx,oldy) = pb;    
        }
        else
        {
          //TODO: this is some bug.  testcase: 188005.jpg
          fprintf(stderr,"no move!! %d %d %3.3f %3.3f\n",dx,dy,xincr,yincr);
        }
        oldx = xi;
        oldy = yi;
      }

      assert(fabs(x-x2)<0.01);
      assert(fabs(y-y2)<0.01);
  }
                                                                                                                 

  //
  // given a set of parabolas, intersect them all with the dual lattice
  //
  void intersectPixels (const ParabolaSet& parabolas, DualLattice& L)
  {
      const float edgelLength = 1.1f;
      const int width = parabolas.edgemap.size(0);
      const int height = parabolas.edgemap.size(1);
      L.width = width;
      L.height = height;
      L.H.resize(width,height+1);
      L.V.resize(width+1,height);
      L.H.init(0);
      L.V.init(0);

      for (int x = 0; x < width; x++)
      {
        for (int y = 0; y < height; y++)
        {
          if (parabolas.edgemap(x,y) != 1) 
          {
            //nearest edgel is more than 1 pixel away
            //from my center
            continue; 
          }

          // edgel parameters
          //const float dval = parabolas.dist(x,y);
          const float energy = parabolas.energy(x,y);
          const float theta = parabolas.orientation(x,y);
          const float dx = 0.5*edgelLength*cos(theta);
          const float dy = 0.5*edgelLength*sin(theta);
          float x1 = (float)x + parabolas.X(x,y) + dx;
          float x2 = (float)x + parabolas.X(x,y) - dx;
          float y1 = (float)y + parabolas.Y(x,y) + dy;
          float y2 = (float)y + parabolas.Y(x,y) - dy;
          traceBoundary(energy,x1,y1,x2,y2,L);
        }
      }
  }


  //
  // given the subpixel parabolas and the orientation fits, this function
  // intersects edges with the pixels: one pixel per parabola
  //
  void intersectPixels (const ParabolaSet& parabolas, Util::Image& pb)
  {
      const int width = parabolas.edgemap.size(0);
      const int height = parabolas.edgemap.size(1);
      pb.resize(width,height);
      pb.init(0);
      for (int x = 0; x < width; x++)
      {
          for (int y = 0; y < height; y++)
          {
              if (parabolas.edgemap(x,y) != 1) 
              {
                continue; 
              }
              float dval = parabolas.dist(x,y);
              const float dx = parabolas.X(x,y);
              const float dy = parabolas.Y(x,y);
              int xi = (int) rint (x + dx);
              int yi = (int) rint (y + dy);
              if (xi < 0 || xi >= width) 
              {
                continue; 
              }
              if (yi < 0 || yi >= height) 
              { 
                continue; 
              }

              if (dval < (sqrt(2)/2))
              {
                float energy = parabolas.energy(x,y);
                energy = Util::max(energy,0.0f);
                energy = Util::min(energy,1.0f);
                pb(xi,yi) = Util::max(pb(xi,yi),energy);
              }
          }
      }
  }
  

  //
  // estimate smoothed derivatives using LOWESS style smoothing
  //
  void lowess(const float theta, const float r, const Util::Image& signal, 
              Util::Image& d0, Util::Image& d1, Util::Image& d2)

  {
    const int width = signal.size(0);
    const int height = signal.size(1);
    d0.resize(width,height);
    d1.resize(width,height);
    d2.resize(width,height);
    for(int x = 0; x < width; x++)
    {
      for(int y = 0; y < height; y++)
      {
        float a, b, c, error;
        fitParabolaLeastSqr(x, y, theta, r, signal, a, b, c, error);
        d2(x,y) = 2*a;
        d1(x,y) = b;
        d0(x,y) = c;
      }
    }
   
  }


  //
  // estimate smoothed 0th,1st,2nd derivatives using gaussian derivative convolution
  // gaussian sigma = r / 3
  //
  void derivs(const float theta, const float r, const Util::Image& signal, 
               Util::Image& d0, Util::Image& d1, Util::Image& d2)
  {
    const float sigma = (sqrt(2)*r) / 3.0;

    Util::Image d0Filt,d1Filt,d2Filt;
    FilterBank::createGaussianKernel(sigma, sigma, 3.0, theta, 0, false, d0Filt);
    FilterBank::createGaussianKernel(sigma, sigma, 3.0, theta, 1, false, d1Filt);
    FilterBank::createGaussianKernel(sigma, sigma, 3.0, theta, 2, false, d2Filt);

    getFiltered(signal,d0Filt,d0);
    getFiltered(signal,d1Filt,d1);
    getFiltered(signal,d2Filt,d2);
  }

} //namespace Group