Opencv HOG行人检测 源码分析(二)

前一篇博客大体讲了下思路，对比较难理解的关系有些图示 http://blog.youkuaiyun.com/soidnhp/article/details/11874285  

/*M///
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/
#include <stdio.h>
#include "precomp.hpp"	//包含了 objdetect.hpp 
#include <iterator>
#ifdef HAVE_IPP
#include "ipp.h"
#endif
/****************************************************************************************\
      The code below is implementation of HOG (Histogram-of-Oriented Gradients)
      descriptor and object detection, introduced by Navneet Dalal and Bill Triggs.

      The computed feature vectors are compatible with the
      INRIA Object Detection and Localization Toolkit
      (http://pascal.inrialpes.fr/soft/olt/)
\****************************************************************************************/

namespace cv
{

size_t HOGDescriptor::getDescriptorSize() const
{
    CV_Assert(blockSize.width % cellSize.width == 0 &&
        blockSize.height % cellSize.height == 0);
    CV_Assert((winSize.width - blockSize.width) % blockStride.width == 0 &&
        (winSize.height - blockSize.height) % blockStride.height == 0 );
    return (size_t)nbins*
        (blockSize.width/cellSize.width)*
        (blockSize.height/cellSize.height)*
        ((winSize.width - blockSize.width)/blockStride.width + 1)*
        ((winSize.height - blockSize.height)/blockStride.height + 1);//描述向量总长度
}

double HOGDescriptor::getWinSigma() const
{
    return winSigma >= 0 ? winSigma : (blockSize.width + blockSize.height)/8.; //默认-1
}

bool HOGDescriptor::checkDetectorSize() const
{
    size_t detectorSize = svmDetector.size(), descriptorSize = getDescriptorSize();
    return detectorSize == 0 ||
        detectorSize == descriptorSize ||
        detectorSize == descriptorSize + 1;
}

void HOGDescriptor::setSVMDetector(InputArray _svmDetector)
{
    _svmDetector.getMat().convertTo(svmDetector, CV_32F);
    CV_Assert( checkDetectorSize() );
}

#define CV_TYPE_NAME_HOG_DESCRIPTOR "opencv-object-detector-hog"

bool HOGDescriptor::read(FileNode& obj)
{
    if( !obj.isMap() )
        return false;
    FileNodeIterator it = obj["winSize"].begin();
    it >> winSize.width >> winSize.height;
    it = obj["blockSize"].begin();
    it >> blockSize.width >> blockSize.height;
    it = obj["blockStride"].begin();
    it >> blockStride.width >> blockStride.height;
    it = obj["cellSize"].begin();
    it >> cellSize.width >> cellSize.height;
    obj["nbins"] >> nbins;
    obj["derivAperture"] >> derivAperture;
    obj["winSigma"] >> winSigma;
    obj["histogramNormType"] >> histogramNormType;
    obj["L2HysThreshold"] >> L2HysThreshold;
    obj["gammaCorrection"] >> gammaCorrection;
    obj["nlevels"] >> nlevels;

    FileNode vecNode = obj["SVMDetector"];
    if( vecNode.isSeq() )
    {
        vecNode >> svmDetector;
        CV_Assert(checkDetectorSize());
    }
    return true;
}

void HOGDescriptor::write(FileStorage& fs, const String& objName) const
{
    if( !objName.empty() )
        fs << objName;

    fs << "{" CV_TYPE_NAME_HOG_DESCRIPTOR
    << "winSize" << winSize
    << "blockSize" << blockSize
    << "blockStride" << blockStride
    << "cellSize" << cellSize
    << "nbins" << nbins
    << "derivAperture" << derivAperture
    << "winSigma" << getWinSigma()
    << "histogramNormType" << histogramNormType
    << "L2HysThreshold" << L2HysThreshold
    << "gammaCorrection" << gammaCorrection
    << "nlevels" << nlevels;
    if( !svmDetector.empty() )
        fs << "SVMDetector" << svmDetector;
    fs << "}";
}

bool HOGDescriptor::load(const String& filename, const String& objname)
{
    FileStorage fs(filename, FileStorage::READ);
    FileNode obj = !objname.empty() ? fs[objname] : fs.getFirstTopLevelNode();
    return read(obj);
}

void HOGDescriptor::save(const String& filename, const String& objName) const
{
    FileStorage fs(filename, FileStorage::WRITE);
    write(fs, !objName.empty() ? objName : FileStorage::getDefaultObjectName(filename));
}

void HOGDescriptor::copyTo(HOGDescriptor& c) const
{
    c.winSize = winSize;
    c.blockSize = blockSize;
    c.blockStride = blockStride;
    c.cellSize = cellSize;
    c.nbins = nbins;
    c.derivAperture = derivAperture;
    c.winSigma = winSigma;
    c.histogramNormType = histogramNormType;
    c.L2HysThreshold = L2HysThreshold;
    c.gammaCorrection = gammaCorrection;
    c.svmDetector = svmDetector;
    c.nlevels = nlevels;
}
//返回 grad：梯度的模在与梯度方向相邻的两个bin的插值值，qangle：与梯度方向相邻的两个bin的编号
void HOGDescriptor::computeGradient(const Mat& img, Mat& grad, Mat& qangle,
                                    Size paddingTL, Size paddingBR) const
{
    CV_Assert( img.type() == CV_8U || img.type() == CV_8UC3 );

    Size gradsize(img.cols + paddingTL.width + paddingBR.width,
                  img.rows + paddingTL.height + paddingBR.height);
    grad.create(gradsize, CV_32FC2);  // <magnitude*(1-alpha), magnitude*alpha>，与该点梯度方向相邻两个bin的梯度模值，由该点线性插值得到
    qangle.create(gradsize, CV_8UC2); // [0..nbins-1] - quantized gradient orientation，与该点梯度方向相邻两个bin的编号
    Size wholeSize;
    Point roiofs;
    img.locateROI(wholeSize, roiofs);	//img如果是一个大图像IMG的Region of interesting,那么IMG和img共享内存
										//比如IMG(120x120),img取自IMG的一部分TL坐标（10,10），BR坐标（109,109）那么尺寸为（100x100）
										//这个函数就返回父矩阵IMG的size（120x120），以及img在IMG中的坐标偏移（roiofs.x=10,roiofs.y=10）
	/*
	Locates the matrix header within a parent matrix.
	wholeSize – Output parameter that contains the size of the whole matrix containing *this as a part
	ofs – Output parameter that contains an offset of *this inside the whole matrix.
	*/
    int i, x, y;
    int cn = img.channels();

    Mat_<float> _lut(1, 256); //gamma 校正Look up table,Mat_ 简化版的 Mat,元素访问直接 用（x,y），无需 .at,但是速度是一样的
    const float* lut = &_lut(0,0);	//只能读

    if( gammaCorrection )
        for( i = 0; i < 256; i++ )
            _lut(0,i) = std::sqrt((float)i);	//gammma 校正 r=0.5,暗区对比度提高，亮区对比度下降
    else
        for( i = 0; i < 256; i++ )
            _lut(0,i) = (float)i;

    AutoBuffer<int> mapbuf(gradsize.width + gradsize.height + 4);	//自动buffer，就不需要我们malloc，free
    int* xmap = (int*)mapbuf + 1;
    int* ymap = xmap + gradsize.width + 2;

    const int borderType = (int)BORDER_REFLECT_101;
	//一种很奇怪的插值方式，扩展出来的边缘用原图像中的像素值，并没有真正扩展存储空间
	//比如说原图为 100x100，现在要访问（-10，-10）的值，但是内存里面不不存在这个值，这种插值方法就是在原图像中找个像素点（比如（5,6））的值作为（-10,-10）的值
	//也就是将扩展后的坐标范围比如（120x120）映射到（100x100）。x,y坐标分别映射，映射表存在xmap,ymap。上面的例子中xmap[-10]=5,ymap[-10]=6
    for( x = -1; x < gradsize.width + 1; x++ )
        xmap[x] = borderInterpolate(x - paddingTL.width + roiofs.x,wholeSize.width, borderType) - roiofs.x;
    for( y = -1; y < gradsize.height + 1; y++ )
        ymap[y] = borderInterpolate(y - paddingTL.height + roiofs.y,wholeSize.height, borderType) - roiofs.y;

    // x- & y- derivatives for the whole row
    int width = gradsize.width;
    AutoBuffer<float> _dbuf(width*4);
    float* dbuf = _dbuf;
    Mat Dx(1, width, CV_32F, dbuf);
    Mat Dy(1, width, CV_32F, dbuf + width);
    Mat Mag(1, width, CV_32F, dbuf + width*2);
    Mat Angle(1, width, CV_32F, dbuf + width*3);

    int _nbins = nbins;
    float angleScale = (float)(_nbins/CV_PI);	//算某一弧度，对应落在哪一个bin的scale
#ifdef HAVE_IPP	//intel的ipp库，优化
    Mat lutimg(img.rows,img.cols,CV_MAKETYPE(CV_32F,cn));	//cn ,为1/3，对于类型 CV_32FC1、CV_32FC3
    Mat hidxs(1, width, CV_32F);
    Ipp32f* pHidxs  = (Ipp32f*)hidxs.data;
    Ipp32f* pAngles = (Ipp32f*)Angle.data;

    IppiSize roiSize;
    roiSize.width = img.cols;
    roiSize.height = img.rows;

	//对原始图像，进行gamma校正，结果保存在 imglutPtr
    for( y = 0; y < roiSize.height; y++ )
    {
       const uchar* imgPtr = img.data + y*img.step;
       float* imglutPtr = (float*)(lutimg.data + y*lutimg.step);

       for( x = 0; x < roiSize.width*cn; x++ )
       {
          imglutPtr[x] = lut[imgPtr[x]];	//查表 gamma校正
       }
    }

#endif
	//好长的循环体，计算了四个梯度的四个量 Dx，Dy, Angle,Mag,最终是保存了Angle,Mag两个量给后续的工作用
    for( y = 0; y < gradsize.height; y++ )
    {
		//行指针(加上了补丁)
#ifdef HAVE_IPP
        const float* imgPtr  = (float*)(lutimg.data + lutimg.step*ymap[y]);
        const float* prevPtr = (float*)(lutimg.data + lutimg.step*ymap[y-1]);
        const float* nextPtr = (float*)(lutimg.data + lutimg.step*ymap[y+1]);
#else
        const uchar* imgPtr  = img.data + img.step*ymap[y];
        const uchar* prevPtr = img.data + img.step*ymap[y-1];
        const uchar* nextPtr = img.data + img.step*ymap[y+1];
#endif
        float* gradPtr = (float*)grad.ptr(y);	//Returns a pointer to the specified matrix row.
        uchar* qanglePtr = (uchar*)qangle.ptr(y);
		
		//计算 水平和垂直梯度 保存在 dbuf 的前两段
        if( cn == 1 )
        {
            for( x = 0; x < width; x++ )
            {
                int x1 = xmap[x];
#ifdef HAVE_IPP
                dbuf[x] = (float)(imgPtr[xmap[x+1]] - imgPtr[xmap[x-1]]);	//水平微分模板 [-1 0 1]
                dbuf[width + x] = (float)(nextPtr[x1] - prevPtr[x1]);		//垂直微分模板 [-1 0 1]'
#else
                dbuf[x] = (float)(lut[imgPtr[xmap[x+1]]] - lut[imgPtr[xmap[x-1]]]);		//dbuf length: width*4
                dbuf[width + x] = (float)(lut[nextPtr[x1]] - lut[prevPtr[x1]]);		//使用了IPP优化，就已经gamm校正了，这里是先gamm校正然后在计算梯度
#endif
            }
        }
        else	取B,G,R通道中梯度模最大的梯度作为该点的梯度，
        {
            for( x = 0; x < width; x++ )
            {
                int x1 = xmap[x]*3; //height*width*element,element:8UC3/32FC3
                float dx0, dy0, dx, dy, mag0, mag;
#ifdef HAVE_IPP
                const float* p2 = imgPtr + xmap[x+1]*3;	
                const float* p0 = imgPtr + xmap[x-1]*3;	
				//R通道的梯度
                dx0 = p2[2] - p0[2];
                dy0 = nextPtr[x1+2] - prevPtr[x1+2];	
                mag0 = dx0*dx0 + dy0*dy0;
				//G通道的梯度
                dx = p2[1] - p0[1];
                dy = nextPtr[x1+1] - prevPtr[x1+1];
                mag = dx*dx + dy*dy;
		
                if( mag0 < mag )	//取G,R通道中梯度模最大的
                {
                    dx0 = dx;
                    dy0 = dy;
                    mag0 = mag;
                }
				//B通道的梯度
                dx = p2[0] - p0[0];
                dy = nextPtr[x1] - prevPtr[x1];
                mag = dx*dx + dy*dy;
#else
                const uchar* p2 = imgPtr + xmap[x+1]*3;
                const uchar* p0 = imgPtr + xmap[x-1]*3;

                dx0 = lut[p2[2]] - lut[p0[2]];
                dy0 = lut[nextPtr[x1+2]] - lut[prevPtr[x1+2]];
                mag0 = dx0*dx0 + dy0*dy0;

                dx = lut[p2[1]] - lut[p0[1]];
                dy = lut[nextPtr[x1+1]] - lut[prevPtr[x1+1]];
                mag = dx*dx + dy*dy;

                if( mag0 < mag )
                {
                    dx0 = dx;
                    dy0 = dy;
                    mag0 = mag;
                }

                dx = lut[p2[0]] - lut[p0[0]];
                dy = lut[nextPtr[x1]] - lut[prevPtr[x1]];
                mag = dx*dx + dy*dy;
 #endif
                if( mag0 < mag )	//取B,G,R通道中梯度模最大的
                {
                    dx0 = dx;
                    dy0 = dy;
                    mag0 = mag;
                }

                dbuf[x] = dx0;
                dbuf[x+width] = dy0;
            }
        }
#ifdef HAVE_IPP
        ippsCartToPolar_32f((const Ipp32f*)Dx.data, (const Ipp32f*)Dy.data, (Ipp32f*)Mag.data, pAngles, width);
        for( x = 0; x < width; x++ )
        {
           if(pAngles[x] < 0.f)
             pAngles[x] += (Ipp32f)(CV_PI*2.);
        }

        ippsNormalize_32f(pAngles, pAngles, width, 0.5f/angleScale, 1.f/angleScale);
        ippsFloor_32f(pAngles,(Ipp32f*)hidxs.data,width);
        ippsSub_32f_I((Ipp32f*)hidxs.data,pAngles,width);
        ippsMul_32f_I((Ipp32f*)Mag.data,pAngles,width);

        ippsSub_32f_I(pAngles,(Ipp32f*)Mag.data,width);
        ippsRealToCplx_32f((Ipp32f*)Mag.data,pAngles,(Ipp32fc*)gradPtr,width);
#else
		//计算梯度的模和角度，默认结果为弧度
        cartToPolar( Dx, Dy, Mag, Angle, false );	//Calculates the magnitude and angle of 2D vectors. angle(I) = atan2(y(I); x(I))
#endif
        for( x = 0; x < width; x++ )
        {
#ifdef HAVE_IPP
            int hidx = (int)pHidxs[x];
#else
			//保存该梯度方向在左右相邻的bin的模，本来只有一个模何来的两个？插值！
			//线性插值，比如某点算出来应该属于 bin 7.6,但是我们的bin都是整数的，四舍五入，把他划分到bin 8又太粗糙了
			//那就按该点到bin7,bin8的距离分配，这样部分属于8，部分属于7。
            float mag = dbuf[x+width*2], angle = dbuf[x+width*3]*angleScale - 0.5f;	// 每一格 pi/9, 那现在算 t落在哪一格自然是 t/(pi/9)
            int hidx = cvFloor(angle);	//向下取整
            angle -= hidx;
            gradPtr[x*2] = mag*(1.f - angle);	//binx的大小是梯度方向和模的共同体现
            gradPtr[x*2+1] = mag*angle;
#endif
            if( hidx < 0 )
                hidx += _nbins;
            else if( hidx >= _nbins )
                hidx -= _nbins;
            assert( (unsigned)hidx < (unsigned)_nbins );
			
			//保存与该梯度方向相邻的左右两个bin编号
            qanglePtr[x*2] = (uchar)hidx;	//也是向下取整
            hidx++;
            hidx &= hidx < _nbins ? -1 : 0;	// hidx &= ( (hidx < _nbins ) ? -1 : 0;)，如果hidx < nbins good;如果超过了，就算子bin 0 ；-1的补码是全1
            qanglePtr[x*2+1] = (uchar)hidx;
        }
    }
}


struct HOGCache
{
    struct BlockData
    {
        BlockData() : histOfs(0), imgOffset() {}
        int histOfs;
        Point imgOffset;
    };

