Is the OpenCV phase() function working correctly in this edge case?

I'm using OpenCV 3.2.0 to do some Fourier space calculations. To obtain a phase image after an inverse DFT, I tried using cv::phase() but I noticed that in some cases, it returned values close to 2*Pi where it should (in my eyes) return a value close to zero. I wonder if this function is implemented badly or if I'm using it wrong.

This is my example data, a 7x8 FFT where the imaginary part is zero or, due to rounding errors, very close to zero (value pairs in the form of real, imag):

 0.75686288, 0, 0.74509817, -3.6017641e-19, 0.74117655, -4.8023428e-19, 0.76078451, -1.3206505e-18, 0.77647072, 0, 0.74509817, -3.6017641e-19, 0.72549027, 4.8023428e-19, 0.70588243, 2.0410032e-18;
 0.70980388, 0, 0.66666675, -6.6032515e-19, 0.69803929, -3.8418834e-18, 0.73725492, -5.3426161e-18, 0.69803923, 0, 0.6549021, -6.6032515e-19, 0.5725491, 3.8418834e-18, 0.5411765, 6.6632662e-18;
 0.63529414, 0, 0.6352942, -1.7408535e-18, 0.63921577, -5.1625314e-18, 0.61960787, -3.1815585e-18, 0.60784316, 0, 0.55686277, -1.7408535e-18, 0.4705883, 5.1625314e-18, 0.45882356, 6.6632657e-18;
 0.58039224, 0, 0.58431381, -6.6032412e-19, 0.63921583, -7.8038246e-18, 0.63921577, -7.9839117e-18, 0.50196087, 0, 0.45490205, -6.6032412e-19, 0.38431379, 7.8038246e-18, 0.35686284, 9.3045593e-18;
 0.54117656, 0, 0.58431375, -9.0044183e-19, 0.68627465, -9.1244722e-18, 0.6156863, -6.7833236e-18, 0.48627454, 0, 0.45490202, -9.0044183e-19, 0.38823539, 9.1244722e-18, 0.36470592, 8.5842074e-18;
 0.50980395, 0, 0.56470597, -6.0029469e-19, 0.57254916, -4.8023546e-18, 0.54901963, -3.9619416e-18, 0.4784314, 0, 0.42352945, -6.0029469e-19, 0.41568634, 4.8023546e-18, 0.39999998, 5.162531e-18;
 0.49411768, 0, 0.50588238, 4.8023392e-19, 0.54509813, -1.6808249e-18, 0.56078434, -3.241587e-18, 0.49803928, 0, 0.49411774, 4.8023392e-19, 0.49019611, 1.6808249e-18, 0.47058827, 2.2811191e-18

I then applied cv::phase() like this:

Mat planes[2];
split(output,planes);
Mat ph;
phase(planes[0],planes[1],ph);

Then, cout<<ph yields:

 0, 6.2831855, 6.2831855, 6.2831855, 0, 6.2831855, 6.6180405e-19, 2.8908079e-18;
 0, 6.2831855, 6.2831855, 6.2831855, 0, 6.2831855, 6.7087144e-18, 1.2309944e-17;
 0, 6.2831855, 6.2831855, 6.2831855, 0, 6.2831855, 1.096805e-17, 1.451942e-17;
 0, 6.2831855, 6.2831855, 6.2831855, 0, 6.2831855, 2.0301558e-17, 2.6067677e-17;
 0, 6.2831855, 6.2831855, 6.2831855, 0, 6.2831855, 2.3497438e-17, 2.3532349e-17;
 0, 6.2831855, 6.2831855, 6.2831855, 0, 6.2831855, 1.1550381e-17, 1.2903592e-17;
 0, 9.4909814e-19, 6.2831855, 6.2831855, 0, 9.7169579e-19, 3.4281555e-18, 4.8463495e-18

So the output is sort of oscillating between the lowest and highest value. I was awaiting a matrix of (near) zeros though, because a non-existing phase shift would be in line with the underlying physics application. I then tried computing the phase image pixel by pixel:

Mat_<double> myPhase = Mat_<double>(8,7);
for(int i = 0; i < fftReal.rows; i++) {
  for(int j = 0; j < fftReal.cols; j++) {
      float fftRealVal = planes[0].at<float>(i,j);
      float fftImagVal = planes[1].at<float>(i,j);
      double angle = atan2(fftImagVal, fftRealVal);
      myPhase(i,j) = angle;
  }

Here, the output of cout<<myPhase is what I expected to see, a matrix of near zeros:

 0, -4.833945789050036e-19, -6.479350716073673e-19, -1.735906137457605e-18, 0, -4.833945789050036e-19, 6.619444609555068e-19;
 0, -9.904875721669217e-19, -5.503821154321125e-18, -7.246633215917781e-18, 0, -1.00828074413082e-18, 6.710137932686301e-18;
 0, -2.740232027682232e-18, -8.076351618590122e-18, -5.13479354918468e-18, 0, -3.126180439429062e-18, 1.097037782204674e-17;
 0, -1.130084743690479e-18, -1.220843476128649e-17, -1.249016668765776e-17, 0, -1.451574279023501e-18, 2.030586625060691e-17;
 0, -1.541024556489219e-18, -1.329565697018843e-17, -1.101749982000204e-17, 0, -1.979419217601631e-18, 2.350242300975683e-17;
 0, -1.063021695913201e-18, -8.387671795472417e-18, -7.216393147084068e-18, 0, -1.417362295683461e-18, 1.155283208729227e-17;
 0, 9.492995611309157e-19, -3.083527284733071e-18, -5.78045227144509e-18, 0, 9.719018577786209e-19, 3.428882635099473e-18;
 4.847377651234249e-18, 6.937607420147441e-310, 6.937607420153765e-310, 6.93760742011582e-310, 6.93760742011503e-310, 6.937607420163251e-310, 6.937607420188547e-310

So does cv::phase() yield a wrong result here due to some rounding errors or does it work as it should and I'm missing some pre-processing or anything?