ic.cc

用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵

// Copyright (C) 2002 Charless C. Fowlkes <[e-mail]>
// Copyright (C) 2002 David R. Martin <[e-mail]>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
// 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html.


#include <assert.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <stdio.h>
#include <ctype.h>
#include <memory.h>

#include "util.hh"
#include "image.hh"
#include "array.hh"

#include "pb.hh"
#include "ic.hh"

namespace Group
{
  //
  // given a pb image and window radius, computes a support map for each
  // pixel out to the given radius.
  //
  void computeSupport(const DualLattice& boundaries, const int wr,
                       const float thresh, SupportMap& support)
  {
    support.resize(boundaries.width,boundaries.height);
    Util::Array1D<PointIC> adj;
    int count = 0;

    Util::Message::startBlock(boundaries.width,"computing support");
    for (int x = 0; x < boundaries.width; x++)
    {
      Util::Message::stepBlock();
      for (int y = 0; y < boundaries.height; y++)
      {
        interveningContour(boundaries,thresh,x,y,wr,adj,count);    
        Util::Array1D<PointIC> map(count);
        for (int i = 0; i < count; i++)
        {
          map(i) = adj(i);    
          const int ix = map(i).x;
          const int iy = map(i).y;
          assert(ix >= 0);
          assert(ix < boundaries.width);
          assert(iy >= 0);
          assert(iy < boundaries.height);
        }
        support(x,y) = map;
      }
    }
    Util::Message::endBlock();
  }

