void GaussSmoothSep(imatrix& image, double sigma)
{
	int	i, j;
	double g, max_g, min_g;
	int s, t;
	int x, y;
	double weight, w_sum;

	int image_x = image.getRow();
	int image_y = image.getCol();

	myvec GAU1;
	MakeGaussianVector(sigma, GAU1); 
	int half = GAU1.getMax()-1;

	mymatrix tmp(image_x, image_y);
		
	max_g = -1;
	min_g = 10000000;
	for (j = 0; j < image_y; j++) {
		for (i = 0; i < image_x; i++) {
			g = 0.0;
			weight = w_sum = 0.0;
			for (s = -half; s <= half; s++) {
				x = i+s; y = j;
				if (x > image_x-1) x = image_x-1;
				else if (x < 0) x = 0;
				if (y > image_y-1) y = image_y-1;
				else if (y < 0) y = 0;
				weight = GAU1[ABS(s)];
				g += weight * image[x][y];
				w_sum += weight;
			}
			g /= w_sum;
			if (g > max_g) max_g = g;
			if (g < min_g) min_g = g;
			tmp[i][j] = g;
		}
	}
	for (j = 0; j < image_y; j++) {
		for (i = 0; i < image_x; i++) {
			g = 0.0;
			weight = w_sum = 0.0;
			for (t = -half; t <= half; t++) {
					x = i; y = j+t;
					if (x > image_x-1) x = image_x-1;
					else if (x < 0) x = 0;
					if (y > image_y-1) y = image_y-1;
					else if (y < 0) y = 0;
					weight = GAU1[ABS(t)];
					g += weight * tmp[x][y];
					w_sum += weight;
			}
			g /= w_sum;
			if (g > max_g) max_g = g;
			if (g < min_g) min_g = g;
			image[i][j] = round(g);
		}
	}
}

  //get most(best) available label from histogram
  int kNearestNeighFilter::getMostLabel(const ivector& histogram,
			       const imatrix& src,
			       const int& row, const int& col) const{

    int numOfMax = 0;
    int maxIndex = -1; // first index, which is max
    int max = 0; //
    for(int i=0;i<histoSize;++i) {
      if(histogram.at(i) < max); // for speed up (probability)
      else if(histogram.at(i) > max) {
	max = histogram.at(i);
	numOfMax = 1;
	maxIndex = i;
      }
      else //if(histogram.at(i) == max)
	++numOfMax;
    }

    //is there more than one possibility ?
    if (numOfMax == 1)
      return maxIndex;
    // is the kernel center one of the max's?
    else if(histogram.at(src.at(row,col)) == max)
      return src.at(row,col);
    else
      return getMedian(histogram,max,numOfMax);
  };

  // applies the histogramMethod for the type boundary NoBoundary
  bool kNearestNeighFilter::histogramMethodNoBoundary(const imatrix& src,
					     imatrix& dest) const {
    const int limit = sizeOfKernel/2; //half size of the kernel
    int i,j,row,col;

    ivector histogram(histoSize,0);
    dest.resize(src.size(),0,false,false);

    //runs through the src's columns,inside the image
    const int lastRow = src.lastRow()-limit;
    for(row=limit;row<=lastRow;++row) {
      histogram.fill(0);

      // first block per row
      const int r = row+limit;
      for(i=row-limit;i<=r;++i)
	for(j=0;j<sizeOfKernel;++j)
          ++histogram.at(src.at(i,j));
      dest.at(row,limit)=getMostLabel(histogram,src,row,limit);

      // runs inside the image
      histogramMethodMiddle(src,dest,histogram,row,col);
    }
    return true;
 };

int writeIntVecs(imatrix & im, string fname, int buffsize) {
  int i, j, nrows, fmt, ncols, nnz, v;
  ostream *ofstr = open_out_buf(fname.c_str(), buffsize);
  ncols = 2;
  nnz = 0;
  for (i = 0, nrows = 0; i < im.size(); i++) nrows += im[i].size();
  fmt = 110;
  ofstr->write((const char *)&fmt, 4);
  ofstr->write((const char *)&nrows, 4);
  ofstr->write((const char *)&ncols, 4);
  ofstr->write((const char *)&nnz, 4);
  for (i = 0; i < im.size(); i++) {
    for (j = 0; j < im[i].size(); j++) {
      ofstr->write((const char *)&i, 4);
    }
  }
  for (i = 0; i < im.size(); i++) {
    ivector & ivv = im[i];
    for (j = 0; j < ivv.size(); j++) {
      v = ivv[j];
      ofstr->write((const char *)&v, 4);
    }
  }
  closeos(ofstr);
  return 0;
}