    struct PixData
    {
        size_t gradOfs, qangleOfs;
        int histOfs[4];
        float histWeights[4];
        float gradWeight;
    };

    HOGCache();
    HOGCache(const HOGDescriptor* descriptor,
        const Mat& img, Size paddingTL, Size paddingBR,
        bool useCache, Size cacheStride);
    virtual ~HOGCache() {};
    virtual void init(const HOGDescriptor* descriptor,
        const Mat& img, Size paddingTL, Size paddingBR,
        bool useCache, Size cacheStride);

    Size windowsInImage(Size imageSize, Size winStride) const;
    Rect getWindow(Size imageSize, Size winStride, int idx) const;

    const float* getBlock(Point pt, float* buf);
    virtual void normalizeBlockHistogram(float* histogram) const;

    vector<PixData> pixData;
    vector<BlockData> blockData;

    bool useCache;
    vector<int> ymaxCached;
    Size winSize, cacheStride;
    Size nblocks, ncells;
    int blockHistogramSize;
    int count1, count2, count4;
    Point imgoffset;
    Mat_<float> blockCache;
    Mat_<uchar> blockCacheFlags;

    Mat grad, qangle;
    const HOGDescriptor* descriptor;
};


HOGCache::HOGCache()
{
    useCache = false;
    blockHistogramSize = count1 = count2 = count4 = 0;
    descriptor = 0;
}

HOGCache::HOGCache(const HOGDescriptor* _descriptor,
        const Mat& _img, Size _paddingTL, Size _paddingBR,
        bool _useCache, Size _cacheStride)
{
    init(_descriptor, _img, _paddingTL, _paddingBR, _useCache, _cacheStride);
}

void HOGCache::init(const HOGDescriptor* _descriptor,
        const Mat& _img, Size _paddingTL, Size _paddingBR,
        bool _useCache, Size _cacheStride)
{
    descriptor = _descriptor;
    cacheStride = _cacheStride;
    useCache = _useCache;
	/*--------------------------------------计算梯度----------------------------------------------*/
	//返回值
	//size：img.cols + paddingTL.width + paddingBR.width,img.rows + paddingTL.height + paddingBR.height,类型 CV_32FC2
	//grad：梯度的模在与梯度方向相邻的两个bin的插值值
	//qangle：与梯度方向相邻的两个bin的编号
    descriptor->computeGradient(_img, grad, qangle, _paddingTL, _paddingBR);
    imgoffset = _paddingTL;

    winSize = descriptor->winSize;	//默认值：winSize(64,128)
    Size blockSize = descriptor->blockSize;//blockSize(16,16)
    Size blockStride = descriptor->blockStride;//lockStride(8,8)
    Size cellSize = descriptor->cellSize;//cellSize(8,8)
    int i, j, nbins = descriptor->nbins;//nbins(9)
    int rawBlockSize = blockSize.width*blockSize.height;

    nblocks = Size((winSize.width - blockSize.width)/blockStride.width + 1,
                   (winSize.height - blockSize.height)/blockStride.height + 1);
		   //这种算法非常直观，也许你会觉得可以和下面一样直接除，但是当(winSize.height - blockSize.height) % blockStride.height 不为0时，就不一定
		   //比如 blockSize=4,blockStride=3,winSize.width =9,那么直接除9/3=3，但是只能有两个block, 4|3|2,只能移动一次
    ncells = Size(blockSize.width/cellSize.width, blockSize.height/cellSize.height);
    blockHistogramSize = ncells.width*ncells.height*nbins;//默认2*2*9

    if( useCache )	//
    {
        Size cacheSize((grad.cols - blockSize.width)/cacheStride.width+1,
                       (winSize.height/cacheStride.height)+1);
        blockCache.create(cacheSize.height, cacheSize.width*blockHistogramSize);
        blockCacheFlags.create(cacheSize);
        size_t cacheRows = blockCache.rows;
        ymaxCached.resize(cacheRows);
        for(size_t ii = 0; ii < cacheRows; ii++ )
            ymaxCached[ii] = -1;
    }

    Mat_<float> weights(blockSize);//16*16 高斯模板
    float sigma = (float)descriptor->getWinSigma();//-1
    float scale = 1.f/(sigma*sigma*2);

    for(i = 0; i < blockSize.height; i++)
        for(j = 0; j < blockSize.width; j++)
        {
            float di = i - blockSize.height*0.5f;
            float dj = j - blockSize.width*0.5f;//中心
            weights(i,j) = std::exp(-(di*di + dj*dj)*scale);//weights(i,j)=exp(-(distance/sigma)^2)
        }

    blockData.resize(nblocks.width*nblocks.height);
    pixData.resize(rawBlockSize*3);// vector::resize(newsize,value),不是Mat::resize,16*16*3个结构体
	/* 
	vector<PixData> pixData;
	struct PixData{
        size_t gradOfs, qangleOfs;
        int histOfs[4];
        float histWeights[4];
        float gradWeight;
    };
	*/
    // Initialize 2 lookup tables, pixData & blockData.
    // Here is why:
    //
    // The detection algorithm runs in 4 nested loops (at each pyramid layer):
    //  loop over the windows within the input image
    //    loop over the blocks within each window
    //      loop over the cells within each block
    //        loop over the pixels in each cell
    //
    // As each of the loops runs over a 2-dimensional array,
    // we could get 8(!) nested loops in total, which is very-very slow.
    //
    // To speed the things up, we do the following:
    //   1. loop over windows is unrolled in the HOGDescriptor::{compute|detect} methods;
    //         inside we compute the current search window using getWindow() method.
    //         Yes, it involves some overhead (function call + couple of divisions),
    //         but it's tiny in fact.
    //   2. loop over the blocks is also unrolled. Inside we use pre-computed blockData[j]
    //         to set up gradient and histogram pointers.
    //   3. loops over cells and pixels in each cell are merged
    //       (since there is no overlap between cells, each pixel in the block is processed once)
    //      and also unrolled. Inside we use PixData[k] to access the gradient values and
    //      update the histogram
    //
    count1 = count2 = count4 = 0;
    for( j = 0; j < blockSize.width; j++ )//16,先水平，再垂直
        for( i = 0; i < blockSize.height; i++ )//16
        {
            PixData* data = 0;
            float cellX = (j+0.5f)/cellSize.width - 0.5f;	//这是干什么 ？？？
            float cellY = (i+0.5f)/cellSize.height - 0.5f;
            int icellX0 = cvFloor(cellX);	//-1（j=0..3）,0(j=4..11),1(j=12..15)
            int icellY0 = cvFloor(cellY);	
            int icellX1 = icellX0 + 1, icellY1 = icellY0 + 1;//0 1 2 
            cellX -= icellX0;
            cellY -= icellY0;
			
            if( (unsigned)icellX0 < (unsigned)ncells.width &&	// icellX0 == 0
                (unsigned)icellX1 < (unsigned)ncells.width )	//判断条件时特别小心，int 转成了 unsigned,(unsigned)(-1)=2^32-1，真对这作者无语
            {
				//  icellX0 == 0，icellY0 == 0 对相邻的四个cell都有贡献，即F,J,G,K区域
                if( (unsigned)icellY0 < (unsigned)ncells.height && // cellX，cellY 范围（0,1）
                    (unsigned)icellY1 < (unsigned)ncells.height )
                {
                    data = &pixData[rawBlockSize*2 + (count4++)];
                    data->histOfs[0] = (icellX0*ncells.height + icellY0)*nbins;//cell 0 在整个block的bin中的偏移
                    data->histWeights[0] = (1.f - cellX)*(1.f - cellY);	//到对称中心的“距离”即cell 3
                    data->histOfs[1] = (icellX1*ncells.height + icellY0)*nbins;//cell 1的偏移 2*9
                    data->histWeights[1] = cellX*(1.f - cellY);         //到对称中心的“距离”即 cell 2
                    data->histOfs[2] = (icellX0*ncells.height + icellY1)*nbins;//cell 2的偏移 1*9
                    data->histWeights[2] = (1.f - cellX)*cellY;         //到对称中心的“距离”即 cell 1
                    data->histOfs[3] = (icellX1*ncells.height + icellY1)*nbins;//cell 3的偏移3*9
                    data->histWeights[3] = cellX*cellY;                 //到对称中心的“距离”即 cell 0
                }
                else	// icellX0 == 0，icellY0 == -1/1，对左右相邻的两个cell有贡献，即B,C,N,O
                {
					// cellX 范围（0,1），cellY 范围 （0.5,1）/（0,0.5）
                    data = &pixData[rawBlockSize + (count2++)];
					//下部分的cellY范围也落在（0.5,1），icellY1==icellY0 == 1
                    if( (unsigned)icellY0 < (unsigned)ncells.height )//icellY0 == 1
                    {
                        icellY1 = icellY0;
                        cellY = 1.f - cellY;
                    }
                    data->histOfs[0] = (icellX0*ncells.height + icellY1)*nbins;// 上部分0；下部分1
                    data->histWeights[0] = (1.f - cellX)*cellY;
                    data->histOfs[1] = (icellX1*ncells.height + icellY1)*nbins;// 上部分2；下部分3
                    data->histWeights[1] = cellX*cellY;
                    data->histOfs[2] = data->histOfs[3] = 0;// 均为0
                    data->histWeights[2] = data->histWeights[3] = 0;
                }
            }
            else //icellX0 == -1/1,cellX范围（0.5,1）/（0,0.5）
            {
				//右部分的cellX范围也落在（0.5,1），icellX1==icellX0 == 1
                if( (unsigned)icellX0 < (unsigned)ncells.width )
                {
                    icellX1 = icellX0;
                    cellX = 1.f - cellX;
                }
				//E,H,I,L
                if( (unsigned)icellY0 < (unsigned)ncells.height &&
                    (unsigned)icellY1 < (unsigned)ncells.height )
                {
                    data = &pixData[rawBlockSize + (count2++)];
                    data->histOfs[0] = (icellX1*ncells.height + icellY0)*nbins;//左：0，右：2*9
                    data->histWeights[0] = cellX*(1.f - cellY);
                    data->histOfs[1] = (icellX1*ncells.height + icellY1)*nbins;//左：1*9 右：3*9
                    data->histWeights[1] = cellX*cellY;
                    data->histOfs[2] = data->histOfs[3] = 0;
                    data->histWeights[2] = data->histWeights[3] = 0;
                }
				// A,D,M,P
                else
                {
					data = &pixData[count1++];
                    if( (unsigned)icellY0 < (unsigned)ncells.height )
                    {
                        icellY1 = icellY0;
                        cellY = 1.f - cellY;
                    }
                    data->histOfs[0] = (icellX1*ncells.height + icellY1)*nbins;
                    data->histWeights[0] = cellX*cellY;
                    data->histOfs[1] = data->histOfs[2] = data->histOfs[3] = 0;
                    data->histWeights[1] = data->histWeights[2] = data->histWeights[3] = 0;
                }
            }
            data->gradOfs = (grad.cols*i + j)*2;	//block窗口的（0,0）位置有相对于整个图像的偏移，此偏移为相对于block(0,0)的偏移
            data->qangleOfs = (qangle.cols*i + j)*2;//计算方式很古怪，但是你画张图就明白了（grad.cols*i多算的==+j少算的），实际上 block窗口的（0,0）的offset加上此offset就可以直接在grad中找到对应的梯度
            data->gradWeight = weights(i,j);	//该点的高斯权值，大小与到block中心的距离成反比
        }

    assert( count1 + count2 + count4 == rawBlockSize );//16*16
    // defragment pixData,整理碎片.
	//数据合并  xxx.........yyy.........zzz.........->xxxyyyzzz..................
	//（.表示未赋值空间，x为count1存储的数据，y为count2存储的数据...）
    for( j &#