  //
  // Walk the bresenham line from (x0,y0) to (x2,y2) ignoring any points outside
  // a circular window of radius wr. 
  // (x1,y1) and (x3,y3) should be on either side of (x2,y2)
  //
  // For each line, stop if the max-accumulated pb is > thresh. 
  //
  // Append any points on the line (for which the line is the best approximant
  // and distance from x0 is less than wr) to scanlines array.
  //
  // points is preallocated scratch space.
  // scanCount and scanLines store the results
  //
  void
  ic_walk (const DualLattice& boundaries, const float thresh, 
           const int x0, const int y0, 
           const int x1, const int y1,
           const int x2, const int y2,
           const int x3, const int y3,
           const int wr,
           Util::Array1D <PointIC> &points, 
           Util::Array1D <int>&scanCount,
           Util::Array2D <PointIC> &scanLines)
  {
      const int width = boundaries.width;
      const int height = boundaries.height;

      // the predicate that uses long longs will overflow if the image
      // is too large
      assert (2*wr+1 < 1000);

      // make sure points array is big enough
      assert ((int) points.size () >= 4*wr+2);

      // make sure scan arrays are the right size
      assert ((int)scanCount.size() == 2*wr+1);
      assert ((int)scanLines.size(0) == 2*wr+1);
      assert ((int)scanLines.size(1) == 2*wr+1);

      //sanity check the points
      assert (x0 >= 0 && x0 < width);
      assert (y0 >= 0 && y0 < height);
      assert (x1 >= 0 && x1 < width);
      assert (y1 >= 0 && y1 < height);
      assert (x2 >= 0 && x2 < width);
      assert (y2 >= 0 && y2 < height);
      assert (x3 >= 0 && x3 < width);
      assert (y3 >= 0 && y3 < height);

      // make sure points are all distinct
      assert (x0 != x1 || y0 != y1);
      assert (x0 != x2 || y0 != y2);
      assert (x0 != x3 || y0 != y3);
      assert (x1 != x2 || y1 != y2);
      assert (x1 != x3 || y1 != y3);
      assert (x2 != x3 || y2 != y3);

      //constants used in testing whether this is
      //the best path
      const long long dx1 = x1 - x0; 
      const long long dy1 = y1 - y0;
      const long long dx2 = x2 - x0;
      const long long dy2 = y2 - y0;
      const long long dx3 = x3 - x0;
      const long long dy3 = y3 - y0;
      const long long dot11 = dx1 * dx1 + dy1 * dy1;
      const long long dot22 = dx2 * dx2 + dy2 * dy2;
      const long long dot33 = dx3 * dx3 + dy3 * dy3;

      // compute dx,dy for the bresenham line
      const int dx = x2 - x0;
      const int dy = y2 - y0;
      const int adx = abs (dx);
      const int ady = abs (dy);

      // figure out what octant we're in for the bresenham algorithm;
      // octant i covers pi/4 * [i,i+1)
      int octant = -1;
      if (dx > 0 && dy >= 0) {           // quadrant 0
        octant = (adx > ady) ? 0 : 1;
      } else if (dx <= 0 && dy > 0) {    // quadrant 1
        octant = (adx < ady) ? 2 : 3;
      } else if (dy <= 0 && dx < 0) {    // quadrant 2
        octant = (adx > ady) ? 4 : 5;
      } else if (dx >= 0 && dy < 0) {    // quadrant 3
        octant = (adx < ady) ? 6 : 7;
      } else {
        assert (0);
      }

      // t is our bresenham counter
      int t = 0;
      switch (octant)
      {
        case 0: t = -adx; break;
        case 1: t = -ady; break;
        case 2: t = -ady; break;
        case 3: t = -adx; break;
        case 4: t = -adx; break;
        case 5: t = -ady; break;
        case 6: t = -ady; break;
        case 7: t = -adx; break;
        default: assert (0);
      }

      // maxpb contains the max-accumulation of pb from (x0,y0) to (x,y)
      // on the bresenham line. 
      float maxpb = 0.0f; 

      // (xi,yi) is our current location on the bresenham line
      int xi = x0;
      int yi = y0;

      // accumulate the points in the order we find them
      int count = 0;
      int oldx = xi;
      int oldy = yi;

      //walk the line
      while (xi != x2 || yi != y2)
      {
        // step one pixel on the bresenham line
        switch (octant)
        {
          case 0:
            xi++; t += (ady << 1);
            if (t > 0) { yi++; t -= (adx << 1); }
            break;
          case 1:
            yi++; t += (adx << 1);
            if (t > 0) { xi++; t -= (ady << 1); }
            break;
          case 2:
            yi++; t += (adx << 1);
            if (t > 0) { xi--; t -= (ady << 1); }
            break;
          case 3:
            xi--; t += (ady << 1);
            if (t > 0) { yi++; t -= (adx << 1); }
            break;
          case 4:
            xi--; t += (ady << 1);
            if (t > 0) { yi--; t -= (adx << 1); }
            break;
          case 5:
            yi--; t += (adx << 1);
            if (t > 0) { xi--; t -= (ady << 1); }
            break;
          case 6:
            yi--; t += (adx << 1);
            if (t > 0) { xi++; t -= (ady << 1); }
            break;
          case 7:
            xi++; t += (ady << 1);
            if (t > 0) { yi--; t -= (adx << 1); }
            break;
          default:
            assert (0);
        }

        // Figure out if the bresenham line from (x0,y0) to (x2,y2) is the
        // best approximant we will see for the line from (x0,y0) to (xi,yi).
        // We need:
        //              T(i,2) < T(i,1) && T(i,2) <= T(i,3)
        // Where T(a,b) is the angle between the two lines (x0,y0)-(xa,ya)
        // and (x0,y0)-(xb,yb).  
        // We can compute an exact integer predicate; let C be the square
        // of the cosine of T:
        //              C(i,2) > C(i,1) && C(i,2) >= C(i,3)
        // Use the identity:
        //              cos(t) = a.b/|a||b|
        // Square and cross-multiply to get rid of the divides and square
        // roots.
        // Note that we first check to see if T(i,2) == 0, in which case 
        // the line is a perfect approximant.

        const long long dxi = xi - x0;
        const long long dyi = yi - y0;
        const long long dotii = dxi * dxi + dyi * dyi;
        const long long doti1 = dxi * dx1 + dyi * dy1;
        const long long doti2 = dxi * dx2 + dyi * dy2;
        const long long doti3 = dxi * dx3 + dyi * dy3;


        const bool good = (doti2*doti2 == dotii*dot22)
                          || (dot11*doti2*doti2 > dot22*doti1*doti1
                              && dot33*doti2*doti2 >= dot22*doti3*doti3);


        // otherwise accumulate the pb value if we've crossed an edge
        float intersected = 0.0f;
        if (oldx == xi)
        {
          if (yi > oldy)
          {
            intersected = boundaries.H(xi,yi);  
          }
          else if (yi < oldy)
          {
            intersected = boundaries.H(oldx,oldy);  
          }
        } 
        else if (oldy == yi)
        {
          if (xi > oldx)
          {
            intersected = boundaries.V(xi,yi);  
          }
          else if (xi < oldx)
          {
            intersected = boundaries.V(oldx,oldy);  
          }
        }
        else
        {
          if ((xi > oldx) && (yi > oldy))  //down to right
          {
            intersected = Util::max(boundaries.H(oldx,yi),intersected); 
            intersected = Util::max(boundaries.H(xi,yi),intersected); 
            intersected = Util::max(boundaries.V(xi,oldy),intersected); 
            intersected = Util::max(boundaries.V(xi,yi),intersected); 
          }
          else if ((xi > oldx) && (yi < oldy)) //up to right
          {
            intersected = Util::max(boundaries.H(oldx,oldy),intersected); 
            intersected = Util::max(boundaries.H(xi,oldy),intersected); 
            intersected = Util::max(boundaries.V(xi,oldy),intersected); 
            intersected = Util::max(boundaries.V(xi,yi),intersected);