void GetFDoG(imatrix& image, ETF& e, double sigma, double sigma3, double tau) 
{
	int	i, j;

    int image_x = image.getRow();
	int image_y = image.getCol();

	myvec GAU1, GAU2, GAU3;
	MakeGaussianVector(sigma, GAU1); 
	MakeGaussianVector(sigma*1.6, GAU2); 

	int half_w1, half_w2, half_l;
	half_w1 = GAU1.getMax()-1;
	half_w2 = GAU2.getMax()-1;
	
	MakeGaussianVector(sigma3, GAU3); 
	half_l = GAU3.getMax()-1;
	
	mymatrix tmp(image_x, image_y);
	mymatrix dog(image_x, image_y);

	GetDirectionalDoG(image, e, dog, GAU1, GAU2, tau);
	GetFlowDoG(e, dog, tmp, GAU3);

	for (i = 0; i < image_x; i++) { 
		for (j = 0; j < image_y; j++) {
			image[i][j] = round(tmp[i][j] * 255.);
		}
	}
}

/**
 * Description not yet available.
 * \param
 */
int min(const imatrix & M)
{
  int mmin=M.indexmin();
  int mmax=M.indexmax();
  int n1=min(M(mmin));
  for (int i=mmin+1;i<=mmax;i++)
  {
    int n=min(M(i));
    if (n1>n) n1=n;
  }
  return n1;
}

void ConstructMergedImage(imatrix& image, imatrix& gray, imatrix& merged) 
{
	int x, y;

	int image_x = image.getRow();
	int image_y = image.getCol();

	for (y = 0; y < image_y; y++) {
		for (x = 0; x < image_x; x++) {
			if (gray[x][y] == 0) merged[x][y] = 0;
			else merged[x][y] = image[x][y];
		}
	}
}

  // the kernel runs inside the image
  void kNearestNeighFilter::histogramMethodMiddle(const imatrix& src,
                                              imatrix& dest,
                                              ivector& histogram,
                                              const int& row,int& col) const {
    int i,j;//index
    int numOfMax, maxIndex;
    int max=0;
    const int maxChange = sizeOfKernel+1;//max change for "max"

    const int limit = sizeOfKernel/2;   //half size of the kernel
    const int lastCol = src.lastColumn()-limit;
    const int r = row+limit;
    col = limit;

    int v; //del test
    while(col <= (lastCol-1)) {
      j = col-limit;
      // sub labels left form the kernel
      for(i=row-limit;i<=r;++i) {
        --histogram.at(src.at(i,j));
      }
      // add labels right from the kernel
      ++col;
      j = col+limit;
      for(i=row-limit;i<=r;++i) {
	v = src.at(i,j);
        ++histogram.at(src.at(i,j));
      }

      //get most(best) available label
      numOfMax = 0;
      maxIndex = -1;
      max -= maxChange; //=0;
      for(i=0;i<histoSize;++i) {
	if(histogram.at(i) < max);// for speed up (probability)
	else if(histogram.at(i) > max) {
	  max = histogram.at(i);
	  numOfMax = 1;
	  maxIndex = i;
	}
	else //if(histogram.at(i) == max)
	  ++numOfMax;
      }

      //is there more than one possibility ?
      if(numOfMax == 1)
	dest.at(row,col) =  maxIndex;
      // is the kernel center one of the max's?
      else if(histogram.at(src.at(row,col)) == max)
	dest.at(row,col) = src.at(row,col);
      else
	dest.at(row,col) = getMedian(histogram,max,numOfMax);
    }//while
  };

