You can calculate the circulation using the bounds of a contour in 2018b.
I completed this on cylindrical slices for my purposes, but I'm sure you can adjust this. You will need velocity data, not vorticity to use the code I have provided but the results will be equivalent.
R = 1;
Nt = floor(2.1*50*R);
[Xc,Zc,Yc] = cylinder(R*ones(1,50),Nt);
surf(Xc,Yc,Zc)
hsli = contourslice(Xw,Yw,Zw,D_ta,Xc,Yc,Zc,[-3,-3],'linear');
% D_ta is the name of my vorticity field in this case.
% You can always use other terms to define your contour, or just choose an arbitrary shape.
if ~isempty(hsli)
hsli.LineWidth = 2;
hsli.EdgeColor = [144,144,144]/255;
axis equal
Xcs = hsli.XData; Xcs = Xcs(1:end-1);
Ycs = hsli.YData; Ycs = Ycs(1:end-1);
Zcs = hsli.ZData; Zcs = Zcs(1:end-1);
ux_loc = interp3(Xw,Yw,Zw,Ux,Xcs,Ycs,Zcs);
uy_loc = interp3(Xw,Yw,Zw,Uy,Xcs,Ycs,Zcs);
uz_loc = interp3(Xw,Yw,Zw,Uz,Xcs,Ycs,Zcs);
GammaLoc = zeros(1,length(Xcs));
for g = 1:length(Xcs)-1
dl = [Xcs(g+1)-Xcs(g),...
Ycs(g+1)-Ycs(g),...
Zcs(g+1)-Zcs(g)];
GammaLoc(g) = dot([ux_loc(g),uy_loc(g),uz_loc(g)],...
dl);
end
Gamma = sum(GammaLoc);
end