docs: public void features(double c_x, double c_y,
  						 double c , double s,
  						 double scale,  SparseImageGradient gradient , double[] features)
  	{
  		int regionSize = widthLargeGrid*widthSubRegion;
  
  		int totalSampleWidth = widthSubRegion+overLap*2;
  
  		int regionR = regionSize/2;
  		int regionEnd = regionSize-regionR;
  
  		int sampleGridWidth = regionSize+2*overLap;
  		int regionIndex = 0;
  
  		// when computing the pixel coordinates it is more precise to round to the nearest integer
  		// since pixels are always positive round() is equivalent to adding 0.5 and then converting
  		// to an int, which floors the variable.
  		c_x += 0.5;
  		c_y += 0.5;
  
  		// first sample the whole grid at once to avoid sampling overlapping regions twice
  		int index = 0;
  		for( int rY = -regionR-overLap; rY < regionEnd+overLap; rY++) {
  			double regionY = rY*scale;
  			for( int rX = -regionR-overLap; rX < regionEnd+overLap; rX++,index++ ) {
  				double regionX = rX*scale;
  
  				// rotate the pixel along the feature's direction
  				int pixelX = (int)(c_x + c*regionX - s*regionY);
  				int pixelY = (int)(c_y + s*regionX + c*regionY);
  
  				GradientValue g = gradient.compute(pixelX,pixelY);
  				samplesX[index] = g.getX();
  				samplesY[index] = g.getY();
  			}
  		}
  
  		// compute descriptor using precomputed samples
  		int indexGridWeight = 0;
  		for( int rY = -regionR; rY < regionEnd; rY += widthSubRegion ) {
  			for( int rX = -regionR; rX < regionEnd; rX += widthSubRegion ) {
  				double sum_dx = 0, sum_dy=0, sum_adx=0, sum_ady=0;
  
  				// compute and sum up the response  inside the sub-region
  				for( int i = 0; i < totalSampleWidth; i++ ) {
  					index = (rY+regionR+i)*sampleGridWidth + rX+regionR;
  					for( int j = 0; j < totalSampleWidth; j++ , index++ ) {
  						double w = weightSub.get(j,i);
  
  						double dx = w*samplesX[index];
  						double dy = w*samplesY[index];
  
  						// align the gradient along image patch
  						// note the transform is transposed
  						double pdx =  c*dx + s*dy;
  						double pdy = -s*dx + c*dy;
  
  						sum_dx += pdx;
  						sum_adx += Math.abs(pdx);
  						sum_dy += pdy;
  						sum_ady += Math.abs(pdy);
  					}
  				}
  
  				double w = weightGrid.data[indexGridWeight++];
  				features[regionIndex++] = w*sum_dx;
  				features[regionIndex++] = w*sum_adx;
  				features[regionIndex++] = w*sum_dy;
  				features[regionIndex++] = w*sum_ady;
  			}
  		}
  	}