void GrayThresholding(imatrix& image, double thres) 
{
	int	i, j;
	double val;

	int image_x = image.getRow();
	int image_y = image.getCol();
	
	for (i = 0; i < image_x; i++) { 
		for (j = 0; j < image_y; j++) {
			val = image[i][j] / 255.0; 
			if (val < thres)
				image[i][j] = round(val * 255.0); 
			else image[i][j] = 255; 
		}
	}
}

void Binarize(imatrix& image, double thres) 
{
	int	i, j;
	double val;

	int image_x = image.getRow();
	int image_y = image.getCol();
	
	for (i = 0; i < image_x; i++) { 
		for (j = 0; j < image_y; j++) {
			val = image[i][j] / 255.0; 
			if (val < thres)
				image[i][j] = 0; 
			else image[i][j] = 255; 
		}
	}
}

void ConstructMergedImageMult(imatrix& image, imatrix& gray, imatrix& merged) 
// using multiplication
{
	int x, y;
	double gray_val, line_darkness;

	int image_x = image.getRow();
	int image_y = image.getCol();

	for (y = 0; y < image_y; y++) {
		for (x = 0; x < image_x; x++) {
			gray_val = image[x][y] / 255.0; 
			line_darkness = gray[x][y] / 255.0; 
			gray_val *= line_darkness; 
			merged[x][y] = round(gray_val * 255.0);
		}
	}
}

