docs: public static void hessianInner(GrayS32 integral, int skip , int size ,
  									GrayF32 intensity)
  	{
  		final int w = intensity.width;
  		final int h = intensity.height;
  
  		float norm = 1.0f/(size*size);
  
  		int blockSmall = size/3;
  		int blockLarge = size-blockSmall-1;
  		int radiusFeature = size/2;
  		int radiusSkinny = blockLarge/2;
  
  		int blockW2 = 2*blockSmall;
  		int blockW3 = 3*blockSmall;
  
  
  		int rowOff1 = blockSmall*integral.stride;
  		int rowOff2 = 2*rowOff1;
  		int rowOff3 = 3*rowOff1;
  
  		// make sure it starts on the correct pixel
  		final int borderOrig = radiusFeature+ 1 + (skip-(radiusFeature+1)%skip);
  		final int border = borderOrig/skip;
  		final int lostPixel = borderOrig - radiusFeature-1;
  		final int endY = h - border;
  		final int endX = w - border;
  
  		for( int y = border; y < endY; y++ ) {
  
  			// pixel location in original input image
  			int yy = y*skip;
  
  			// index for output
  			int indexDst = intensity.startIndex + y*intensity.stride+border;
  
  			// indexes for Dxx
  			int indexTop = integral.startIndex + (yy-radiusSkinny-1)*integral.stride+lostPixel;
  			int indexBottom = indexTop + (blockLarge)*integral.stride;
  
  			// indexes for Dyy
  			int indexL = integral.startIndex + (yy-radiusFeature-1)*integral.stride + (radiusFeature-radiusSkinny)+lostPixel;
  			int indexR = indexL + blockLarge;
  
  			// indexes for Dxy
  			int indexY1 = integral.startIndex + (yy-blockSmall-1)*integral.stride + (radiusFeature-blockSmall)+lostPixel;
  			int indexY2 = indexY1 + blockSmall*integral.stride;
  			int indexY3 = indexY2 + integral.stride;
  			int indexY4 = indexY3 + blockSmall*integral.stride;
  
  			for( int x = border; x < endX; x++ , indexDst++) {
  				float Dxx = integral.data[indexBottom+blockW3] - integral.data[indexTop+blockW3] - integral.data[indexBottom] + integral.data[indexTop];
  				Dxx -= 3*(integral.data[indexBottom+blockW2] - integral.data[indexTop+blockW2] - integral.data[indexBottom+blockSmall] + integral.data[indexTop+blockSmall]);
  
  				float Dyy = integral.data[indexR+rowOff3] - integral.data[indexL+rowOff3] - integral.data[indexR] + integral.data[indexL];
  				Dyy -= 3*(integral.data[indexR+rowOff2] - integral.data[indexL+rowOff2] - integral.data[indexR+rowOff1] + integral.data[indexL+rowOff1]);
  
  				int x3 = blockSmall+1;
  				int x4 = x3+blockSmall;
  
  				float Dxy = integral.data[indexY2+blockSmall] - integral.data[indexY1+blockSmall] - integral.data[indexY2] + integral.data[indexY1];
  				Dxy -= integral.data[indexY2+x4] - integral.data[indexY1+x4] - integral.data[indexY2+x3] + integral.data[indexY1+x3];
  				Dxy += integral.data[indexY4+x4] - integral.data[indexY3+x4] - integral.data[indexY4+x3] + integral.data[indexY3+x3];
  				Dxy -= integral.data[indexY4+blockSmall] - integral.data[indexY3+blockSmall] - integral.data[indexY4] + integral.data[indexY3];
  
  				Dxx *= norm;
  				Dxy *= norm;
  				Dyy *= norm;
  
  				intensity.data[indexDst] = Dxx*Dyy-0.81f*Dxy*Dxy;
  
  				indexTop += skip;
  				indexBottom += skip;
  				indexL += skip;
  				indexR += skip;
  				indexY1 += skip;
  				indexY2 += skip;
  				indexY3 += skip;
  				indexY4 += skip;
  			}
  		}
  	}