C++ OutputArray类代码示例

  //  Reconstruction function for API
  void
  reconstruct(InputArrayOfArrays points2d, OutputArray Ps, OutputArray points3d, InputOutputArray K,
              bool is_projective)
  {
    const int nviews = points2d.total();
    CV_Assert( nviews >= 2 );

    // OpenCV data types
    std::vector<Mat> pts2d;
    points2d.getMatVector(pts2d);
    const int depth = pts2d[0].depth();

    Matx33d Ka = K.getMat();

    // Projective reconstruction

    if (is_projective)
    {

      if ( nviews == 2 )
      {
        // Get Projection matrices
        Matx33d F;
        Matx34d P, Pp;

        normalizedEightPointSolver(pts2d[0], pts2d[1], F);
        projectionsFromFundamental(F, P, Pp);
        Ps.create(2, 1, depth);
        Mat(P).copyTo(Ps.getMatRef(0));
        Mat(Pp).copyTo(Ps.getMatRef(1));

        // Triangulate and find 3D points using inliers
        triangulatePoints(points2d, Ps, points3d);
      }
      else
      {
        std::vector<Mat> Rs, Ts;
        reconstruct(points2d, Rs, Ts, Ka, points3d, is_projective);

        // From Rs and Ts, extract Ps
        const int nviews = Rs.size();
        Ps.create(nviews, 1, depth);

        Matx34d P;
        for (size_t i = 0; i < nviews; ++i)
        {
          projectionFromKRt(Ka, Rs[i], Vec3d(Ts[i]), P);
          Mat(P).copyTo(Ps.getMatRef(i));
        }

        Mat(Ka).copyTo(K.getMat());
      }

    }


    // Affine reconstruction

    else
    {
      // TODO: implement me
      CV_Error(Error::StsNotImplemented, "Affine reconstruction not yet implemented");
    }

  }

    void BTVL1_Base::process(InputArrayOfArrays _src, OutputArray _dst, InputArrayOfArrays _forwardMotions,
                             InputArrayOfArrays _backwardMotions, int baseIdx)
    {
        CV_Assert( scale_ > 1 );
        CV_Assert( iterations_ > 0 );
        CV_Assert( tau_ > 0.0 );
        CV_Assert( alpha_ > 0.0 );
        CV_Assert( btvKernelSize_ > 0 );
        CV_Assert( blurKernelSize_ > 0 );
        CV_Assert( blurSigma_ >= 0.0 );

        CV_OCL_RUN(_src.isUMatVector() && _dst.isUMat() && _forwardMotions.isUMatVector() &&
                   _backwardMotions.isUMatVector(),
                   ocl_process(_src, _dst, _forwardMotions, _backwardMotions, baseIdx))

        std::vector<Mat> & src = *(std::vector<Mat> *)_src.getObj(),
                & forwardMotions = *(std::vector<Mat> *)_forwardMotions.getObj(),
                & backwardMotions = *(std::vector<Mat> *)_backwardMotions.getObj();

        // update blur filter and btv weights
        if (blurKernelSize_ != curBlurKernelSize_ || blurSigma_ != curBlurSigma_ || src[0].type() != curSrcType_)
        {
            //filter_ = createGaussianFilter(src[0].type(), Size(blurKernelSize_, blurKernelSize_), blurSigma_);
            curBlurKernelSize_ = blurKernelSize_;
            curBlurSigma_ = blurSigma_;
            curSrcType_ = src[0].type();
        }

        if (btvWeights_.empty() || btvKernelSize_ != curBtvKernelSize_ || alpha_ != curAlpha_)
        {
            calcBtvWeights(btvKernelSize_, alpha_, btvWeights_);
            curBtvKernelSize_ = btvKernelSize_;
            curAlpha_ = alpha_;
        }

        // calc high res motions
        calcRelativeMotions(forwardMotions, backwardMotions, lowResForwardMotions_, lowResBackwardMotions_, baseIdx, src[0].size());

        upscaleMotions(lowResForwardMotions_, highResForwardMotions_, scale_);
        upscaleMotions(lowResBackwardMotions_, highResBackwardMotions_, scale_);

        forwardMaps_.resize(highResForwardMotions_.size());
        backwardMaps_.resize(highResForwardMotions_.size());
        for (size_t i = 0; i < highResForwardMotions_.size(); ++i)
            buildMotionMaps(highResForwardMotions_[i], highResBackwardMotions_[i], forwardMaps_[i], backwardMaps_[i]);

        // initial estimation
        const Size lowResSize = src[0].size();
        const Size highResSize(lowResSize.width * scale_, lowResSize.height * scale_);

        resize(src[baseIdx], highRes_, highResSize, 0, 0, INTER_CUBIC);

        // iterations
        diffTerm_.create(highResSize, highRes_.type());
        a_.create(highResSize, highRes_.type());
        b_.create(highResSize, highRes_.type());
        c_.create(lowResSize, highRes_.type());

        for (int i = 0; i < iterations_; ++i)
        {
            diffTerm_.setTo(Scalar::all(0));

            for (size_t k = 0; k < src.size(); ++k)
            {
                // a = M * Ih
                remap(highRes_, a_, backwardMaps_[k], noArray(), INTER_NEAREST);
                // b = HM * Ih
                GaussianBlur(a_, b_, Size(blurKernelSize_, blurKernelSize_), blurSigma_);
                // c = DHM * Ih
                resize(b_, c_, lowResSize, 0, 0, INTER_NEAREST);

                diffSign(src[k], c_, c_);

                // a = Dt * diff
                upscale(c_, a_, scale_);
                // b = HtDt * diff
                GaussianBlur(a_, b_, Size(blurKernelSize_, blurKernelSize_), blurSigma_);
                // a = MtHtDt * diff
                remap(b_, a_, forwardMaps_[k], noArray(), INTER_NEAREST);

                add(diffTerm_, a_, diffTerm_);
            }

            if (lambda_ > 0)
            {
                calcBtvRegularization(highRes_, regTerm_, btvKernelSize_, btvWeights_, ubtvWeights_);
                addWeighted(diffTerm_, 1.0, regTerm_, -lambda_, 0.0, diffTerm_);
            }

            addWeighted(highRes_, 1.0, diffTerm_, tau_, 0.0, highRes_);
        }

        Rect inner(btvKernelSize_, btvKernelSize_, highRes_.cols - 2 * btvKernelSize_, highRes_.rows - 2 * btvKernelSize_);
        highRes_(inner).copyTo(_dst);
    }