/**
 * Description not yet available.
 * \param
 */
  int sub_unallocated(const imatrix& m)
  {
    int iflag=0;
    int mmin=m.indexmin();
    int mmax=m.indexmax();
    if (!allocated(m))
    {
      iflag=1;
      return iflag;
    }
    for (int i=mmin;i<=mmax;i++)
    {
      if (!allocated(m(i)))
      {
        iflag=1;
        break;
      }
    }
    return iflag;
  }

  bool kNearestNeighFilter::apply(const imatrix& src,imatrix& dest) {

    if (src.empty()) {
      dest.clear();
      return true;
    }  

    const parameters& param = getParameters();
    if (param.kernelSize <= 1) {
      dest.copy(src);
      return true;
    }
    sizeOfKernel = param.kernelSize + ((param.kernelSize%2 == 0) ? 1 : 0);
    histoSize = src.maximum()+1;
    bool control = false;   // return variable

    // choose border behaviour
    switch(param.boundaryType) {
      case lti::Zero:
        control = histogramMethodZero(src,dest);
        break;
      case lti::Mirror:
        control = histogramMethodMirror(src,dest);
        break;
      case lti::Constant:
        control = histogramMethodConstant(src,dest);
        break;
      case lti::Periodic:
        control = histogramMethodPeriodic(src,dest);
        break;
      case lti::NoBoundary:
        control = histogramMethodNoBoundary(src,dest);
        break;
      default:
        setStatusString("Unknown boundary type");
        break;
    }

    return control;
  };

  bool regionMerge::apply(const imatrix& srcmask,
                          const dmatrix& simMat,
                          const dvector& thresholds,
                          imatrix& destmask) const {
    int i,j;
    const dvector& thrs = thresholds;

    ivector eLab(simMat.rows());
    for (i=0;i<eLab.size();++i)
      eLab.at(i) = i;

    double a;
    for (j=0;j<simMat.rows();++j)
      for (i=0;i<simMat.columns();++i) {
        if (simMat.at(j,i) > (a=max(thrs.at(j),thrs.at(i))) ) {
//          cout<<j<<" "<<i<<" "<<simMat.at(j,i)<<" "<<a<<endl;
          if (eLab.at(j)>eLab.at(i)) {
            eLab.at(j)=eLab.at(i);
          } else {
            eLab.at(i)=eLab.at(j);
          }

        }
      }


    // now correct the labels
    for (j=eLab.size()-1;j>=0;--j) { //todo
      i = j;
      while (eLab.at(i) != i) {
        i=eLab.at(i);
      }
      eLab.at(j) = i;
    }

    destmask.resize(srcmask.size(),0,false,false);
    for (j=0;j<srcmask.rows();++j)
      for (i=0;i<srcmask.columns();++i) {
        destmask.at(j,i) = eLab.at(srcmask.at(j,i));
      }

    return true;
 };

   // applies the histogramMethod for the type boundary Constant
  bool kNearestNeighFilter::histogramMethodConstant(const imatrix& src,
					   imatrix& dest) const {

    // image is divided in 9 areas, which are calc one by one
    //  1|  2  |3
    //  ---------
    //   |     |
    //  4|  5  |6
    //   |     |
    // ----------
    //  7|  8  |9

    int i,j,row,col;//index

    ivector histogram(histoSize,0);
    const int rowSize = src.rows();
    const int columnSize = src.columns();
    const int limit = sizeOfKernel/2;   //half size of the kernel
    const int lastCol = src.lastColumn()-limit;
    const int lastRow = src.lastRow()-limit;
    dest.resize(src.size(),0,false,false);

    //runs through the src's columns
    // (area:4,5,6) only kernels with full kernel-height
    for(row=limit;row<=lastRow;++row) {
      histogram.fill(0);

      // first full kernel (area4)
      col=0;
      const int r = row+limit;
      const int c = col+limit;
      for(i=row-limit;i<=r;++i)
	for(j=col-limit;j<=c;++j) { //todo better
	  if(j<0)
            ++histogram.at(src.at(i,0));
	  else// if(j>=0)
            ++histogram.at(src.at(i,j));
        }
      dest.at(row,col)=getMostLabel(histogram,src,row,col);

      // rest (area 4)
      while(col < limit) {
        for(i=row-limit;i<=r;++i)
          --histogram.at(src.at(i,0));
        ++col;
        j=col+limit;
        for(i=row-limit;i<=r;++i)
          ++histogram.at(src.at(i,j));
	dest.at(row,col)=getMostLabel(histogram,src,row,col);
      }

      // area 5
      histogramMethodMiddle(src,dest,histogram,row,col);

      // area 6
      col=lastCol;
      while(col < (columnSize-1)) {
        j=col-limit;
        for(i=row-limit;i<=r;++i)
          --histogram.at(src.at(i,j));
        ++col;
	j = columnSize-1;
	for(i=row-limit;i<=r;++i)
          ++histogram.at(src.at(i,j));
	dest.at(row,col)=getMostLabel(histogram,src,row,col);
      }
    } // area 4,5,6

    // areas 1,2,3,7,8,9
    for(row=0;row<limit;++row) {     //runs top rows (1,2,3)
      const int r=row+limit;

      // runs middle top rows (area 2)
      for(col=limit;col<=lastCol;++col) {
	histogram.fill(0);
        const int c=col+limit;
        for(i=row-limit;i<=r;++i)
	  for(j=col-limit;j<=c;++j) {
	    if(i<0)
              ++histogram.at(src.at(0,j));
	    else//if(i>=0) {
              ++histogram.at(src.at(i,j));
          }
	dest.at(row,col)=getMostLabel(histogram,src,row,col);
      }
      // runs left top corner (area 1)
      for(col=0;col<limit;++col) {
	histogram.fill(0);
        const int c=col+limit;
        for(i=row-limit;i<=r;++i)
	  for(j=col-limit;j<=c;++j) {
            if(i<0 && j<0)
              ++histogram.at(src.at(0,0));
            else if(i>=0 && j<0)
              ++histogram.at(src.at(i,0));
            else if(i<0 && j>=0)
              ++histogram.at(src.at(0,j));
            else //if(i>=0 && j>=0)
              ++histogram.at(src.at(i,j));
//.........这里部分代码省略.........

  // applies the histogramMethod for the type boundary Zero
  bool kNearestNeighFilter::histogramMethodZero(const imatrix& src,
					imatrix& dest) const {
    int i,j,row,col;//index

    const int limit = sizeOfKernel/2;   //half size of the kernel
    const int rowSize = src.rows();
    const int columnSize = src.columns();
    const int lastCol = src.lastColumn()-limit;
    const int lastRow = src.lastRow()-limit;

    ivector histogram(histoSize,0);  //the histogram for a channel8
    dest.resize(src.size(),ubyte(0),false,false);

    //runs area 4,5,6
    for(row=limit;row<=lastRow;++row) {
      const int r = row+limit;
      histogram.fill(0);
      col=0;

      //number of 0's are known
      histogram.at(0) = sizeOfKernel*(sizeOfKernel-limit-1);
      for(i=row-limit;i<=r;++i)
        for(j=0;j<=limit;++j)
          ++histogram.at(src.at(i,j));
      dest.at(row,col) = getMostLabel(histogram,src,row,col);

      // the kernel at the position between the border and the image
      while(col < limit) {
        histogram.at(0) -= sizeOfKernel; // cut all the 0 in the leftmost column
	++col;
        j = col+limit;        // for each pixel in the rightmost column
        for(i=row-limit;i<=r;++i)
          ++histogram.at(src.at(i,j));
	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }

      //runs area 5
      histogramMethodMiddle(src,dest,histogram,row,col);

      // area 6
      col = lastCol;
      while(col < (columnSize-1)) {
        j = col-limit;
        for(i=row-limit;i<=r;++i)
          --histogram.at(src.at(i,j));
        ++col;
        histogram.at(0) += sizeOfKernel;

	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
    }

    // runs area 1,2,3
    for(row=0;row<limit;++row) {
      const int r=row+limit;
      // runs middle top rows (area 2)
      for(col=limit;col<=lastCol;++col) {
	histogram.fill(0);
        const int c=col+limit;
        histogram.at(0) = sizeOfKernel*(sizeOfKernel-(limit+1)-row);
	for(i=0;i<=r;++i)
          for(j=col-limit;j<=c;++j)
	    ++histogram.at(src.at(i,j));
        dest.at(row,col) = getMostLabel(histogram,src,row,col) ;
      }
      // runs left top corner (area 1);
      for(col=0;col<limit;++col) {
	histogram.fill(0);
	const int c=col+limit;
        histogram.at(0) = sizeOfKernel*sizeOfKernel-(limit+1)*(limit+1+row+col)-col*row;
        for(i=0;i<=r;++i)
	  for(j=0;j<=c;++j)
	    ++histogram.at(src.at(i,j));
	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
      // runs top right corner (area 3)
      for(col=lastCol+1;col<=columnSize-1;++col) {
	histogram.fill(0);
	const int c=columnSize-1;
	histogram.at(0) = sizeOfKernel*sizeOfKernel
	  -(limit+1)*(limit+1+row+(columnSize-1-col))-row*(columnSize-1-col);
        for(i=0;i<=r;++i)
          for(j=col-limit;j<=c;++j)
            ++histogram.at(src.at(i,j));
	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
    }

    //runs the rows at the bottom (area 7,8,9)
    for(row=lastRow+1;row<=rowSize-1;++row) {
      //runs middle bottom rows (area 8)
      for(col=limit;col<=lastCol;++col) {
	histogram.fill(0);
        const int c=col+limit;
        histogram.at(0) = sizeOfKernel*(sizeOfKernel-(limit+1)-(rowSize-1-row));
        for(i=row-limit;i<=rowSize-1;++i)
          for(j=col-limit;j<=c;++j)
            ++histogram.at(src.at(i,j));
	dest.at(row,col) = getMostLabel(histogram,src,row,col);
//.........这里部分代码省略.........

  // applies the histogramMethod for the type boundary Mirror
  bool kNearestNeighFilter::histogramMethodMirror(const imatrix& src,
					 imatrix& dest) const {
    int i,j,row,col;//index

    const int limit = sizeOfKernel/2;   //half size of the kernel
    const int rowSize = src.rows();
    const int columnSize = src.columns();
    const int lastCol = src.lastColumn()-limit;
    const int lastRow = src.lastRow()-limit;

    ivector histogram(histoSize,0);
    dest.resize(src.size(),0,false,false);

    //runs through the src's columns,inside the image
    for(row=limit;row<=lastRow;++row) {
      histogram.fill(0);

      col=0;
      const int c=col+limit;
      const int r = row+limit;
      for(i=row-limit;i<=r;++i)
        for(j=col-limit;j<=c;++j) {
          if(j<0)
            ++histogram.at(src.at(i,-j-1));
          else // if(j>=0)
            ++histogram.at(src.at(i,j));
	}
      dest.at(row,col) = getMostLabel(histogram,src,row,col);

      while(col < limit) {
        j = col-limit;
        for(i=row-limit;i<=r;++i)
          --histogram.at(src.at(i,-j-1));
        ++col;
        j=col+limit;
        for(i=row-limit;i<=r;++i)
	  ++histogram.at(src.at(i,j));

	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
      // runs inside the image
      histogramMethodMiddle(src,dest,histogram,row,col);

      col=lastCol;
      while(col < (columnSize-1)) {
        j = col-limit;

        for(i=row-limit;i<=r;++i)
          --histogram.at(src.at(i,j));
        ++col;
        j=col+limit;
        for(i=row-limit;i<=r;++i)
          ++histogram.at(src.at(i,columnSize-1+(columnSize-j)));
	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }

    }

    // areas 1,2,3
    for(row=0;row<limit;++row) {
      const int r=row+limit;
      // runs middle top rows (area 2)
      for(col=limit;col<=lastCol;++col) {
	histogram.fill(0);
	const int c=col+limit;
        for(i=row-limit;i<=r;++i)
	  for(j=col-limit;j<=c;++j) {
	    if(i<0)
              ++histogram.at(src.at(-i-1,j));
            else //if(i>=0)
              ++histogram.at(src.at(i,j));
	  }
	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
      // runs top left corner (area 1)
      for(col=0;col<limit;++col) {
	histogram.fill(0);
	const int c=col+limit;
        for(i=row-limit;i<=r;++i)
	  for(j=col-limit;j<=c;++j) {
            if(i<0 && j<0)
              ++histogram.at(src.at(-i-1,-j-1));
            else if(i>=0 && j<0)
              ++histogram.at(src.at(i,-j-1));
            else if(i<0 && j>=0)
              ++histogram.at(src.at(-i-1,j));
            else //if(i>=0 && j>=0)
              ++histogram.at(src.at(i,j));
          }
        dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
      //runs top right corner (area 3)
      for(col=lastCol+1;col<columnSize;++col) {
	histogram.fill(0);
	const int c=col+limit;
        for(i=row-limit;i<=r;++i)
          for(j=col-limit;j<=c;++j) {
            if(i<0 && j<=columnSize-1)
              ++histogram.at(src.at(-i-1,j));
//.........这里部分代码省略.........

void GetDirectionalDoG(imatrix& image, ETF& e, mymatrix& dog, myvec& GAU1, myvec& GAU2, double tau)
{
	myvec vn(2);
	double x, y, d_x, d_y;
	double weight1, weight2, w_sum1, sum1, sum2, w_sum2;

	int s;
	int x1, y1;
	int i, j;
	int dd;
	double val;

	int half_w1, half_w2;

	half_w1 = GAU1.getMax()-1;
	half_w2 = GAU2.getMax()-1;

	int image_x, image_y;

	image_x = image.getRow();
	image_y = image.getCol();
	
	for (i = 0; i < image_x; i++) {
		for (j = 0; j < image_y; j++) {
			sum1 = sum2 = 0.0;
			w_sum1 = w_sum2 = 0.0;
			weight1 = weight2 = 0.0;
	
			vn[0] = -e[i][j].ty;
			vn[1] = e[i][j].tx;

			if (vn[0] == 0.0 && vn[1] == 0.0) {
				sum1 = 255.0;
				sum2 = 255.0;
				dog[i][j] = sum1 - tau * sum2;
				continue;
			}
			d_x = i; d_y = j;
			////////////////////////////////////////
			for (s = -half_w2; s <= half_w2; s++) { 
				////////////////////////
				x = d_x + vn[0] * s;
				y = d_y + vn[1] * s;
				/////////////////////////////////////////////////////
				if (x > (double)image_x-1 || x < 0.0 || y > (double)image_y-1 || y < 0.0) 
					continue;
				x1 = round(x);	if (x1 < 0) x1 = 0; if (x1 > image_x-1) x1 = image_x-1;
				y1 = round(y);	if (y1 < 0) y1 = 0; if (y1 > image_y-1) y1 = image_y-1;
				val = image[x1][y1];
				/////////////////////////////////////////////////////////
				dd = ABS(s);
				if (dd > half_w1) weight1 = 0.0;
				else weight1 = GAU1[dd];
				//////////////////////////////////
				sum1 += val * weight1;
				w_sum1 += weight1;
				/////////////////////////////////////////////////////
				weight2 = GAU2[dd];
				sum2 += val * weight2;
				w_sum2 += weight2;
			}
			/////////////////////////
			sum1 /= w_sum1; 
			sum2 /= w_sum2; 
			//////////////////////////////////////
			dog[i][j] = sum1 - tau * sum2;
		}
	}

}

  // applies the histogramMethod for the type boundary Periodic
  bool kNearestNeighFilter::histogramMethodPeriodic(const imatrix& src,
					   imatrix& dest) const {
    int i,j,row,col;//index
    ivector histogram(histoSize,0);

    const int limit = sizeOfKernel/2;   //half size of the kernel
    const int rowSize = src.rows();
    const int columnSize = src.columns();
    const int lastCol = src.lastColumn()-limit;
    const int lastRow = src.lastRow()-limit;

    dest.resize(src.size(),0,false,false);

    //runs through the src's columns area 4,5,6
    for(row=limit;row<=lastRow;++row) {
      histogram.fill(0);

      col=0;
      const int c = col+limit;
      const int r = row+limit;
      for(i=row-limit;i<=r;++i)
	for(j=col-limit;j<=c;++j) {
          if(j<0)
            ++histogram.at(src.at(i,j+columnSize));
          else // if(j>=0)
	    ++histogram.at(src.at(i,j));
        }
      dest.at(row,col)=getMostLabel(histogram,src,row,col);

      while(col < limit) {
	j = col-limit;
        for(i=row-limit;i<=r;++i)
          --histogram.at(src.at(i,j+columnSize));
	++col;
        j=col+limit;
        for(i=row-limit;i<=r;++i)
	  ++histogram.at(src.at(i,j));

	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }

      //runs inside the image
      histogramMethodMiddle(src,dest,histogram,row,col);

      col=lastCol;
      while(col < (columnSize-1)) {
        j = col-limit;
	for(i=row-limit;i<=r;++i)
          --histogram.at(src.at(i,j));  // column of intensity,take off in the histogram
        ++col;
        j=col+limit;
        for(i=row-limit;i<=r;++i)
          ++histogram.at(src.at(i,j-columnSize));

	dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
    }

    //runs top rows (area 1,2,3)
    for(row=0;row<limit;++row) {
      const int r=row+limit;
      for(col=limit;col<=lastCol;++col) {  // runs top middle rows (area 2)
	histogram.fill(0);
        const int c=col+limit;
        for(i=row-limit;i<=r;++i)
	  for(j=col-limit;j<=c;++j) {
	    if(i<0)
              ++histogram.at(src.at(i+rowSize,j));
            else //if(i>=0)
              ++histogram.at(src.at(i,j));
          }
        dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
      for(col=0;col<limit;++col) {   // runs top left corner (area 1)
	histogram.fill(0);
        const int c=col+limit;
        for(i=row-limit;i<=r;++i)
	  for(j=col-limit;j<=c;++j) {
            if(i<0 && j<0)
              ++histogram.at(src.at(rowSize+i,columnSize+j));
            else if(i>=0 && j<0)
              ++histogram.at(src.at(i,columnSize+j));
            else if(i<0 && j>=0)
              ++histogram.at(src.at(rowSize+i,j));
            else //if(i>=0 && j>=0)
              ++histogram.at(src.at(i,j));
          }
        dest.at(row,col) = getMostLabel(histogram,src,row,col);
      }
      for(col=lastCol+1;col<columnSize;++col) { //runs top right corner area 3
      	histogram.fill(0);
	const int c=col+limit;
        for(i=row-limit;i<=r;++i)
          for(j=col-limit;j<=c;++j) {
            if(i<0 && j<=columnSize-1)
              ++histogram.at(src.at(rowSize+i,j));
            else if(i<0 && j>columnSize-1)
              ++histogram.at(src.at(rowSize+i,j-columnSize));
            else if(i>=0 && j>columnSize-1)
//.........这里部分代码省略.........

 void allocate(int lb, int ub, int l1, ivector & iv, int l2, int u2)
 {
     D.allocate(lb, ub, l1, iv, l2, u2);
     count.allocate(lb, ub, l1, iv);
     count = 1;
 }