docs: private float iterationSor(GrayF32 image1, GrayF32 deriv1X, GrayF32 deriv1Y,
  							   int i, int ipx, int imx, int ipy, int imy) {
  		float w = SOR_RELAXATION;
  
  		// these variables could be precomputed once.  See equation 11
  		float psid = psiData.data[i];
  		float psig = gamma*psiGradient.data[i];
  
  		float di = warpImage2.data[i] - image1.data[i];
  		float dx2 = warpDeriv2X.data[i];
  		float dy2 = warpDeriv2Y.data[i];
  		float dxx2 = warpDeriv2XX.data[i];
  		float dyy2 = warpDeriv2YY.data[i];
  		float dxy2 = warpDeriv2XY.data[i];
  
  		float Au = -psid*di*dx2 + alpha*divU.data[i]
  				- psig*((dx2 - deriv1X.data[i])*warpDeriv2XX.data[i] + (dy2 - deriv1Y.data[i])*warpDeriv2XY.data[i]);
  
  		float Av = -psid*di*dy2 + alpha*divV.data[i]
  				- psig*((dx2 - deriv1X.data[i])*warpDeriv2XY.data[i] + (dy2 - deriv1Y.data[i])*warpDeriv2YY.data[i]);
  
  		float Du = psid*dx2*dx2 + psig*(dxx2*dxx2 + dxy2*dxy2) + alpha*divD.data[i];
  		float Dv = psid*dy2*dy2 + psig*(dyy2*dyy2 + dxy2*dxy2) + alpha*divD.data[i];
  		float D = psid*dx2*dy2 + psig*(dxx2 + dyy2)*dxy2;
  
  		// update the change in flow
  		float psi_index = psiSmooth.data[i];
  		float coef0 = 0.5f*(psiSmooth.data[ipx] + psi_index);
  		float coef1 = 0.5f*(psiSmooth.data[imx] + psi_index);
  		float coef2 = 0.5f*(psiSmooth.data[ipy] + psi_index);
  		float coef3 = 0.5f*(psiSmooth.data[imy] + psi_index);
  
  		float div_du = coef0 * du.data[ipx] + coef1 * du.data[imx] + coef2 * du.data[ipy] + coef3 * du.data[imy];
  		float div_dv = coef0 * dv.data[ipx] + coef1 * dv.data[imx] + coef2 * dv.data[ipy] + coef3 * dv.data[imy] ;
  
  		final float dui = du.data[i];
  		final float dvi = dv.data[i];
  
  		du.data[i] = (1f-w)*dui + w*(Au - D*dvi + alpha*div_du)/Du;
  		dv.data[i] = (1f-w)*dvi + w*(Av - D*du.data[i] + alpha*div_dv)/Dv;
  
  		return (du.data[i] - dui) * (du.data[i] - dui) + (dv.data[i] - dvi) * (dv.data[i] - dvi);
  	}