/*******************************************************************************
* Function:      subtractBGOpenDiagonal  
* Description:   BG subtraction via opening with diagonal structuring elements
* Arguments:
	inImg           -   input image
	bgsImg          -   BG subtracted image
	threshVal       -   threshold value for converting to binary image
	seLength        -   length of structuring elements
	
* Returns:       void
* Comments:
* Revision: 
*******************************************************************************/
int
FGExtraction::subtractBGOpenDiagonal(InputArray src, OutputArray dst, int threshVal, int seLength)
{
    // generate binary image by thresholding
	Mat bin;
	double thresh = threshold(src, bin, threshVal, 255, THRESH_BINARY);

	// opening by horizontal structuring element
	//Mat structElemHorizontal = Mat::ones(1, seLength, CV_8U);
	//morphologyEx(bin, dst, MORPH_OPEN, structElemHorizontal);

	// opening by vertical structuring element
	//Mat structElemVertical = Mat::ones(seLength, 1, CV_8U);
	//morphologyEx(dst, dst, MORPH_OPEN, structElemVertical);

	//imshow("src", src);
	//imshow("bin", bin);
	//waitKey(0);

    // opening by first diagonal structuring element
	Mat structElemBackSlash = Mat::eye(seLength, seLength, CV_8U);
	morphologyEx(bin, dst, MORPH_OPEN, structElemBackSlash);

	//imshow("dst1", dst);
	//waitKey(0);

    // opening by second diagonal structuring element
	Mat structElemSlash;
	flip(structElemBackSlash, structElemSlash, 0);
	morphologyEx(dst, dst, MORPH_OPEN, structElemSlash);

	//imshow("dst2", dst);
	//waitKey(0);

	// eliminate small noise
	Mat structElemEllip = getStructuringElement(MORPH_ELLIPSE, Size(seLength, seLength));
	morphologyEx(dst, dst, MORPH_OPEN, structElemEllip);

	//imshow("dst3", dst);
	//waitKey(0);


	// get object size
	Mat dstImg = dst.getMat();
	vector<vector<Point>> contours = extractContours(dstImg);
	if (contours.size()==0)
		return 1;

    Mat mask = Mat::zeros(_bgsImg.size(), CV_8U); 
	vector<int> areas(contours.size());
	int cnt = 0;
	int argMax = 0;
	int max_area = 0;
    for(vector<vector<Point> >::const_iterator it = contours.begin(); it != contours.end(); ++it){
        Rect uprightBox = boundingRect(*it);
		areas[cnt] = uprightBox.height*uprightBox.width;
		if (areas[cnt]>max_area) {
			max_area = areas[cnt];
			argMax = cnt;
		}
		cnt++;
	}
	vector<Point> largestContour = contours[argMax];	//***** only use the largest contour
	RotatedRect orientedBox = orientedBoundingBox(largestContour);

	int updateSeL = int(min(orientedBox.size.width, orientedBox.size.height)/5.0+0.5);

	// opening by first diagonal structuring element
	structElemBackSlash = Mat::eye(updateSeL, updateSeL, CV_8U);
	morphologyEx(bin, dst, MORPH_OPEN, structElemBackSlash);
	
	//imshow("dst1", dst);
	//waitKey(0);

    // opening by second diagonal structuring element
	flip(structElemBackSlash, structElemSlash, 0);
	morphologyEx(dst, dst, MORPH_OPEN, structElemSlash);

	//imshow("dst2", dst);
	//waitKey(0);

	// eliminate small noise
	structElemEllip = getStructuringElement(MORPH_ELLIPSE, Size(updateSeL, updateSeL));
	morphologyEx(dst, dst, MORPH_OPEN, structElemEllip);

	//imshow("dst3", dst);
	//waitKey(0);
//.........这里部分代码省略.........

static void divSpectrums( InputArray _srcA, InputArray _srcB, OutputArray _dst, int flags, bool conjB)
{
    Mat srcA = _srcA.getMat(), srcB = _srcB.getMat();
    int depth = srcA.depth(), cn = srcA.channels(), type = srcA.type();
    int rows = srcA.rows, cols = srcA.cols;
    int j, k;

    CV_Assert( type == srcB.type() && srcA.size() == srcB.size() );
    CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 );

    _dst.create( srcA.rows, srcA.cols, type );
    Mat dst = _dst.getMat();

    bool is_1d = (flags & DFT_ROWS) || (rows == 1 || (cols == 1 &&
             srcA.isContinuous() && srcB.isContinuous() && dst.isContinuous()));

    if( is_1d && !(flags & DFT_ROWS) )
        cols = cols + rows - 1, rows = 1;

    int ncols = cols*cn;
    int j0 = cn == 1;
    int j1 = ncols - (cols % 2 == 0 && cn == 1);

    if( depth == CV_32F )
    {
        const float* dataA = srcA.ptr<float>();
        const float* dataB = srcB.ptr<float>();
        float* dataC = dst.ptr<float>();
        float eps = FLT_EPSILON; // prevent div0 problems

        size_t stepA = srcA.step/sizeof(dataA[0]);
        size_t stepB = srcB.step/sizeof(dataB[0]);
        size_t stepC = dst.step/sizeof(dataC[0]);

        if( !is_1d && cn == 1 )
        {
            for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
            {
                if( k == 1 )
                    dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
                dataC[0] = dataA[0] / (dataB[0] + eps);
                if( rows % 2 == 0 )
                    dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / (dataB[(rows-1)*stepB] + eps);
                if( !conjB )
                    for( j = 1; j <= rows - 2; j += 2 )
                    {
                        double denom = (double)dataB[j*stepB]*dataB[j*stepB] +
                                       (double)dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + (double)eps;

                        double re = (double)dataA[j*stepA]*dataB[j*stepB] +
                                    (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];

                        double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] -
                                    (double)dataA[j*stepA]*dataB[(j+1)*stepB];

                        dataC[j*stepC] = (float)(re / denom);
                        dataC[(j+1)*stepC] = (float)(im / denom);
                    }
                else
                    for( j = 1; j <= rows - 2; j += 2 )
                    {

                        double denom = (double)dataB[j*stepB]*dataB[j*stepB] +
                                       (double)dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + (double)eps;

                        double re = (double)dataA[j*stepA]*dataB[j*stepB] -
                                    (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];

                        double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] +
                                    (double)dataA[j*stepA]*dataB[(j+1)*stepB];

                        dataC[j*stepC] = (float)(re / denom);
                        dataC[(j+1)*stepC] = (float)(im / denom);
                    }
                if( k == 1 )
                    dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1;
            }
        }

        for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC )
        {
            if( is_1d && cn == 1 )
            {
                dataC[0] = dataA[0] / (dataB[0] + eps);
                if( cols % 2 == 0 )
                    dataC[j1] = dataA[j1] / (dataB[j1] + eps);
            }

            if( !conjB )
                for( j = j0; j < j1; j += 2 )
                {
                    double denom = (double)(dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps);
                    double re = (double)(dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1]);
                    double im = (double)(dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1]);
                    dataC[j] = (float)(re / denom);
                    dataC[j+1] = (float)(im / denom);
                }
            else
                for( j = j0; j < j1; j += 2 )
                {
//.........这里部分代码省略.........

    void FeatureShiCorner::computeRawCornerMat( InputArray _image, OutputArray _corner )
    {
        // TODO check: _corner must be CV_32SC1

        const Mat image  = _image.getMat();
        const int height = image.rows;
        const int width  = image.cols;
        const int radius = 1;
        Mat derX( height, width, CV_32SC1, Scalar( 0 ) );
        Mat derY( height, width, CV_32SC1, Scalar( 0 ) );
        Mat Mx2( height, width, CV_32SC1, Scalar( 0 ) );
        Mat My2( height, width, CV_32SC1, Scalar( 0 ) );
        Mat Mxy( height, width, CV_32SC1, Scalar( 0 ) );


        applyFilter< uchar, int32_t >( _image, derX, &filter_derX[0][0], 3, 1, 0, true );
        applyFilter< uchar, int32_t >( _image, derY, &filter_derY[0][0], 1, 3, 0, true );


        int normDivisor     = 0;
        const int * pGauss  = &FeatureShiCorner::filter_gauss[0][0];
        int const * pGaussE = pGauss + 9;
        for(; pGauss != pGaussE; pGauss++ )
        {
            normDivisor += abs( *pGauss );
        }

        int32_t maxVal = 0;
        for( int y = 0; y < height; y++ )
        {
            for( int x = 0; x < width; x++ )
            {
                for( int dy = -radius; dy <= radius; dy++ )
                {
                    for( int dx = -radius; dx <= radius; dx++ )
                    {
                        int fx = x + dx;
                        if( (fx < 0) || (fx >= width) ) { continue; }
                        int fy = y + dy;
                        if( (fy < 0) || (fy >= height) ) { continue; }

                        int f = FeatureShiCorner::filter_gauss[(radius + dx)][(radius + dy)];

                        Mx2.at< int32_t >( y, x ) += int32_t( f * pow( derX.at< int32_t >( fy, fx ), 2 ) );
                        My2.at< int32_t >( y, x ) += int32_t( f * pow( derY.at< int32_t >( fy, fx ), 2 ) );
                        Mxy.at< int32_t >( y, x ) += int32_t( f * derX.at< int32_t >( fy, fx ) * derY.at< int >( fy, fx ) );
                    }
                }

                Mx2.at< int32_t >( y, x ) /= normDivisor;
                My2.at< int32_t >( y, x ) /= normDivisor;
                Mxy.at< int32_t >( y, x ) /= normDivisor;
                maxVal = max( Mx2.at< int32_t >( y, x ), maxVal );
                maxVal = max( My2.at< int32_t >( y, x ), maxVal );
                maxVal = max( Mxy.at< int32_t >( y, x ), maxVal );
            }
        }

        Mat corners      = _corner.getMat();
        const auto it_cE = corners.end< int32_t >();
        auto it_cS       = corners.begin< int32_t >();
        auto it_Mx2S     = Mx2.begin< int32_t >();
        auto it_My2S     = My2.begin< int32_t >();
        auto it_MxyS     = Mxy.begin< int32_t >();

        // reduce to high values if necessary
        // maxval: 0..1 * 255^2, maxval^2 should not overflow for the next step
        // reduce to sqrt( 2^31-1 (signed int)  ) -> 46340
        const int maxValC = 46340;
        if( maxVal > maxValC )
        {
            cout << "maxVal > maxValC | maxVal: " << maxVal << endl;

            const double scaleFac = maxValC / (double) maxVal; // scaleFac = 0.xxxx

            while( it_cS != it_cE )
            {
                *it_cS   *= int32_t( scaleFac );
                *it_Mx2S *= int32_t( scaleFac );
                *it_My2S *= int32_t( scaleFac );
                *it_MxyS *= int32_t( scaleFac );

                it_cS++;
                it_Mx2S++;
                it_My2S++;
                it_MxyS++;
            }

            // reset iterators
            it_cS   = corners.begin< int32_t >();
            it_Mx2S = Mx2.begin< int32_t >();
            it_My2S = My2.begin< int32_t >();
            it_MxyS = Mxy.begin< int32_t >();
        }

        maxVal = 0;

        // calc eigenvalues
        int32_t trc, det;
        double ev_sqrt, trc_halve, eigVal1, eigVal2;
//.........这里部分代码省略.........

bool GrayCodePattern_Impl::decode( InputArrayOfArrays patternImages, OutputArray disparityMap,
                                   InputArrayOfArrays blackImages, InputArrayOfArrays whitheImages, int flags ) const
{
  std::vector<std::vector<Mat> >& acquired_pattern = *( std::vector<std::vector<Mat> >* ) patternImages.getObj();

  if( flags == DECODE_3D_UNDERWORLD )
  {
    // Computing shadows mask
    std::vector<Mat> shadowMasks;
    computeShadowMasks( blackImages, whitheImages, shadowMasks );

    int cam_width = acquired_pattern[0][0].cols;
    int cam_height = acquired_pattern[0][0].rows;

    Point projPixel;

    // Storage for the pixels of the two cams that correspond to the same pixel of the projector
    std::vector<std::vector<std::vector<Point> > > camsPixels;
    camsPixels.resize( acquired_pattern.size() );

    // TODO: parallelize for (k and j)
    for( size_t k = 0; k < acquired_pattern.size(); k++ )
    {
      camsPixels[k].resize( params.height * params.width );
      for( int i = 0; i < cam_width; i++ )
      {
        for( int j = 0; j < cam_height; j++ )
        {
          //if the pixel is not shadowed, reconstruct
          if( shadowMasks[k].at<uchar>( j, i ) )
          {
            //for a (x,y) pixel of the camera returns the corresponding projector pixel by calculating the decimal number
            bool error = getProjPixel( acquired_pattern[k], i, j, projPixel );

            if( error )
            {
              continue;
            }

            camsPixels[k][projPixel.x * params.height + projPixel.y].push_back( Point( i, j ) );
          }
        }
      }
    }

    std::vector<Point> cam1Pixs, cam2Pixs;

    Mat& disparityMap_ = *( Mat* ) disparityMap.getObj();
    disparityMap_ = Mat( cam_height, cam_width, CV_64F, double( 0 ) );

    double number_of_pixels_cam1 = 0;
    double number_of_pixels_cam2 = 0;

    for( int i = 0; i < params.width; i++ )
    {
      for( int j = 0; j < params.height; j++ )
      {
        cam1Pixs = camsPixels[0][i * params.height + j];
        cam2Pixs = camsPixels[1][i * params.height + j];

        if( cam1Pixs.size() == 0 || cam2Pixs.size() == 0 )
          continue;

        Point p1;
        Point p2;

        double sump1x = 0;
        double sump2x = 0;

        number_of_pixels_cam1 += cam1Pixs.size();
        number_of_pixels_cam2 += cam2Pixs.size();
        for( int c1 = 0; c1 < (int) cam1Pixs.size(); c1++ )
        {
          p1 = cam1Pixs[c1];
          sump1x += p1.x;
        }
        for( int c2 = 0; c2 < (int) cam2Pixs.size(); c2++ )
        {
          p2 = cam2Pixs[c2];
          sump2x += p2.x;
        }

        sump2x /= cam2Pixs.size();
        sump1x /= cam1Pixs.size();
        for( int c1 = 0; c1 < (int) cam1Pixs.size(); c1++ )
        {
          p1 = cam1Pixs[c1];
          disparityMap_.at<double>( p1.y, p1.x ) = ( double ) (sump2x - sump1x);
        }

        sump2x = 0;
        sump1x = 0;
      }
    }

    return true;
  }  // end if flags

  return false;
}

void convexityDefects( InputArray _points, InputArray _hull, OutputArray _defects )
{
    CV_INSTRUMENT_REGION()

    Mat points = _points.getMat();
    int i, j = 0, npoints = points.checkVector(2, CV_32S);
    CV_Assert( npoints >= 0 );

    if( npoints <= 3 )
    {
        _defects.release();
        return;
    }

    Mat hull = _hull.getMat();
    int hpoints = hull.checkVector(1, CV_32S);
    CV_Assert( hpoints > 2 );

    const Point* ptr = points.ptr<Point>();
    const int* hptr = hull.ptr<int>();
    std::vector<Vec4i> defects;

    // 1. recognize co-orientation of the contour and its hull
    bool rev_orientation = ((hptr[1] > hptr[0]) + (hptr[2] > hptr[1]) + (hptr[0] > hptr[2])) != 2;

    // 2. cycle through points and hull, compute defects
    int hcurr = hptr[rev_orientation ? 0 : hpoints-1];
    CV_Assert( 0 <= hcurr && hcurr < npoints );

    for( i = 0; i < hpoints; i++ )
    {
        int hnext = hptr[rev_orientation ? hpoints - i - 1 : i];
        CV_Assert( 0 <= hnext && hnext < npoints );

        Point pt0 = ptr[hcurr], pt1 = ptr[hnext];
        double dx0 = pt1.x - pt0.x;
        double dy0 = pt1.y - pt0.y;
        double scale = dx0 == 0 && dy0 == 0 ? 0. : 1./std::sqrt(dx0*dx0 + dy0*dy0);

        int defect_deepest_point = -1;
        double defect_depth = 0;
        bool is_defect = false;
        j=hcurr;
        for(;;)
        {
            // go through points to achieve next hull point
            j++;
            j &= j >= npoints ? 0 : -1;
            if( j == hnext )
                break;

            // compute distance from current point to hull edge
            double dx = ptr[j].x - pt0.x;
            double dy = ptr[j].y - pt0.y;
            double dist = fabs(-dy0*dx + dx0*dy) * scale;

            if( dist > defect_depth )
            {
                defect_depth = dist;
                defect_deepest_point = j;
                is_defect = true;
            }
        }

        if( is_defect )
        {
            int idepth = cvRound(defect_depth*256);
            defects.push_back(Vec4i(hcurr, hnext, defect_deepest_point, idepth));
        }

        hcurr = hnext;
    }

    Mat(defects).copyTo(_defects);
}

template <typename _Tp> static
inline void ulbp_(InputArray _src, OutputArray _dst, int radius, int neighbors, vector<int> m_uniform)
{
    if (neighbors != 8 || m_uniform.size() != 58)
    {
        cout << "neighbor must be 8! and uniform size be 58!\n";
        system("pause");
        exit(-1);
    }


    //get matrices
    Mat src = _src.getMat();
    // allocate memory for result
    _dst.create(src.rows-2*radius, src.cols-2*radius, CV_32SC1);
    Mat dst = _dst.getMat();
    // zero
    dst.setTo(0);
    for(int n=0; n<neighbors; n++)
    {
        // sample points
        float x = static_cast<float>(-radius) * sin(2.0*CV_PI*n/static_cast<float>(neighbors));
        float y = static_cast<float>(radius) * cos(2.0*CV_PI*n/static_cast<float>(neighbors));
        // relative indices
        int fx = static_cast<int>(floor(x));
        int fy = static_cast<int>(floor(y));
        int cx = static_cast<int>(ceil(x));
        int cy = static_cast<int>(ceil(y));
        // fractional part
        float ty = y - fy;
        float tx = x - fx;
        // set interpolation weights
        float w1 = (1 - tx) * (1 - ty);
        float w2 =      tx  * (1 - ty);
        float w3 = (1 - tx) *      ty;
        float w4 =      tx  *      ty;
        // iterate through your data
        for(int i=radius; i < src.rows-radius; i++) {
            for(int j=radius; j < src.cols-radius; j++) {
                // calculate interpolated value
                float t = w1*src.at<_Tp>(i+fy,j+fx) + w2*src.at<_Tp>(i+fy,j+cx) + w3*src.at<_Tp>(i+cy,j+fx) + w4*src.at<_Tp>(i+cy,j+cx);
                // floating point precision, so check some machine-dependent epsilon
                dst.at<int>(i-radius,j-radius) += ((t > src.at<_Tp>(i,j)) ||
                                                   (std::abs(t-src.at<_Tp>(i,j)) < std::numeric_limits<float>::epsilon()))
                                                  << n;
            }
        }
    }


    for (int i = 0; i < dst.rows; ++i)
    {
        for (int j = 0; j < dst.cols; ++j)
        {
            int data = dst.at<int>(i, j);
            vector<int>::iterator iter = find(m_uniform.begin(), m_uniform.end(), data);
            if (iter == m_uniform.end())
            {
                dst.at<int>(i, j) = 0;
            }
            else
            {
                int new_data = iter - m_uniform.begin() ;
                dst.at<int>(i, j) = new_data + 1;
            }
        }
    }
}

void cv::goodFeaturesToTrack( InputArray _image, OutputArray _corners,
                              int maxCorners, double qualityLevel, double minDistance,
                              InputArray _mask, int blockSize,
                              bool useHarrisDetector, double harrisK )
{
    Mat image = _image.getMat(), mask = _mask.getMat();
    
    CV_Assert( qualityLevel > 0 && minDistance >= 0 && maxCorners >= 0 );
    CV_Assert( mask.empty() || (mask.type() == CV_8UC1 && mask.size() == image.size()) );

    Mat eig, tmp;
    if( useHarrisDetector )
        cornerHarris( image, eig, blockSize, 3, harrisK );
    else
        cornerMinEigenVal( image, eig, blockSize, 3 );

    double maxVal = 0;
    minMaxLoc( eig, 0, &maxVal, 0, 0, mask );
    threshold( eig, eig, maxVal*qualityLevel, 0, THRESH_TOZERO );
    dilate( eig, tmp, Mat());

    Size imgsize = image.size();

    vector<const float*> tmpCorners;

    // collect list of pointers to features - put them into temporary image
    for( int y = 1; y < imgsize.height - 1; y++ )
    {
        const float* eig_data = (const float*)eig.ptr(y);
        const float* tmp_data = (const float*)tmp.ptr(y);
        const uchar* mask_data = mask.data ? mask.ptr(y) : 0;

        for( int x = 1; x < imgsize.width - 1; x++ )
        {
            float val = eig_data[x];
            if( val != 0 && val == tmp_data[x] && (!mask_data || mask_data[x]) )
                tmpCorners.push_back(eig_data + x);
        }
    }

    sort( tmpCorners, greaterThanPtr<float>() );
    vector<Point2f> corners;
    size_t i, j, total = tmpCorners.size(), ncorners = 0;

    if(minDistance >= 1)
    {
         // Partition the image into larger grids
        int w = image.cols;
        int h = image.rows;

        const int cell_size = cvRound(minDistance);
        const int grid_width = (w + cell_size - 1) / cell_size;
        const int grid_height = (h + cell_size - 1) / cell_size;

        std::vector<std::vector<Point2f> > grid(grid_width*grid_height);

        minDistance *= minDistance;

        for( i = 0; i < total; i++ )
        {
            int ofs = (int)((const uchar*)tmpCorners[i] - eig.data);
            int y = (int)(ofs / eig.step);
            int x = (int)((ofs - y*eig.step)/sizeof(float));

	        bool good = true;

            int x_cell = x / cell_size;
            int y_cell = y / cell_size;

            int x1 = x_cell - 1;
            int y1 = y_cell - 1;
            int x2 = x_cell + 1;
            int y2 = y_cell + 1;

            // boundary check
            x1 = std::max(0, x1);
            y1 = std::max(0, y1);
            x2 = std::min(grid_width-1, x2);
            y2 = std::min(grid_height-1, y2);

            for( int yy = y1; yy <= y2; yy++ )
            {
                for( int xx = x1; xx <= x2; xx++ )
                {   
                    vector <Point2f> &m = grid[yy*grid_width + xx];

                    if( m.size() )
                    {
                        for(j = 0; j < m.size(); j++)
                        {
                            float dx = x - m[j].x;
                            float dy = y - m[j].y;

                            if( dx*dx + dy*dy < minDistance )
                            {
                                good = false;
                                goto break_out;
                            }
                        }
                    }                
//.........这里部分代码省略.........

void convexHull( InputArray _points, OutputArray _hull, bool clockwise, bool returnPoints )
{
    CV_INSTRUMENT_REGION()

    Mat points = _points.getMat();
    int i, total = points.checkVector(2), depth = points.depth(), nout = 0;
    int miny_ind = 0, maxy_ind = 0;
    CV_Assert(total >= 0 && (depth == CV_32F || depth == CV_32S));

    if( total == 0 )
    {
        _hull.release();
        return;
    }

    returnPoints = !_hull.fixedType() ? returnPoints : _hull.type() != CV_32S;

    bool is_float = depth == CV_32F;
    AutoBuffer<Point*> _pointer(total);
    AutoBuffer<int> _stack(total + 2), _hullbuf(total);
    Point** pointer = _pointer;
    Point2f** pointerf = (Point2f**)pointer;
    Point* data0 = points.ptr<Point>();
    int* stack = _stack;
    int* hullbuf = _hullbuf;

    CV_Assert(points.isContinuous());

    for( i = 0; i < total; i++ )
        pointer[i] = &data0[i];

    // sort the point set by x-coordinate, find min and max y
    if( !is_float )
    {
        std::sort(pointer, pointer + total, CHullCmpPoints<int>());
        for( i = 1; i < total; i++ )
        {
            int y = pointer[i]->y;
            if( pointer[miny_ind]->y > y )
                miny_ind = i;
            if( pointer[maxy_ind]->y < y )
                maxy_ind = i;
        }
    }
    else
    {
        std::sort(pointerf, pointerf + total, CHullCmpPoints<float>());
        for( i = 1; i < total; i++ )
        {
            float y = pointerf[i]->y;
            if( pointerf[miny_ind]->y > y )
                miny_ind = i;
            if( pointerf[maxy_ind]->y < y )
                maxy_ind = i;
        }
    }

    if( pointer[0]->x == pointer[total-1]->x &&
        pointer[0]->y == pointer[total-1]->y )
    {
        hullbuf[nout++] = 0;
    }
    else
    {
        // upper half
        int *tl_stack = stack;
        int tl_count = !is_float ?
            Sklansky_( pointer, 0, maxy_ind, tl_stack, -1, 1) :
            Sklansky_( pointerf, 0, maxy_ind, tl_stack, -1, 1);
        int *tr_stack = stack + tl_count;
        int tr_count = !is_float ?
            Sklansky_( pointer, total-1, maxy_ind, tr_stack, -1, -1) :
            Sklansky_( pointerf, total-1, maxy_ind, tr_stack, -1, -1);

        // gather upper part of convex hull to output
        if( !clockwise )
        {
            std::swap( tl_stack, tr_stack );
            std::swap( tl_count, tr_count );
        }

        for( i = 0; i < tl_count-1; i++ )
            hullbuf[nout++] = int(pointer[tl_stack[i]] - data0);
        for( i = tr_count - 1; i > 0; i-- )
            hullbuf[nout++] = int(pointer[tr_stack[i]] - data0);
        int stop_idx = tr_count > 2 ? tr_stack[1] : tl_count > 2 ? tl_stack[tl_count - 2] : -1;

        // lower half
        int *bl_stack = stack;
        int bl_count = !is_float ?
            Sklansky_( pointer, 0, miny_ind, bl_stack, 1, -1) :
            Sklansky_( pointerf, 0, miny_ind, bl_stack, 1, -1);
        int *br_stack = stack + bl_count;
        int br_count = !is_float ?
            Sklansky_( pointer, total-1, miny_ind, br_stack, 1, 1) :
            Sklansky_( pointerf, total-1, miny_ind, br_stack, 1, 1);

        if( clockwise )
        {
            std::swap( bl_stack, br_stack );
//.........这里部分代码省略.........

void SIFT::operator()(InputArray _image, InputArray _mask,
                      std::vector<KeyPoint>& keypoints,
                      OutputArray _descriptors,
                      bool useProvidedKeypoints) const
{
    int firstOctave = -1, actualNOctaves = 0, actualNLayers = 0;
    Mat image = _image.getMat(), mask = _mask.getMat();

    if( image.empty() || image.depth() != CV_8U )
        CV_Error( Error::StsBadArg, "image is empty or has incorrect depth (!=CV_8U)" );

    if( !mask.empty() && mask.type() != CV_8UC1 )
        CV_Error( Error::StsBadArg, "mask has incorrect type (!=CV_8UC1)" );

    if( useProvidedKeypoints )
    {
        firstOctave = 0;
        int maxOctave = INT_MIN;
        for( size_t i = 0; i < keypoints.size(); i++ )
        {
            int octave, layer;
            float scale;
            unpackOctave(keypoints[i], octave, layer, scale);
            firstOctave = std::min(firstOctave, octave);
            maxOctave = std::max(maxOctave, octave);
            actualNLayers = std::max(actualNLayers, layer-2);
        }

        firstOctave = std::min(firstOctave, 0);
        CV_Assert( firstOctave >= -1 && actualNLayers <= nOctaveLayers );
        actualNOctaves = maxOctave - firstOctave + 1;
    }

    Mat base = createInitialImage(image, firstOctave < 0, (float)sigma);
    std::vector<Mat> gpyr, dogpyr;
    int nOctaves = actualNOctaves > 0 ? actualNOctaves : cvRound(std::log( (double)std::min( base.cols, base.rows ) ) / std::log(2.) - 2) - firstOctave;

    //double t, tf = getTickFrequency();
    //t = (double)getTickCount();
    buildGaussianPyramid(base, gpyr, nOctaves);
    buildDoGPyramid(gpyr, dogpyr);

    //t = (double)getTickCount() - t;
    //printf("pyramid construction time: %g\n", t*1000./tf);

    if( !useProvidedKeypoints )
    {
        //t = (double)getTickCount();
        findScaleSpaceExtrema(gpyr, dogpyr, keypoints);
        KeyPointsFilter::removeDuplicated( keypoints );

        if( !mask.empty() )
            KeyPointsFilter::runByPixelsMask( keypoints, mask );

        if( nfeatures > 0 )
            KeyPointsFilter::retainBest(keypoints, nfeatures);
        //t = (double)getTickCount() - t;
        //printf("keypoint detection time: %g\n", t*1000./tf);

        if( firstOctave < 0 )
            for( size_t i = 0; i < keypoints.size(); i++ )
            {
                KeyPoint& kpt = keypoints[i];
                float scale = 1.f/(float)(1 << -firstOctave);
                kpt.octave = (kpt.octave & ~255) | ((kpt.octave + firstOctave) & 255);
                kpt.pt *= scale;
                kpt.size *= scale;
            }
    }
    else
    {
        // filter keypoints by mask
        //KeyPointsFilter::runByPixelsMask( keypoints, mask );
    }

    if( _descriptors.needed() )
    {
        //t = (double)getTickCount();
        int dsize = descriptorSize();
        _descriptors.create((int)keypoints.size(), dsize, CV_32F);
        Mat descriptors = _descriptors.getMat();

        calcDescriptors(gpyr, keypoints, descriptors, nOctaveLayers, firstOctave);
        //t = (double)getTickCount() - t;
        //printf("descriptor extraction time: %g\n", t*1000./tf);
    }
}

static bool matchTemplate_CCOEFF_NORMED(InputArray _image, InputArray _templ, OutputArray _result)
{
    matchTemplate(_image, _templ, _result, CV_TM_CCORR);

    UMat temp, image_sums, image_sqsums;
    integral(_image, image_sums, image_sqsums, CV_32F, CV_32F);

    int type = image_sums.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);

    ocl::Kernel k("matchTemplate_CCOEFF_NORMED", ocl::imgproc::match_template_oclsrc,
        format("-D CCOEFF_NORMED -D type=%s -D elem_type=%s -D cn=%d", ocl::typeToStr(type), ocl::typeToStr(depth), cn));
    if (k.empty())
        return false;

    UMat templ = _templ.getUMat();
    Size size = _image.size(), tsize = templ.size();
    _result.create(size.height - templ.rows + 1, size.width - templ.cols + 1, CV_32F);
    UMat result = _result.getUMat();

    float scale = 1.f / tsize.area();

    if (cn == 1)
    {
        float templ_sum = (float)sum(templ)[0];

        multiply(templ, templ, temp, 1, CV_32F);
        float templ_sqsum = (float)sum(temp)[0];

        templ_sqsum -= scale * templ_sum * templ_sum;
        templ_sum   *= scale;

        if (templ_sqsum < DBL_EPSILON)
        {
            result = Scalar::all(1);
            return true;
        }

        k.args(ocl::KernelArg::ReadOnlyNoSize(image_sums), ocl::KernelArg::ReadOnlyNoSize(image_sqsums),
                      ocl::KernelArg::ReadWrite(result), templ.rows, templ.cols, scale, templ_sum, templ_sqsum);
    }
    else
    {
        Vec4f templ_sum = Vec4f::all(0), templ_sqsum = Vec4f::all(0);
        templ_sum = sum(templ);

        multiply(templ, templ, temp, 1, CV_32F);
        templ_sqsum = sum(temp);

        float templ_sqsum_sum = 0;
        for (int i = 0; i < cn; i ++)
            templ_sqsum_sum += templ_sqsum[i] - scale * templ_sum[i] * templ_sum[i];

        templ_sum *= scale;

        if (templ_sqsum_sum < DBL_EPSILON)
        {
            result = Scalar::all(1);
            return true;
        }

        if (cn == 2)
            k.args(ocl::KernelArg::ReadOnlyNoSize(image_sums), ocl::KernelArg::ReadOnlyNoSize(image_sqsums),
                   ocl::KernelArg::ReadWrite(result), templ.rows, templ.cols, scale,
                   templ_sum[0], templ_sum[1], templ_sqsum_sum);
        else
            k.args(ocl::KernelArg::ReadOnlyNoSize(image_sums), ocl::KernelArg::ReadOnlyNoSize(image_sqsums),
                   ocl::KernelArg::ReadWrite(result), templ.rows, templ.cols, scale,
                   templ_sum[0], templ_sum[1], templ_sum[2], templ_sum[3], templ_sqsum_sum);
    }

    size_t globalsize[2] = { result.cols, result.rows };
    return k.run(2, globalsize, NULL, false);
}

/**
 * Function that computes the Harris responses in a
 * 2*r x 2*r patch at given points in the image
 */
static void
HarrisResponses(InputArray _img, InputArray _diff_x, InputArray _diff_y,
                std::vector<KeyPoint>& pts,
                OutputArray _response, int r, float harris_k)
{
    size_t ptidx, ptsize = pts.size();

    // Get mats
    Mat img = _img.getMat(), diff_x = _diff_x.getMat(),
        diff_y = _diff_y.getMat(), response;
    CV_Assert( img.type() == CV_8UC1 );


    bool compute_response = _response.needed();
    if (compute_response) response = _response.getMat();

    const int* dx00 = diff_x.ptr<int>();
    const int* dy00 = diff_y.ptr<int>();
    float* r00 = response.ptr<float>();
    int step = diff_x.step1();
    int r_step = response.step1();

    for( ptidx = 0; ptidx < ptsize; ptidx++ )
    {
        float kp_x = pts[ptidx].pt.x;
        float kp_y = pts[ptidx].pt.y;
        int x0 = (int)kp_x;
        int y0 = (int)kp_y;

        float xd = 2;
        float yd = 2;
        //float xd = 0.5;
        //float yd = 0.5;

        const int* dx0 = dx00 + (y0)*step + x0;
        const int* dy0 = dy00 + (y0)*step + x0;
        int a = 0, b = 0, c = 0, d = 0;
        float* r0 = r00 + ptidx*r_step;

        for( int i = -r; i < r; i++ )
        {
            for( int j = -r; j < r; j++ )
            {
                const int ofs = i*step + j;
                const int* dx = dx0 + ofs;
                const int* dy = dy0 + ofs;
                const int Ix = (float)dx[-1]*(xd) + (float)dx[0] + (float)dx[1]*xd + (float)dx[-step]*(yd) + (float)dx[step]*yd;
                const int Iy = (float)dy[-1]*(xd) + (float)dy[0] + (float)dy[1]*xd + (float)dy[-step]*(yd) + (float)dy[step]*yd;
                a += (Ix*Ix);
                b += (Iy*Iy);
                c += (Ix*Iy);
                d += Ix;
            }
        }
        if (compute_response) {
            r0[0] = (float)a;
            r0[1] = (float)b;
            r0[2] = (float)c;
            r0[3] = (float)d;
        }
        else
            pts[ptidx].response = ((float)a * b - (float)c * c -
                                   harris_k * ((float)a + b) * ((float)a + b));
    }
}

void spatialGradient( InputArray _src, OutputArray _dx, OutputArray _dy,
                      int ksize, int borderType )
{
    CV_INSTRUMENT_REGION()

    // Prepare InputArray src
    Mat src = _src.getMat();
    CV_Assert( !src.empty() );
    CV_Assert( src.type() == CV_8UC1 );
    CV_Assert( borderType == BORDER_DEFAULT || borderType == BORDER_REPLICATE );

    // Prepare OutputArrays dx, dy
    _dx.create( src.size(), CV_16SC1 );
    _dy.create( src.size(), CV_16SC1 );
    Mat dx = _dx.getMat(),
        dy = _dy.getMat();

    // TODO: Allow for other kernel sizes
    CV_Assert(ksize == 3);

    // Get dimensions
    const int H = src.rows,
              W = src.cols;

    // Row, column indices
    int i = 0,
        j = 0;

    // Handle border types
    int i_top    = 0,     // Case for H == 1 && W == 1 && BORDER_REPLICATE
        i_bottom = H - 1,
        j_offl   = 0,     // j offset from 0th   pixel to reach -1st pixel
        j_offr   = 0;     // j offset from W-1th pixel to reach Wth  pixel

    if ( borderType == BORDER_DEFAULT ) // Equiv. to BORDER_REFLECT_101
    {
        if ( H > 1 )
        {
            i_top    = 1;
            i_bottom = H - 2;
        }
        if ( W > 1 )
        {
            j_offl = 1;
            j_offr = -1;
        }
    }

    // Pointer to row vectors
    uchar *p_src, *c_src, *n_src; // previous, current, next row
    short *c_dx,  *c_dy;

    int i_start = 0;
    int j_start = 0;
#if CV_SIMD128 && CV_SSE2
    if(hasSIMD128())
    {
        uchar *m_src;
        short *n_dx, *n_dy;

        // Characters in variable names have the following meanings:
        // u: unsigned char
        // s: signed int
        //
        // [row][column]
        // m: offset -1
        // n: offset  0
        // p: offset  1
        // Example: umn is offset -1 in row and offset 0 in column
        for ( i = 0; i < H - 1; i += 2 )
        {
            if   ( i == 0 ) p_src = src.ptr<uchar>(i_top);
            else            p_src = src.ptr<uchar>(i-1);

            c_src = src.ptr<uchar>(i);
            n_src = src.ptr<uchar>(i+1);

            if ( i == H - 2 ) m_src = src.ptr<uchar>(i_bottom);
            else              m_src = src.ptr<uchar>(i+2);

            c_dx = dx.ptr<short>(i);
            c_dy = dy.ptr<short>(i);
            n_dx = dx.ptr<short>(i+1);
            n_dy = dy.ptr<short>(i+1);

            v_uint8x16 v_select_m = v_uint8x16(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                               0, 0, 0, 0xFF);

            // Process rest of columns 16-column chunks at a time
            for ( j = 1; j < W - 16; j += 16 )
            {
                // Load top row for 3x3 Sobel filter
                v_uint8x16 v_um = v_load(&p_src[j-1]);
                v_uint8x16 v_up = v_load(&p_src[j+1]);
                // TODO: Replace _mm_slli_si128 with hal method
                v_uint8x16 v_un = v_select(v_select_m, v_uint8x16(_mm_slli_si128(v_up.val, 1)),
                                                       v_uint8x16(_mm_srli_si128(v_um.val, 1)));
                v_uint16x8 v_um1, v_um2, v_un1, v_un2, v_up1, v_up2;
                v_expand(v_um, v_um1, v_um2);
                v_expand(v_un, v_un1, v_un2);
//.........这里部分代码省略.........

double cv::kmeans( InputArray _data, int K,
                   InputOutputArray _bestLabels,
                   TermCriteria criteria, int attempts,
                   int flags, OutputArray _centers )
{
    const int SPP_TRIALS = 3;
    Mat data0 = _data.getMat();
    bool isrow = data0.rows == 1;
    int N = isrow ? data0.cols : data0.rows;
    int dims = (isrow ? 1 : data0.cols)*data0.channels();
    int type = data0.depth();

    attempts = std::max(attempts, 1);
    CV_Assert( data0.dims <= 2 && type == CV_32F && K > 0 );
    CV_Assert( N >= K );

    Mat data(N, dims, CV_32F, data0.ptr(), isrow ? dims * sizeof(float) : static_cast<size_t>(data0.step));

    _bestLabels.create(N, 1, CV_32S, -1, true);

    Mat _labels, best_labels = _bestLabels.getMat();
    if( flags & CV_KMEANS_USE_INITIAL_LABELS )
    {
        CV_Assert( (best_labels.cols == 1 || best_labels.rows == 1) &&
                   best_labels.cols*best_labels.rows == N &&
                   best_labels.type() == CV_32S &&
                   best_labels.isContinuous());
        best_labels.copyTo(_labels);
    }
    else
    {
        if( !((best_labels.cols == 1 || best_labels.rows == 1) &&
                best_labels.cols*best_labels.rows == N &&
                best_labels.type() == CV_32S &&
                best_labels.isContinuous()))
            best_labels.create(N, 1, CV_32S);
        _labels.create(best_labels.size(), best_labels.type());
    }
    int* labels = _labels.ptr<int>();

    Mat centers(K, dims, type), old_centers(K, dims, type), temp(1, dims, type);
    std::vector<int> counters(K);
    std::vector<Vec2f> _box(dims);
    Vec2f* box = &_box[0];
    double best_compactness = DBL_MAX, compactness = 0;
    RNG& rng = theRNG();
    int a, iter, i, j, k;

    if( criteria.type & TermCriteria::EPS )
        criteria.epsilon = std::max(criteria.epsilon, 0.);
    else
        criteria.epsilon = FLT_EPSILON;
    criteria.epsilon *= criteria.epsilon;

    if( criteria.type & TermCriteria::COUNT )
        criteria.maxCount = std::min(std::max(criteria.maxCount, 2), 100);
    else
        criteria.maxCount = 100;

    if( K == 1 )
    {
        attempts = 1;
        criteria.maxCount = 2;
    }

    const float* sample = data.ptr<float>(0);
    for( j = 0; j < dims; j++ )
        box[j] = Vec2f(sample[j], sample[j]);

    for( i = 1; i < N; i++ )
    {
        sample = data.ptr<float>(i);
        for( j = 0; j < dims; j++ )
        {
            float v = sample[j];
            box[j][0] = std::min(box[j][0], v);
            box[j][1] = std::max(box[j][1], v);
        }
    }

    for( a = 0; a < attempts; a++ )
    {
        double max_center_shift = DBL_MAX;
        for( iter = 0;; )
        {
            swap(centers, old_centers);

            if( iter == 0 && (a > 0 || !(flags & KMEANS_USE_INITIAL_LABELS)) )
            {
                if( flags & KMEANS_PP_CENTERS )
                    generateCentersPP(data, centers, K, rng, SPP_TRIALS);
                else
                {
                    for( k = 0; k < K; k++ )
                        generateRandomCenter(_box, centers.ptr<float>(k), rng);
                }
            }
            else
            {
                if( iter == 0 && a == 0 && (flags & KMEANS_USE_INITIAL_LABELS) )
//.........这里部分代码省略.........