docs: private void computeDivUVD(GrayF32 u , GrayF32 v , GrayF32 psi ,
  							   GrayF32 divU , GrayF32 divV , GrayF32 divD ) {
  
  		final int stride = psi.stride;
  
  		// compute the inside pixel
  		for (int y = 1; y < psi.height-1; y++) {
  
  			// index of the current pixel
  			int index = y*stride + 1;
  
  			for (int x = 1; x < psi.width-1; x++ , index++) {
  
  				float psi_index = psi.data[index];
  
  				float coef0 = 0.5f*(psi.data[index+1] + psi_index);
  				float coef1 = 0.5f*(psi.data[index-1] + psi_index);
  				float coef2 = 0.5f*(psi.data[index+stride] + psi_index);
  				float coef3 = 0.5f*(psi.data[index-stride] + psi_index);
  
  				float u_index = u.data[index];
  
  				divU.data[index] = coef0*(u.data[index+1] - u_index) + coef1*(u.data[index-1] - u_index) +
  								   coef2*(u.data[index+stride] - u_index) + coef3*(u.data[index-stride] - u_index);
  
  				float v_index = v.data[index];
  
  				divV.data[index] = coef0*(v.data[index+1] - v_index) + coef1*(v.data[index-1] - v_index) +
  								   coef2*(v.data[index+stride] - v_index) + coef3*(v.data[index-stride] - v_index);
  
  				divD.data[index] = coef0 + coef1 + coef2 + coef3;
  			}
  		}
  
  		// handle the image borders
  		for( int x = 0; x < psi.width; x++ ) {
  			computeDivUVD_safe(x,0,u,v,psi,divU,divV,divD);
  			computeDivUVD_safe(x,psi.height-1,u,v,psi,divU,divV,divD);
  		}
  		for( int y = 1; y < psi.height-1; y++ ) {
  			computeDivUVD_safe(0,y,u,v,psi,divU,divV,divD);
  			computeDivUVD_safe(psi.width-1,y,u,v,psi,divU,divV,divD);
  		}
  	}