Volume in dclab<=0.36.1 is overestimated by on average 2µm³

[B-Hartmann]

The method counter_clockwise() in the file dclab/features/volume is supposed to return contour coordinates in a counter-clockwise order, yet the returned coordinates are in clockwise order.

minimal working example:

import numpy as np
cx = np.array([+1, +1, -1, -1])
cy = np.array([+1, -1, -1, +1])

angles = np.unwrap(np.arctan2(cy, cx))
grad = np.gradient(angles)
if np.average(grad) > 0:
    x = cx[::-1]
    y = cy[::-1]
else:
    x = cx
    y = cy
print(x)
print(y)

gives the output:
[ 1 1 -1 -1]
[ 1 -1 -1 1]

[paulmueller]

This is an unsettling find. According to https://de.mathworks.com/matlabcentral/fileexchange/36525-volrevolve (which this implementation is based on), this is what happens if you pass a clock-wise contour:

%  - R and Z vectors must be in order, counter-clockwise around the area
%    being defined. If not, this will give the volume of the
%    counter-clockwise parts, minus the volume of the clockwise parts.

@maikherbig Do you have any insights?

[paulmueller]

After comparing with the original matlab script, It also appears as if this line

https://github.com/ZELLMECHANIK-DRESDEN/dclab/blob/aade02014f4807053b2c2a668474dcdcdca9588a/dclab/features/volume.py#L152

is actually supposed to read

v3 = -1 * d_r**2 * z_vec_m1

I'm afraid this calls for thorough testing, starting with these tests in the original script:

clear all
R0 = 5 ;
a = 1 ;
npoints = 100 ;
theta = 2*pi*[0:1:npoints-1]'/double(npoints-1) ;
R = R0 + a*cos(theta) ;
Z =      a*sin(theta) ;
vol_analytic = 2 * pi^2 * R0 * a^2 ;
% >> 98.6960
vol = volRevolve(R,Z) ;
%  >> 98.6298 (6.7e-04 relative error)
% Do it again with npoints = 1000, get:
%  >> 98.6954 (6.6e-06 relative error)
% As expected, it's always slightly small because the polygon inscribes the
% circle.

% Verification for washer (rectangular toroid), with the radius of the
% 'hole' in the washer being a, and the outer radius of the washer being b.
% (Thus the width of the metal cross section is b-a.) The height of the
% washer is h. Then the volume is pi * (b^2 - a^2) * h. Run this code:
clear all
a = 1 ;
b = 2 ;
h = 10 ;
R = [a; b; b; a; a] ;
Z = [0; 0; h; h; 0] ;
vol_analytic = pi * (b^2 - a^2) * h ;
% >> 94.2478
vol = volRevolve(R,Z) ;
% >> 94.2478

[maikherbig]

Hi,
assuming ShapeOut and dclab are limited (so)RT(F)DC data, this is not a problem since the contours, returned by cv2.findContours (in ShapeIn) will always be counterclockwise (because only the outer contour is returned -> cv2.retr_external). In case, ShapeIn gets an update where it starts to return inner contours too (e.g. for assessing the shape of bubbles or the nucleus inside cells) it would actually cause a problem as OpenCV returns clockwise inner contours. However, computing the volume for inner contours does not really make sense as rotational symmetry cannot be assumed.
In conclusion, one could either simply delete line 85
https://github.com/ZELLMECHANIK-DRESDEN/dclab/blob/97f9f4ebeb84fd6f7cf7d239874a8675b5d90feb/dclab/features/volume.py#L85
or fix the function def counter_clockwise(cx, cy), or insert a test if a contour is counter-clockwise (return np.nan if not).

[paulmueller]

Another issue is this line:
https://github.com/ZELLMECHANIK-DRESDEN/dclab/blob/aade02014f4807053b2c2a668474dcdcdca9588a/dclab/features/volume.py#L146

After taking a look at the test function, I believe that this is a misinterpretation of what "R" is. It is not the length of the vector (x,z). It is, in fact, x.

[paulmueller]

What's also fishy is the abs part when returning the computed volume:

https://github.com/ZELLMECHANIK-DRESDEN/dclab/blob/aade02014f4807053b2c2a668474dcdcdca9588a/dclab/features/volume.py#L157

[paulmueller]

To answer the original issue, the counter_clockwise function indeed computed a clockwise contour. But that was correct, because the x- and y- coordinates were flipped when passing it to the function that computed the volume (channel axis is x, but contour revolution takes place around z). If you switch the coordinates of the contour afterwards, you get the counter-clockwise contour.

But, there was still an issue: The lower part of the contour had negative values. As a result that part always returned negative values for the volume - hence the fishy abs I mentioned above.

[paulmueller]

Here are the first results using these reference datasets: https://dcor.mpl.mpg.de/dataset/figshare-7771184-v2

This means that on average, the previous implementation overestimated the volume by roughly 2µm³ on average.

Here the same plot, normalized by the new volume (roughly 1-1.5 % overestimation in previous implementation on average):