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);
		}
	}
private void computeDivUVD ( GrayF32 u , GrayF32 v , GrayF32 psi , GrayF32 divU , GrayF32 divV , GrayF32 divD ) { final int stride = psi . stride ; for ( int y = 1 ; y < psi . height - 1 ; y ++ ) { 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 ; } } 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 ) ; } }
Computes the divergence for u,v, and d. Equation 8 and Equation 10.