Reproducing Raicevic et al. (2023) from dandiset 000491¶
Paper. Raicevic N, Forer JM, Ladrón-de-Guevara A, Du T, Nedergaard M, Kelley DH, Boster K (2023). Sizes and shapes of perivascular spaces surrounding murine pial arteries. Research Square preprint, version 1, doi:10.21203/rs.3.rs-2587250/v1. The claims file was written from the full text of this preprint as served by PubMed Central (PMC9949243), with its figure legends and Tables 1 to 3; the Additional Files were not in the text. The paper links no code.
Dandiset. DANDI:000491, version 0.230602.1307 (doi:10.48324/dandi.000491/0.230602.1307). The paper cites no version. This is the most recent published version, and the only one that can support the analysis: the earlier version 0.230427.1809 holds no PVS mask for locations L, M and N (the data check below lists the masks each version holds). The draft has the same 14 assets and size.
What the dandiset holds. For each of the 14 pial artery segments (locations A to N, 9 mice), the 3D two-photon stack and the final binary segmentations of the vessel and of the PVS. The 2D time series used to validate the segmentation are not in it. Everything here starts from the deposited segmentations, as the claims file says; the paper's segmentation includes manual steps that cannot be redone.
Method. The pipeline P1 to P9 of the claims file's opening comment, in pvs.py: isotropic
sampling, a 3D skeleton of the vessel cleaned of erroneous points, second-order polynomial fits of
the centerline, normal cross sections, the PVS lobes beside the vessel in each, the edge exclusion
from the vessel center, areas, the circle fit of the vessel, the ellipse, spline and
polynomial fits, and the hydraulic resistance by a finite-difference solution of the Poisson
equation (in place of MATLAB's finite elements). Every reading the claims file names is run, with
one exception stated in the attempt record: reading (b) of P2 (all skeleton points, 10^4 to 10^5
per segment) is run for the counts only. No random numbers are used.
Revision after the first verification. The verifier found the longest path through the skeleton much longer than the vessel (it runs into branches and doubles back) and the PVS of a cross section taken over the whole plane. The path is now cleaned of erroneous points, the PVS of a cross section is the two largest PVS regions beside the vessel, planes that the PVS mask fills edge to edge (E and F) are cleared, and the edge distance is measured from the vessel center in the slice. The data check shows the path before and after cleaning, the PVS content of the cross sections and one cross section per location with the fits drawn over it.
Revision after the second verification. The data check shows how far the cleaned path still meanders, cross-section-counts adds the counts a path along the fitted centerline would give, and a cell before the summary computes every claim's table again without location E, whose cleaned path follows a side branch. Statements the verifier found unsupported are corrected.
Third verification. The verifier reviewed the proposed change of cross-section-counts to reproduced and did not agree; the outcome stays not reproduced.
Answer from a reader. A question on fig8-ellipse-parameters asked whether the ellipse fit is the PVS's alone or that of the PVS and vessel together. The reader amended the claims file (questions/000491.yaml): the ellipse is now fit to the region the PVS's outer boundary encloses, the PVS and vessel segmentations together with holes filled, for every claim that uses it (fig5d, table2, fig5b, ellipse-annulus-underestimate, fig7 and fig8). This notebook uses that ellipse; the ellipse of the PVS alone is shown as a variant under fig8.
Run time. About 29 minutes on four cores, almost all in the cross sections and the resistance solves.
| Claim | Tested |
|---|---|
| cross-section-counts | tested |
| segmentation-validation-overlap | blocked: data missing (no 2D time series) |
| fig5a-pvs-area-scales-with-vessel-area | tested |
| fig5e-area-ratio-distribution | tested |
| fig5bc-area-intersegment-variation | tested (blocked by its criterion: methods underspecified) |
| area-ratio-per-location | tested |
| fig5d-polynomial-closest-area-ratio | tested |
| table2-fit-area-ratio-peaks | tested |
| fig6a-resistance-power-law | tested |
| fig6b-resistance-segmentation-peak | tested |
| resistance-per-location | tested |
| fig5b-polynomial-closest-resistance | tested |
| ellipse-annulus-underestimate | tested |
| spline-resistance-five-times | tested |
| fig7-velocity-profiles | tested |
| vessel-radius-peak | tested |
| fig8-ellipse-parameters | tested |
| fig9c-pvs-two-lobes | tested |
| fig9e-pvs-width-asymmetry | tested |
| left-right-no-difference | tested (blocked by its criterion: methods underspecified) |
| fig9a-heights-taper | tested |
| fig9f-pvs-concave | tested |
In [1]:
import contextlib
import io
import os
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import curve_fit
import pvs
pd.set_option("display.width", 250)
pd.set_option("display.max_columns", 60)
pd.set_option("display.max_rows", 400)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=pd.errors.PerformanceWarning)
os.environ["OMP_NUM_THREADS"] = "1"
assets = pvs.list_assets(pvs.VERSION)
started = time.time()
TARGET = {"a": 0.648, "b": 1.0} # claims P1 readings: isotropic voxel size in um
DATA, META = {}, {}
for p1, target in TARGET.items():
t0 = time.time()
DATA[p1], META[p1] = pvs.run_all(assets, target)
print(f"P1 reading {p1} ({target} um): {len(DATA[p1])} cross sections from {DATA[p1].location.nunique()} "
f"files in {time.time() - t0:.0f} s")
COUNTS_B = {}
for p1, target in TARGET.items():
t0 = time.time()
COUNTS_B[p1] = pvs.run_counts_b(assets, target)
print(f"P2 reading b at P1 {p1}: {sum(r['n_points'] for r in COUNTS_B[p1])} skeleton points in {time.time() - t0:.0f} s")
print(f"dandiset {pvs.DANDISET} version {pvs.VERSION}, {len(assets)} assets, {time.time() - started:.0f} s in all")
P1 reading a (0.648 um): 7534 cross sections from 14 files in 999 s
P1 reading b (1.0 um): 4974 cross sections from 14 files in 321 s
P2 reading b at P1 a: 438100 skeleton points in 74 s
P2 reading b at P1 b: 179969 skeleton points in 35 s dandiset 000491 version 0.230602.1307, 14 assets, 1429 s in all
Data check¶
Each check of the metadata against the values is stated under its output. The stacks are static 3D images, with no time base and no trial table, so there is no clock to align; the check that matters here is the geometry (voxel size, depth direction, which mask is which).
In [2]:
def depth_offset(ves, pv):
"""Per (y, x) column holding both masks, the PVS's mean plane minus the vessel's mean plane."""
z = np.arange(ves.shape[0])[:, None, None]
vz = (z * ves).sum(0) / np.maximum(ves.sum(0), 1)
pz = (z * pv).sum(0) / np.maximum(pv.sum(0), 1)
both = ves.any(0) & pv.any(0)
return (pz - vz)[both]
rows = []
for a in assets:
f = pvs.open_asset(a)
ves, pv = pvs.load_masks(f)
vx = pvs.voxel_size(f)
zv = np.nonzero(ves.any((1, 2)))[0]
zp = np.nonzero(pv.any((1, 2)))[0]
acq = f["acquisition"]
# vessel intensity against depth, on 8 planes through the vessel (ChanA, the vessel channel)
d = acq["TwoPhotonSeriesChanA/data"]
zs = np.linspace(zv[0], zv[-1], 8).astype(int)
prof = [float(d[z, :, :, 0][ves[z]].mean()) for z in zs]
rho = stats.spearmanr(zs, prof).statistic
rows.append({
"location": pvs.location_of(f), "subject": f["general/subject/subject_id"][()].decode(),
"strain": f["general/subject/strain"][()].decode(), "sex": f["general/subject/sex"][()].decode(),
"age": f["general/subject/age"][()].decode(), "planes": ves.shape[0], "shape_yx": ves.shape[1:],
"voxel_um": vx, "channels": ",".join(k[-1] for k in acq), "image_dtype": str(d.dtype),
"masks": ",".join(k[-1] for k in f["processing/ophys/ImageSegmentation"]),
"vessel_vox": int(ves.sum()), "pvs_vox": int(pv.sum()), "overlap_vox": int((ves & pv).sum()),
"pvs_frac_volume": round(pv.mean(), 3), "vessel_z": (int(zv[0]), int(zv[-1])),
"pvs_z": (int(zp[0]), int(zp[-1])), "deepest_plane_pvs_frac": round(pv[zp[-1]].mean(), 2),
"intensity_vs_depth_rho": round(rho, 2),
# whole-plane mean intensity of every channel against depth (every 4th pixel), same planes
"plane_rho_by_channel": {k[-1]: round(stats.spearmanr(zs, [float(np.asarray(acq[k]["data"][z, ::4, ::4, 0]).mean())
for z in zs]).statistic, 2) for k in sorted(acq)},
"planes_over_90pct_pvs": int((pv.mean((1, 2)) > pvs.FILLED_PLANE).sum()),
# depth of the PVS against the vessel, in columns holding both, after clearing filled planes
"pvs_minus_vessel_depth_planes": round(float(np.nanmedian(depth_offset(ves, pvs.clear_filled_planes(pv)))), 1),
"emission_nm": [float(f[f"general/optophysiology/{p}/optical_channel/emission_lambda"][()])
for p in f["general/optophysiology"]],
})
info = pd.DataFrame(rows).sort_values("location").set_index("location")
print(info.to_string())
print("\nDescription: 14 datasets from 9 subjects; locations per subject:")
print(info.groupby("subject").apply(lambda d: "".join(d.index), include_groups=False).to_string())
print("\nChanA description:", f["acquisition/TwoPhotonSeriesChanA"].attrs["description"])
print("Mask descriptions:", {k: f[f"processing/ophys/ImageSegmentation/{k}/image_mask"].attrs["description"]
for k in f["processing/ophys/ImageSegmentation"]})
for a in assets:
f = pvs.open_asset(a)
if pvs.voxel_size(f)[2] != 1.0:
print(f"Session description ({pvs.location_of(f)}, stated voxel depth {pvs.voxel_size(f)[2]} um):",
f["session_description"][()].decode())
subject strain sex age planes shape_yx voxel_um channels image_dtype masks vessel_vox pvs_vox overlap_vox pvs_frac_volume vessel_z pvs_z deepest_plane_pvs_frac intensity_vs_depth_rho plane_rho_by_channel planes_over_90pct_pvs pvs_minus_vessel_depth_planes emission_nm
location
A BPN-OLD-M2 BPN/3J M P603D 111 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 2743873 4080587 770 0.140 (16, 100) (0, 99) 0.00 -0.62 {'A': -0.74, 'B': -0.98, 'C': -0.98} 0 -33.5 [698.0, 607.0, 525.0]
B BPN-OLD-M2 BPN/3J M P573D 131 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 2597689 4852230 0 0.141 (29, 112) (0, 116) 0.00 -0.93 {'A': -0.57, 'B': -0.93, 'C': -0.98} 0 -44.4 [698.0, 607.0, 525.0]
C BPN-OLD-M3 BPN/3J M P575D 131 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 813216 583622 24 0.017 (22, 108) (21, 92) 0.00 -0.19 {'A': -0.76, 'B': -0.93, 'C': -0.9} 0 -14.5 [698.0, 607.0, 525.0]
D BPN-OLD-M3 BPN/3J M P575D 131 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 1026398 1817420 1654 0.053 (39, 91) (28, 92) 0.00 0.90 {'A': -0.33, 'B': -0.48, 'C': -0.55} 0 -17.0 [698.0, 607.0, 525.0]
E BPN-M4 BPN/3J M P97D 181 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 3593270 28381243 344565 0.598 (15, 101) (32, 180) 0.99 -0.43 {'A': -0.83, 'B': 0.19, 'C': 0.19} 87 10.2 [698.0, 607.0, 525.0]
F BPN-M4 BPN/3J M P97D 181 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 3534910 23848232 18716 0.503 (28, 116) (5, 180) 0.98 -0.81 {'A': -0.67, 'B': -0.79, 'C': -0.93} 58 -16.5 [698.0, 607.0, 525.0]
G BPN-M6 BPN/3J M P99D 151 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 767570 7239025 0 0.183 (40, 95) (0, 108) 0.03 -0.62 {'A': -0.93, 'B': -0.88, 'C': -0.74} 0 -22.6 [698.0, 607.0, 525.0]
H BPN-M6 BPN/3J M P99D 201 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 4020768 6252765 25657 0.119 (38, 179) (0, 172) 0.00 -0.48 {'A': -1.0, 'B': -0.69, 'C': -0.6} 0 -51.2 [698.0, 607.0, 525.0]
I BPN-M7 BPN/3J M P100D 151 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 2229002 2185133 1871 0.055 (28, 139) (29, 104) 0.00 -0.60 {'A': -0.62, 'B': -0.93, 'C': -0.83} 0 -27.4 [698.0, 607.0, 525.0]
J BPN-M7 BPN/3J M P100D 201 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 1604412 3632352 177417 0.069 (51, 130) (0, 200) 0.00 0.02 {'A': -0.26, 'B': -0.76, 'C': -0.69} 0 -23.0 [698.0, 607.0, 525.0]
K 21-07-19_b-act β-act-EGFP M P91D 234 (512, 512) (0.648, 0.648, 1.0) A,B,C float64 A,B,C 2111778 4653193 15087 0.076 (36, 112) (0, 127) 0.00 -0.12 {'A': -0.98, 'B': -1.0, 'C': -1.0} 0 -49.0 [698.0, 607.0, 525.0]
L 21-10-08_b-act β-act-EGFP F P113D 309 (512, 512) (0.648, 0.648, 0.648) A float64 A,C 3914381 4233377 2242 0.052 (50, 211) (30, 200) 0.00 -0.93 {'A': -0.83} 0 -48.5 [698.0]
M 21-09-01_b-act β-act-EGFP F P76D 232 (512, 512) (0.648, 0.648, 0.648) A float64 A,C 1644982 2036897 0 0.033 (1, 152) (1, 167) 0.00 -0.24 {'A': 0.24} 0 -40.5 [698.0]
N 21-09-20_b-act β-act-EGFP F P95D 599 (512, 512) (0.648, 0.648, 0.648) A float64 A,C 6872151 6793881 12102 0.043 (36, 253) (12, 280) 0.00 -0.95 {'A': -0.79} 0 -72.5 [698.0]
Description: 14 datasets from 9 subjects; locations per subject:
subject
21-07-19_b-act K
21-09-01_b-act M
21-09-20_b-act N
21-10-08_b-act L
BPN-M4 EF
BPN-M6 GH
BPN-M7 IJ
BPN-OLD-M2 AB
BPN-OLD-M3 CD
ChanA description: Two-photon series for FITC Dextran (tracer in blood vessel).
Mask descriptions: {'PlaneSegmentationChanA': 'Segmentation masks for blood vessel.', 'PlaneSegmentationChanB': 'Segmentation masks for ChanB.', 'PlaneSegmentationChanC': 'Segmentation masks for perivascular space.'}
Session description (M, stated voxel depth 0.648 um): Zstack (series of 2D microscopy images acquired with focal planes at various cortical depths separated by 1 micron) in the pial artery from a bAct-GFP mouse under ketamine/xylazine anesthesia and 3 min after atipamezole injection
Session description (N, stated voxel depth 0.648 um): Zstack (series of 2D microscopy images acquired with focal planes at various cortical depths separated by 1 micron) in the pial artery from a bAct-GFP mouse under ketamine/xylazine anesthesia and 31 min after atipamezole injection
Session description (L, stated voxel depth 0.648 um): Zstack (series of 2D microscopy images acquired with focal planes at various cortical depths separated by 1 micron) in the pial artery from a bAct-GFP mouse under ketamine/xylazine anesthesia
Checks of the metadata against the values:
- Nine subjects and 14 files, one per location, as the description maps them. Every file has a vessel mask (ChanA) and a PVS mask (ChanC). Image channels are float64 in "arbitrary unit".
- Depth direction. The mean vessel-channel intensity inside the vessel mask falls with plane index in most files (negative rank correlation above), as the Methods expect from attenuation with depth, so plane 0 is the surface side. This is consistent with, not proof of, the orientation.
- Location D, the one file where the intensity inside the vessel rises with plane index. The mean
intensity over the whole plane falls with plane index in all its channels (
plane_rho_by_channel), as in the other files, and its PVS lies above the vessel (negativepvs_minus_vessel_depth_planes), as in every file but E and F. So plane 0 is taken as the surface side in D too; the rise inside its vessel mask does not by itself show a reversed depth axis. - The PVS mask fills whole planes in E and F (
planes_over_90pct_pvs: every plane below a depth); no plane of any other file is more than half PVS. These planes are cleared before the cross sections are taken (pvs.clear_filled_planes, an assumption stated in the attempt record). - Overlap: the Methods say overlaps between the vessel and PVS segmentations were removed by hand;
the deposited masks still overlap in most files (
overlap_vox). - Emission wavelengths: the file gives 698 nm for ChanA, described as FITC-dextran, and 525 nm for ChanC, the PVS tracer; the Methods say the green (525 nm) channel holds the FITC-dextran.
In [3]:
# Voxel size in depth. L, M and N state 0.648 um per plane, the others 1 um, while every session
# description says the planes are 1 um apart. A pial artery is close to round in cross section,
# so the height-to-width ratio of the vessel in the normal cross sections tests the stated depth.
d = DATA["a"]
d = d[d.keep_b & (d.ves_area > 0)]
rnd = d.assign(hw=d.ves_h / d.ves_w).groupby("location").hw.median()
print("Median vessel height / width in normal cross sections, with each file's stated voxel depth:")
print(pd.DataFrame({"height_over_width": rnd.round(2), "stated_dz_um": info.voxel_um.map(lambda v: v[2])}).to_string())
print("\nDigests of the masks (duplicate recordings would share one):")
dig = META["a"].digest
print("distinct:", dig.nunique(), "of", len(dig))
# Versions: what each published version holds
for v in ("0.230427.1809", pvs.VERSION):
held = {}
for a in pvs.list_assets(v):
f = pvs.open_asset(a)
held[pvs.location_of(f)] = ",".join(k[-1] for k in f["processing/ophys/ImageSegmentation"])
print(f"version {v}: masks per location", dict(sorted(held.items())))
Median vessel height / width in normal cross sections, with each file's stated voxel depth:
height_over_width stated_dz_um
location
A 1.10 1.000
B 0.93 1.000
C 0.78 1.000
D 0.76 1.000
E 1.16 1.000
F 0.92 1.000
G 0.73 1.000
H 1.34 1.000
I 1.24 1.000
J 0.82 1.000
K 1.49 1.000
L 0.91 0.648
M 1.08 0.648
N 1.13 0.648
Digests of the masks (duplicate recordings would share one):
distinct: 14 of 14
version 0.230427.1809: masks per location {'A': 'A,B,C', 'B': 'A,B,C', 'C': 'A,B,C', 'D': 'A,B,C', 'E': 'A,B,C', 'F': 'A,B,C', 'G': 'A,B,C', 'H': 'A,B,C', 'I': 'A,B,C', 'J': 'A,B,C', 'K': 'A,B,C', 'L': 'A', 'M': 'A', 'N': 'A'}
version 0.230602.1307: masks per location {'A': 'A,B,C', 'B': 'A,B,C', 'C': 'A,B,C', 'D': 'A,B,C', 'E': 'A,B,C', 'F': 'A,B,C', 'G': 'A,B,C', 'H': 'A,B,C', 'I': 'A,B,C', 'J': 'A,B,C', 'K': 'A,B,C', 'L': 'A,C', 'M': 'A,C', 'N': 'A,C'}
If the depth were mis-stated by the ratio 1 / 0.648, the vessel in those files would look 1.54 times taller or shorter than wide. The ratios above say whether the stated voxel sizes are consistent with round vessels. The earlier published version lacks the PVS masks of L, M and N, so it cannot support the analysis; the pinned version holds all 14.
In [4]:
# Cross sections per location and what the exclusion of P4 keeps
P4 = ("a", "b")
P6 = pvs.EROSIONS
P7 = ("a", "b")
summ = []
for p1 in TARGET:
g = DATA[p1].groupby("location")
summ.append(pd.DataFrame({"P1": p1, "points": g.size(), "kept_P4a": g.keep_a.sum(), "kept_P4b": g.keep_b.sum(),
"med_width_a_um": META[p1].med_a.round(1), "med_width_b_um": META[p1].med_b.round(1),
"no_vessel_fit": g.R1_e3.apply(lambda s: int(s.isna().sum())),
"no_poly_fit": g.HC_e3.apply(lambda s: int(s.isna().sum()))}))
print(pd.concat(summ).to_string())
def sel(p1, p4):
d = DATA[p1]
return d[d["keep_" + p4]]
def ratio_seg(d, p5, k=3):
"""Segmentation area ratio: P5 reading a, segmented vessel area; b, the fitted circle."""
return d.pvs_area / (d.ves_area if p5 == "a" else np.pi * d[f"R1_e{k}"] ** 2)
def within(x, value, tol):
return bool(np.isfinite(x) and abs(x - value) <= tol + 1e-12)
d = sel("a", "b")
share = d.overlap_area / (d.pvs_area + d.overlap_area)
print("Pixels in both masks, as a share of the PVS mask in each slice (P1 a, P4 b): median",
round(share.median(), 4), "; slices above 5%:", int((share > 0.05).sum()), "of", len(d), "; by location above 5%:",
share.gt(0.05).groupby(d.location).sum().to_dict())
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
ax[0].hist(d.pvs_area, bins=60); ax[0].set_xlabel("PVS area (um^2)")
ax[1].hist(d.ves_area, bins=60); ax[1].set_xlabel("vessel area (um^2)")
ax[2].hist(d.R1_e3.dropna(), bins=60); ax[2].set_xlabel("R1 (um), erosion 3")
plt.tight_layout(); plt.show()
print("Ranges (P1 a, P4 b): PVS area", d.pvs_area.min().round(1), "to", d.pvs_area.max().round(1),
"; vessel area", d.ves_area.min().round(1), "to", d.ves_area.max().round(1),
"; R1", d.R1_e3.min().round(1), "to", d.R1_e3.max().round(1), "um")
P1 points kept_P4a kept_P4b med_width_a_um med_width_b_um no_vessel_fit no_poly_fit
location
A a 508 0 460 233.9 86.5 16 28
B a 711 0 430 226.2 108.8 8 235
C a 514 491 514 85.5 26.1 114 189
D a 544 267 492 146.4 58.0 12 41
E a 361 22 202 182.1 95.6 76 164
F a 619 0 472 261.8 130.3 20 116
G a 445 0 409 225.5 84.1 1 19
H a 527 309 508 219.7 99.5 1 240
I a 552 453 537 140.0 39.1 25 64
J a 494 458 485 174.3 51.4 4 110
K a 516 477 516 153.6 44.0 0 33
L a 540 449 522 154.9 49.2 0 67
M a 571 341 532 103.7 38.0 10 139
N a 632 435 604 195.0 59.7 10 120
A b 375 0 325 233.0 86.6 10 53
B b 450 0 283 226.0 107.7 5 145
C b 351 320 346 85.0 26.4 84 140
D b 361 176 324 146.0 57.8 13 30
E b 240 1 124 202.5 97.7 48 94
F b 424 0 305 258.5 123.8 2 87
G b 271 0 264 227.0 82.4 0 0
H b 339 232 319 185.0 95.9 32 193
I b 372 292 357 136.0 39.4 13 46
J b 338 292 331 173.0 48.7 14 69
K b 333 302 325 154.0 44.2 0 20
L b 351 289 328 153.0 49.4 10 54
M b 362 215 334 104.0 38.5 7 91
N b 407 276 394 198.0 60.3 5 77
Pixels in both masks, as a share of the PVS mask in each slice (P1 a, P4 b): median 0.0 ; slices above 5%: 381 of 6683 ; by location above 5%: {'A': 0, 'B': 0, 'C': 0, 'D': 0, 'E': 131, 'F': 0, 'G': 0, 'H': 0, 'I': 0, 'J': 247, 'K': 3, 'L': 0, 'M': 0, 'N': 0}
Ranges (P1 a, P4 b): PVS area 0.0 to 14744.5 ; vessel area 0.4 to 10686.1 ; R1 2.4 to 109.9 um
Centerline and PVS content of the cross sections¶
The verifier of round 1 found the longest path through the skeleton longer than any straight line
through the volume, and the PVS of a slice taken over the whole 332 um plane. The path is now
cleaned of erroneous points (pvs.clean_path: points farther from the polynomial fit than the
vessel radius, and points that revisit a stretch already passed), the PVS of a slice is the two
largest PVS regions within one vessel radius of the vessel (pvs.pvs_lobes), and the edge
distance of P4 is measured from the vessel center found in the slice, to which the paper shifts
the slice. The tables below show the path before and after cleaning, the PVS content of the kept
cross sections, and one cross section per location with the fits drawn over it.
In [5]:
# Path length against the straight distance between its ends, before and after cleaning
pl = pd.concat({p1: META[p1][["path_points", "path_length_um", "path_chord_um", "kept_points", "kept_chord_um",
"tube_radius_um", "gaps", "spacing", "spacing_all"]] for p1 in TARGET}, axis=1)
print("Centerline per location (lengths in um; spacing: mean step between neighbouring kept points, spacing_all with gaps):")
print(pl.round(2).to_string())
for p1 in TARGET:
m = META[p1]
print(f"P1 {p1}: path length over end-to-end distance {(m.path_length_um / m.path_chord_um).min():.2f} to "
f"{(m.path_length_um / m.path_chord_um).max():.2f}; points kept {m.kept_points.sum()} of {m.path_points.sum()}")
ex = META["a"].attrs["examples"]
fig, axs = plt.subplots(2, 7, figsize=(21, 6.5))
for ax, loc in zip(axs.flat, sorted(ex)):
e = ex[loc]
ax.imshow(e["mip_vessel"].astype(float) + 0.4 * e["mip_pvs"], cmap="gray", vmin=0, vmax=1.4,
extent=(0, 512 * pvs.XY_UM, 512 * pvs.XY_UM, 0))
ax.plot(*e["path"][~e["keep"], :2].T, "r.", ms=1)
ax.plot(*e["path"][e["keep"], :2].T, "c.", ms=1)
ax.set_title(loc, fontsize=9); ax.set_xticks([]); ax.set_yticks([])
fig.suptitle("Vessel (white) and PVS (grey) projected over depth; path points kept (cyan) and removed (red), P1 a")
plt.tight_layout(); plt.show()
# How much the cleaned path still meanders: its length (kept points times the mean step) against
# the straight distance between its kept ends and against the arc length of the second-order
# polynomial fit of the centerline over the kept points (P3), along which the slices are oriented
MEANDER = {}
for p1 in TARGET:
m = META[p1]
ex1 = m.attrs["examples"]
arc = {}
for loc in m.index:
_, s = pvs._curve(ex1[loc]["path"], ex1[loc]["keep"])
arc[loc] = float(np.ptp(s[ex1[loc]["keep"]]))
MEANDER[p1] = pd.DataFrame({"kept_points": m.kept_points, "kept_length_um": (m.kept_points * m.spacing).round(1),
"kept_chord_um": m.kept_chord_um.round(1), "fit_arc_um": pd.Series(arc).round(1)})
MEANDER[p1]["length_over_chord"] = (MEANDER[p1].kept_length_um / MEANDER[p1].kept_chord_um).round(2)
MEANDER[p1]["length_over_fit_arc"] = (MEANDER[p1].kept_length_um / MEANDER[p1].fit_arc_um).round(2)
print("\nMeander of the cleaned path (kept length = kept points x mean step):")
print(pd.concat(MEANDER, axis=1).to_string())
for p1 in TARGET:
print(f"P1 {p1}: kept length over end-to-end distance {MEANDER[p1].length_over_chord.min():.2f} to "
f"{MEANDER[p1].length_over_chord.max():.2f}; over the fitted arc {MEANDER[p1].length_over_fit_arc.min():.2f} to "
f"{MEANDER[p1].length_over_fit_arc.max():.2f}")
# PVS content of the kept cross sections (P4 b)
rows = []
for p1 in TARGET:
d = sel(p1, "b")
g = d.groupby("location")
rows.append(pd.DataFrame({"P1": p1, "n": g.size(), "regions_in_slice_median": g.pvs_regions.median(),
"pvs_area_lobes_median": g.pvs_area.median().round(0), "pvs_area_all_median": g.pvs_area_all.median().round(0),
"share_lobes_within_3R1": g.pvs_near_share.median().round(3),
"share_all_within_3R1": g.pvs_all_near_share.median().round(3),
"hext_lobes_um": g.pvs_hext.median().round(1), "hext_all_um": g.pvs_hext_all.median().round(1),
"edge_shift_um": (d.edge_dist - d.edge_dist_skel).abs().groupby(d.location).median().round(2)}))
print("\nPVS content per location (medians over the cross sections kept under P4 b):")
print(pd.concat(rows).to_string())
for p1 in TARGET:
d = DATA[p1]
print(f"P1 {p1}: kept under P4 b with the edge distance from the slice's vessel center {int(d.keep_b.sum())}; "
f"from the skeleton point {int((d.edge_dist_skel >= d.location.map(META[p1].med_b) / 2.5).sum())}")
fig, axs = plt.subplots(2, 7, figsize=(21, 6.5))
for ax, loc in zip(axs.flat, sorted(ex)):
e = ex[loc]
fm = pvs.fit_masks(e["v"], e["p"])
img = np.zeros(e["v"].shape + (3,))
img[e["p"] & ~e["v"]] = (0.35, 0.35, 0.35) # all PVS pixels
img[fm["pvs"]] = (0.2, 0.6, 1.0) # the lobes kept
img[fm["vessel"]] = (1.0, 0.3, 0.3)
ax.imshow(img)
for name, col in (("ellipse", "y"), ("spline", "m"), ("polynomial", "lime")):
if name in fm:
ax.contour(fm[name], levels=[0.5], colors=col, linewidths=0.7)
if "circle" in fm:
x0, y0, r = fm["circle"]
ax.add_patch(plt.Circle((x0, y0), r, fill=False, color="w", lw=0.7))
ax.set_title(f"{loc}, slice {e['i']}", fontsize=9); ax.set_xticks([]); ax.set_yticks([])
fig.suptitle("One kept cross section per location (P1 a): vessel red, PVS lobes blue, other PVS grey; circle white, "
"ellipse yellow, spline magenta, polynomial green")
plt.tight_layout(); plt.show()
Centerline per location (lengths in um; spacing: mean step between neighbouring kept points, spacing_all with gaps):
a b
path_points path_length_um path_chord_um kept_points kept_chord_um tube_radius_um gaps spacing spacing_all path_points path_length_um path_chord_um kept_points kept_chord_um tube_radius_um gaps spacing spacing_all
location
K 695 570.55 370.23 516 330.76 22.05 11 0.83 0.87 432 543.52 357.85 333 331.48 22.07 4 1.27 1.29
M 584 538.60 422.12 571 422.12 16.11 5 0.93 0.94 364 518.47 426.38 362 426.38 16.21 1 1.43 1.44
N 770 685.61 383.12 632 401.41 27.34 20 0.89 0.96 476 651.16 389.97 407 404.27 28.12 8 1.38 1.44
L 600 496.10 352.47 540 355.12 27.70 9 0.83 0.86 353 447.89 357.70 351 357.70 27.51 1 1.27 1.28
F 717 642.82 394.29 619 394.29 29.11 19 0.91 0.96 431 611.23 439.10 424 439.10 29.02 5 1.42 1.44
E 939 813.68 363.27 361 343.54 17.24 19 0.86 1.27 592 798.99 358.87 240 344.58 17.34 11 1.32 1.95
G 667 578.05 359.10 445 343.50 12.55 8 0.87 1.09 424 566.19 400.61 271 323.33 13.01 6 1.32 1.54
H 716 624.48 283.85 527 290.79 25.95 17 0.88 1.00 460 625.17 360.21 339 326.42 25.83 10 1.38 1.50
I 644 573.42 319.40 552 319.14 24.03 20 0.89 0.92 407 553.28 325.54 372 325.32 23.96 12 1.36 1.44
J 554 494.36 323.68 494 323.96 21.09 7 0.89 0.96 352 482.17 339.85 338 339.85 21.07 4 1.37 1.41
B 740 658.56 356.51 711 353.89 27.90 21 0.89 0.91 468 633.72 356.07 450 348.68 28.00 9 1.35 1.38
A 617 542.94 328.35 508 328.35 31.02 9 0.87 0.97 397 536.41 370.40 375 370.40 31.04 11 1.35 1.41
D 714 625.63 460.09 544 380.84 15.01 14 0.88 0.90 461 613.25 454.72 361 384.37 15.07 1 1.34 1.34
C 607 510.92 326.90 514 318.52 13.25 6 0.83 0.85 393 499.84 325.05 351 324.17 13.24 6 1.27 1.29
P1 a: path length over end-to-end distance 1.28 to 2.24; points kept 7534 of 9564
P1 b: path length over end-to-end distance 1.22 to 2.23; points kept 4974 of 6010
Meander of the cleaned path (kept length = kept points x mean step):
a b
kept_points kept_length_um kept_chord_um fit_arc_um length_over_chord length_over_fit_arc kept_points kept_length_um kept_chord_um fit_arc_um length_over_chord length_over_fit_arc
K 516 428.3 330.8 328.4 1.29 1.30 333 422.4 331.5 330.8 1.27 1.28
M 571 529.0 422.1 422.2 1.25 1.25 362 517.6 426.4 426.2 1.21 1.21
N 632 564.1 401.4 398.9 1.41 1.41 407 560.1 404.3 400.2 1.39 1.40
L 540 447.8 355.1 353.7 1.26 1.27 351 446.7 357.7 365.2 1.25 1.22
F 619 560.2 394.3 415.1 1.42 1.35 424 603.5 439.1 448.3 1.37 1.35
E 361 311.8 343.5 344.2 0.91 0.91 240 316.8 344.6 343.1 0.92 0.92
G 445 387.9 343.5 360.2 1.13 1.08 271 357.5 323.3 337.8 1.11 1.06
H 527 463.5 290.8 333.1 1.59 1.39 339 467.6 326.4 334.0 1.43 1.40
I 552 492.8 319.1 342.4 1.54 1.44 372 505.9 325.3 349.2 1.56 1.45
J 494 441.5 324.0 322.8 1.36 1.37 338 464.7 339.8 340.0 1.37 1.37
B 711 630.8 353.9 475.8 1.78 1.33 450 608.8 348.7 466.4 1.75 1.31
A 508 443.9 328.3 329.2 1.35 1.35 375 506.4 370.4 371.3 1.37 1.36
D 544 480.1 380.8 384.9 1.26 1.25 361 483.5 384.4 387.9 1.26 1.25
C 514 428.0 318.5 329.6 1.34 1.30 351 444.0 324.2 338.8 1.37 1.31
P1 a: kept length over end-to-end distance 0.91 to 1.78; over the fitted arc 0.91 to 1.44
P1 b: kept length over end-to-end distance 0.92 to 1.75; over the fitted arc 0.92 to 1.45
PVS content per location (medians over the cross sections kept under P4 b):
P1 n regions_in_slice_median pvs_area_lobes_median pvs_area_all_median share_lobes_within_3R1 share_all_within_3R1 hext_lobes_um hext_all_um edge_shift_um
location
A a 460 6.0 4088.0 4455.0 0.848 0.833 237.8 266.7 1.53
B a 430 7.0 4918.0 5419.0 0.644 0.620 241.1 274.1 2.37
C a 514 2.0 692.0 716.0 1.000 0.978 85.5 88.8 1.70
D a 492 4.0 1959.0 2020.0 0.702 0.696 149.7 163.9 1.75
E a 202 17.0 3067.0 5852.0 0.696 0.674 208.7 295.8 1.28
F a 472 40.0 9154.0 9413.0 0.772 0.768 279.3 279.3 2.32
G a 409 53.0 7574.0 7717.0 0.526 0.521 227.4 243.6 1.79
H a 508 81.0 7442.0 8159.0 0.997 0.873 224.9 241.4 2.99
I a 537 22.0 1961.0 1996.0 0.990 0.982 141.3 156.2 3.31
J a 485 55.0 4317.0 4359.0 0.801 0.791 175.0 187.9 1.54
K a 516 4.0 5499.0 5619.0 0.992 0.951 153.6 206.7 3.06
L a 522 3.0 2986.0 3046.0 1.000 1.000 157.1 162.0 1.27
M a 532 8.0 1188.0 1248.0 0.925 0.901 105.0 112.8 0.75
N a 604 29.0 5413.0 5567.0 0.996 0.991 195.7 208.0 3.81
A b 325 5.0 3939.0 4216.0 0.830 0.829 239.0 265.0 1.46
B b 283 6.0 4879.0 5384.0 0.637 0.611 241.0 274.0 2.78
C b 346 2.0 686.0 702.0 0.993 0.975 85.5 88.0 1.65
D b 324 3.0 1948.0 2011.0 0.670 0.664 150.0 164.0 1.63
E b 124 15.0 1595.0 5465.0 0.723 0.687 205.0 297.5 1.43
F b 305 29.0 9289.0 9548.0 0.775 0.772 281.0 281.0 2.38
G b 264 43.5 7186.0 7280.0 0.480 0.476 227.0 295.0 1.62
H b 319 60.0 6247.0 8894.0 0.997 0.885 187.0 266.0 2.48
I b 357 16.0 2074.0 2123.0 0.987 0.979 140.0 153.0 3.43
J b 331 39.0 4143.0 4251.0 0.796 0.785 173.0 188.0 1.57
K b 325 4.0 5489.0 5621.0 0.986 0.951 154.0 206.0 2.89
L b 328 2.0 3026.0 3090.0 1.000 1.000 154.0 160.0 1.26
M b 334 6.0 1208.0 1264.0 0.916 0.890 105.0 113.0 1.09
N b 394 28.0 5364.0 5550.0 0.997 0.992 198.5 209.0 3.85
P1 a: kept under P4 b with the edge distance from the slice's vessel center 6683; from the skeleton point 6665
P1 b: kept under P4 b with the edge distance from the slice's vessel center 4359; from the skeleton point 4369
Before cleaning, the longest path is 1.28 to 2.24 times the straight distance between its ends at 0.648 um (1.22 to 2.23 at 1 um), as the verifier found. Cleaning keeps 7534 of 9564 points (4974 of 6010 at 1 um), and the mean step between neighbouring kept points is 0.83 to 0.93 um at 0.648 um. The cleaned path still meanders: at 0.648 um its length (kept points times the mean step) is 1.13 to 1.78 times the distance between its kept ends and 1.08 to 1.44 times the arc of the fitted centerline at every location but E (0.91 there, where the gaps between kept points make that product shorter than the path). A 26-connected path along a straight line is at most about 8% longer than the line, so most of the excess is the path meandering inside the vessel, not the zigzag of the voxels. The vessel mask of E holds a network of branching vessels, which the Methods say were removed by hand; its cleaned path keeps 361 of 939 points along a side branch (right of the projection), while the points removed run along the vessel with PVS beside it (left). E's cross sections are therefore of the branch. Its median R1 is the largest of the 14 locations (40.9 um, under vessel-radius-peak); the paper gives no per-location radii to compare with. By median area, the PVS lobes kept are 91 to 99% of the slice's PVS pixels at every location but E (52% at 0.648 um, 29% at 1 um; the rest is the band left above the cleared planes and the PVS of other vessels) and, at 1 um, H (70%). The lobes are wide: by median, 48 to 72% of their area lies within 3 R1 of the vessel center at B, D, E and G, and 77 to 100% elsewhere. At F the lobes still join a band of PVS along the bottom of the slice (the example above). Measuring the edge distance from the vessel center moves it by 0.75 to 3.85 um by median and keeps 6683 cross sections in place of 6665 (P1 a, P4 b). The effect of E on each claim is shown before the summary.
cross-section-counts¶
After excluding cross sections near the edge of the imaging volume, each of the 14 PVS segments yields between 54 and 451 normal cross sections, about 0.8 um apart.
Every combination of P1 (voxel 0.648 or 1 um), P2 (longest path or all skeleton points) and P4 (median width from the horizontal extent of the PVS, or from the polynomial fit). Holds when the fewest and most per segment lie in 48.6 to 496.1 and the mean spacing is 0.72 to 0.88 um.
In [6]:
rows = []
for p1 in TARGET:
for p4 in P4:
g = DATA[p1].groupby("location")["keep_" + p4].sum()
sp = META[p1].spacing.mean()
rows.append({"P1": p1, "P2": "a", "P4": p4, "fewest": int(g.min()), "most": int(g.max()),
"segments_with_0": int((g == 0).sum()), "mean_spacing_um": round(sp, 3)})
for p4 in P4:
n = []
for r in COUNTS_B[p1]:
med = r["med_hext"] if p4 == "a" else r["med_wpol"]
n.append(int((r["edge_dist"] >= med / 2.5).sum()))
sp = np.mean([r["spacing"] for r in COUNTS_B[p1]])
rows.append({"P1": p1, "P2": "b", "P4": p4, "fewest": min(n), "most": max(n),
"segments_with_0": sum(x == 0 for x in n), "mean_spacing_um": round(sp, 3)})
cc = pd.DataFrame(rows)
cc["counts_in_range"] = (cc.fewest >= 54 - 5.4) & (cc.most <= 451 + 45.1)
cc["spacing_ok"] = (cc.mean_spacing_um - 0.8).abs() <= 0.08
cc["holds"] = cc.counts_in_range & cc.spacing_ok
print(cc.to_string(index=False))
print("\nPer location, P2 a (path order spacing in um):")
print(pd.concat({p1: META[p1][["points", "spacing"]].round(3) for p1 in TARGET}, axis=1).to_string())
print("\nP2 b points per location:", {p1: {r["location"]: r["n_points"] for r in COUNTS_B[p1]} for p1 in TARGET})
# Path meander (data check): the counts a path without meander would give, at the same mean step
# along the fitted centerline, with the same share of cross sections kept under P4 b. A diagnostic,
# not a reading of the claims file, which takes one slice per centerline point.
est = []
for p1 in TARGET:
g = DATA[p1].groupby("location").keep_b.sum()
n_fit = (MEANDER[p1].fit_arc_um / META[p1].spacing).reindex(g.index)
kept = (g / META[p1].points * n_fit).round().astype(int)
est.append({"P1": p1, "P4": "b", "step_um": "path mean", "fewest": int(kept.min()), "most": int(kept.max()),
"mean_spacing_um": round(META[p1].spacing.mean(), 3),
"counts_in_range": bool(kept.min() >= 54 - 5.4 and kept.max() <= 451 + 45.1), "per_location": kept.to_dict()})
# the same at the paper's 0.8 um and at the bounds of the spacing tolerance
for st in (0.72, 0.8, 0.88):
k_st = (g / META[p1].points * MEANDER[p1].fit_arc_um.reindex(g.index) / st).round().astype(int)
est.append({"P1": p1, "P4": "b", "step_um": st, "fewest": int(k_st.min()), "most": int(k_st.max()), "mean_spacing_um": st,
"counts_in_range": bool(k_st.min() >= 54 - 5.4 and k_st.max() <= 451 + 45.1), "per_location": ""})
print("\nDiagnostic, counts along the fitted centerline at a given step (P4 b):")
print(pd.DataFrame(est).to_string(index=False))
P1 P2 P4 fewest most segments_with_0 mean_spacing_um counts_in_range spacing_ok holds
a a a 0 491 4 0.876 False True False
a a b 202 604 0 0.876 False True False
a b a 0 30121 5 25.358 False False False
a b b 1122 78442 0 25.358 False False False
b a a 0 320 4 1.345 False False False
b a b 124 394 0 1.345 True False False
b b a 0 19532 5 24.496 False False False
b b b 439 36049 0 24.496 False False False
Per location, P2 a (path order spacing in um):
a b
points spacing points spacing
location
K 516 0.830 333 1.268
M 571 0.926 362 1.430
N 632 0.893 407 1.376
L 540 0.829 351 1.273
F 619 0.905 424 1.423
E 361 0.864 240 1.320
G 445 0.872 271 1.319
H 527 0.880 339 1.379
I 552 0.893 372 1.360
J 494 0.894 338 1.375
B 711 0.887 450 1.353
A 508 0.874 375 1.350
D 544 0.882 361 1.339
C 514 0.833 351 1.265
P2 b points per location: {'a': {'K': 3653, 'M': 4288, 'N': 37285, 'L': 1170, 'F': 49864, 'E': 61131, 'G': 14817, 'H': 91207, 'I': 48886, 'J': 38368, 'B': 32835, 'A': 31666, 'D': 12642, 'C': 10288}, 'b': {'K': 1163, 'M': 745, 'N': 8359, 'L': 455, 'F': 20007, 'E': 25103, 'G': 6377, 'H': 42010, 'I': 22053, 'J': 16662, 'B': 14453, 'A': 12833, 'D': 5134, 'C': 4615}}
Diagnostic, counts along the fitted centerline at a given step (P4 b):
P1 P4 step_um fewest most mean_spacing_um counts_in_range per_location
a b path mean 223 427 0.876 True {'A': 341, 'B': 324, 'C': 396, 'D': 394, 'E': 223, 'F': 350, 'G': 380, 'H': 365, 'I': 373, 'J': 355, 'K': 396, 'L': 412, 'M': 425, 'N': 427}
a b 0.72 267 546 0.720 False
a b 0.8 241 492 0.800 True
a b 0.88 219 447 0.880 True
b b path mean 134 282 1.345 True {'A': 238, 'B': 217, 'C': 264, 'D': 260, 'E': 134, 'F': 227, 'G': 249, 'H': 228, 'I': 246, 'J': 242, 'K': 255, 'L': 268, 'M': 275, 'N': 282}
b b 0.72 246 546 0.720 False
b b 0.8 222 492 0.800 True
b b 0.88 201 447 0.880 True
Result: not reproduced. No combination of P1, P2 and P4 holds. With the cleaned path (P2 a) and the polynomial width (P4 b), the segments give 202 to 604 cross sections at 0.648 um, with a mean spacing of 0.876 um (within tolerance) but the most above 496.1, and 124 to 394 at 1 um, inside the range widened by 10%, but with a mean spacing of 1.345 um. With the horizontal extent of the PVS (P4 a), 4 segments keep no cross section at either voxel size. All skeleton points (P2 b) give up to 78712 per segment at 0.648 um and 36106 at 1 um.
The cleaned path still meanders (data check), so it holds more points than a path along the vessel would. As a diagnostic, not a reading of the claims file, counted along the fitted centerline at each segment's own mean step and with the share kept under P4 (b), the segments would give 223 to 427 cross sections at 0.648 um, with the mean spacing of 0.876 um, which meets the criterion (134 to 282 at 1 um, where the spacing stays 1.345 um). The count depends on the step a path along the vessel would have: at the paper's 0.8 um it is 241 to 492, inside 48.6 to 496.1; at 0.88 um 219 to 447; at 0.72 um 267 to 546, above it. So at 0.648 um the claim would hold with a path along the vessel at a step of 0.8 um or more, and the outcome turns on how the path is cleaned, which is an assumption. The attempt record proposed a change to reproduced. The third verification did not agree: these counts are scaled from the kept share, not taken from slices along a path that follows the vessel, that path is a reading outside the claims file, and the step of 0.876 um is the mean step of the meandering path, not one measured along the vessel (an average 26-connected step of about 0.73 um would give a most of about 539). The outcome stays not reproduced.
fig5a-pvs-area-scales-with-vessel-area¶
PVS area increases with vessel area only approximately; the least-squares line of PVS area on vessel area over all cross sections has slope 1.05 (95% CI 0.973 to 1.12) and R2 of 0.158.
Ordinary least squares on the segmented areas of every retained cross section, (a) with an intercept, (b) through the origin. Slope within 0.077 of 1.05 and R2 within 0.0158 of 0.158. R2 through the origin is computed against the mean (centred), so that it is comparable with (a); the uncentred value is also shown.
In [7]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
x, y = d.ves_area.values, d.pvs_area.values
r = stats.linregress(x, y)
tq = stats.t.ppf(0.975, len(x) - 2)
b0 = (x @ y) / (x @ x)
res = y - b0 * x
se0 = np.sqrt(res @ res / (len(x) - 1) / (x @ x))
r2c = 1 - (res @ res) / ((y - y.mean()) @ (y - y.mean()))
r2u = 1 - (res @ res) / (y @ y)
for rd, slope, lo, hi, r2 in (("a", r.slope, r.slope - tq * r.stderr, r.slope + tq * r.stderr, r.rvalue ** 2),
("b", b0, b0 - tq * se0, b0 + tq * se0, r2c)):
rows.append({"P1": p1, "P4": p4, "reading": rd, "n": len(x), "slope": round(slope, 3),
"ci_low": round(lo, 3), "ci_high": round(hi, 3), "R2": round(r2, 3),
"R2_uncentred": round(r2u, 3) if rd == "b" else np.nan,
"intercept_um2": round(r.intercept, 1) if rd == "a" else 0.0,
"match": within(slope, 1.05, 0.077) and within(r2, 0.158, 0.0158)})
f5a = pd.DataFrame(rows)
print(f5a.to_string(index=False))
d = sel("a", "b")
print("\nPer location (P1 a, P4 b), slope with intercept:")
print(d.groupby("location").apply(lambda g: pd.Series({"n": len(g), "slope": stats.linregress(g.ves_area, g.pvs_area).slope,
"R2": stats.linregress(g.ves_area, g.pvs_area).rvalue ** 2}), include_groups=False).round(3).to_string())
fig, ax = plt.subplots(figsize=(5, 4))
for loc, g in d.groupby("location"):
ax.scatter(g.ves_area, g.pvs_area, s=2, label=loc)
ax.set_xlabel("vessel area (um^2)"); ax.set_ylabel("PVS area (um^2)"); ax.legend(markerscale=4, fontsize=6, ncol=2)
plt.show()
print("Plotted ranges: vessel", d.ves_area.min().round(), "to", d.ves_area.max().round(), "; PVS", d.pvs_area.min().round(), "to", d.pvs_area.max().round(), "um^2; no points outside the axes")
P1 P4 reading n slope ci_low ci_high R2 R2_uncentred intercept_um2 match
a a a 3702 1.252 1.203 1.302 0.398 NaN 603.1 False
a a b 3702 1.442 1.418 1.466 0.386 0.791 0.0 False
a b a 6683 0.912 0.865 0.959 0.176 NaN 2013.3 False
a b b 6683 1.526 1.501 1.550 0.072 0.691 0.0 False
b a a 2395 1.543 1.484 1.601 0.527 NaN -25.3 False
b a b 2395 1.535 1.508 1.562 0.527 0.840 0.0 False
b b a 4359 0.951 0.894 1.007 0.199 NaN 1842.7 False
b b b 4359 1.510 1.480 1.540 0.108 0.697 0.0 False
Per location (P1 a, P4 b), slope with intercept:
n slope R2
location
A 460.0 -1.456 0.232
B 430.0 0.555 0.066
C 514.0 1.177 0.178
D 492.0 2.462 0.705
E 202.0 0.265 0.046
F 472.0 0.076 0.001
G 409.0 -6.911 0.289
H 508.0 1.917 0.214
I 537.0 1.623 0.658
J 485.0 1.201 0.145
K 516.0 0.896 0.074
L 522.0 -0.693 0.093
M 532.0 -1.172 0.101
N 604.0 -1.433 0.260
Plotted ranges: vessel 0.0 to 10686.0 ; PVS 0.0 to 14745.0 um^2; no points outside the axes
Result: not reproduced. Over the eight combinations the slope is 0.912 to 1.543, against 1.05 within 0.077, and R2 is 0.072 to 0.527, against 0.158 within 0.0158. The nearest are P1 (a), P4 (b) with an intercept, slope 0.912 (95% CI 0.865 to 0.959) and R2 0.176, and P1 (b), P4 (b), slope 0.951 (0.894 to 1.007) and R2 0.199; neither puts both values within tolerance. The slopes within single segments run from -6.911 (G) to 2.462 (D) (P1 a, P4 b).
fig5e-area-ratio-distribution¶
The PVS-to-vessel area ratio of the segmentation varies widely across cross sections; its distribution peaks at 1.12, its median is 1.88 and its IQR is 1.79 (25th percentile 1.21, 75th 3.00), so the ratio varies by more than 100%.
Every combination of P1, P4, P5 (segmented vessel area, or the fitted circle for each erosion of P6) and P7 (bins over all values, or over the values without outliers).
In [8]:
EXP_5E = [("peak", 1.12, 0.112), ("median", 1.88, 0.188), ("IQR", 1.79, 0.179), ("q25", 1.21, 0.121), ("q75", 3.0, 0.3)]
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for p5, k in [("a", 3)] + [("b", k) for k in P6]:
r = ratio_seg(d, p5, k).replace([np.inf], np.nan).dropna()
q25, q50, q75 = np.percentile(r, [25, 50, 75])
for p7 in P7:
v = {"peak": pvs.peak(r, p7), "median": q50, "IQR": q75 - q25, "q25": q25, "q75": q75}
ok = {n: within(v[n], e, t) for n, e, t in EXP_5E}
rows.append({"P1": p1, "P4": p4, "P5": p5, "P6": k if p5 == "b" else "-", "P7": p7, "n": len(r),
**{n: round(v[n], 3) for n in v}, "all_within": all(ok.values()),
"median_IQR_within": ok["median"] and ok["IQR"]})
f5e = pd.DataFrame(rows)
print(f5e.to_string(index=False))
print("\nRatio range (P1 a, P4 b, P5 a):", ratio_seg(sel("a", "b"), "a").min().round(3), "to", ratio_seg(sel("a", "b"), "a").max().round(2))
print("Per location median ratio (P1 a, P4 b, P5 a):", ratio_seg(sel("a", "b"), "a").groupby(sel("a", "b").location).median().round(2).to_dict())
P1 P4 P5 P6 P7 n peak median IQR q25 q75 all_within median_IQR_within
a a a - a 3702 5.273 1.402 1.040 0.943 1.983 False False
a a a - b 3702 1.324 1.402 1.040 0.943 1.983 False False
a a b 1 a 3590 3.845 1.226 0.874 0.875 1.749 False False
a a b 1 b 3590 1.152 1.226 0.874 0.875 1.749 False False
a a b 3 a 3590 3.848 1.277 0.924 0.940 1.864 False False
a a b 3 b 3590 1.227 1.277 0.924 0.940 1.864 False False
a a b 5 a 3590 3.849 1.317 0.942 0.966 1.908 False False
a a b 5 b 3590 1.256 1.317 0.942 0.966 1.908 False False
a b a - a 6683 7.142 1.565 1.287 0.995 2.282 False False
a b a - b 6683 1.568 1.565 1.287 0.995 2.282 False False
a b b 1 a 6511 3.845 1.348 1.093 0.908 2.001 False False
a b b 1 b 6511 1.194 1.348 1.093 0.908 2.001 False False
a b b 3 a 6511 3.848 1.426 1.148 0.967 2.115 False False
a b b 3 b 6511 1.257 1.426 1.148 0.967 2.115 False False
a b b 5 a 6511 3.849 1.472 1.209 0.989 2.197 False False
a b b 5 b 6511 1.285 1.472 1.209 0.989 2.197 False False
b a a - a 2395 1.260 1.400 1.009 0.971 1.979 False False
b a a - b 2395 1.260 1.400 1.009 0.971 1.979 False False
b a b 1 a 2328 1.133 1.253 0.892 0.881 1.774 False False
b a b 1 b 2328 1.133 1.253 0.892 0.881 1.774 False False
b a b 3 a 2328 1.230 1.334 0.940 0.960 1.900 False False
b a b 3 b 2328 1.230 1.334 0.940 0.960 1.900 False False
b a b 5 a 2328 1.432 1.376 0.968 0.989 1.956 False False
b a b 5 b 2328 1.432 1.376 0.968 0.989 1.956 False False
b b a - a 4359 1.257 1.533 1.214 0.999 2.213 False False
b b a - b 4359 1.510 1.533 1.214 0.999 2.213 False False
b b b 1 a 4213 1.251 1.377 1.064 0.923 1.988 False False
b b b 1 b 4213 1.178 1.377 1.064 0.923 1.988 False False
b b b 3 a 4213 1.342 1.461 1.173 0.993 2.166 False False
b b b 3 b 4213 1.287 1.461 1.173 0.993 2.166 False False
b b b 5 a 4213 1.342 1.510 1.249 1.016 2.265 False False
b b b 5 b 4213 1.346 1.510 1.249 1.016 2.265 False False
Ratio range (P1 a, P4 b, P5 a): 0.0 to 285.68
Per location median ratio (P1 a, P4 b, P5 a): {'A': 1.31, 'B': 1.64, 'C': 0.71, 'D': 2.18, 'E': 1.21, 'F': 2.66, 'G': 8.79, 'H': 1.58, 'I': 0.77, 'J': 2.1, 'K': 2.04, 'L': 1.05, 'M': 1.21, 'N': 1.29}
Result: not reproduced. No combination puts the median and IQR within tolerance, which the partial rule needs: the median is 1.226 to 1.565 (expected 1.88 within 0.188) and the IQR 0.874 to 1.287 (expected 1.79 within 0.179); the 25th percentile is 0.875 to 1.016 and the 75th 1.749 to 2.282 (expected 1.21 and 3.00). The peak is within tolerance in eight combinations, all with P5 (b), at 1.133 to 1.230 (for example 1.152 with P1 a, P4 a, erosion 1 and outliers out). With P1 (a), P4 (b) and P5 (a) the median is 1.565 and the IQR 1.287, so the ratio varies by about 80% of its median.
fig5bc-area-intersegment-variation¶
Variation in PVS area is much larger between PVS segments than along the length of a single segment.
Per segment, the median and IQR of segmented PVS area. Reading (a): the IQR of the 14 segment medians is larger than the median of the 14 within-segment IQRs; (b): the range of the 14 medians is larger than the largest within-segment IQR.
In [9]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
g = d.groupby("location").pvs_area
med = g.median(); iqr = g.quantile(0.75) - g.quantile(0.25)
med, iqr = med.dropna(), iqr.dropna()
a = (med.quantile(0.75) - med.quantile(0.25)) > iqr.median()
b = (med.max() - med.min()) > iqr.max()
rows.append({"P1": p1, "P4": p4, "segments": len(med), "IQR_of_medians": round(med.quantile(0.75) - med.quantile(0.25), 1),
"median_within_IQR": round(iqr.median(), 1), "range_of_medians": round(med.max() - med.min(), 1),
"largest_within_IQR": round(iqr.max(), 1), "reading_a": a, "reading_b": b})
if p1 == "a":
print(f"P1 {p1}, P4 {p4}: per segment median and IQR of PVS area (um^2)")
print(pd.DataFrame({"n": g.size(), "median": med.round(1), "IQR": iqr.round(1)}).T.to_string())
f5bc = pd.DataFrame(rows)
print(f5bc.to_string(index=False))
P1 a, P4 a: per segment median and IQR of PVS area (um^2) location C D E H I J K L M N n 491.0 267.0 22.0 309.0 453.0 458.0 477.0 449.0 341.0 435.0 median 706.3 2661.4 5444.5 9593.5 1803.9 4444.7 5546.5 3214.8 1219.8 5542.3 IQR 1214.8 929.5 4424.6 7661.6 1337.0 1285.5 1433.6 1228.2 500.1 433.8
P1 a, P4 b: per segment median and IQR of PVS area (um^2) location A B C D E F G H I J K L M N n 460.0 430.0 514.0 492.0 202.0 472.0 409.0 508.0 537.0 485.0 516.0 522.0 532.0 604.0 median 4087.6 4917.9 691.8 1959.1 3067.4 9153.7 7574.2 7442.4 1961.4 4317.0 5499.5 2986.4 1188.3 5413.4 IQR 1957.1 1423.5 1182.0 1388.0 4748.7 3768.5 2403.1 9476.5 1618.3 1872.4 1443.1 1211.0 675.8 1624.6 P1 P4 segments IQR_of_medians median_within_IQR range_of_medians largest_within_IQR reading_a reading_b a a 10 3499.6 1256.9 8887.3 7661.6 True True a b 14 3260.3 1621.5 8461.9 9476.5 True False b a 10 3382.8 1242.4 6178.5 6886.5 True False b b 14 3478.0 1708.5 8603.5 10464.5 True False
Result: blocked (methods underspecified). Reading (a) holds in all four combinations of P1 and P4: with P1 (a) and P4 (b) the IQR of the 14 segment medians is 3260.3 um^2 against a median within-segment IQR of 1621.5. Reading (b) holds only with P1 (a) and P4 (a), where 4 of the 14 segments keep no cross section (range of 10 medians 8887.3 against 7661.6 at H). With all 14 segments (P4 b) the range of the medians, 8461.9 at P1 (a) and 8603.5 at P1 (b), is below the largest within-segment IQR, 9476.5 (H) and 10464.5. The readings split, which the criterion sends to methods_underspecified.
area-ratio-per-location¶
The peak of the segmented area ratio distribution ranges from 0.6 to 11.43 across the 14 locations, and in all but one location the ratio's variation along the vessel lies within 30% of the peak, so variation between segments dominates.
Per location, P5 and P7 as above. "Within 30% of the peak": (a) the 25th and 75th percentiles both lie within 30% of the peak; (b) IQR at most 30% of the peak; (c) IQR at most 60% of the peak; (d) all non-outlier values within 30% of the peak. Expected: smallest peak 0.54 to 0.66, largest 10.29 to 12.57, count of locations 11.7 to 14.3.
In [10]:
def per_location_peaks(d, values, p7, frac, readings):
out = []
for loc, v in values.groupby(d.location):
v = v.replace([np.inf], np.nan).dropna()
if len(v) < 5:
continue
pk = pvs.peak(v, p7)
q25, q75 = np.percentile(v, [25, 75])
nof = pvs.outlier_free(v)
r = {"location": loc, "n": len(v), "peak": pk, "q25": q25, "q75": q75,
"a": (abs(q25 - pk) <= frac * pk) and (abs(q75 - pk) <= frac * pk),
"b": (q75 - q25) <= frac * pk, "c": (q75 - q25) <= 2 * frac * pk,
"d": bool(np.all(np.abs(nof - pk) <= frac * pk))}
out.append({k: r[k] for k in ["location", "n", "peak", "q25", "q75"] + list(readings)})
return pd.DataFrame(out)
rows, shown = [], None
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for p5, k in [("a", 3)] + [("b", k) for k in P6]:
for p7 in P7:
t = per_location_peaks(d, ratio_seg(d, p5, k), p7, 0.30, "abcd")
if (p1, p4, p5, p7) == ("a", "b", "a", "b"):
shown = t
for rd in "abcd":
cnt = int(t[rd].sum())
rows.append({"P1": p1, "P4": p4, "P5": p5, "P6": k if p5 == "b" else "-", "P7": p7, "reading": rd,
"locations": len(t), "min_peak": round(t.peak.min(), 3), "max_peak": round(t.peak.max(), 2),
"count_within": cnt,
"match": within(t.peak.min(), 0.6, 0.06) and within(t.peak.max(), 11.43, 1.143) and within(cnt, 13, 1.3)})
arl = pd.DataFrame(rows)
print("Per location, P1 a, P4 b, P5 a, P7 b:")
print(shown.round(3).to_string(index=False))
print(arl.to_string(index=False))
print("\nCombinations matching:", int(arl.match.sum()), "of", len(arl))
print("Ranges over all combinations: smallest peak", arl.min_peak.min(), "to", arl.min_peak.max(), "; largest peak", arl.max_peak.min(),
"to", arl.max_peak.max(), "; locations within, at most", arl.count_within.max())
print("\nThe 8 combinations with the most locations within (compact):")
print(arl.sort_values("count_within", ascending=False).head(8).to_string(index=False))
Per location, P1 a, P4 b, P5 a, P7 b:
location n peak q25 q75 a b c d
A 460 1.015 1.019 1.764 False False False False
B 430 1.508 1.467 2.001 False False True False
C 514 0.053 0.075 1.040 False False False False
D 492 2.373 1.802 2.438 True True True False
E 202 1.067 0.449 17.556 False False False False
F 472 2.583 2.261 3.616 False False True False
G 409 9.045 7.143 10.348 True False True False
H 508 0.230 0.388 2.017 False False False False
I 537 0.717 0.644 1.086 False False False False
J 485 2.017 1.689 2.718 False False True False
K 516 1.841 1.834 2.605 False False True False
L 522 0.741 0.775 1.407 False False False False
M 532 0.970 0.916 1.580 False False False False
N 604 1.302 0.897 1.511 False False True False
P1 P4 P5 P6 P7 reading locations min_peak max_peak count_within match
a a a - a a 10 0.053 6.41 2 False
a a a - a b 10 0.053 6.41 3 False
a a a - a c 10 0.053 6.41 6 False
a a a - a d 10 0.053 6.41 0 False
a a a - b a 10 0.053 2.75 1 False
a a a - b b 10 0.053 2.75 2 False
a a a - b c 10 0.053 2.75 5 False
a a a - b d 10 0.053 2.75 0 False
a a b 1 a a 10 0.141 4.76 3 False
a a b 1 a b 10 0.141 4.76 2 False
a a b 1 a c 10 0.141 4.76 6 False
a a b 1 a d 10 0.141 4.76 0 False
a a b 1 b a 10 0.141 1.70 3 False
a a b 1 b b 10 0.141 1.70 1 False
a a b 1 b c 10 0.141 1.70 5 False
a a b 1 b d 10 0.141 1.70 0 False
a a b 3 a a 10 0.172 4.82 3 False
a a b 3 a b 10 0.172 4.82 2 False
a a b 3 a c 10 0.172 4.82 5 False
a a b 3 a d 10 0.172 4.82 0 False
a a b 3 b a 10 0.172 2.29 4 False
a a b 3 b b 10 0.172 2.29 1 False
a a b 3 b c 10 0.172 2.29 5 False
a a b 3 b d 10 0.172 2.29 0 False
a a b 5 a a 10 0.186 4.88 2 False
a a b 5 a b 10 0.186 4.88 2 False
a a b 5 a c 10 0.186 4.88 4 False
a a b 5 a d 10 0.186 4.88 0 False
a a b 5 b a 10 0.186 2.51 2 False
a a b 5 b b 10 0.186 2.51 1 False
a a b 5 b c 10 0.186 2.51 4 False
a a b 5 b d 10 0.186 2.51 0 False
a b a - a a 14 0.053 8.83 2 False
a b a - a b 14 0.053 8.83 3 False
a b a - a c 14 0.053 8.83 7 False
a b a - a d 14 0.053 8.83 0 False
a b a - b a 14 0.053 9.04 2 False
a b a - b b 14 0.053 9.04 1 False
a b a - b c 14 0.053 9.04 7 False
a b a - b d 14 0.053 9.04 0 False
a b b 1 a a 14 0.161 7.20 4 False
a b b 1 a b 14 0.161 7.20 3 False
a b b 1 a c 14 0.161 7.20 7 False
a b b 1 a d 14 0.161 7.20 0 False
a b b 1 b a 14 0.158 6.91 1 False
a b b 1 b b 14 0.158 6.91 2 False
a b b 1 b c 14 0.158 6.91 7 False
a b b 1 b d 14 0.158 6.91 0 False
a b b 3 a a 14 0.172 7.42 3 False
a b b 3 a b 14 0.172 7.42 2 False
a b b 3 a c 14 0.172 7.42 8 False
a b b 3 a d 14 0.172 7.42 0 False
a b b 3 b a 14 0.172 7.42 3 False
a b b 3 b b 14 0.172 7.42 1 False
a b b 3 b c 14 0.172 7.42 8 False
a b b 3 b d 14 0.172 7.42 0 False
a b b 5 a a 14 0.186 7.98 4 False
a b b 5 a b 14 0.186 7.98 1 False
a b b 5 a c 14 0.186 7.98 8 False
a b b 5 a d 14 0.186 7.98 0 False
a b b 5 b a 14 0.186 7.83 3 False
a b b 5 b b 14 0.186 7.83 0 False
a b b 5 b c 14 0.186 7.83 9 False
a b b 5 b d 14 0.186 7.83 0 False
b a a - a a 9 0.052 2.10 2 False
b a a - a b 9 0.052 2.10 1 False
b a a - a c 9 0.052 2.10 6 False
b a a - a d 9 0.052 2.10 0 False
b a a - b a 9 0.052 2.10 1 False
b a a - b b 9 0.052 2.10 1 False
b a a - b c 9 0.052 2.10 6 False
b a a - b d 9 0.052 2.10 0 False
b a b 1 a a 9 0.653 2.06 3 False
b a b 1 a b 9 0.653 2.06 1 False
b a b 1 a c 9 0.653 2.06 6 False
b a b 1 a d 9 0.653 2.06 0 False
b a b 1 b a 9 0.629 2.06 3 False
b a b 1 b b 9 0.629 2.06 1 False
b a b 1 b c 9 0.629 2.06 6 False
b a b 1 b d 9 0.629 2.06 0 False
b a b 3 a a 9 0.722 1.99 3 False
b a b 3 a b 9 0.722 1.99 1 False
b a b 3 a c 9 0.722 1.99 5 False
b a b 3 a d 9 0.722 1.99 0 False
b a b 3 b a 9 0.703 1.99 3 False
b a b 3 b b 9 0.703 1.99 1 False
b a b 3 b c 9 0.703 1.99 5 False
b a b 3 b d 9 0.703 1.99 0 False
b a b 5 a a 9 0.740 2.73 2 False
b a b 5 a b 9 0.740 2.73 1 False
b a b 5 a c 9 0.740 2.73 6 False
b a b 5 a d 9 0.740 2.73 0 False
b a b 5 b a 9 0.722 2.73 2 False
b a b 5 b b 9 0.722 2.73 1 False
b a b 5 b c 9 0.722 2.73 6 False
b a b 5 b d 9 0.722 2.73 0 False
b b a - a a 14 0.052 8.25 2 False
b b a - a b 14 0.052 8.25 1 False
b b a - a c 14 0.052 8.25 7 False
b b a - a d 14 0.052 8.25 0 False
b b a - b a 14 0.052 8.28 2 False
b b a - b b 14 0.052 8.28 1 False
b b a - b c 14 0.052 8.28 7 False
b b a - b d 14 0.052 8.28 0 False
b b b 1 a a 14 0.190 7.21 2 False
b b b 1 a b 14 0.190 7.21 1 False
b b b 1 a c 14 0.190 7.21 6 False
b b b 1 a d 14 0.190 7.21 0 False
b b b 1 b a 14 0.190 6.93 1 False
b b b 1 b b 14 0.190 6.93 0 False
b b b 1 b c 14 0.190 6.93 7 False
b b b 1 b d 14 0.190 6.93 0 False
b b b 3 a a 14 0.204 7.94 4 False
b b b 3 a b 14 0.204 7.94 1 False
b b b 3 a c 14 0.204 7.94 9 False
b b b 3 a d 14 0.204 7.94 0 False
b b b 3 b a 14 0.204 8.80 3 False
b b b 3 b b 14 0.204 8.80 1 False
b b b 3 b c 14 0.204 8.80 9 False
b b b 3 b d 14 0.204 8.80 0 False
b b b 5 a a 14 0.212 8.71 3 False
b b b 5 a b 14 0.212 8.71 1 False
b b b 5 a c 14 0.212 8.71 8 False
b b b 5 a d 14 0.212 8.71 0 False
b b b 5 b a 14 0.212 7.92 2 False
b b b 5 b b 14 0.212 7.92 1 False
b b b 5 b c 14 0.212 7.92 9 False
b b b 5 b d 14 0.212 7.92 0 False
Combinations matching: 0 of 128
Ranges over all combinations: smallest peak 0.052 to 0.74 ; largest peak 1.7 to 9.04 ; locations within, at most 9
The 8 combinations with the most locations within (compact):
P1 P4 P5 P6 P7 reading locations min_peak max_peak count_within match
a b b 5 b c 14 0.186 7.83 9 False
b b b 3 b c 14 0.204 8.80 9 False
b b b 3 a c 14 0.204 7.94 9 False
b b b 5 b c 14 0.212 7.92 9 False
b b b 5 a c 14 0.212 8.71 8 False
a b b 5 a c 14 0.186 7.98 8 False
a b b 3 b c 14 0.172 7.42 8 False
a b b 3 a c 14 0.172 7.42 8 False
Result: not reproduced. None of the 128 combinations matches. The smallest per-location peak is 0.052 to 0.74 and the largest 1.7 to 9.04 (expected 0.54 to 0.66 and 10.29 to 12.57); with P1 (a), P4 (b), P5 (a) and P7 (b) they are 0.053 (C) and 9.045 (G). At most 9 locations meet a reading of "within 30% of the peak" (reading c, the loosest), against 11.7 to 14.3, and none meets reading (d) in any combination.
fig5d-polynomial-closest-area-ratio¶
Of the three fits, the polynomial fit matches the segmentation's area ratio most closely in 73% of cross sections, the spline fit in 19% and the ellipse fit in 8%, and its distribution matches the segmentation's best.
Per cross section, the fit whose area ratio is closest to the segmentation's. Every combination of P5 (for the segmentation and, separately, for the fits), P6 and P9 (ellipse region a, ellipse minus its overlap with the circle; b, ellipse area minus circle area), under each P1 and P4. Only cross sections where all three fits exist are counted. Shares within 7.3, 1.9 and 0.8 points of 73, 19 and 8%.
In [11]:
def fit_ratios(d, k, p5_fit, p9):
den = np.pi * d[f"R1_e{k}"] ** 2 if p5_fit == "a" else d.ves_area
ell = d[f"ell_area_a_e{k}"] if p9 == "a" else d[f"ell_area_b_e{k}"]
return pd.DataFrame({"polynomial": d[f"pol_area_e{k}"] / den, "spline": d[f"spl_area_e{k}"] / den,
"ellipse": ell / den})
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
for p5s in "ab":
seg = ratio_seg(d, p5s, k)
for p5f in "ab":
for p9 in "ab":
fr = fit_ratios(d, k, p5f, p9)
ok = fr.notna().all(axis=1) & np.isfinite(seg)
best = (fr[ok].sub(seg[ok], axis=0)).abs().idxmin(axis=1)
sh = best.value_counts(normalize=True).reindex(["polynomial", "spline", "ellipse"]).fillna(0) * 100
rows.append({"P1": p1, "P4": p4, "P6": k, "P5_seg": p5s, "P5_fit": p5f, "P9": p9, "n": int(ok.sum()),
"dropped": int((~ok).sum()), **sh.round(1).to_dict(),
"match": within(sh.polynomial, 73, 7.3) and within(sh.spline, 19, 1.9) and within(sh.ellipse, 8, 0.8),
"partial": sh.polynomial == sh.max() and sh.polynomial > 50})
f5d = pd.DataFrame(rows)
print(f5d.to_string(index=False))
print("\nCombinations matching:", int(f5d.match.sum()), "; meeting the partial rule:", int(f5d.partial.sum()), "of", len(f5d))
P1 P4 P6 P5_seg P5_fit P9 n dropped polynomial spline ellipse match partial a a 1 a a a 3306 396 50.0 20.2 29.8 False False a a 1 a a b 3306 396 47.4 20.4 32.2 False False a a 1 a b a 3306 396 78.4 20.0 1.6 False True a a 1 a b b 3306 396 75.5 19.8 4.7 False True a a 1 b a a 3306 396 78.4 20.0 1.6 False True a a 1 b a b 3306 396 75.5 19.8 4.7 False True a a 1 b b a 3306 396 63.8 36.2 0.1 False True a a 1 b b b 3306 396 60.0 36.2 3.8 False True a a 3 a a a 3309 393 66.9 19.7 13.4 False True a a 3 a a b 3309 393 65.1 20.7 14.3 False True a a 3 a b a 3309 393 77.7 22.0 0.4 False True a a 3 a b b 3309 393 76.8 21.5 1.7 False True a a 3 b a a 3309 393 77.6 22.0 0.4 False True a a 3 b a b 3309 393 76.8 21.5 1.7 False True a a 3 b b a 3309 393 71.7 28.3 0.1 False True a a 3 b b b 3309 393 70.6 28.3 1.1 False True a a 5 a a a 3317 385 74.9 20.3 4.8 False True a a 5 a a b 3317 385 74.0 20.4 5.5 False True a a 5 a b a 3317 385 77.0 22.8 0.2 False True a a 5 a b b 3317 385 76.9 22.8 0.3 False True a a 5 b a a 3317 385 77.0 22.8 0.2 False True a a 5 b a b 3317 385 76.9 22.8 0.3 False True a a 5 b b a 3317 385 74.6 25.3 0.1 False True a a 5 b b b 3317 385 74.4 25.3 0.3 False True a b 1 a a a 5778 905 49.3 21.1 29.6 False False a b 1 a a b 5778 905 47.8 21.4 30.8 False False a b 1 a b a 5778 905 77.8 20.1 2.1 False True a b 1 a b b 5778 905 76.0 19.9 4.2 False True a b 1 b a a 5778 905 77.8 20.1 2.1 False True a b 1 b a b 5778 905 76.0 19.9 4.2 False True a b 1 b b a 5778 905 65.0 34.9 0.1 False True a b 1 b b b 5778 905 61.7 34.7 3.6 False True a b 3 a a a 5786 897 63.2 20.2 16.6 False True a b 3 a a b 5786 897 62.5 20.8 16.7 False True a b 3 a b a 5786 897 76.7 22.0 1.3 False True a b 3 a b b 5786 897 76.4 21.8 1.8 False True a b 3 b a a 5786 897 76.7 22.0 1.3 False True a b 3 b a b 5786 897 76.4 21.8 1.8 False True a b 3 b b a 5786 897 70.1 29.8 0.1 False True a b 3 b b b 5786 897 68.8 29.7 1.5 False True a b 5 a a a 5796 887 70.7 20.6 8.7 True True a b 5 a a b 5796 887 70.4 20.8 8.8 True True a b 5 a b a 5796 887 76.8 21.8 1.4 False True a b 5 a b b 5796 887 77.1 21.8 1.1 False True a b 5 b a a 5796 887 76.8 21.8 1.4 False True a b 5 b a b 5796 887 77.1 21.8 1.1 False True a b 5 b b a 5796 887 73.6 26.2 0.2 False True a b 5 b b b 5796 887 72.9 26.1 1.0 False True b a 1 a a a 2135 260 55.8 18.3 25.9 False True b a 1 a a b 2135 260 53.2 19.3 27.5 False True b a 1 a b a 2135 260 79.1 20.4 0.5 False True b a 1 a b b 2135 260 76.5 20.2 3.2 False True b a 1 b a a 2135 260 79.1 20.4 0.5 False True b a 1 b a b 2135 260 76.5 20.2 3.2 False True b a 1 b b a 2135 260 67.0 32.6 0.4 False True b a 1 b b b 2135 260 64.2 32.6 3.2 False True b a 3 a a a 2135 260 75.4 20.3 4.3 False True b a 3 a a b 2135 260 74.3 20.3 5.4 False True b a 3 a b a 2135 260 78.1 21.2 0.7 False True b a 3 a b b 2135 260 77.8 21.2 1.0 False True b a 3 b a a 2135 260 78.1 21.2 0.7 False True b a 3 b a b 2135 260 77.8 21.2 1.0 False True b a 3 b b a 2135 260 75.7 24.1 0.1 False True b a 3 b b b 2135 260 75.5 24.0 0.5 False True b a 5 a a a 2133 262 77.0 21.3 1.7 False True b a 5 a a b 2133 262 76.8 21.0 2.2 False True b a 5 a b a 2133 262 77.4 21.8 0.9 False True b a 5 a b b 2133 262 77.2 21.8 1.1 False True b a 5 b a a 2133 262 77.4 21.8 0.9 False True b a 5 b a b 2133 262 77.2 21.8 1.1 False True b a 5 b b a 2133 262 76.3 23.6 0.1 False True b a 5 b b b 2133 262 76.1 23.6 0.2 False True b b 1 a a a 3716 643 54.0 19.8 26.2 False True b b 1 a a b 3716 643 52.8 20.5 26.8 False True b b 1 a b a 3716 643 78.6 20.1 1.3 False True b b 1 a b b 3716 643 77.3 20.1 2.6 False True b b 1 b a a 3716 643 78.6 20.1 1.3 False True b b 1 b a b 3716 643 77.3 20.1 2.6 False True b b 1 b b a 3716 643 67.2 32.5 0.3 False True b b 1 b b b 3716 643 65.1 32.3 2.6 False True b b 3 a a a 3729 630 71.4 20.8 7.9 True True b b 3 a a b 3729 630 71.0 20.8 8.2 True True b b 3 a b a 3729 630 78.1 20.6 1.2 False True b b 3 a b b 3729 630 78.3 20.7 1.0 False True b b 3 b a a 3729 630 78.1 20.6 1.2 False True b b 3 b a b 3729 630 78.3 20.7 1.0 False True b b 3 b b a 3729 630 75.4 24.0 0.6 False True b b 3 b b b 3729 630 74.7 23.9 1.3 False True b b 5 a a a 3732 627 74.7 21.6 3.8 False True b b 5 a a b 3732 627 75.0 21.6 3.4 False True b b 5 a b a 3732 627 78.0 20.7 1.3 False True b b 5 a b b 3732 627 78.0 20.6 1.4 False True b b 5 b a a 3732 627 78.0 20.7 1.3 False True b b 5 b a b 3732 627 78.0 20.6 1.4 False True b b 5 b b a 3732 627 76.4 22.2 1.4 False True b b 5 b b b 3732 627 75.6 22.2 2.2 False True Combinations matching: 4 ; meeting the partial rule: 92 of 96
Result: reproduced. Four of the 96 combinations put the three shares within tolerance, all with P4 (b), the segmented vessel area for the segmentation's ratio and the fitted circle for the fits: P1 (a) with erosion 5 (polynomial 70.7%, spline 20.6%, ellipse 8.8% with P9 a; 70.4, 20.8 and 8.8% with P9 b) and P1 (b) with erosion 3 (71.4, 20.8 and 7.9%; 71.0, 20.8 and 8.3%). The polynomial share is the largest and above 50% in 92 of the 96 combinations, and runs from 47.4 to 79.1% over all. The outcome rests on those four; under the other readings of P5 the ellipse is closest in only 0.1 to 4.7% of cross sections.
table2-fit-area-ratio-peaks¶
The area ratio distribution peaks at 1.21 for the polynomial fit (close to the segmentation's 1.12), 2.08 for the ellipse fit, which does not narrow at the ends, and 1.48 for the spline fit, whose ratio varies little because it scales with vessel size alone.
Peaks, medians and the spline IQR of the fits' area ratios under each reading of P5 for the fits, P6, P7, P9 (ellipse region), P1 and P4.
In [12]:
EXP_T2 = [("pk_pol", 1.21, 0.121), ("pk_ell", 2.08, 0.208), ("pk_spl", 1.48, 0.148), ("md_pol", 2.15, 0.215),
("md_ell", 2.89, 0.289), ("md_spl", 1.58, 0.158), ("iqr_spl", 0.3, 0.03)]
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
for p5f in "ab":
for p9 in "ab":
fr = fit_ratios(d, k, p5f, p9).replace([np.inf, -np.inf], np.nan)
iqr = {c: np.nanpercentile(fr[c], 75) - np.nanpercentile(fr[c], 25) for c in fr}
for p7 in P7:
v = {"pk_pol": pvs.peak(fr.polynomial.dropna(), p7), "pk_ell": pvs.peak(fr.ellipse.dropna(), p7),
"pk_spl": pvs.peak(fr.spline.dropna(), p7), "md_pol": fr.polynomial.median(),
"md_ell": fr.ellipse.median(), "md_spl": fr.spline.median(), "iqr_spl": iqr["spline"]}
rows.append({"P1": p1, "P4": p4, "P6": k, "P5_fit": p5f, "P9": p9, "P7": p7,
**{n: round(x, 3) for n, x in v.items()},
"iqr_pol": round(iqr["polynomial"], 3), "iqr_ell": round(iqr["ellipse"], 3),
"n_within": sum(within(v[n], e, t) for n, e, t in EXP_T2),
"match": all(within(v[n], e, t) for n, e, t in EXP_T2),
"partial": (v["pk_ell"] > v["pk_spl"] > v["pk_pol"]) and iqr["spline"] < min(iqr["polynomial"], iqr["ellipse"])})
t2 = pd.DataFrame(rows)
print(t2.to_string(index=False))
print("\nCombinations meeting the partial rule:")
print(t2[t2.partial].to_string(index=False))
print("\nCombinations matching:", int(t2.match.sum()), "; meeting the partial rule:", int(t2.partial.sum()), "of", len(t2),
"; most values within tolerance in one combination:", int(t2.n_within.max()), "of 7")
print("\nThe 8 combinations with the most values within tolerance (compact):")
print(t2.sort_values("n_within", ascending=False).head(8)[["P1", "P4", "P6", "P5_fit", "P9", "P7"] + [n for n, _, _ in EXP_T2] + ["n_within"]].to_string(index=False))
print("Ranges over all combinations:", {n: (t2[n].min(), t2[n].max()) for n, _, _ in EXP_T2})
P1 P4 P6 P5_fit P9 P7 pk_pol pk_ell pk_spl md_pol md_ell md_spl iqr_spl iqr_pol iqr_ell n_within match partial
a a 1 a a a 4.369 5.672 1.494 1.467 1.934 1.506 0.022 0.915 1.340 2 False False
a a 1 a a b 1.718 1.957 1.509 1.467 1.934 1.506 0.022 0.915 1.340 3 False False
a a 1 a b a 4.369 4.901 1.494 1.467 1.848 1.506 0.022 0.915 1.361 2 False False
a a 1 a b b 1.718 1.999 1.509 1.467 1.848 1.506 0.022 0.915 1.361 3 False False
a a 1 b a a 6.109 8.570 1.624 1.653 2.192 1.710 0.225 1.011 1.497 2 False False
a a 1 b a b 1.627 2.123 1.568 1.653 2.192 1.710 0.225 1.011 1.497 3 False False
a a 1 b b a 6.109 4.752 1.624 1.653 2.089 1.710 0.225 1.011 1.539 2 False False
a a 1 b b b 1.627 2.252 1.568 1.653 2.089 1.710 0.225 1.011 1.539 3 False False
a a 3 a a a 4.381 5.705 1.495 1.567 2.082 1.508 0.019 0.978 1.409 2 False False
a a 3 a a b 1.847 2.062 1.510 1.567 2.082 1.508 0.019 0.978 1.409 3 False False
a a 3 a b a 4.381 5.041 1.495 1.567 2.001 1.508 0.019 0.978 1.427 2 False False
a a 3 a b b 1.847 2.169 1.510 1.567 2.001 1.508 0.019 0.978 1.427 3 False False
a a 3 b a a 6.089 8.593 1.586 1.664 2.227 1.613 0.148 1.015 1.503 2 False False
a a 3 b a b 0.915 2.150 1.494 1.664 2.227 1.613 0.148 1.015 1.503 3 False True
a a 3 b b a 6.089 6.469 1.586 1.664 2.148 1.613 0.148 1.015 1.528 2 False False
a a 3 b b b 0.915 2.135 1.494 1.664 2.148 1.613 0.148 1.015 1.528 3 False True
a a 5 a a a 4.391 5.743 1.499 1.602 2.148 1.508 0.020 1.009 1.390 2 False False
a a 5 a a b 0.852 2.127 1.511 1.602 2.148 1.508 0.020 1.009 1.390 3 False True
a a 5 a b a 4.391 5.137 1.499 1.602 2.080 1.508 0.020 1.009 1.417 2 False False
a a 5 a b b 0.852 2.003 1.511 1.602 2.080 1.508 0.020 1.009 1.417 3 False True
a a 5 b a a 6.089 8.645 1.599 1.671 2.250 1.571 0.114 1.020 1.505 2 False False
a a 5 b a b 0.922 2.184 1.484 1.671 2.250 1.571 0.114 1.020 1.505 3 False True
a a 5 b b a 6.089 6.910 1.599 1.671 2.178 1.571 0.114 1.020 1.511 2 False False
a a 5 b b b 0.922 2.139 1.484 1.671 2.178 1.571 0.114 1.020 1.511 3 False True
a b 1 a a a 4.369 5.672 1.490 1.608 2.120 1.505 0.031 1.195 1.712 2 False False
a b 1 a a b 1.748 2.113 1.511 1.608 2.120 1.505 0.031 1.195 1.712 3 False False
a b 1 a b a 4.369 4.901 1.490 1.608 2.023 1.505 0.031 1.195 1.737 2 False False
a b 1 a b b 1.748 1.911 1.511 1.608 2.023 1.505 0.031 1.195 1.737 3 False False
a b 1 b a a 12.237 14.500 1.802 1.828 2.385 1.679 0.260 1.375 1.756 1 False False
a b 1 b a b 1.821 2.304 1.543 1.828 2.385 1.679 0.260 1.375 1.756 2 False False
a b 1 b b a 12.237 10.687 1.802 1.828 2.251 1.679 0.260 1.375 1.776 1 False False
a b 1 b b b 1.821 2.005 1.543 1.828 2.251 1.679 0.260 1.375 1.776 3 False False
a b 3 a a a 4.381 5.705 1.531 1.710 2.250 1.506 0.027 1.308 1.768 2 False False
a b 3 a a b 1.680 2.222 1.508 1.710 2.250 1.506 0.027 1.308 1.768 3 False False
a b 3 a b a 4.381 5.041 1.531 1.710 2.158 1.506 0.027 1.308 1.814 2 False False
a b 3 a b b 1.680 2.119 1.508 1.710 2.158 1.506 0.027 1.308 1.814 3 False False
a b 3 b a a 12.239 14.524 1.897 1.845 2.425 1.601 0.157 1.374 1.753 1 False False
a b 3 b a b 1.599 2.328 1.516 1.845 2.425 1.601 0.157 1.374 1.753 2 False False
a b 3 b b a 12.239 12.404 1.897 1.845 2.301 1.601 0.157 1.374 1.767 1 False False
a b 3 b b b 1.599 2.298 1.516 1.845 2.301 1.601 0.157 1.374 1.767 2 False False
a b 5 a a a 4.391 5.743 1.522 1.757 2.313 1.506 0.025 1.360 1.772 2 False False
a b 5 a a b 1.752 2.290 1.509 1.757 2.313 1.506 0.025 1.360 1.772 2 False False
a b 5 a b a 4.391 5.137 1.522 1.757 2.234 1.506 0.025 1.360 1.776 2 False False
a b 5 a b b 1.752 2.193 1.509 1.757 2.234 1.506 0.025 1.360 1.776 3 False False
a b 5 b a a 12.242 14.544 1.837 1.850 2.445 1.559 0.114 1.377 1.747 1 False False
a b 5 b a b 1.837 2.343 1.490 1.850 2.445 1.559 0.114 1.377 1.747 2 False False
a b 5 b b a 12.242 12.846 1.837 1.850 2.327 1.559 0.114 1.377 1.748 1 False False
a b 5 b b b 1.837 2.300 1.490 1.850 2.327 1.559 0.114 1.377 1.748 2 False False
b a 1 a a a 1.634 2.136 1.514 1.478 1.983 1.505 0.036 0.915 1.320 3 False False
b a 1 a a b 1.690 1.987 1.512 1.478 1.983 1.505 0.036 0.915 1.320 3 False False
b a 1 a b a 1.634 2.065 1.514 1.478 1.902 1.505 0.036 0.915 1.370 3 False False
b a 1 a b b 1.690 1.891 1.512 1.478 1.902 1.505 0.036 0.915 1.370 3 False False
b a 1 b a a 1.470 2.131 1.681 1.629 2.212 1.682 0.205 0.983 1.458 2 False True
b a 1 b a b 2.003 2.131 1.690 1.629 2.212 1.682 0.205 0.983 1.458 2 False False
b a 1 b b a 1.470 2.216 1.681 1.629 2.117 1.682 0.205 0.983 1.488 2 False True
b a 1 b b b 2.003 2.256 1.690 1.629 2.117 1.682 0.205 0.983 1.488 2 False False
b a 3 a a a 0.887 2.302 1.520 1.588 2.165 1.506 0.035 0.977 1.419 2 False True
b a 3 a a b 0.840 2.200 1.510 1.588 2.165 1.506 0.035 0.977 1.419 3 False True
b a 3 a b a 0.887 2.247 1.520 1.588 2.099 1.506 0.035 0.977 1.445 3 False True
b a 3 a b b 0.840 2.120 1.510 1.588 2.099 1.506 0.035 0.977 1.445 3 False True
b a 3 b a a 1.474 2.418 1.580 1.647 2.269 1.563 0.121 0.997 1.454 2 False True
b a 3 b a b 1.534 2.418 1.574 1.647 2.269 1.563 0.121 0.997 1.454 2 False True
b a 3 b b a 1.474 2.121 1.580 1.647 2.200 1.563 0.121 0.997 1.473 3 False True
b a 3 b b b 1.534 1.179 1.574 1.647 2.200 1.563 0.121 0.997 1.473 2 False False
b a 5 a a a 1.624 2.302 1.495 1.633 2.240 1.506 0.036 1.016 1.399 2 False False
b a 5 a a b 0.863 2.458 1.505 1.633 2.240 1.506 0.036 1.016 1.399 2 False True
b a 5 a b a 1.624 2.331 1.495 1.633 2.186 1.506 0.036 1.016 1.427 2 False False
b a 5 a b b 0.863 2.215 1.505 1.633 2.186 1.506 0.036 1.016 1.427 3 False True
b a 5 b a a 1.476 2.418 1.489 1.651 2.287 1.515 0.093 0.999 1.455 2 False True
b a 5 b a b 1.545 2.418 1.470 1.651 2.287 1.515 0.093 0.999 1.455 2 False False
b a 5 b b a 1.476 2.361 1.489 1.651 2.235 1.515 0.093 0.999 1.473 2 False True
b a 5 b b b 1.545 2.361 1.470 1.651 2.235 1.515 0.093 0.999 1.473 2 False False
b b 1 a a a 1.341 1.885 1.515 1.610 2.151 1.504 0.039 1.157 1.667 3 False True
b b 1 a a b 1.519 2.097 1.509 1.610 2.151 1.504 0.039 1.157 1.667 3 False False
b b 1 a b a 1.341 1.083 1.515 1.610 2.066 1.504 0.039 1.157 1.691 2 False False
b b 1 a b b 1.519 2.196 1.509 1.610 2.066 1.504 0.039 1.157 1.691 3 False False
b b 1 b a a 1.673 2.251 1.688 1.811 2.381 1.654 0.236 1.290 1.650 2 False True
b b 1 b a b 1.524 2.222 1.532 1.811 2.381 1.654 0.236 1.290 1.650 3 False True
b b 1 b b a 1.673 2.631 1.688 1.811 2.260 1.654 0.236 1.290 1.669 1 False True
b b 1 b b b 1.524 2.223 1.532 1.811 2.260 1.654 0.236 1.290 1.669 3 False True
b b 3 a a a 1.685 2.377 1.516 1.744 2.322 1.505 0.035 1.297 1.742 2 False False
b b 3 a a b 1.691 2.285 1.511 1.744 2.322 1.505 0.035 1.297 1.742 3 False False
b b 3 a b a 1.685 1.741 1.516 1.744 2.245 1.505 0.035 1.297 1.764 2 False False
b b 3 a b b 1.691 2.192 1.511 1.744 2.245 1.505 0.035 1.297 1.764 3 False False
b b 3 b a a 1.683 2.291 1.543 1.823 2.433 1.550 0.121 1.308 1.639 2 False False
b b 3 b a b 1.543 2.259 1.499 1.823 2.433 1.550 0.121 1.308 1.639 3 False False
b b 3 b b a 1.683 0.682 1.543 1.823 2.325 1.550 0.121 1.308 1.664 2 False False
b b 3 b b b 1.543 2.278 1.499 1.823 2.325 1.550 0.121 1.308 1.664 3 False False
b b 5 a a a 1.690 2.383 1.491 1.792 2.410 1.506 0.035 1.367 1.746 2 False False
b b 5 a a b 1.776 2.322 1.511 1.792 2.410 1.506 0.035 1.367 1.746 2 False False
b b 5 a b a 1.690 1.774 1.491 1.792 2.333 1.506 0.035 1.367 1.729 2 False False
b b 5 a b b 1.776 2.227 1.511 1.792 2.333 1.506 0.035 1.367 1.729 3 False False
b b 5 b a a 1.685 2.300 1.504 1.829 2.455 1.508 0.087 1.317 1.640 2 False False
b b 5 b a b 1.554 2.272 1.474 1.829 2.455 1.508 0.087 1.317 1.640 3 False False
b b 5 b b a 1.685 0.926 1.504 1.829 2.358 1.508 0.087 1.317 1.651 2 False False
b b 5 b b b 1.554 2.325 1.474 1.829 2.358 1.508 0.087 1.317 1.651 2 False False
Combinations meeting the partial rule:
P1 P4 P6 P5_fit P9 P7 pk_pol pk_ell pk_spl md_pol md_ell md_spl iqr_spl iqr_pol iqr_ell n_within match partial
a a 3 b a b 0.915 2.150 1.494 1.664 2.227 1.613 0.148 1.015 1.503 3 False True
a a 3 b b b 0.915 2.135 1.494 1.664 2.148 1.613 0.148 1.015 1.528 3 False True
a a 5 a a b 0.852 2.127 1.511 1.602 2.148 1.508 0.020 1.009 1.390 3 False True
a a 5 a b b 0.852 2.003 1.511 1.602 2.080 1.508 0.020 1.009 1.417 3 False True
a a 5 b a b 0.922 2.184 1.484 1.671 2.250 1.571 0.114 1.020 1.505 3 False True
a a 5 b b b 0.922 2.139 1.484 1.671 2.178 1.571 0.114 1.020 1.511 3 False True
b a 1 b a a 1.470 2.131 1.681 1.629 2.212 1.682 0.205 0.983 1.458 2 False True
b a 1 b b a 1.470 2.216 1.681 1.629 2.117 1.682 0.205 0.983 1.488 2 False True
b a 3 a a a 0.887 2.302 1.520 1.588 2.165 1.506 0.035 0.977 1.419 2 False True
b a 3 a a b 0.840 2.200 1.510 1.588 2.165 1.506 0.035 0.977 1.419 3 False True
b a 3 a b a 0.887 2.247 1.520 1.588 2.099 1.506 0.035 0.977 1.445 3 False True
b a 3 a b b 0.840 2.120 1.510 1.588 2.099 1.506 0.035 0.977 1.445 3 False True
b a 3 b a a 1.474 2.418 1.580 1.647 2.269 1.563 0.121 0.997 1.454 2 False True
b a 3 b a b 1.534 2.418 1.574 1.647 2.269 1.563 0.121 0.997 1.454 2 False True
b a 3 b b a 1.474 2.121 1.580 1.647 2.200 1.563 0.121 0.997 1.473 3 False True
b a 5 a a b 0.863 2.458 1.505 1.633 2.240 1.506 0.036 1.016 1.399 2 False True
b a 5 a b b 0.863 2.215 1.505 1.633 2.186 1.506 0.036 1.016 1.427 3 False True
b a 5 b a a 1.476 2.418 1.489 1.651 2.287 1.515 0.093 0.999 1.455 2 False True
b a 5 b b a 1.476 2.361 1.489 1.651 2.235 1.515 0.093 0.999 1.473 2 False True
b b 1 a a a 1.341 1.885 1.515 1.610 2.151 1.504 0.039 1.157 1.667 3 False True
b b 1 b a a 1.673 2.251 1.688 1.811 2.381 1.654 0.236 1.290 1.650 2 False True
b b 1 b a b 1.524 2.222 1.532 1.811 2.381 1.654 0.236 1.290 1.650 3 False True
b b 1 b b a 1.673 2.631 1.688 1.811 2.260 1.654 0.236 1.290 1.669 1 False True
b b 1 b b b 1.524 2.223 1.532 1.811 2.260 1.654 0.236 1.290 1.669 3 False True
Combinations matching: 0 ; meeting the partial rule: 24 of 96 ; most values within tolerance in one combination: 3 of 7
The 8 combinations with the most values within tolerance (compact):
P1 P4 P6 P5_fit P9 P7 pk_pol pk_ell pk_spl md_pol md_ell md_spl iqr_spl n_within
a a 1 a a b 1.718 1.957 1.509 1.467 1.934 1.506 0.022 3
a a 1 a b b 1.718 1.999 1.509 1.467 1.848 1.506 0.022 3
a a 1 b b b 1.627 2.252 1.568 1.653 2.089 1.710 0.225 3
a a 1 b a b 1.627 2.123 1.568 1.653 2.192 1.710 0.225 3
a a 3 b a b 0.915 2.150 1.494 1.664 2.227 1.613 0.148 3
a a 3 b b b 0.915 2.135 1.494 1.664 2.148 1.613 0.148 3
a a 3 a b b 1.847 2.169 1.510 1.567 2.001 1.508 0.019 3
a a 3 a a b 1.847 2.062 1.510 1.567 2.082 1.508 0.019 3
Ranges over all combinations: {'pk_pol': (np.float64(0.84), np.float64(12.242)), 'pk_ell': (np.float64(0.682), np.float64(14.544)), 'pk_spl': (np.float64(1.47), np.float64(1.897)), 'md_pol': (np.float64(1.467), np.float64(1.85)), 'md_ell': (np.float64(1.848), np.float64(2.455)), 'md_spl': (np.float64(1.504), np.float64(1.71)), 'iqr_spl': (np.float64(0.019), np.float64(0.26))}
Result: partially reproduced. No combination puts all seven values within tolerance; the most in one combination is 3 of 7. The order of the peaks, ellipse above spline above polynomial, with the spline IQR the smallest, holds in 24 of 96 combinations, 19 of them with P4 (a). Over all combinations the polynomial peak is 0.842 to 12.237 (expected 1.21), the ellipse peak 0.682 to 14.544 (2.08) and the spline peak 1.47 to 1.897 (1.48); the medians are 1.467 to 1.85 (2.15), 1.848 to 2.455 (2.89) and 1.504 to 1.71 (1.58), and the spline IQR is 0.019 to 0.26 (0.3). The spline's ratio is nearly constant, as the paper says.
fig6a-resistance-power-law¶
The hydraulic resistance per unit length of the segmented PVS is well predicted by a power law in PVS area (R2 = 0.944), with exponent -1.75 (95% CI -1.763 to -1.737), close to the -2 of a shape scaled in size alone.
Readings: (a) least squares of log resistance on log area, R2 on the log scale; (b) nonlinear least squares of resistance on area, R2 on the linear scale. Exponent within 0.013 of -1.75, R2 within 0.0944 of 0.944. The solver check of P8 (a circle of the same area within 1% of the exact resistance) is shown first.
In [13]:
# Solver check (Methods, Hydraulic resistance): circles of the areas met here
chk = []
for p1, px in TARGET.items():
for area_um2 in (300, 1000, 3000, 10000):
rpx = np.sqrt(area_um2 / np.pi) / px
n = int(2 * rpx + 20)
m = pvs.disk((n, n), n / 2, n / 2, rpx)
R, _, _ = pvs.poiseuille(m, px)
a = (np.sqrt(m.sum() / np.pi) - 1.5) * px * 1e-6 # the domain: 1.5 pixels inside the mask
chk.append({"P1": p1, "area_um2": area_um2, "error_vs_exact": round(R / (8 * pvs.MU / (np.pi * a ** 4)) - 1, 4)})
print(pd.DataFrame(chk).to_string(index=False))
d = sel("a", "b")
eq_r = np.sqrt(d.pvs_area / np.pi)
print("Median PVS area", round(d.pvs_area.median()), "um^2; median radius of the circle of the same area", round(eq_r.median(), 1), "um")
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
ok = np.isfinite(d.R_seg) & (d.R_seg > 0) & (d.pvs_area > 0)
A, R = d.pvs_area[ok].values * 1e-12, d.R_seg[ok].values
r = stats.linregress(np.log(A), np.log(R))
tq = stats.t.ppf(0.975, len(A) - 2)
rows.append({"P1": p1, "P4": p4, "reading": "a", "n": len(A), "exponent": round(r.slope, 4),
"ci_low": round(r.slope - tq * r.stderr, 4), "ci_high": round(r.slope + tq * r.stderr, 4),
"R2": round(r.rvalue ** 2, 4), "match": within(r.slope, -1.75, 0.013) and within(r.rvalue ** 2, 0.944, 0.0944),
"partial": within(r.rvalue ** 2, 0.944, 0.0944) and -2 <= r.slope <= -1.5})
pp, cov = curve_fit(lambda a, c, k: c * a ** k, A / A.mean(), R / R.mean(), p0=(1, r.slope), maxfev=20000)
pred = pp[0] * (A / A.mean()) ** pp[1] * R.mean()
r2 = 1 - np.sum((R - pred) ** 2) / np.sum((R - R.mean()) ** 2)
se = np.sqrt(cov[1, 1])
rows.append({"P1": p1, "P4": p4, "reading": "b", "n": len(A), "exponent": round(pp[1], 4),
"ci_low": round(pp[1] - 1.96 * se, 4), "ci_high": round(pp[1] + 1.96 * se, 4), "R2": round(r2, 4),
"match": within(pp[1], -1.75, 0.013) and within(r2, 0.944, 0.0944),
"partial": within(r2, 0.944, 0.0944) and -2 <= pp[1] <= -1.5})
f6a = pd.DataFrame(rows)
print(f6a.to_string(index=False))
d = sel("a", "b")
print("\nPer location log-log exponent (P1 a, P4 b):")
print(d.groupby("location").apply(lambda g: stats.linregress(np.log(g.pvs_area), np.log(g.R_seg)).slope, include_groups=False).round(3).to_dict())
fig, ax = plt.subplots(figsize=(5, 4))
fin = d[np.isfinite(d.R_seg) & (d.pvs_area > 0)]
ax.loglog(fin.pvs_area, fin.R_seg, ".", ms=2); ax.set_xlabel("PVS area (um^2)"); ax.set_ylabel("R (Pa s/m^4)"); plt.show()
print("Plotted ranges: area", fin.pvs_area.min().round(1), "to", fin.pvs_area.max().round(), "um^2; R", f"{fin.R_seg.min():.3g} to {fin.R_seg.max():.3g};",
len(d) - len(fin), "cross sections with no flow domain (infinite resistance) not plotted or fitted")
P1 area_um2 error_vs_exact
a 300 0.0023
a 1000 -0.0078
a 3000 -0.0039
a 10000 -0.0022
b 300 -0.0246
b 1000 -0.0027
b 3000 -0.0083
b 10000 -0.0047
Median PVS area 3901 um^2; median radius of the circle of the same area 35.2 um
P1 P4 reading n exponent ci_low ci_high R2 match partial
a a a 3590 -1.9767 -1.9925 -1.9610 0.9441 False True
a a b 3590 -3.9541 -3.9641 -3.9441 0.9989 False False
a b a 6511 -1.9256 -1.9392 -1.9119 0.9216 False True
a b b 6511 -2.9877 -3.0196 -2.9557 0.9131 False False
b a a 2327 -2.0393 -2.0606 -2.0180 0.9382 False False
b a b 2327 -4.6871 -4.7301 -4.6440 0.9802 False False
b b a 4208 -2.0022 -2.0211 -1.9833 0.9114 False False
b b b 4208 -4.4499 -4.5334 -4.3665 0.8538 False False
Per location log-log exponent (P1 a, P4 b):
{'A': -1.995, 'B': -2.89, 'C': nan, 'D': -1.315, 'E': nan, 'F': -1.712, 'G': -1.968, 'H': -2.051, 'I': nan, 'J': nan, 'K': -1.808, 'L': -1.804, 'M': -2.214, 'N': -2.316}
Plotted ranges: area 16.4 to 14745.0 um^2; R 4.08e+14 to 1.88e+21; 172 cross sections with no flow domain (infinite resistance) not plotted or fitted
Result: partially reproduced. No combination gives the exponent within 0.013 of -1.75. The log-log fit (reading a) gives -1.926 to -2.039 with R2 0.911 to 0.944, and two combinations meet the partial rule (R2 within tolerance, exponent between -2 and -1.5): P1 (a), P4 (a), exponent -1.977 (95% CI -1.993 to -1.961) and R2 0.944, and P1 (a), P4 (b), -1.926 (-1.939 to -1.912) and R2 0.922. The nonlinear fit (reading b) gives -2.99 to -4.69; a few cross sections with very large resistance (up to 1.88e21 Pa s/m^4) dominate it, so it says little about the power law. The solver is within 2.5% of the exact resistance for circles of 300 to 10000 um^2 (within 1% from 1000 um^2). Per-location exponents (P1 a, P4 b) run from -2.89 (B) to -1.315 (D); they are not computed for C, E, I and J, which hold cross sections with no flow domain.
fig6b-resistance-segmentation-peak¶
The distribution of segmented PVS hydraulic resistance per unit length over all cross sections peaks at 1.73e15 Pa s/m^4, with a median of 3.84e15 Pa s/m^4.
P7 readings, with the 20 bins even on a linear (a) or logarithmic (b) scale. Peak within 1.73e14 of 1.73e15, median within 3.84e14 of 3.84e15.
In [14]:
def rpeak(x, p7, scale):
return pvs.peak(np.asarray(x, float), p7, log=(scale == "log"))
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
R = d.R_seg[np.isfinite(d.R_seg)]
for p7 in P7:
for sc in ("linear", "log"):
pk, md_ = rpeak(R, p7, sc), R.median()
rows.append({"P1": p1, "P4": p4, "P7": p7, "bins": sc, "n": len(R), "peak": f"{pk:.3g}", "median": f"{md_:.3g}",
"match": within(pk, 1.73e15, 1.73e14) and within(md_, 3.84e15, 3.84e14),
"median_within": within(md_, 3.84e15, 3.84e14)})
f6b = pd.DataFrame(rows)
print(f6b.to_string(index=False))
d = sel("a", "b")
print("\nPer location median R_seg (P1 a, P4 b):", {k: f"{v:.3g}" for k, v in d.groupby("location").R_seg.median().items()})
P1 P4 P7 bins n peak median match median_within
a a a linear 3590 4.69e+19 5.95e+15 False False
a a a log 3590 1.86e+15 5.95e+15 False False
a a b linear 3590 2.15e+15 5.95e+15 False False
a a b log 3590 2.16e+15 5.95e+15 False False
a b a linear 6511 4.69e+19 5.28e+15 False False
a b a log 6511 2.78e+15 5.28e+15 False False
a b b linear 6511 1.68e+15 5.28e+15 False False
a b b log 6511 2.5e+15 5.28e+15 False False
b a a linear 2327 2.65e+19 6.87e+15 False False
b a a log 2327 2.18e+15 6.87e+15 False False
b a b linear 2327 2.04e+15 6.87e+15 False False
b a b log 2327 2.47e+15 6.87e+15 False False
b b a linear 4208 2.65e+19 6.1e+15 False False
b b a log 4208 2.15e+15 6.1e+15 False False
b b b linear 4208 1.83e+15 6.1e+15 False False
b b b log 4208 2.33e+15 6.1e+15 False False
Per location median R_seg (P1 a, P4 b): {'A': '5.72e+15', 'B': '3.53e+15', 'C': '1.11e+17', 'D': '2.43e+16', 'E': '8.09e+15', 'F': '2.92e+15', 'G': '1.72e+15', 'H': '2.23e+15', 'I': '1.27e+16', 'J': '2.69e+15', 'K': '2.02e+15', 'L': '6.28e+15', 'M': '6.19e+16', 'N': '2.6e+15'}
Result: not reproduced. The peak is within tolerance (1.557e15 to 1.903e15) in three combinations: 1.86e15 Pa s/m^4 with P1 (a), P4 (a), all values and logarithmic bins; 1.68e15 with P1 (a), P4 (b), outliers out and linear bins; and 1.83e15 with P1 (b), P4 (b), outliers out and linear bins. The median is 5.28e15 to 6.87e15 in every combination, against 3.84e15 within 3.84e14. The partial rule needs the median within tolerance, so neither the claim nor its partial rule holds.
resistance-per-location¶
Variation in resistance is much larger between locations than along each PVS; the per-location peaks range from 3.23e14 to 1.98e16, and in 65% of vessels all the variation lies within 37% of the peak.
Per location, P7 readings on linear bins. "All the variation within 37% of the peak": (a) 25th and 75th percentiles within 37% of the peak; (b) all non-outlier values within 37%; (c) IQR at most 37% of the peak. The readings are labelled as in the claims file.
In [15]:
rows, shown = [], None
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for p7 in P7:
t = per_location_peaks(d, d.R_seg, p7, 0.37, "abd")
t = t.rename(columns={"b": "c_iqr", "d": "b_all"}).rename(columns={"c_iqr": "c", "b_all": "b"})
if (p1, p4, p7) == ("a", "b", "b"):
shown = t
for rd in "abc":
share = 100 * t[rd].mean()
rows.append({"P1": p1, "P4": p4, "P7": p7, "reading": rd, "locations": len(t),
"min_peak": f"{t.peak.min():.3g}", "max_peak": f"{t.peak.max():.3g}", "share_pct": round(share, 1),
"match": within(t.peak.min(), 3.23e14, 3.23e13) and within(t.peak.max(), 1.98e16, 1.98e15) and within(share, 65, 6.5)})
rpl = pd.DataFrame(rows)
print("Per location (P1 a, P4 b, P7 b):")
print(shown.assign(peak=shown.peak.map("{:.3g}".format), q25=shown.q25.map("{:.3g}".format), q75=shown.q75.map("{:.3g}".format)).to_string(index=False))
print(rpl.to_string(index=False))
Per location (P1 a, P4 b, P7 b):
location n peak q25 q75 a c b
A 460 5.59e+15 3.41e+15 8.15e+15 False False False
B 430 1.66e+15 1.81e+15 5.94e+15 False False False
C 400 3.78e+16 3.7e+16 2.54e+17 False False False
D 492 1.73e+16 1.79e+16 3.74e+16 False False False
E 160 2.97e+15 1.6e+15 4.71e+16 False False False
F 472 1.52e+15 1.56e+15 3.78e+15 False False False
G 409 1.14e+15 1.14e+15 2.83e+15 False False False
H 508 3.83e+15 1.01e+15 5.61e+16 False False False
I 522 5.21e+15 5.18e+15 2.26e+16 False False False
J 484 2.53e+15 2.17e+15 6.1e+15 False False False
K 516 2.03e+15 1.71e+15 2.69e+15 True False False
L 522 5.81e+15 4.66e+15 9.86e+15 False False False
M 532 2.49e+16 3.32e+16 9.2e+16 False False False
N 604 2.48e+15 2.36e+15 6.34e+15 False False False
P1 P4 P7 reading locations min_peak max_peak share_pct match
a a a a 10 2.13e+15 4.69e+19 20.0 False
a a a b 10 2.13e+15 4.69e+19 0.0 False
a a a c 10 2.13e+15 4.69e+19 30.0 False
a a b a 10 2.01e+15 4.62e+16 20.0 False
a a b b 10 2.01e+15 4.62e+16 10.0 False
a a b c 10 2.01e+15 4.62e+16 10.0 False
a b a a 14 1.34e+15 4.69e+19 0.0 False
a b a b 14 1.34e+15 4.69e+19 0.0 False
a b a c 14 1.34e+15 4.69e+19 28.6 False
a b b a 14 1.14e+15 3.78e+16 7.1 False
a b b b 14 1.14e+15 3.78e+16 0.0 False
a b b c 14 1.14e+15 3.78e+16 0.0 False
b a a a 9 2.35e+15 2.65e+19 22.2 False
b a a b 9 2.35e+15 2.65e+19 11.1 False
b a a c 9 2.35e+15 2.65e+19 22.2 False
b a b a 9 1.92e+15 5.46e+16 22.2 False
b a b b 9 1.92e+15 5.46e+16 11.1 False
b a b c 9 1.92e+15 5.46e+16 11.1 False
b b a a 14 8.93e+14 2.65e+19 0.0 False
b b a b 14 8.93e+14 2.65e+19 0.0 False
b b a c 14 8.93e+14 2.65e+19 35.7 False
b b b a 14 1.26e+15 5.36e+16 0.0 False
b b b b 14 1.26e+15 5.36e+16 0.0 False
b b b c 14 1.26e+15 5.36e+16 0.0 False
Result: not reproduced. No combination matches. The smallest per-location peak is 8.93e14 to 2.35e15, against 3.23e14; the largest 3.78e16 to 4.69e19, against 1.98e16; and the share of locations within 37% of their peak is 0 to 35.7%, against 65%. The per-location medians (P1 a, P4 b) run from 1.72e15 (G) to 1.11e17 (C), so between-location variation is large, as the paper says.
fig5b-polynomial-closest-resistance¶
The distribution of hydraulic resistance of the polynomial fit matches that of the segmentation most closely, followed by the ellipse fit, and the polynomial fit is the best match for the majority of cross sections; the distribution peaks are 1.44e15 (polynomial), 7.06e14 (ellipse), 1.07e16 (spline) and 1.32e15 (optimal annulus) Pa s/m^4.
The fits' resistances are solved with the circle fit of each erosion of P6. The annulus follows Equation 3 as the claims file reads it, with K the segmentation's area ratio (P5 a or b) and R1 of the same erosion. Competing models: (a) the three fits, (b) the fits and the annulus. Peaks on linear bins (P7); the logarithmic bins are shown as a variant.
In [16]:
def models(d, p5, k=3):
K = ratio_seg(d, p5, k)
return pd.DataFrame({"segmentation": d.R_seg, "polynomial": d[f"R_pol_e{k}"], "ellipse": d[f"R_ell_e{k}"],
"spline": d[f"R_spl_e{k}"], "annulus": pvs.optimal_annulus(K, d[f"R1_e{k}"] * 1e-6)}
).replace([np.inf, -np.inf], np.nan)
EXP_PK = {"polynomial": (1.44e15, 1.44e14), "ellipse": (7.06e14, 7.06e13), "spline": (1.07e16, 1.07e15), "annulus": (1.32e15, 1.32e14)}
rows, rows_v = [], []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k, p5 in [(k, p5) for k in P6 for p5 in "ab"]:
M = models(d, p5, k)
ok = M.notna().all(axis=1)
M = M[ok]
for comp, names in (("a", ["polynomial", "ellipse", "spline"]), ("b", ["polynomial", "ellipse", "spline", "annulus"])):
best = M[names].sub(M.segmentation, axis=0).abs().idxmin(axis=1).value_counts(normalize=True).reindex(names).fillna(0) * 100
for p7 in P7:
for sc in ("linear", "log"):
pk = {n: rpeak(M[n], p7, sc) for n in ["segmentation"] + list(EXP_PK)}
r = {"P1": p1, "P4": p4, "P6": k, "P5": p5, "models": comp, "P7": p7, "bins": sc, "n": len(M), "dropped": int((~ok).sum()),
**{f"best_{n[:4]}": round(best.get(n, np.nan), 1) for n in ["polynomial", "ellipse", "spline", "annulus"]},
**{f"pk_{n[:4]}": f"{pk[n]:.3g}" for n in pk},
"peaks_within": sum(within(pk[n], *EXP_PK[n]) for n in EXP_PK)}
r["match"] = best["polynomial"] > 50 and r["peaks_within"] == 4
r["partial"] = best["polynomial"] == best.max()
(rows if sc == "linear" else rows_v).append(r)
f5b = pd.DataFrame(rows)
print(f5b.to_string(index=False))
print("\nVariant, logarithmic bins (not part of the criterion):")
print(pd.DataFrame(rows_v).to_string(index=False))
print("\nSummary by competing models (linear bins): polynomial best-match share and peaks within tolerance")
print(f5b.groupby("models").agg(best_poly_min=("best_poly", "min"), best_poly_max=("best_poly", "max"),
partial=("partial", "sum"), combos=("partial", "size"), peaks_within_max=("peaks_within", "max")).to_string())
M = models(sel("a", "b"), "a")
print("\nMedian resistance by model (P1 a, P4 b, P5 a):", {c: f"{M[c].median():.3g}" for c in M})
print("Distribution distance to the segmentation (two-sample KS statistic), P1 a, P4 b:",
{c: round(stats.ks_2samp(M[c].dropna(), M.segmentation.dropna()).statistic, 3) for c in M if c != "segmentation"})
P1 P4 P6 P5 models P7 bins n dropped best_poly best_elli best_spli best_annu pk_segm pk_poly pk_elli pk_spli pk_annu peaks_within match partial
a a 1 a a a linear 3303 399 81.4 0.2 18.5 NaN 1.78e+16 1.72e+16 7.03e+15 5.12e+15 1.29e+16 0 False True
a a 1 a a b linear 3303 399 81.4 0.2 18.5 NaN 1.89e+15 1.27e+15 5.45e+14 7.74e+15 7.29e+14 0 False True
a a 1 a b a linear 3303 399 69.3 0.1 17.5 13.1 1.78e+16 1.72e+16 7.03e+15 5.12e+15 1.29e+16 0 False True
a a 1 a b b linear 3303 399 69.3 0.1 17.5 13.1 1.89e+15 1.27e+15 5.45e+14 7.74e+15 7.29e+14 0 False True
a a 1 b a a linear 3303 399 81.4 0.2 18.5 NaN 1.78e+16 1.72e+16 7.03e+15 5.12e+15 1.85e+16 0 False True
a a 1 b a b linear 3303 399 81.4 0.2 18.5 NaN 1.89e+15 1.27e+15 5.45e+14 7.74e+15 1.12e+15 0 False True
a a 1 b b a linear 3303 399 53.7 0.1 15.4 30.9 1.78e+16 1.72e+16 7.03e+15 5.12e+15 1.85e+16 0 False True
a a 1 b b b linear 3303 399 53.7 0.1 15.4 30.9 1.89e+15 1.27e+15 5.45e+14 7.74e+15 1.12e+15 0 False True
a a 3 a a a linear 3306 396 82.3 0.0 17.7 NaN 1.78e+16 1.77e+16 5.95e+15 5.13e+15 1.62e+16 0 False True
a a 3 a a b linear 3306 396 82.3 0.0 17.7 NaN 1.9e+15 1.22e+15 5.35e+14 9.2e+15 8.85e+14 0 False True
a a 3 a b a linear 3306 396 61.3 0.0 16.1 22.6 1.78e+16 1.77e+16 5.95e+15 5.13e+15 1.62e+16 0 False True
a a 3 a b b linear 3306 396 61.3 0.0 16.1 22.6 1.9e+15 1.22e+15 5.35e+14 9.2e+15 8.85e+14 0 False True
a a 3 b a a linear 3306 396 82.3 0.0 17.7 NaN 1.78e+16 1.77e+16 5.95e+15 5.13e+15 1.86e+16 0 False True
a a 3 b a b linear 3306 396 82.3 0.0 17.7 NaN 1.9e+15 1.22e+15 5.35e+14 9.2e+15 1.13e+15 0 False True
a a 3 b b a linear 3306 396 53.9 0.0 13.3 32.8 1.78e+16 1.77e+16 5.95e+15 5.13e+15 1.86e+16 0 False True
a a 3 b b b linear 3306 396 53.9 0.0 13.3 32.8 1.9e+15 1.22e+15 5.35e+14 9.2e+15 1.13e+15 0 False True
a a 5 a a a linear 3314 388 82.6 0.1 17.3 NaN 1.78e+16 1.78e+16 5.5e+15 5.31e+15 1.76e+16 0 False True
a a 5 a a b linear 3314 388 82.6 0.1 17.3 NaN 1.9e+15 1.19e+15 5.19e+14 1.03e+16 9.38e+14 1 False True
a a 5 a b a linear 3314 388 57.5 0.1 15.0 27.4 1.78e+16 1.78e+16 5.5e+15 5.31e+15 1.76e+16 0 False True
a a 5 a b b linear 3314 388 57.5 0.1 15.0 27.4 1.9e+15 1.19e+15 5.19e+14 1.03e+16 9.38e+14 1 False True
a a 5 b a a linear 3314 388 82.6 0.1 17.3 NaN 1.78e+16 1.78e+16 5.5e+15 5.31e+15 1.87e+16 0 False True
a a 5 b a b linear 3314 388 82.6 0.1 17.3 NaN 1.9e+15 1.19e+15 5.19e+14 1.03e+16 1.13e+15 1 False True
a a 5 b b a linear 3314 388 52.8 0.0 13.6 33.6 1.78e+16 1.78e+16 5.5e+15 5.31e+15 1.87e+16 0 False True
a a 5 b b b linear 3314 388 52.8 0.0 13.6 33.6 1.9e+15 1.19e+15 5.19e+14 1.03e+16 1.13e+15 1 False True
a b 1 a a a linear 5775 908 81.4 2.6 16.0 NaN 4e+18 2.19e+17 3.02e+16 1.06e+18 5.51e+16 0 False True
a b 1 a a b linear 5775 908 81.4 2.6 16.0 NaN 1.35e+15 8.64e+14 3.98e+14 9.64e+15 1.22e+15 2 False True
a b 1 a b a linear 5775 908 71.4 0.9 15.2 12.6 4e+18 2.19e+17 3.02e+16 1.06e+18 5.51e+16 0 False True
a b 1 a b b linear 5775 908 71.4 0.9 15.2 12.6 1.35e+15 8.64e+14 3.98e+14 9.64e+15 1.22e+15 2 False True
a b 1 b a a linear 5775 908 81.4 2.6 16.0 NaN 4e+18 2.19e+17 3.02e+16 1.06e+18 1.48e+17 0 False True
a b 1 b a b linear 5775 908 81.4 2.6 16.0 NaN 1.35e+15 8.64e+14 3.98e+14 9.64e+15 1.68e+15 1 False True
a b 1 b b a linear 5775 908 59.4 0.5 13.8 26.3 4e+18 2.19e+17 3.02e+16 1.06e+18 1.48e+17 0 False True
a b 1 b b b linear 5775 908 59.4 0.5 13.8 26.3 1.35e+15 8.64e+14 3.98e+14 9.64e+15 1.68e+15 1 False True
a b 3 a a a linear 5783 900 82.1 2.6 15.4 NaN 4e+18 2.19e+17 1.65e+16 1.93e+18 8.48e+16 0 False True
a b 3 a a b linear 5783 900 82.1 2.6 15.4 NaN 1.36e+15 8.48e+14 3.9e+14 1.15e+16 1.35e+15 2 False True
a b 3 a b a linear 5783 900 65.9 0.8 14.1 19.2 4e+18 2.19e+17 1.65e+16 1.93e+18 8.48e+16 0 False True
a b 3 a b b linear 5783 900 65.9 0.8 14.1 19.2 1.36e+15 8.48e+14 3.9e+14 1.15e+16 1.35e+15 2 False True
a b 3 b a a linear 5783 900 82.1 2.6 15.4 NaN 4e+18 2.19e+17 1.65e+16 1.93e+18 1.49e+17 0 False True
a b 3 b a b linear 5783 900 82.1 2.6 15.4 NaN 1.36e+15 8.48e+14 3.9e+14 1.15e+16 1.69e+15 1 False True
a b 3 b b a linear 5783 900 58.9 0.4 12.2 28.5 4e+18 2.19e+17 1.65e+16 1.93e+18 1.49e+17 0 False True
a b 3 b b b linear 5783 900 58.9 0.4 12.2 28.5 1.36e+15 8.48e+14 3.9e+14 1.15e+16 1.69e+15 1 False True
a b 5 a a a linear 5793 890 82.3 2.7 15.0 NaN 4e+18 2.19e+17 1.32e+16 1.85e+18 1.12e+17 0 False True
a b 5 a a b linear 5793 890 82.3 2.7 15.0 NaN 1.36e+15 8.37e+14 3.82e+14 1.25e+16 1.43e+15 1 False True
a b 5 a b a linear 5793 890 62.8 0.8 12.8 23.6 4e+18 2.19e+17 1.32e+16 1.85e+18 1.12e+17 0 False True
a b 5 a b b linear 5793 890 62.8 0.8 12.8 23.6 1.36e+15 8.37e+14 3.82e+14 1.25e+16 1.43e+15 1 False True
a b 5 b a a linear 5793 890 82.3 2.7 15.0 NaN 4e+18 2.19e+17 1.32e+16 1.85e+18 1.5e+17 0 False True
a b 5 b a b linear 5793 890 82.3 2.7 15.0 NaN 1.36e+15 8.37e+14 3.82e+14 1.25e+16 1.69e+15 0 False True
a b 5 b b a linear 5793 890 57.9 0.1 11.8 30.1 4e+18 2.19e+17 1.32e+16 1.85e+18 1.5e+17 0 False True
a b 5 b b b linear 5793 890 57.9 0.1 11.8 30.1 1.36e+15 8.37e+14 3.82e+14 1.25e+16 1.69e+15 0 False True
b a 1 a a a linear 2135 260 82.6 1.2 16.2 NaN 2.13e+16 2.45e+16 6.57e+15 8e+15 1.38e+16 0 False True
b a 1 a a b linear 2135 260 82.6 1.2 16.2 NaN 1.9e+15 1.58e+15 5.89e+14 1.09e+16 7.3e+14 2 False True
b a 1 a b a linear 2135 260 76.5 0.7 15.9 6.9 2.13e+16 2.45e+16 6.57e+15 8e+15 1.38e+16 0 False True
b a 1 a b b linear 2135 260 76.5 0.7 15.9 6.9 1.9e+15 1.58e+15 5.89e+14 1.09e+16 7.3e+14 2 False True
b a 1 b a a linear 2135 260 82.6 1.2 16.2 NaN 2.13e+16 2.45e+16 6.57e+15 8e+15 1.59e+16 0 False True
b a 1 b a b linear 2135 260 82.6 1.2 16.2 NaN 1.9e+15 1.58e+15 5.89e+14 1.09e+16 1.14e+15 2 False True
b a 1 b b a linear 2135 260 64.4 0.6 14.2 20.8 2.13e+16 2.45e+16 6.57e+15 8e+15 1.59e+16 0 False True
b a 1 b b b linear 2135 260 64.4 0.6 14.2 20.8 1.9e+15 1.58e+15 5.89e+14 1.09e+16 1.14e+15 2 False True
b a 3 a a a linear 2135 260 84.7 0.9 14.4 NaN 2.13e+16 2.45e+16 6.12e+15 7.83e+15 1.58e+16 0 False True
b a 3 a a b linear 2135 260 84.7 0.9 14.4 NaN 1.9e+15 1.5e+15 5.59e+14 1.38e+16 9.22e+14 1 False True
b a 3 a b a linear 2135 260 70.6 0.7 13.7 15.0 2.13e+16 2.45e+16 6.12e+15 7.83e+15 1.58e+16 0 False True
b a 3 a b b linear 2135 260 70.6 0.7 13.7 15.0 1.9e+15 1.5e+15 5.59e+14 1.38e+16 9.22e+14 1 False True
b a 3 b a a linear 2135 260 84.7 0.9 14.4 NaN 2.13e+16 2.45e+16 6.12e+15 7.83e+15 1.59e+16 0 False True
b a 3 b a b linear 2135 260 84.7 0.9 14.4 NaN 1.9e+15 1.5e+15 5.59e+14 1.38e+16 1.14e+15 1 False True
b a 3 b b a linear 2135 260 64.3 0.7 13.1 21.9 2.13e+16 2.45e+16 6.12e+15 7.83e+15 1.59e+16 0 False True
b a 3 b b b linear 2135 260 64.3 0.7 13.1 21.9 1.9e+15 1.5e+15 5.59e+14 1.38e+16 1.14e+15 1 False True
b a 5 a a a linear 2133 262 85.7 0.8 13.5 NaN 2.13e+16 2.45e+16 6.11e+15 9.38e+15 1.58e+16 0 False True
b a 5 a a b linear 2133 262 85.7 0.8 13.5 NaN 1.9e+15 1.48e+15 5.43e+14 1.51e+16 1.01e+15 1 False True
b a 5 a b a linear 2133 262 66.9 0.6 12.7 19.8 2.13e+16 2.45e+16 6.11e+15 9.38e+15 1.58e+16 0 False True
b a 5 a b b linear 2133 262 66.9 0.6 12.7 19.8 1.9e+15 1.48e+15 5.43e+14 1.51e+16 1.01e+15 1 False True
b a 5 b a a linear 2133 262 85.7 0.8 13.5 NaN 2.13e+16 2.45e+16 6.11e+15 9.38e+15 1.59e+16 0 False True
b a 5 b a b linear 2133 262 85.7 0.8 13.5 NaN 1.9e+15 1.48e+15 5.43e+14 1.51e+16 1.15e+15 1 False True
b a 5 b b a linear 2133 262 64.0 0.6 12.3 23.1 2.13e+16 2.45e+16 6.11e+15 9.38e+15 1.59e+16 0 False True
b a 5 b b b linear 2133 262 64.0 0.6 12.3 23.1 1.9e+15 1.48e+15 5.43e+14 1.51e+16 1.15e+15 1 False True
b b 1 a a a linear 3716 643 81.6 3.4 15.0 NaN 8.26e+18 3.34e+17 9.52e+16 6.76e+18 7.35e+16 0 False True
b b 1 a a b linear 3716 643 81.6 3.4 15.0 NaN 1.53e+15 1.12e+15 4.4e+14 1.31e+16 1.28e+15 1 False True
b b 1 a b a linear 3716 643 75.6 1.3 14.7 8.4 8.26e+18 3.34e+17 9.52e+16 6.76e+18 7.35e+16 0 False True
b b 1 a b b linear 3716 643 75.6 1.3 14.7 8.4 1.53e+15 1.12e+15 4.4e+14 1.31e+16 1.28e+15 1 False True
b b 1 b a a linear 3716 643 81.6 3.4 15.0 NaN 8.26e+18 3.34e+17 9.52e+16 6.76e+18 1.14e+17 0 False True
b b 1 b a b linear 3716 643 81.6 3.4 15.0 NaN 1.53e+15 1.12e+15 4.4e+14 1.31e+16 1.76e+15 0 False True
b b 1 b b a linear 3716 643 67.3 0.8 13.6 18.2 8.26e+18 3.34e+17 9.52e+16 6.76e+18 1.14e+17 0 False True
b b 1 b b b linear 3716 643 67.3 0.8 13.6 18.2 1.53e+15 1.12e+15 4.4e+14 1.31e+16 1.76e+15 0 False True
b b 3 a a a linear 3729 630 83.1 3.1 13.8 NaN 8.26e+18 3.31e+17 4.2e+16 1.51e+19 8.51e+16 0 False True
b b 3 a a b linear 3729 630 83.1 3.1 13.8 NaN 1.54e+15 1.06e+15 4.23e+14 5.74e+15 1.49e+15 0 False True
b b 3 a b a linear 3729 630 71.0 1.1 12.7 15.2 8.26e+18 3.31e+17 4.2e+16 1.51e+19 8.51e+16 0 False True
b b 3 a b b linear 3729 630 71.0 1.1 12.7 15.2 1.54e+15 1.06e+15 4.23e+14 5.74e+15 1.49e+15 0 False True
b b 3 b a a linear 3729 630 83.1 3.1 13.8 NaN 8.26e+18 3.31e+17 4.2e+16 1.51e+19 1.15e+17 0 False True
b b 3 b a b linear 3729 630 83.1 3.1 13.8 NaN 1.54e+15 1.06e+15 4.23e+14 5.74e+15 1.77e+15 0 False True
b b 3 b b a linear 3729 630 66.9 0.7 12.4 20.0 8.26e+18 3.31e+17 4.2e+16 1.51e+19 1.15e+17 0 False True
b b 3 b b b linear 3729 630 66.9 0.7 12.4 20.0 1.54e+15 1.06e+15 4.23e+14 5.74e+15 1.77e+15 0 False True
b b 5 a a a linear 3732 627 84.3 3.1 12.6 NaN 8.26e+18 3.31e+17 2.32e+16 2.35e+19 1.01e+17 0 False True
b b 5 a a b linear 3732 627 84.3 3.1 12.6 NaN 1.54e+15 1.04e+15 4.13e+14 6.31e+15 1.58e+15 0 False True
b b 5 a b a linear 3732 627 68.2 1.0 11.5 19.3 8.26e+18 3.31e+17 2.32e+16 2.35e+19 1.01e+17 0 False True
b b 5 a b b linear 3732 627 68.2 1.0 11.5 19.3 1.54e+15 1.04e+15 4.13e+14 6.31e+15 1.58e+15 0 False True
b b 5 b a a linear 3732 627 84.3 3.1 12.6 NaN 8.26e+18 3.31e+17 2.32e+16 2.35e+19 1.15e+17 0 False True
b b 5 b a b linear 3732 627 84.3 3.1 12.6 NaN 1.54e+15 1.04e+15 4.13e+14 6.31e+15 1.78e+15 0 False True
b b 5 b b a linear 3732 627 66.5 0.6 11.3 21.5 8.26e+18 3.31e+17 2.32e+16 2.35e+19 1.15e+17 0 False True
b b 5 b b b linear 3732 627 66.5 0.6 11.3 21.5 1.54e+15 1.04e+15 4.13e+14 6.31e+15 1.78e+15 0 False True
Variant, logarithmic bins (not part of the criterion):
P1 P4 P6 P5 models P7 bins n dropped best_poly best_elli best_spli best_annu pk_segm pk_poly pk_elli pk_spli pk_annu peaks_within match partial
a a 1 a a a log 3303 399 81.4 0.2 18.5 NaN 2.07e+15 1.72e+15 1.45e+15 9.32e+15 1.07e+15 0 False True
a a 1 a a b log 3303 399 81.4 0.2 18.5 NaN 2.57e+15 1.4e+15 5.89e+14 1.22e+16 1.4e+15 2 False True
a a 1 a b a log 3303 399 69.3 0.1 17.5 13.1 2.07e+15 1.72e+15 1.45e+15 9.32e+15 1.07e+15 0 False True
a a 1 a b b log 3303 399 69.3 0.1 17.5 13.1 2.57e+15 1.4e+15 5.89e+14 1.22e+16 1.4e+15 2 False True
a a 1 b a a log 3303 399 81.4 0.2 18.5 NaN 2.07e+15 1.72e+15 1.45e+15 9.32e+15 1.38e+15 1 False True
a a 1 b a b log 3303 399 81.4 0.2 18.5 NaN 2.57e+15 1.4e+15 5.89e+14 1.22e+16 1.5e+15 1 False True
a a 1 b b a log 3303 399 53.7 0.1 15.4 30.9 2.07e+15 1.72e+15 1.45e+15 9.32e+15 1.38e+15 1 False True
a a 1 b b b log 3303 399 53.7 0.1 15.4 30.9 2.57e+15 1.4e+15 5.89e+14 1.22e+16 1.5e+15 1 False True
a a 3 a a a log 3306 396 82.3 0.0 17.7 NaN 2.07e+15 1.73e+15 6.34e+14 1.03e+16 1.35e+15 2 False True
a a 3 a a b log 3306 396 82.3 0.0 17.7 NaN 2.58e+15 1.37e+15 5.79e+14 1.07e+16 1.32e+15 3 False True
a a 3 a b a log 3306 396 61.3 0.0 16.1 22.6 2.07e+15 1.73e+15 6.34e+14 1.03e+16 1.35e+15 2 False True
a a 3 a b b log 3306 396 61.3 0.0 16.1 22.6 2.58e+15 1.37e+15 5.79e+14 1.07e+16 1.32e+15 3 False True
a a 3 b a a log 3306 396 82.3 0.0 17.7 NaN 2.07e+15 1.73e+15 6.34e+14 1.03e+16 1.38e+15 2 False True
a a 3 b a b log 3306 396 82.3 0.0 17.7 NaN 2.58e+15 1.37e+15 5.79e+14 1.07e+16 1.51e+15 2 False True
a a 3 b b a log 3306 396 53.9 0.0 13.3 32.8 2.07e+15 1.73e+15 6.34e+14 1.03e+16 1.38e+15 2 False True
a a 3 b b b log 3306 396 53.9 0.0 13.3 32.8 2.58e+15 1.37e+15 5.79e+14 1.07e+16 1.51e+15 2 False True
a a 5 a a a log 3314 388 82.6 0.1 17.3 NaN 2.07e+15 1.71e+15 6.03e+14 1.39e+16 1.25e+15 1 False True
a a 5 a a b log 3314 388 82.6 0.1 17.3 NaN 2.58e+15 1.35e+15 5.58e+14 1.07e+16 1.26e+15 3 False True
a a 5 a b a log 3314 388 57.5 0.1 15.0 27.4 2.07e+15 1.71e+15 6.03e+14 1.39e+16 1.25e+15 1 False True
a a 5 a b b log 3314 388 57.5 0.1 15.0 27.4 2.58e+15 1.35e+15 5.58e+14 1.07e+16 1.26e+15 3 False True
a a 5 b a a log 3314 388 82.6 0.1 17.3 NaN 2.07e+15 1.71e+15 6.03e+14 1.39e+16 1.39e+15 1 False True
a a 5 b a b log 3314 388 82.6 0.1 17.3 NaN 2.58e+15 1.35e+15 5.58e+14 1.07e+16 1.51e+15 2 False True
a a 5 b b a log 3314 388 52.8 0.0 13.6 33.6 2.07e+15 1.71e+15 6.03e+14 1.39e+16 1.39e+15 1 False True
a a 5 b b b log 3314 388 52.8 0.0 13.6 33.6 2.58e+15 1.35e+15 5.58e+14 1.07e+16 1.51e+15 2 False True
a b 1 a a a log 5775 908 81.4 2.6 16.0 NaN 2.04e+15 1.38e+15 9.28e+14 8.81e+15 9.28e+14 1 False True
a b 1 a a b log 5775 908 81.4 2.6 16.0 NaN 2.23e+15 2.09e+15 8.15e+14 7.19e+15 1.52e+15 0 False True
a b 1 a b a log 5775 908 71.4 0.9 15.2 12.6 2.04e+15 1.38e+15 9.28e+14 8.81e+15 9.28e+14 1 False True
a b 1 a b b log 5775 908 71.4 0.9 15.2 12.6 2.23e+15 2.09e+15 8.15e+14 7.19e+15 1.52e+15 0 False True
a b 1 b a a log 5775 908 81.4 2.6 16.0 NaN 2.04e+15 1.38e+15 9.28e+14 8.81e+15 1.32e+15 2 False True
a b 1 b a b log 5775 908 81.4 2.6 16.0 NaN 2.23e+15 2.09e+15 8.15e+14 7.19e+15 1.49e+15 0 False True
a b 1 b b a log 5775 908 59.4 0.5 13.8 26.3 2.04e+15 1.38e+15 9.28e+14 8.81e+15 1.32e+15 2 False True
a b 1 b b b log 5775 908 59.4 0.5 13.8 26.3 2.23e+15 2.09e+15 8.15e+14 7.19e+15 1.49e+15 0 False True
a b 3 a a a log 5783 900 82.1 2.6 15.4 NaN 2.04e+15 1.38e+15 7.99e+14 6.61e+15 1.31e+15 2 False True
a b 3 a a b log 5783 900 82.1 2.6 15.4 NaN 2.24e+15 1.56e+15 8e+14 8.81e+15 1.16e+15 1 False True
a b 3 a b a log 5783 900 65.9 0.8 14.1 19.2 2.04e+15 1.38e+15 7.99e+14 6.61e+15 1.31e+15 2 False True
a b 3 a b b log 5783 900 65.9 0.8 14.1 19.2 2.24e+15 1.56e+15 8e+14 8.81e+15 1.16e+15 1 False True
a b 3 b a a log 5783 900 82.1 2.6 15.4 NaN 2.04e+15 1.38e+15 7.99e+14 6.61e+15 1.33e+15 2 False True
a b 3 b a b log 5783 900 82.1 2.6 15.4 NaN 2.24e+15 1.56e+15 8e+14 8.81e+15 1.5e+15 1 False True
a b 3 b b a log 5783 900 58.9 0.4 12.2 28.5 2.04e+15 1.38e+15 7.99e+14 6.61e+15 1.33e+15 2 False True
a b 3 b b b log 5783 900 58.9 0.4 12.2 28.5 2.24e+15 1.56e+15 8e+14 8.81e+15 1.5e+15 1 False True
a b 5 a a a log 5793 890 82.3 2.7 15.0 NaN 2.04e+15 1.37e+15 7.36e+14 1.08e+16 1.69e+15 3 False True
a b 5 a a b log 5793 890 82.3 2.7 15.0 NaN 2.24e+15 1.54e+15 9.78e+14 8.53e+15 1.3e+15 2 False True
a b 5 a b a log 5793 890 62.8 0.8 12.8 23.6 2.04e+15 1.37e+15 7.36e+14 1.08e+16 1.69e+15 3 False True
a b 5 a b b log 5793 890 62.8 0.8 12.8 23.6 2.24e+15 1.54e+15 9.78e+14 8.53e+15 1.3e+15 2 False True
a b 5 b a a log 5793 890 82.3 2.7 15.0 NaN 2.04e+15 1.37e+15 7.36e+14 1.08e+16 1.33e+15 4 True True
a b 5 b a b log 5793 890 82.3 2.7 15.0 NaN 2.24e+15 1.54e+15 9.78e+14 8.53e+15 1.5e+15 1 False True
a b 5 b b a log 5793 890 57.9 0.1 11.8 30.1 2.04e+15 1.37e+15 7.36e+14 1.08e+16 1.33e+15 4 True True
a b 5 b b b log 5793 890 57.9 0.1 11.8 30.1 2.24e+15 1.54e+15 9.78e+14 8.53e+15 1.5e+15 1 False True
b a 1 a a a log 2135 260 82.6 1.2 16.2 NaN 2.88e+15 1.47e+15 6.58e+14 1.39e+16 1.06e+15 2 False True
b a 1 a a b log 2135 260 82.6 1.2 16.2 NaN 2.39e+15 1.34e+15 6.09e+14 1.24e+16 1.39e+15 2 False True
b a 1 a b a log 2135 260 76.5 0.7 15.9 6.9 2.88e+15 1.47e+15 6.58e+14 1.39e+16 1.06e+15 2 False True
b a 1 a b b log 2135 260 76.5 0.7 15.9 6.9 2.39e+15 1.34e+15 6.09e+14 1.24e+16 1.39e+15 2 False True
b a 1 b a a log 2135 260 82.6 1.2 16.2 NaN 2.88e+15 1.47e+15 6.58e+14 1.39e+16 1.31e+15 3 False True
b a 1 b a b log 2135 260 82.6 1.2 16.2 NaN 2.39e+15 1.34e+15 6.09e+14 1.24e+16 1.49e+15 1 False True
b a 1 b b a log 2135 260 64.4 0.6 14.2 20.8 2.88e+15 1.47e+15 6.58e+14 1.39e+16 1.31e+15 3 False True
b a 1 b b b log 2135 260 64.4 0.6 14.2 20.8 2.39e+15 1.34e+15 6.09e+14 1.24e+16 1.49e+15 1 False True
b a 3 a a a log 2135 260 84.7 0.9 14.4 NaN 2.88e+15 1.42e+15 6.11e+14 1.51e+16 1.37e+15 2 False True
b a 3 a a b log 2135 260 84.7 0.9 14.4 NaN 2.39e+15 1.67e+15 5.7e+14 1.16e+16 1.34e+15 2 False True
b a 3 a b a log 2135 260 70.6 0.7 13.7 15.0 2.88e+15 1.42e+15 6.11e+14 1.51e+16 1.37e+15 2 False True
b a 3 a b b log 2135 260 70.6 0.7 13.7 15.0 2.39e+15 1.67e+15 5.7e+14 1.16e+16 1.34e+15 2 False True
b a 3 b a a log 2135 260 84.7 0.9 14.4 NaN 2.88e+15 1.42e+15 6.11e+14 1.51e+16 1.31e+15 2 False True
b a 3 b a b log 2135 260 84.7 0.9 14.4 NaN 2.39e+15 1.67e+15 5.7e+14 1.16e+16 1.5e+15 1 False True
b a 3 b b a log 2135 260 64.3 0.7 13.1 21.9 2.88e+15 1.42e+15 6.11e+14 1.51e+16 1.31e+15 2 False True
b a 3 b b b log 2135 260 64.3 0.7 13.1 21.9 2.39e+15 1.67e+15 5.7e+14 1.16e+16 1.5e+15 1 False True
b a 5 a a a log 2133 262 85.7 0.8 13.5 NaN 2.88e+15 1.4e+15 5.81e+14 1.34e+16 1.64e+15 1 False True
b a 5 a a b log 2133 262 85.7 0.8 13.5 NaN 2.39e+15 1.65e+15 5.42e+14 1.34e+16 1.59e+15 0 False True
b a 5 a b a log 2133 262 66.9 0.6 12.7 19.8 2.88e+15 1.4e+15 5.81e+14 1.34e+16 1.64e+15 1 False True
b a 5 a b b log 2133 262 66.9 0.6 12.7 19.8 2.39e+15 1.65e+15 5.42e+14 1.34e+16 1.59e+15 0 False True
b a 5 b a a log 2133 262 85.7 0.8 13.5 NaN 2.88e+15 1.4e+15 5.81e+14 1.34e+16 1.32e+15 2 False True
b a 5 b a b log 2133 262 85.7 0.8 13.5 NaN 2.39e+15 1.65e+15 5.42e+14 1.34e+16 1.51e+15 0 False True
b a 5 b b a log 2133 262 64.0 0.6 12.3 23.1 2.88e+15 1.4e+15 5.81e+14 1.34e+16 1.32e+15 2 False True
b a 5 b b b log 2133 262 64.0 0.6 12.3 23.1 2.39e+15 1.65e+15 5.42e+14 1.34e+16 1.51e+15 0 False True
b b 1 a a a log 3716 643 81.6 3.4 15.0 NaN 1.86e+15 1.64e+15 1.2e+15 7.96e+15 1.34e+15 1 False True
b b 1 a a b log 3716 643 81.6 3.4 15.0 NaN 2.75e+15 1.9e+15 8.64e+14 8.69e+15 1.18e+15 0 False True
b b 1 a b a log 3716 643 75.6 1.3 14.7 8.4 1.86e+15 1.64e+15 1.2e+15 7.96e+15 1.34e+15 1 False True
b b 1 a b b log 3716 643 75.6 1.3 14.7 8.4 2.75e+15 1.9e+15 8.64e+14 8.69e+15 1.18e+15 0 False True
b b 1 b a a log 3716 643 81.6 3.4 15.0 NaN 1.86e+15 1.64e+15 1.2e+15 7.96e+15 1.24e+15 1 False True
b b 1 b a b log 3716 643 81.6 3.4 15.0 NaN 2.75e+15 1.9e+15 8.64e+14 8.69e+15 1.51e+15 0 False True
b b 1 b b a log 3716 643 67.3 0.8 13.6 18.2 1.86e+15 1.64e+15 1.2e+15 7.96e+15 1.24e+15 1 False True
b b 1 b b b log 3716 643 67.3 0.8 13.6 18.2 2.75e+15 1.9e+15 8.64e+14 8.69e+15 1.51e+15 0 False True
b b 3 a a a log 3729 630 83.1 3.1 13.8 NaN 1.86e+15 1.62e+15 9.45e+14 9.6e+15 1.56e+15 0 False True
b b 3 a a b log 3729 630 83.1 3.1 13.8 NaN 2.76e+15 1.83e+15 1.04e+15 9.89e+15 1.37e+15 2 False True
b b 3 a b a log 3729 630 71.0 1.1 12.7 15.2 1.86e+15 1.62e+15 9.45e+14 9.6e+15 1.56e+15 0 False True
b b 3 a b b log 3729 630 71.0 1.1 12.7 15.2 2.76e+15 1.83e+15 1.04e+15 9.89e+15 1.37e+15 2 False True
b b 3 b a a log 3729 630 83.1 3.1 13.8 NaN 1.86e+15 1.62e+15 9.45e+14 9.6e+15 1.25e+15 1 False True
b b 3 b a b log 3729 630 83.1 3.1 13.8 NaN 2.76e+15 1.83e+15 1.04e+15 9.89e+15 1.52e+15 1 False True
b b 3 b b a log 3729 630 66.9 0.7 12.4 20.0 1.86e+15 1.62e+15 9.45e+14 9.6e+15 1.25e+15 1 False True
b b 3 b b b log 3729 630 66.9 0.7 12.4 20.0 2.76e+15 1.83e+15 1.04e+15 9.89e+15 1.52e+15 1 False True
b b 5 a a a log 3732 627 84.3 3.1 12.6 NaN 1.86e+15 1.61e+15 7.86e+14 1.34e+16 1.07e+15 0 False True
b b 5 a a b log 3732 627 84.3 3.1 12.6 NaN 2.76e+15 1.81e+15 9.99e+14 1.18e+16 1.6e+15 0 False True
b b 5 a b a log 3732 627 68.2 1.0 11.5 19.3 1.86e+15 1.61e+15 7.86e+14 1.34e+16 1.07e+15 0 False True
b b 5 a b b log 3732 627 68.2 1.0 11.5 19.3 2.76e+15 1.81e+15 9.99e+14 1.18e+16 1.6e+15 0 False True
b b 5 b a a log 3732 627 84.3 3.1 12.6 NaN 1.86e+15 1.61e+15 7.86e+14 1.34e+16 1.25e+15 1 False True
b b 5 b a b log 3732 627 84.3 3.1 12.6 NaN 2.76e+15 1.81e+15 9.99e+14 1.18e+16 1.53e+15 0 False True
b b 5 b b a log 3732 627 66.5 0.6 11.3 21.5 1.86e+15 1.61e+15 7.86e+14 1.34e+16 1.25e+15 1 False True
b b 5 b b b log 3732 627 66.5 0.6 11.3 21.5 2.76e+15 1.81e+15 9.99e+14 1.18e+16 1.53e+15 0 False True
Summary by competing models (linear bins): polynomial best-match share and peaks within tolerance
best_poly_min best_poly_max partial combos peaks_within_max
models
a 81.4 85.7 48 48 2
b 52.8 76.5 48 48 2
Median resistance by model (P1 a, P4 b, P5 a): {'segmentation': '5.28e+15', 'polynomial': '3.52e+15', 'ellipse': '1.33e+15', 'spline': '1.07e+16', 'annulus': '2.56e+15'}
Distribution distance to the segmentation (two-sample KS statistic), P1 a, P4 b: {'polynomial': np.float64(0.207), 'ellipse': np.float64(0.42), 'spline': np.float64(0.346), 'annulus': np.float64(0.272)}
Result: partially reproduced. The polynomial fit is the closest model for the largest share of cross sections in every one of the 96 combinations: 81.4 to 85.7% against the three fits, and 52.8 to 76.5% when the annulus competes. No combination puts all four peaks within tolerance on linear bins, which the criterion uses: at most 2 (for example P1 a, P4 b, erosion 3, P5 a, outliers out: polynomial 8.48e14, ellipse 3.9e14, spline 1.15e16 and annulus 1.35e15, the last two within). On logarithmic bins, a variant the claims file names for fig6b but not for this claim, two combinations put all four within (P1 a, P4 b, erosion 5, P5 b, all values: 1.37e15, 7.36e14, 1.08e16 and 1.33e15, with the polynomial closest in 82.3% and 57.9%). By the two-sample KS statistic (P1 a, P4 b, erosion 3) the polynomial distribution is the closest to the segmentation's (0.207), then the annulus (0.272), the spline (0.346) and the ellipse (0.42), so the ellipse is not second, as the statement has it. Median resistances: segmentation 5.28e15, polynomial 3.52e15, annulus 2.56e15, ellipse 1.33e15, spline 1.07e16.
ellipse-annulus-underestimate¶
The resistance distribution of the optimal elliptical annulus is very similar to that of the ellipse fit, and both tend to underestimate the resistance of the segmented PVS.
Readings of "underestimate": (a) the median of the per-cross-section ratio to the segmentation is below 1 for both; (b) the peak and median of both are below the segmentation's. Similarity: the annulus-over-ellipse ratios of the peaks and of the medians within 10% of 1.87 and 1.61.
In [17]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k, p5 in [(k, p5) for k in P6 for p5 in "ab"]:
M = models(d, p5, k).dropna(subset=["segmentation", "ellipse", "annulus"])
for p7 in P7:
for sc in ("linear", "log"):
pk = {n: rpeak(M[n], p7, sc) for n in ("segmentation", "ellipse", "annulus")}
md_ = {n: M[n].median() for n in pk}
ua = (M.ellipse / M.segmentation).median() < 1 and (M.annulus / M.segmentation).median() < 1
ub = all(pk[n] < pk["segmentation"] and md_[n] < md_["segmentation"] for n in ("ellipse", "annulus"))
rp, rm = pk["annulus"] / pk["ellipse"], md_["annulus"] / md_["ellipse"]
sim = within(rp, 1.87, 0.187) and within(rm, 1.61, 0.161)
rows.append({"P1": p1, "P4": p4, "P6": k, "P5": p5, "P7": p7, "bins": sc + ("" if sc == "linear" else " (variant)"), "n": len(M),
"med_ratio_ell": round((M.ellipse / M.segmentation).median(), 3),
"med_ratio_ann": round((M.annulus / M.segmentation).median(), 3),
"under_a": ua, "under_b": ub, "ann/ell_peak": round(rp, 3), "ann/ell_median": round(rm, 3),
"similar": sim, "match": (ua or ub) and sim if sc == "linear" else np.nan})
eau = pd.DataFrame(rows)
print(eau.to_string(index=False))
lin = eau[eau.bins == "linear"]
lin = eau[eau.bins == "linear"]
print("\nCombinations with the similarity ratios within 10%:")
print(lin[lin.similar].to_string(index=False))
print("Ranges on linear bins: annulus/ellipse peak", lin["ann/ell_peak"].min(), "to", lin["ann/ell_peak"].max(), "; median", lin["ann/ell_median"].min(), "to", lin["ann/ell_median"].max())
print("\nLinear bins, as in fig5b's criterion: combinations with both underestimating (a or b):", int((lin.under_a | lin.under_b).sum()),
"of", len(lin), "; with the similarity ratios within 10%:", int(lin.similar.sum()))
P1 P4 P6 P5 P7 bins n med_ratio_ell med_ratio_ann under_a under_b ann/ell_peak ann/ell_median similar match a a 1 a a linear 3590 0.374 0.518 True True 0.644 1.503 False False a a 1 a a log (variant) 3590 0.374 0.518 True True 1.726 1.503 True NaN a a 1 a b linear 3590 0.374 0.518 True True 1.133 1.503 False False a a 1 a b log (variant) 3590 0.374 0.518 True True 1.530 1.503 False NaN a a 1 b a linear 3590 0.374 0.663 True True 0.695 1.875 False False a a 1 b a log (variant) 3590 0.374 0.663 True True 1.694 1.875 False NaN a a 1 b b linear 3590 0.374 0.663 True True 1.919 1.875 False False a a 1 b b log (variant) 3590 0.374 0.663 True True 1.921 1.875 False NaN a a 3 a a linear 3590 0.351 0.589 True True 1.535 1.784 False False a a 3 a a log (variant) 3590 0.351 0.589 True True 2.640 1.784 False NaN a a 3 a b linear 3590 0.351 0.589 True True 1.437 1.784 False False a a 3 a b log (variant) 3590 0.351 0.589 True True 1.960 1.784 False NaN a a 3 b a linear 3590 0.351 0.665 True True 1.522 1.947 False False a a 3 b a log (variant) 3590 0.351 0.665 True True 1.900 1.947 False NaN a a 3 b b linear 3590 0.351 0.665 True True 2.025 1.947 False False a a 3 b b log (variant) 3590 0.351 0.665 True True 2.571 1.947 False NaN a a 5 a a linear 3590 0.347 0.614 True True 2.034 1.904 False False a a 5 a a log (variant) 3590 0.347 0.614 True False 3.395 1.904 False NaN a a 5 a b linear 3590 0.347 0.614 True True 1.648 1.904 False False a a 5 a b log (variant) 3590 0.347 0.614 True True 2.396 1.904 False NaN a a 5 b a linear 3590 0.347 0.665 True True 1.907 1.980 False False a a 5 b a log (variant) 3590 0.347 0.665 True True 2.027 1.980 False NaN a a 5 b b linear 3590 0.347 0.665 True True 2.091 1.980 False False a a 5 b b log (variant) 3590 0.347 0.665 True True 2.670 1.980 False NaN a b 1 a a linear 6511 0.327 0.445 True True 0.644 1.730 False False a b 1 a a log (variant) 6511 0.327 0.445 True True 1.308 1.730 False NaN a b 1 a b linear 6511 0.327 0.445 True True 1.194 1.730 False False a b 1 a b log (variant) 6511 0.327 0.445 True True 1.367 1.730 False NaN a b 1 b a linear 6511 0.327 0.587 True True 0.695 2.146 False False a b 1 b a log (variant) 6511 0.327 0.587 True True 1.694 2.146 False NaN a b 1 b b linear 6511 0.327 0.587 True True 2.145 2.146 False False a b 1 b b log (variant) 6511 0.327 0.587 True True 1.597 2.146 False NaN a b 3 a a linear 6511 0.311 0.497 True True 0.965 1.919 False False a b 3 a a log (variant) 6511 0.311 0.497 True True 1.719 1.919 False NaN a b 3 a b linear 6511 0.311 0.497 True True 1.372 1.919 False False a b 3 a b log (variant) 6511 0.311 0.497 True True 1.653 1.919 False NaN a b 3 b a linear 6511 0.311 0.588 True True 0.956 2.232 False False a b 3 b a log (variant) 6511 0.311 0.588 True True 1.793 2.232 False NaN a b 3 b b linear 6511 0.311 0.588 True True 2.210 2.232 False False a b 3 b b log (variant) 6511 0.311 0.588 True True 1.638 2.232 False NaN a b 5 a a linear 6511 0.306 0.525 True True 1.484 2.061 False False a b 5 a a log (variant) 6511 0.306 0.525 True True 2.223 2.061 False NaN a b 5 a b linear 6511 0.306 0.525 True True 1.524 2.061 False False a b 5 a b log (variant) 6511 0.306 0.525 True True 1.941 2.061 False NaN a b 5 b a linear 6511 0.306 0.588 True True 1.392 2.264 False False a b 5 b a log (variant) 6511 0.306 0.588 True True 1.949 2.264 False NaN a b 5 b b linear 6511 0.306 0.588 True True 2.259 2.264 False False a b 5 b b log (variant) 6511 0.306 0.588 True True 1.688 2.264 False NaN b a 1 a a linear 2327 0.348 0.474 True True 0.038 1.409 False False b a 1 a a log (variant) 2327 0.348 0.474 True True 1.265 1.409 False NaN b a 1 a b linear 2327 0.348 0.474 True True 1.168 1.409 False False b a 1 a b log (variant) 2327 0.348 0.474 True True 1.027 1.409 False NaN b a 1 b a linear 2327 0.348 0.580 True True 0.048 1.702 False False b a 1 b a log (variant) 2327 0.348 0.580 True True 2.185 1.702 False NaN b a 1 b b linear 2327 0.348 0.580 True True 1.820 1.702 True True b a 1 b b log (variant) 2327 0.348 0.580 True True 1.098 1.702 False NaN b a 3 a a linear 2327 0.321 0.540 True True 0.189 1.754 False False b a 3 a a log (variant) 2327 0.321 0.540 True True 2.176 1.754 False NaN b a 3 a b linear 2327 0.321 0.540 True True 1.551 1.754 False False b a 3 a b log (variant) 2327 0.321 0.540 True True 2.299 1.754 False NaN b a 3 b a linear 2327 0.321 0.582 True True 0.180 1.801 False False b a 3 b a log (variant) 2327 0.321 0.582 True True 2.761 1.801 False NaN b a 3 b b linear 2327 0.321 0.582 True True 1.952 1.801 False False b a 3 b b log (variant) 2327 0.321 0.582 True True 2.587 1.801 False NaN b a 5 a a linear 2327 0.317 0.565 True True 0.189 1.861 False False b a 5 a a log (variant) 2327 0.317 0.565 True True 2.818 1.861 False NaN b a 5 a b linear 2327 0.317 0.565 True True 1.750 1.861 False False b a 5 a b log (variant) 2327 0.317 0.565 True True 2.131 1.861 False NaN b a 5 b a linear 2327 0.317 0.582 True True 0.180 1.834 False False b a 5 b a log (variant) 2327 0.317 0.582 True True 2.934 1.834 False NaN b a 5 b b linear 2327 0.317 0.582 True True 2.029 1.834 False False b a 5 b b log (variant) 2327 0.317 0.582 True True 2.740 1.834 False NaN b b 1 a a linear 4208 0.294 0.389 True True 0.038 1.656 False False b b 1 a a log (variant) 4208 0.294 0.389 True True 2.068 1.656 False NaN b b 1 a b linear 4208 0.294 0.389 True True 1.131 1.656 False False b b 1 a b log (variant) 4208 0.294 0.389 True True 1.530 1.656 False NaN b b 1 b a linear 4208 0.294 0.501 True True 0.048 2.052 False False b b 1 b a log (variant) 4208 0.294 0.501 True True 2.185 2.052 False NaN b b 1 b b linear 4208 0.294 0.501 True True 1.948 2.052 False False b b 1 b b log (variant) 4208 0.294 0.501 True True 1.493 2.052 False NaN b b 3 a a linear 4208 0.277 0.449 True True 0.189 1.941 False False b b 3 a a log (variant) 4208 0.277 0.449 True True 1.432 1.941 False NaN b b 3 a b linear 4208 0.277 0.449 True True 1.474 1.941 False False b b 3 a b log (variant) 4208 0.277 0.449 True True 1.286 1.941 False NaN b b 3 b a linear 4208 0.277 0.502 True True 0.180 2.146 False False b b 3 b a log (variant) 4208 0.277 0.502 True True 2.761 2.146 False NaN b b 3 b b linear 4208 0.277 0.502 True True 2.055 2.146 False False b b 3 b b log (variant) 4208 0.277 0.502 True True 1.592 2.146 False NaN b b 5 a a linear 4208 0.273 0.479 True True 0.189 2.105 False False b b 5 a a log (variant) 4208 0.273 0.479 True True 1.937 2.105 False NaN b b 5 a b linear 4208 0.273 0.479 True True 1.623 2.105 False False b b 5 a b log (variant) 4208 0.273 0.479 True True 1.251 2.105 False NaN b b 5 b a linear 4208 0.273 0.503 True True 0.180 2.178 False False b b 5 b a log (variant) 4208 0.273 0.503 True True 2.934 2.178 False NaN b b 5 b b linear 4208 0.273 0.503 True True 2.106 2.178 False False b b 5 b b log (variant) 4208 0.273 0.503 True True 1.297 2.178 False NaN Combinations with the similarity ratios within 10%: P1 P4 P6 P5 P7 bins n med_ratio_ell med_ratio_ann under_a under_b ann/ell_peak ann/ell_median similar match b a 1 b b linear 2327 0.348 0.58 True True 1.82 1.702 True True Ranges on linear bins: annulus/ellipse peak 0.038 to 2.259 ; median 1.409 to 2.264 Linear bins, as in fig5b's criterion: combinations with both underestimating (a or b): 48 of 48 ; with the similarity ratios within 10%: 1
Result: reproduced. Both the ellipse fit and the annulus underestimate the segmentation in all 48 combinations on linear bins. The annulus-over-ellipse ratios are within 10% of the paper's 1.87 and 1.61 in one of the 48: P1 (b), P4 (a), erosion 1, P5 (b), outliers out, with 1.82 for the peaks and 1.702 for the medians (median per-cross-section ratio 0.348 for the ellipse and 0.58 for the annulus). The criterion's "under one reading" counts it, but the outcome rests on that single combination, under the P4 reading that leaves 4 segments without a cross section. Over all 48 the ratio of the medians is 1.409 to 2.264 and that of the peaks 0.038 to 2.259.
spline-resistance-five-times¶
The spline fit's resistance is considerably larger than the segmentation's, its distribution peak about 5 times larger, despite a larger area ratio, because the spline PVS narrows rapidly away from the vessel.
Ratio of the spline peak to the segmentation peak, with the bins of the claims it depends on: reading (a) 5 within 0.5, (b) 6.2 within 0.62; and the median spline resistance above the median segmented one.
In [18]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k, p7 in [(k, p7) for k in P6 for p7 in P7]:
M = models(d, "a", k).dropna(subset=["segmentation", "spline"])
for sc in ("linear", "log"):
ratio = rpeak(M.spline, p7, sc) / rpeak(M.segmentation, p7, sc)
med_ok = M.spline.median() > M.segmentation.median()
rows.append({"P1": p1, "P4": p4, "P6": k, "P7": p7, "bins": sc, "n": len(M), "peak_ratio": round(ratio, 2),
"reading_a": within(ratio, 5, 0.5), "reading_b": within(ratio, 6.2, 0.62),
"median_spline": f"{M.spline.median():.3g}", "median_seg": f"{M.segmentation.median():.3g}",
"match": (within(ratio, 5, 0.5) or within(ratio, 6.2, 0.62)) and med_ok})
s5 = pd.DataFrame(rows)
print(s5.to_string(index=False))
print("\nCombinations matching:", int(s5.match.sum()), "of", len(s5), "; on linear bins:", int(s5[s5.bins == "linear"].match.sum()))
d = sel("a", "b")
print("\nSpline over segmentation, per cross section (P1 a, P4 b): median", round((d.R_spl_e3 / d.R_seg).median(), 2),
"; share above 1:", round(100 * ((d.R_spl_e3 / d.R_seg) > 1).mean(), 1), "%")
print("Spline area ratio above the segmentation's (P5 a, fits over pi R1^2):",
round(100 * ((d.spl_area_e3 / (np.pi * d.R1_e3 ** 2)) > ratio_seg(d, "a")).mean(), 1), "% of cross sections")
P1 P4 P6 P7 bins n peak_ratio reading_a reading_b median_spline median_seg match a a 1 a linear 3587 0.00 False False 1.15e+16 5.93e+15 False a a 1 a log 3587 6.10 False True 1.15e+16 5.93e+15 True a a 1 b linear 3587 3.61 False False 1.15e+16 5.93e+15 False a a 1 b log 3587 4.67 True False 1.15e+16 5.93e+15 True a a 3 a linear 3587 0.00 False False 1.32e+16 5.93e+15 False a a 3 a log 3587 7.22 False False 1.32e+16 5.93e+15 False a a 3 b linear 3587 4.29 False False 1.32e+16 5.93e+15 False a a 3 b log 3587 5.79 False True 1.32e+16 5.93e+15 True a a 5 a linear 3587 0.00 False False 1.38e+16 5.93e+15 False a a 5 a log 3587 5.49 True False 1.38e+16 5.93e+15 True a a 5 b linear 3587 4.85 True False 1.38e+16 5.93e+15 True a a 5 b log 3587 4.69 True False 1.38e+16 5.93e+15 True a b 1 a linear 6508 0.02 False False 9.39e+15 5.28e+15 False a b 1 a log 6508 2.13 False False 9.39e+15 5.28e+15 False a b 1 b linear 6508 5.57 False False 9.39e+15 5.28e+15 False a b 1 b log 6508 3.18 False False 9.39e+15 5.28e+15 False a b 3 a linear 6508 0.04 False False 1.07e+16 5.28e+15 False a b 3 a log 6508 3.17 False False 1.07e+16 5.28e+15 False a b 3 b linear 6508 6.59 False True 1.07e+16 5.28e+15 True a b 3 b log 6508 3.96 False False 1.07e+16 5.28e+15 False a b 5 a linear 6508 0.04 False False 1.13e+16 5.28e+15 False a b 5 a log 6508 3.51 False False 1.13e+16 5.28e+15 False a b 5 b linear 6508 2.46 False False 1.13e+16 5.28e+15 False a b 5 b log 6508 3.13 False False 1.13e+16 5.28e+15 False b a 1 a linear 2327 0.00 False False 1.44e+16 6.87e+15 False b a 1 a log 2327 5.66 False True 1.44e+16 6.87e+15 True b a 1 b linear 2327 5.25 True False 1.44e+16 6.87e+15 True b a 1 b log 2327 4.40 False False 1.44e+16 6.87e+15 False b a 3 a linear 2327 0.00 False False 1.73e+16 6.87e+15 False b a 3 a log 2327 7.02 False False 1.73e+16 6.87e+15 False b a 3 b linear 2327 6.57 False True 1.73e+16 6.87e+15 True b a 3 b log 2327 5.50 False False 1.73e+16 6.87e+15 False b a 5 a linear 2327 0.00 False False 1.75e+16 6.87e+15 False b a 5 a log 2327 5.63 False True 1.75e+16 6.87e+15 True b a 5 b linear 2327 7.21 False False 1.75e+16 6.87e+15 False b a 5 b log 2327 6.61 False True 1.75e+16 6.87e+15 True b b 1 a linear 4208 0.26 False False 1.17e+16 6.1e+15 False b b 1 a log 4208 5.22 True False 1.17e+16 6.1e+15 True b b 1 b linear 4208 6.93 False False 1.17e+16 6.1e+15 False b b 1 b log 4208 4.23 False False 1.17e+16 6.1e+15 False b b 3 a linear 4208 0.57 False False 1.38e+16 6.1e+15 False b b 3 a log 4208 4.47 False False 1.38e+16 6.1e+15 False b b 3 b linear 4208 3.02 False False 1.38e+16 6.1e+15 False b b 3 b log 4208 4.15 False False 1.38e+16 6.1e+15 False b b 5 a linear 4208 0.89 False False 1.44e+16 6.1e+15 False b b 5 a log 4208 6.24 False True 1.44e+16 6.1e+15 True b b 5 b linear 4208 3.29 False False 1.44e+16 6.1e+15 False b b 5 b log 4208 4.93 True False 1.44e+16 6.1e+15 True Combinations matching: 15 of 48 ; on linear bins: 4 Spline over segmentation, per cross section (P1 a, P4 b): median 2.15 ; share above 1: 74.7 % Spline area ratio above the segmentation's (P5 a, fits over pi R1^2): 43.3 % of cross sections
Result: reproduced. On linear bins the ratio of the spline peak to the segmentation peak is within tolerance in 4 of 24 combinations: 4.85 (P1 a, P4 a, erosion 5) and 5.25 (P1 b, P4 a, erosion
- within 0.5 of the text's 5, and 6.59 (P1 a, P4 b, erosion 3) and 6.57 (P1 b, P4 a, erosion 3) within 0.62 of Table 2's 6.2, all with outliers out. On logarithmic bins, which fig6b's criterion allows for both peaks, 11 more combinations are within. The median spline resistance is above the segmentation's in every combination (1.07e16 against 5.28e15 with P1 a, P4 b, erosion 3), and the spline is the larger in 74.7% of cross sections. The spline's area ratio exceeds the segmentation's in only 43.3% of cross sections, so the "larger area ratio" of the statement does not hold here.
fig7-velocity-profiles¶
In a representative cross section, the polynomial fit's velocity profile resembles the segmentation's, the ellipse fit gives much higher velocities and the spline fit much lower; these trends are generally consistent across cross sections.
Velocity at 500 Pa/m. Readings: (a) maximum velocity, (b) mean velocity. Holds when in more than half the cross sections the ellipse is faster than the segmentation, the spline slower, and the polynomial closer to the segmentation than either.
In [19]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k, (rd, pre) in [(k, r) for k in P6 for r in (("a", "umax"), ("b", "umean"))]:
u = pd.DataFrame({m: d[f"{pre}_{s}"] for m, s in (("seg", "seg"), ("ell", f"ell_e{k}"), ("spl", f"spl_e{k}"), ("pol", f"pol_e{k}"))}).dropna()
fe = 100 * (u.ell > u.seg).mean(); ss = 100 * (u.spl < u.seg).mean()
dp = (u.pol - u.seg).abs()
pc = 100 * ((dp < (u.ell - u.seg).abs()) & (dp < (u.spl - u.seg).abs())).mean()
rows.append({"P1": p1, "P4": p4, "P6": k, "reading": rd, "n": len(u), "ellipse_faster_pct": round(fe, 1),
"spline_slower_pct": round(ss, 1), "polynomial_closest_pct": round(pc, 1),
"median_seg_um_s": round(1e6 * u.seg.median(), 2), "median_ell": round(1e6 * u.ell.median(), 2),
"median_spl": round(1e6 * u.spl.median(), 2), "median_pol": round(1e6 * u.pol.median(), 2),
"match": fe > 50 and ss > 50 and pc > 50})
f7 = pd.DataFrame(rows)
print(f7.to_string(index=False))
P1 P4 P6 reading n ellipse_faster_pct spline_slower_pct polynomial_closest_pct median_seg_um_s median_ell median_spl median_pol match a a 1 a 3306 100.0 89.8 76.6 82.83 132.76 43.49 81.31 True a a 1 b 3306 100.0 93.6 74.5 26.57 47.57 10.33 28.69 True a a 3 a 3309 100.0 91.9 78.4 82.75 136.53 40.78 81.69 True a a 3 b 3309 100.0 95.5 75.6 26.47 48.34 9.62 29.15 True a a 5 a 3317 100.0 92.7 78.9 83.07 138.05 39.76 82.25 True a a 5 b 3317 100.0 96.4 76.7 26.64 48.82 9.35 29.24 True a b 1 a 5778 99.9 85.8 72.1 80.03 168.61 47.14 90.73 True a b 1 b 5778 100.0 91.6 71.5 24.60 65.34 11.26 31.73 True a b 3 a 5786 100.0 88.1 73.5 79.88 171.77 44.12 91.49 True a b 3 b 5786 100.0 93.2 72.5 24.53 66.56 10.48 31.98 True a b 5 a 5796 100.0 89.2 74.2 80.08 172.92 43.12 91.75 True a b 5 b 5796 100.0 94.0 73.2 24.60 66.85 10.22 32.17 True b a 1 a 2135 100.0 90.8 77.3 77.63 127.32 39.11 75.10 True b a 1 b 2135 100.0 95.5 77.0 23.33 44.46 8.35 25.33 True b a 3 a 2135 100.0 93.7 79.7 77.88 130.74 35.61 76.77 True b a 3 b 2135 100.0 97.0 77.9 23.42 45.84 7.51 25.94 True b a 5 a 2133 100.0 94.5 81.0 77.88 131.28 35.33 77.35 True b a 5 b 2133 100.0 97.6 78.2 23.42 45.70 7.41 26.03 True b b 1 a 3716 99.9 85.9 72.6 73.72 157.64 42.62 81.30 True b b 1 b 3716 99.9 92.1 72.9 19.76 58.34 9.17 26.45 True b b 3 a 3729 100.0 88.9 74.0 73.74 162.11 39.29 82.88 True b b 3 b 3729 100.0 93.3 73.5 19.77 59.80 8.38 26.93 True b b 5 a 3732 100.0 90.0 74.9 73.72 162.26 38.52 83.25 True b b 5 b 3732 100.0 93.8 74.0 19.77 60.17 8.17 27.02 True
Result: reproduced. In every combination of P1, P4 and P6 and under both readings the ellipse fit is faster than the segmentation in 99.9 to 100% of cross sections, the spline slower in 85.8 to 97.6%, and the polynomial fit closer to the segmentation than either in 71.5 to 81.0%. The median maximum velocity with P1 (a), P4 (b) and erosion 3 is 79.88 um/s for the segmentation, 171.77 for the ellipse, 44.12 for the spline and 91.49 for the polynomial.
vessel-radius-peak¶
The fitted vessel radius distribution peaks at 28 um.
Peak of R1 under P1, P4, P6 and P7; within 2.8 um of 28 um.
In [20]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
r1 = d[f"R1_e{k}"].dropna()
for p7 in P7:
pk = pvs.peak(r1, p7)
rows.append({"P1": p1, "P4": p4, "P6": k, "P7": p7, "n": len(r1), "peak_um": round(pk, 2),
"median_um": round(r1.median(), 2), "match": within(pk, 28, 2.8)})
vr = pd.DataFrame(rows)
print(vr.to_string(index=False))
d = sel("a", "b")
print("\nPer location median R1 (P1 a, P4 b, erosion 3):", d.groupby("location").R1_e3.median().round(1).to_dict())
print("R1 range:", d.R1_e3.min().round(1), "to", d.R1_e3.max().round(1), "um")
P1 P4 P6 P7 n peak_um median_um match
a a 1 a 3590 29.50 29.83 True
a a 1 b 3590 30.12 29.83 True
a a 3 a 3590 31.96 28.92 False
a a 3 b 3590 31.04 28.92 False
a a 5 a 3590 30.18 28.59 True
a a 5 b 3590 30.90 28.59 False
a b 1 a 6511 29.50 31.26 True
a b 1 b 6511 32.30 31.26 False
a b 3 a 6511 31.96 30.31 False
a b 3 b 6511 32.20 30.31 False
a b 5 a 6511 30.18 29.94 True
a b 5 b 6511 32.11 29.94 False
b a 1 a 2328 30.39 29.55 True
b a 1 b 2328 31.94 29.55 False
b a 3 a 2328 28.75 28.39 True
b a 3 b 2328 30.71 28.39 True
b a 5 a 2328 30.53 28.31 True
b a 5 b 2328 29.20 28.31 True
b b 1 a 4213 32.34 30.95 False
b b 1 b 4213 33.03 30.95 False
b b 3 a 4213 31.79 29.83 False
b b 3 b 4213 30.51 29.83 True
b b 5 a 4213 29.31 29.53 True
b b 5 b 4213 30.21 29.53 True
Per location median R1 (P1 a, P4 b, erosion 3): {'A': 32.9, 'B': 31.8, 'C': 17.8, 'D': 18.2, 'E': 40.9, 'F': 34.1, 'G': 17.9, 'H': 39.9, 'I': 30.5, 'J': 26.8, 'K': 30.6, 'L': 30.4, 'M': 17.8, 'N': 38.5}
R1 range: 2.4 to 109.9 um
Result: reproduced. The R1 peak is within 2.8 um of 28 um in 13 of 24 combinations, for example 29.50 um with P1 (a), P4 (a), erosion 1 and all values, and 28.75 um with P1 (b), P4 (a), erosion 3; over all combinations the peak is 28.75 to 33.03 um and the median 28.31 to 31.26 um. Per location the median R1 (P1 a, P4 b, erosion 3) runs from 17.8 (C and M) to 40.9 um (E).
fig8-ellipse-parameters¶
The ellipse fit shows an elongated, eccentric PVS with the vessel slightly below its center; R1 peaks at 28 um and R3 at 35 um, R3/R1 at 0.80, R2/R1 at 2.89, R3/R2 at 0.398, the vertical offset Ov at 0.24, and the horizontal offset Oh is centered near zero with a much wider distribution than Ov.
R2 and R3 are the semi-axes of the ellipse fit to the PVS and vessel together, holes filled (P9, as amended after the question on this claim). Offsets in units of (a) R1, (b) R3, (c) um; Ov positive when the vessel is below the ellipse center. Oh conditions: median within one Ov peak of zero, IQR of Oh at least twice that of Ov.
In [21]:
EXP_8 = [("R3", 35, 3.5), ("R3/R1", 0.8, 0.08), ("R2/R1", 2.89, 0.289), ("R3/R2", 0.398, 0.0398), ("Ov", 0.24, 0.024)]
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
r1, r2, r3 = d[f"R1_e{k}"], d[f"R2_e{k}"], d[f"R3_e{k}"]
for unit in "abc":
den = r1 if unit == "a" else (r3 if unit == "b" else 1.0)
ov, oh = (d[f"Ov_e{k}"] / den).dropna(), (d[f"Oh_e{k}"] / den).dropna()
for p7 in P7:
v = {"R3": pvs.peak(r3.dropna(), p7), "R3/R1": pvs.peak((r3 / r1).dropna(), p7), "R2/R1": pvs.peak((r2 / r1).dropna(), p7),
"R3/R2": pvs.peak((r3 / r2).dropna(), p7), "Ov": pvs.peak(ov, p7)}
oh_c = abs(oh.median()) <= abs(v["Ov"])
oh_w = (oh.quantile(0.75) - oh.quantile(0.25)) >= 2 * (ov.quantile(0.75) - ov.quantile(0.25))
okv = {n: within(v[n], e, t) for n, e, t in EXP_8}
rows.append({"P1": p1, "P4": p4, "P6": k, "units": unit, "P7": p7, **{n: round(x, 3) for n, x in v.items()},
"Oh_median": round(oh.median(), 3), "Oh_IQR": round(oh.quantile(0.75) - oh.quantile(0.25), 3),
"Ov_IQR": round(ov.quantile(0.75) - ov.quantile(0.25), 3), "Oh_centered": oh_c, "Oh_wider": oh_w,
"n_within": sum(okv.values()), "match": all(okv.values()) and oh_c and oh_w,
"ratios_within": okv["R3/R1"] and okv["R2/R1"] and okv["R3/R2"]})
f8 = pd.DataFrame(rows)
print(f8.to_string(index=False))
print("\nCombinations matching:", int(f8.match.sum()), "; with the three ratios within tolerance:", int(f8.ratios_within.sum()), "of", len(f8))
print("Ranges over all combinations:", {n: (f8[n].min(), f8[n].max()) for n, _, _ in EXP_8},
"; within tolerance in:", {n: int(sum(within(v, e, t) for v in f8[n])) for n, e, t in EXP_8})
print("\nThe 8 combinations with the most values within tolerance (compact):")
print(f8.sort_values("n_within", ascending=False).head(8)[["P1", "P4", "P6", "units", "P7", "R3", "R3/R1", "R2/R1", "R3/R2", "Ov", "Oh_centered", "Oh_wider", "n_within"]].to_string(index=False))
d = sel("a", "b")
print("\nVariant, not used by the criterion since the answer: ellipse of the PVS alone (P1 a, P4 b, P7 b): R3 peak",
round(pvs.peak(d.R3p.dropna(), "b"), 1), "um; R3/R1", round(pvs.peak((d.R3p / d.R1_e3).dropna(), "b"), 3),
"; R2/R1", round(pvs.peak((d.R2p / d.R1_e3).dropna(), "b"), 3), "; R3/R2", round(pvs.peak((d.R3p / d.R2p).dropna(), "b"), 3))
P1 P4 P6 units P7 R3 R3/R1 R2/R1 R3/R2 Ov Oh_median Oh_IQR Ov_IQR Oh_centered Oh_wider n_within match ratios_within
a a 1 a a 30.096 1.037 2.735 0.394 0.211 -0.292 0.739 0.271 False True 2 False False
a a 1 a b 30.096 0.970 2.770 0.369 0.206 -0.292 0.739 0.271 False True 2 False False
a a 1 b a 30.096 1.037 2.735 0.394 0.224 -0.262 0.751 0.221 False True 3 False False
a a 1 b b 30.096 0.970 2.770 0.369 0.212 -0.262 0.751 0.221 False True 2 False False
a a 1 c a 30.096 1.037 2.735 0.394 6.323 -8.466 22.353 7.947 False True 2 False False
a a 1 c b 30.096 0.970 2.770 0.369 5.655 -8.466 22.353 7.947 False True 2 False False
a a 3 a a 30.096 1.078 2.788 0.394 0.223 -0.305 0.752 0.276 False True 3 False False
a a 3 a b 30.096 0.996 2.843 0.369 0.207 -0.305 0.752 0.276 False True 2 False False
a a 3 b a 30.096 1.078 2.788 0.394 0.153 -0.262 0.749 0.225 False True 2 False False
a a 3 b b 30.096 0.996 2.843 0.369 0.198 -0.262 0.749 0.225 False True 2 False False
a a 3 c a 30.096 1.078 2.788 0.394 6.579 -8.426 22.020 8.345 False True 2 False False
a a 3 c b 30.096 0.996 2.843 0.369 5.785 -8.426 22.020 8.345 False True 2 False False
a a 5 a a 30.096 1.107 2.813 0.394 0.224 -0.306 0.769 0.289 False True 3 False False
a a 5 a b 30.096 0.965 2.876 0.369 0.164 -0.306 0.769 0.289 False True 2 False False
a a 5 b a 30.096 1.107 2.813 0.394 0.185 -0.269 0.748 0.227 False True 2 False False
a a 5 b b 30.096 0.965 2.876 0.369 0.177 -0.269 0.748 0.227 False True 2 False False
a a 5 c a 30.096 1.107 2.813 0.394 5.477 -8.296 22.165 8.172 False True 2 False False
a a 5 c b 30.096 0.965 2.876 0.369 5.077 -8.296 22.165 8.172 False True 2 False False
a b 1 a a 30.096 1.037 2.903 0.376 0.092 -0.205 0.899 0.306 False True 2 False False
a b 1 a b 27.918 0.983 2.983 0.377 0.192 -0.205 0.899 0.306 False True 2 False False
a b 1 b a 30.096 1.037 2.903 0.376 0.231 -0.184 0.859 0.279 True True 3 False False
a b 1 b b 27.918 0.983 2.983 0.377 0.189 -0.184 0.859 0.279 True True 2 False False
a b 1 c a 30.096 1.037 2.903 0.376 7.520 -6.015 25.897 9.766 True True 2 False False
a b 1 c b 27.918 0.983 2.983 0.377 5.335 -6.015 25.897 9.766 False True 2 False False
a b 3 a a 30.096 1.078 3.049 0.376 0.423 -0.213 0.919 0.316 True True 2 False False
a b 3 a b 27.918 1.008 3.066 0.377 0.184 -0.213 0.919 0.316 False True 2 False False
a b 3 b a 30.096 1.078 3.049 0.376 0.154 -0.187 0.851 0.284 False True 2 False False
a b 3 b b 27.918 1.008 3.066 0.377 0.166 -0.187 0.851 0.284 False True 2 False False
a b 3 c a 30.096 1.078 3.049 0.376 4.583 -6.022 25.792 9.493 False True 2 False False
a b 3 c b 27.918 1.008 3.066 0.377 4.954 -6.022 25.792 9.493 False True 2 False False
a b 5 a a 30.096 1.107 3.071 0.376 0.407 -0.218 0.928 0.322 True True 2 False False
a b 5 a b 27.918 0.971 2.875 0.377 0.178 -0.218 0.928 0.322 False True 2 False False
a b 5 b a 30.096 1.107 3.071 0.376 0.160 -0.186 0.843 0.282 False True 2 False False
a b 5 b b 27.918 0.971 2.875 0.377 0.169 -0.186 0.843 0.282 False True 2 False False
a b 5 c a 30.096 1.107 3.071 0.376 5.242 -6.062 25.719 9.628 False True 2 False False
a b 5 c b 27.918 0.971 2.875 0.377 4.938 -6.062 25.719 9.628 False True 2 False False
b a 1 a a 27.351 0.989 2.957 0.431 0.171 -0.293 0.740 0.274 False True 2 False False
b a 1 a b 27.581 0.967 2.840 0.437 0.208 -0.293 0.740 0.274 False True 2 False False
b a 1 b a 27.351 0.989 2.957 0.431 0.195 -0.270 0.757 0.226 False True 2 False False
b a 1 b b 27.581 0.967 2.840 0.437 0.202 -0.270 0.757 0.226 False True 2 False False
b a 1 c a 27.351 0.989 2.957 0.431 4.314 -8.735 22.502 7.707 False True 2 False False
b a 1 c b 27.581 0.967 2.840 0.437 5.594 -8.735 22.502 7.707 False True 2 False False
b a 3 a a 27.351 1.024 3.049 0.431 0.190 -0.317 0.767 0.287 False True 2 False False
b a 3 a b 27.581 0.966 3.073 0.437 0.200 -0.317 0.767 0.287 False True 2 False False
b a 3 b a 27.351 1.024 3.049 0.431 0.176 -0.281 0.744 0.230 False True 2 False False
b a 3 b b 27.581 0.966 3.073 0.437 0.187 -0.281 0.744 0.230 False True 2 False False
b a 3 c a 27.351 1.024 3.049 0.431 4.370 -9.028 22.132 8.130 False True 2 False False
b a 3 c b 27.581 0.966 3.073 0.437 5.222 -9.028 22.132 8.130 False True 2 False False
b a 5 a a 27.351 0.977 2.827 0.431 0.216 -0.332 0.777 0.301 False True 2 False False
b a 5 a b 27.581 0.997 2.636 0.437 0.196 -0.332 0.777 0.301 False True 2 False False
b a 5 b a 27.351 0.977 2.827 0.431 0.205 -0.296 0.742 0.232 False True 2 False False
b a 5 b b 27.581 0.997 2.636 0.437 0.180 -0.296 0.742 0.232 False True 2 False False
b a 5 c a 27.351 0.977 2.827 0.431 6.086 -9.029 22.249 8.181 False True 2 False False
b a 5 c b 27.581 0.997 2.636 0.437 5.059 -9.029 22.249 8.181 False True 2 False False
b b 1 a a 27.351 1.003 3.033 0.357 0.130 -0.226 0.881 0.306 False True 1 False False
b b 1 a b 28.055 0.991 3.031 0.347 0.207 -0.226 0.881 0.306 False True 1 False False
b b 1 b a 27.351 1.003 3.033 0.357 0.355 -0.205 0.864 0.276 True True 1 False False
b b 1 b b 28.055 0.991 3.031 0.347 0.195 -0.205 0.864 0.276 False True 1 False False
b b 1 c a 27.351 1.003 3.033 0.357 7.238 -6.430 24.948 9.586 True True 1 False False
b b 1 c b 28.055 0.991 3.031 0.347 5.526 -6.430 24.948 9.586 False True 1 False False
b b 3 a a 27.351 1.000 3.175 0.357 0.273 -0.233 0.922 0.319 True True 1 False False
b b 3 a b 28.055 0.970 3.122 0.347 0.197 -0.233 0.922 0.319 False True 1 False False
b b 3 b a 27.351 1.000 3.175 0.357 0.211 -0.209 0.856 0.280 True True 1 False False
b b 3 b b 28.055 0.970 3.122 0.347 0.185 -0.209 0.856 0.280 False True 1 False False
b b 3 c a 27.351 1.000 3.175 0.357 4.876 -6.618 24.676 9.427 False True 1 False False
b b 3 c b 28.055 0.970 3.122 0.347 4.946 -6.618 24.676 9.427 False True 1 False False
b b 5 a a 27.351 1.016 3.181 0.357 0.272 -0.238 0.943 0.327 True True 0 False False
b b 5 a b 28.055 0.986 2.950 0.347 0.190 -0.238 0.943 0.327 False True 1 False False
b b 5 b a 27.351 1.016 3.181 0.357 0.384 -0.212 0.849 0.282 True True 0 False False
b b 5 b b 28.055 0.986 2.950 0.347 0.166 -0.212 0.849 0.282 False True 1 False False
b b 5 c a 27.351 1.016 3.181 0.357 3.421 -6.574 24.920 9.321 False True 0 False False
b b 5 c b 28.055 0.986 2.950 0.347 4.863 -6.574 24.920 9.321 False True 1 False False
Combinations matching: 0 ; with the three ratios within tolerance: 0 of 72
Ranges over all combinations: {'R3': (np.float64(27.351), np.float64(30.096)), 'R3/R1': (np.float64(0.965), np.float64(1.107)), 'R2/R1': (np.float64(2.636), np.float64(3.181)), 'R3/R2': (np.float64(0.347), np.float64(0.437)), 'Ov': (np.float64(0.092), np.float64(7.52))} ; within tolerance in: {'R3': 0, 'R3/R1': 0, 'R2/R1': 69, 'R3/R2': 54, 'Ov': 5}
The 8 combinations with the most values within tolerance (compact):
P1 P4 P6 units P7 R3 R3/R1 R2/R1 R3/R2 Ov Oh_centered Oh_wider n_within
a a 1 b a 30.096 1.037 2.735 0.394 0.224 False True 3
a a 3 a a 30.096 1.078 2.788 0.394 0.223 False True 3
a b 1 b a 30.096 1.037 2.903 0.376 0.231 True True 3
a a 5 a a 30.096 1.107 2.813 0.394 0.224 False True 3
a a 1 b b 30.096 0.970 2.770 0.369 0.212 False True 2
a a 1 c a 30.096 1.037 2.735 0.394 6.323 False True 2
a a 1 a b 30.096 0.970 2.770 0.369 0.206 False True 2
a a 1 a a 30.096 1.037 2.735 0.394 0.211 False True 2
Variant, not used by the criterion since the answer: ellipse of the PVS alone (P1 a, P4 b, P7 b): R3 peak 19.8 um; R3/R1 0.667 ; R2/R1 3.566 ; R3/R2 0.225
Result: not reproduced. No combination matches, and none puts the three ratios within tolerance, so the partial rule does not hold either. R3 peaks at 27.35 to 30.10 um (expected 35 within 3.5) and R3/R1 at 0.965 to 1.107 (expected 0.80 within 0.08), within tolerance in no combination; R2/R1 (2.636 to 3.181) is within in 69 of 72 combinations and R3/R2 (0.347 to 0.437) in 54 of 72. Ov is within tolerance of 0.24 in 5 of 72. The IQR of Oh is at least twice that of Ov in all 72, and the median of Oh lies within one Ov peak of zero in 11. The ellipse of the PVS alone, the variant above, gives R3 19.8 um and R3/R1 0.667.
fig9c-pvs-two-lobes¶
The PVS usually consists of two mostly or entirely disconnected lobes on either side of the vessel; the polynomial fit's HC (top of PVS above top of vessel) peaks at zero, and when the lobes connect the channel is narrow, usually less than 0.5 R1.
Peak of HC/R1 within 0.1 of zero (the claims file's tolerance), and, among cross sections with HC above zero, more than half with HC below 0.5 R1.
In [22]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
hc = (d[f"HC_e{k}"] / d[f"R1_e{k}"]).dropna()
pos = hc[hc > 0]
for p7 in P7:
pk = pvs.peak(hc, p7)
rows.append({"P1": p1, "P4": p4, "P6": k, "P7": p7, "n": len(hc), "peak_HC/R1": round(pk, 3),
"median": round(hc.median(), 3), "share_HC>0": round(100 * (hc > 0).mean(), 1),
"of_those_below_0.5R1": round(100 * (pos < 0.5).mean(), 1),
"match": within(pk, 0.000689, 0.1) and (pos < 0.5).mean() > 0.5})
f9c = pd.DataFrame(rows)
print(f9c.to_string(index=False))
P1 P4 P6 P7 n peak_HC/R1 median share_HC>0 of_those_below_0.5R1 match a a 1 a 3306 0.376 0.002 50.4 63.6 False a a 1 b 3306 -0.017 0.002 50.4 63.6 True a a 3 a 3309 0.381 0.018 52.5 63.6 False a a 3 b 3309 0.008 0.018 52.5 63.6 True a a 5 a 3317 0.121 0.022 53.0 63.4 False a a 5 b 3317 -0.052 0.022 53.0 63.4 True a b 1 a 5778 -0.401 -0.014 47.1 66.9 False a b 1 b 5778 -0.065 -0.014 47.1 66.9 True a b 3 a 5786 0.471 -0.005 49.3 66.8 False a b 3 b 5786 -0.061 -0.005 49.3 66.8 True a b 5 a 5796 0.446 -0.000 50.0 66.5 False a b 5 b 5796 -0.056 -0.000 50.0 66.5 True b a 1 a 2135 -0.178 0.004 50.8 66.2 False b a 1 b 2135 -0.138 0.004 50.8 66.2 False b a 3 a 2135 -0.037 0.026 53.6 66.6 True b a 3 b 2135 -0.126 0.026 53.6 66.6 False b a 5 a 2133 0.006 0.034 54.6 66.1 True b a 5 b 2133 -0.125 0.034 54.6 66.1 False b b 1 a 3716 -0.077 -0.014 47.4 68.4 True b b 1 b 3716 -0.059 -0.014 47.4 68.4 True b b 3 a 3729 -0.064 -0.001 49.9 68.3 True b b 3 b 3729 -0.053 -0.001 49.9 68.3 True b b 5 a 3732 -0.005 0.003 50.4 67.5 True b b 5 b 3732 -0.058 0.003 50.4 67.5 True
Result: reproduced. The peak of HC/R1 is within 0.1 of zero in 14 of 24 combinations (for example -0.061 with P1 a, P4 b, erosion 3, outliers out; median -0.005), and of the cross sections with HC above zero, 63.4 to 68.4% have HC below 0.5 R1 in every combination. HC is above zero in 47.1 to 54.6% of cross sections.
fig9e-pvs-width-asymmetry¶
The PVS is about twice as wide as the vessel radius and asymmetric about the vessel; the short-side width WS/R1 peaks near 1.1, the long-side width WL/R1 near 1.9, and WS/WL near 0.5.
WS and WL are the shorter and longer of the two side widths of the polynomial fit (vessel edge to PVS end). Peaks within 0.108, 0.185 and 0.0502 of 1.08, 1.85 and 0.502.
In [23]:
def widths(d, k):
w = d[[f"WL_e{k}", f"WR_e{k}"]].values # left and right side widths
return pd.DataFrame({"WS": np.nanmin(w, 1), "WL": np.nanmax(w, 1), "R1": d[f"R1_e{k}"].values,
"left": w[:, 0], "right": w[:, 1], "location": d.location.values}).dropna()
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
w = widths(d, k)
for p7 in P7:
v = {"WS/R1": pvs.peak(w.WS / w.R1, p7), "WL/R1": pvs.peak(w.WL / w.R1, p7), "WS/WL": pvs.peak(w.WS / w.WL, p7)}
rows.append({"P1": p1, "P4": p4, "P6": k, "P7": p7, "n": len(w), **{n: round(x, 3) for n, x in v.items()},
"median_WS/R1": round((w.WS / w.R1).median(), 3), "median_WL/R1": round((w.WL / w.R1).median(), 3),
"median_WS/WL": round((w.WS / w.WL).median(), 3),
"match": within(v["WS/R1"], 1.08, 0.108) and within(v["WL/R1"], 1.85, 0.185) and within(v["WS/WL"], 0.502, 0.0502)})
f9e = pd.DataFrame(rows)
print(f9e.to_string(index=False))
P1 P4 P6 P7 n WS/R1 WL/R1 WS/WL median_WS/R1 median_WL/R1 median_WS/WL match a a 1 a 3306 1.219 1.146 0.728 1.197 2.230 0.587 False a a 1 b 3306 1.117 1.948 0.728 1.197 2.230 0.587 False a a 3 a 3309 1.000 1.333 0.678 1.258 2.339 0.595 False a a 3 b 3309 1.135 2.041 0.678 1.258 2.339 0.595 False a a 5 a 3317 1.211 1.470 0.679 1.290 2.376 0.595 False a a 5 b 3317 1.267 2.063 0.679 1.290 2.376 0.595 False a b 1 a 5778 1.497 3.309 0.629 1.369 2.854 0.564 False a b 1 b 5778 1.198 1.960 0.629 1.369 2.854 0.564 False a b 3 a 5786 1.563 3.474 0.678 1.437 2.997 0.570 False a b 3 b 5786 1.041 2.109 0.678 1.437 2.997 0.570 False a b 5 a 5796 1.582 3.605 0.679 1.464 3.055 0.573 False a b 5 b 5796 1.232 3.905 0.679 1.464 3.055 0.573 False b a 1 a 2135 1.186 1.887 0.729 1.216 2.284 0.581 False b a 1 b 2135 1.246 1.817 0.729 1.216 2.284 0.581 False b a 3 a 2135 1.220 2.116 0.679 1.300 2.397 0.587 False b a 3 b 2135 1.267 2.016 0.679 1.300 2.397 0.587 False b a 5 a 2133 1.256 2.196 0.681 1.327 2.451 0.591 False b a 5 b 2133 1.289 2.090 0.681 1.327 2.451 0.591 False b b 1 a 3716 1.558 2.028 0.631 1.366 2.893 0.564 False b b 1 b 3716 1.199 1.962 0.631 1.366 2.893 0.564 False b b 3 a 3729 1.630 2.459 0.630 1.455 3.055 0.571 False b b 3 b 3729 1.032 1.612 0.630 1.455 3.055 0.571 False b b 5 a 3732 1.624 2.476 0.631 1.484 3.136 0.573 False b b 5 b 3732 1.066 1.641 0.631 1.484 3.136 0.573 False
Result: not reproduced. No combination puts the three peaks within tolerance. WS/R1 peaks at 1.000 to 1.630 (expected 1.08), WL/R1 at 1.146 to 3.905 (expected 1.85) and WS/WL at 0.629 to 0.729 (expected 0.502). With P1 (b), P4 (a), erosion 1 and all values, WS/R1 (1.186) and WL/R1 (1.887) are within tolerance, but WS/WL is 0.729. The medians are WS/R1 1.197 to 1.484, WL/R1 2.230 to 3.136 and WS/WL 0.564 to 0.595, so the two sides are asymmetric.
left-right-no-difference¶
The PVS widths on the left and right sides of the vessel do not differ significantly.
Readings: (a) Wilcoxon signed-rank test on paired left and right widths over all cross sections; (b) two-sample Kolmogorov-Smirnov test; (c) Wilcoxon signed-rank test on the per-segment median widths (14 pairs); alpha 0.05. Left and right are fixed by the direction along the centerline.
In [24]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
w = widths(d, k)
pa = stats.wilcoxon(w.left, w.right).pvalue
pb = stats.ks_2samp(w.left, w.right).pvalue
m = w.groupby("location")[["left", "right"]].median()
pc = stats.wilcoxon(m.left, m.right).pvalue
rows.append({"P1": p1, "P4": p4, "P6": k, "n": len(w), "median_left_minus_right_um": round((w.left - w.right).median(), 2),
"p_a": f"{pa:.3g}", "p_b": f"{pb:.3g}", "segments": len(m),
"median_seg_diff_um": round((m.left - m.right).median(), 2), "p_c": f"{pc:.3g}",
"differences_found": int(pa < 0.05) + int(pb < 0.05) + int(pc < 0.05)})
lr = pd.DataFrame(rows)
print(lr.to_string(index=False))
d = sel("a", "b"); w = widths(d, 3)
print("\nPer location median left and right widths (um), P1 a, P4 b, erosion 3:")
print(w.groupby("location")[["left", "right"]].median().round(1).T.to_string())
P1 P4 P6 n median_left_minus_right_um p_a p_b segments median_seg_diff_um p_c differences_found a a 1 3306 -18.00 3.42e-187 3.99e-137 10 -16.19 0.0645 2 a a 3 3309 -17.61 2.76e-187 2.26e-138 10 -16.87 0.105 2 a a 5 3317 -17.97 8.09e-190 1.81e-142 10 -17.21 0.084 2 a b 1 5778 -9.63 5.98e-42 1.81e-45 14 -12.44 0.67 2 a b 3 5786 -9.75 1.91e-41 2.06e-47 14 -12.68 0.715 2 a b 5 5796 -9.76 8.67e-43 3.9e-49 14 -12.97 0.67 2 b a 1 2135 -17.85 3.08e-122 7.65e-94 9 -24.05 0.0547 2 b a 3 2135 -17.69 8.48e-125 2.08e-96 9 -26.19 0.0547 2 b a 5 2133 -17.96 2.5e-126 1.18e-97 9 -25.83 0.0547 2 b b 1 3716 -9.78 2.09e-30 7.26e-33 14 -12.85 0.583 2 b b 3 3729 -9.83 2.03e-31 5.14e-32 14 -13.59 0.583 2 b b 5 3732 -9.91 2.44e-32 1.26e-31 14 -13.77 0.542 2 Per location median left and right widths (um), P1 a, P4 b, erosion 3: location A B C D E F G H I J K L M N left 76.2 133.9 40.6 37.8 134.1 132.0 76.1 65.3 27.4 21.1 32.3 52.6 38.5 59.8 right 106.4 44.5 24.3 78.5 78.8 80.6 108.0 117.4 67.2 99.5 57.7 43.8 37.8 59.7
Result: blocked (methods underspecified). The readings disagree in every combination: the paired Wilcoxon test (a) and the Kolmogorov-Smirnov test (b) find a difference (p 2e-30 to 8e-190 and 1e-31 to 2e-142; median left minus right -9.6 to -18.0 um), while the Wilcoxon test on the per-segment medians (c) does not (p 0.055 to 0.715). Readings (a) and (b) treat thousands of neighbouring cross sections, which are not independent, as independent samples, so their p-values overstate the evidence; reading (c) has independent units. The criterion sends the split to methods_underspecified. Left and right differ by segment in both directions (left wider in B, C, E, F, L, M and N; right wider in A, D, G, H, I, J and K).
fig9a-heights-taper¶
The PVS is tallest next to the vessel and narrows far from it; the heights H1/R1, H2/R1 and H3/R1 peak at 1.34, 0.792 and 0.307, so H1 is about 4 times H3.
Heights of the polynomial fit over R1, (a) both sides pooled, (b) the mean of the two sides.
In [25]:
def heights(d, k, rd):
r1 = d[f"R1_e{k}"]
if rd == "a":
return pd.DataFrame({h: np.r_[d[f"{h}L_e{k}"] / r1, d[f"{h}R_e{k}"] / r1] for h in ("H1", "H2", "H3")}).dropna()
return pd.DataFrame({h: (d[f"{h}L_e{k}"] + d[f"{h}R_e{k}"]) / 2 / r1 for h in ("H1", "H2", "H3")}).dropna()
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
for rd in "ab":
h = heights(d, k, rd)
for p7 in P7:
pk = {c: pvs.peak(h[c], p7) for c in h}
md_ = h.median()
rows.append({"P1": p1, "P4": p4, "P6": k, "reading": rd, "P7": p7, "n": len(h),
**{f"pk_{c}": round(pk[c], 3) for c in pk}, **{f"md_{c}": round(md_[c], 3) for c in md_.index},
"match": within(pk["H1"], 1.34, 0.134) and within(pk["H2"], 0.792, 0.0792) and within(pk["H3"], 0.307, 0.0307),
"order_holds": pk["H1"] > pk["H2"] > pk["H3"] and md_.H1 > md_.H2 > md_.H3})
f9a = pd.DataFrame(rows)
print(f9a.to_string(index=False))
print("\nShare of side heights below zero (the line under the polynomial), P1 a, P4 b, erosion 3:",
(heights(sel("a", "b"), 3, "a") < 0).mean().round(3).to_dict())
P1 P4 P6 reading P7 n pk_H1 pk_H2 pk_H3 md_H1 md_H2 md_H3 match order_holds
a a 1 a a 6612 0.973 0.724 0.328 1.080 1.112 0.477 False False
a a 1 a b 6612 1.060 1.105 0.388 1.080 1.112 0.477 False False
a a 1 b a 3306 1.445 1.187 0.308 1.054 1.124 0.481 False False
a a 1 b b 3306 0.961 0.879 0.355 1.054 1.124 0.481 False False
a a 3 a a 6618 0.718 0.831 0.315 1.084 1.155 0.486 False False
a a 3 a b 6618 1.052 1.033 0.395 1.084 1.155 0.486 False False
a a 3 b a 3309 1.375 1.343 0.456 1.047 1.160 0.486 False False
a a 3 b b 3309 0.954 0.918 0.368 1.047 1.160 0.486 False False
a a 5 a a 6634 0.343 1.597 0.309 1.087 1.172 0.492 False False
a a 5 a b 6634 1.044 1.062 0.398 1.087 1.172 0.492 False False
a a 5 b a 3317 1.167 1.292 0.451 1.050 1.166 0.486 False False
a a 5 b b 3317 0.963 0.931 0.375 1.050 1.166 0.486 False False
a b 1 a a 11556 0.973 0.724 0.620 1.189 1.108 0.520 False True
a b 1 a b 11556 1.044 1.109 0.400 1.189 1.108 0.520 False False
a b 1 b a 5778 1.445 0.934 0.534 1.162 1.132 0.517 False True
a b 1 b b 5778 0.910 0.815 0.361 1.162 1.132 0.517 False True
a b 3 a a 11572 0.718 0.831 0.615 1.196 1.151 0.531 False False
a b 3 a b 11572 1.052 1.024 0.308 1.196 1.151 0.531 False True
a b 3 b a 5786 1.375 1.117 0.531 1.162 1.171 0.526 False False
a b 3 b b 5786 0.910 0.818 0.470 1.162 1.171 0.526 False False
a b 5 a a 11592 0.343 0.719 0.608 1.199 1.172 0.536 False False
a b 5 a b 11592 1.040 1.150 0.314 1.199 1.172 0.536 False False
a b 5 b a 5796 1.167 1.062 0.522 1.162 1.187 0.532 False False
a b 5 b b 5796 1.080 0.903 0.476 1.162 1.187 0.532 False False
b a 1 a a 4270 1.324 1.015 0.451 1.066 1.119 0.459 False False
b a 1 a b 4270 1.046 1.061 0.365 1.066 1.119 0.459 False False
b a 1 b a 2135 1.159 0.871 0.363 1.037 1.129 0.458 False False
b a 1 b b 2135 0.995 0.842 0.345 1.037 1.129 0.458 False False
b a 3 a a 4270 1.289 1.278 0.447 1.070 1.169 0.470 False False
b a 3 a b 4270 1.022 1.097 0.374 1.070 1.169 0.470 False False
b a 3 b a 2135 1.087 0.982 0.375 1.040 1.169 0.465 False False
b a 3 b b 2135 0.954 0.999 0.351 1.040 1.169 0.465 False False
b a 5 a a 4266 1.371 1.312 0.463 1.082 1.186 0.477 False False
b a 5 a b 4266 1.052 1.120 0.373 1.082 1.186 0.477 False False
b a 5 b a 2133 1.142 0.994 0.392 1.037 1.178 0.471 False False
b a 5 b b 2133 0.939 0.984 0.354 1.037 1.178 0.471 False False
b b 1 a a 7432 0.896 0.773 0.293 1.177 1.116 0.498 False True
b b 1 a b 7432 1.007 0.931 0.398 1.177 1.116 0.498 False True
b b 1 b a 3716 0.882 0.754 0.437 1.147 1.139 0.494 False True
b b 1 b b 3716 0.933 0.862 0.328 1.147 1.139 0.494 False True
b b 3 a a 7458 0.871 0.871 0.371 1.191 1.167 0.509 False True
b b 3 a b 7458 1.228 0.820 0.305 1.191 1.167 0.509 True True
b b 3 b a 3729 1.138 0.807 0.518 1.138 1.187 0.501 False False
b b 3 b b 3729 0.925 0.834 0.350 1.138 1.187 0.501 False False
b b 5 a a 7464 0.951 0.864 0.359 1.199 1.185 0.515 False True
b b 5 a b 7464 1.220 1.051 0.306 1.199 1.185 0.515 False True
b b 5 b a 3732 1.127 0.804 0.515 1.141 1.205 0.507 False False
b b 5 b b 3732 1.078 0.883 0.359 1.141 1.205 0.507 False False
Share of side heights below zero (the line under the polynomial), P1 a, P4 b, erosion 3: {'H1': 0.03, 'H2': 0.003, 'H3': 0.029}
Result: reproduced. One of the 48 combinations puts the three peaks within tolerance: P1 (b), P4 (b), erosion 3, both sides pooled and outliers out, with H1/R1 1.228, H2/R1 0.820 and H3/R1 0.305 (expected 1.34, 0.792 and 0.307). The outcome rests on that one combination; in the others the H1/R1 peak runs from 0.343 to 1.445. The order H1 above H2 above H3 holds for the peaks and the medians in 12 of 48 combinations. The medians are H1/R1 1.037 to 1.199, H2/R1 1.108 to 1.205 and H3/R1 0.458 to 0.536, so by median H1 is about 2.3 times H3, not 4, and with P4 (a) H2 is above H1.
fig9f-pvs-concave¶
The bottom edge of the PVS is most often concave; the distribution of (H2 - H3)/(H1 - H2) peaks at 0.108, below the value of 1 that separates concave from convex.
As fig9a, with the same readings: (a) the ratio per side, both sides pooled; (b) the ratio of the mean heights of the two sides. Peak within 0.0108 of 0.108 and more than half the values below 1.
In [26]:
rows = []
for p1 in TARGET:
for p4 in P4:
d = sel(p1, p4)
for k in P6:
for rd in "ab":
if rd == "a":
c = pd.Series(np.r_[(d[f"H2L_e{k}"] - d[f"H3L_e{k}"]) / (d[f"H1L_e{k}"] - d[f"H2L_e{k}"]),
(d[f"H2R_e{k}"] - d[f"H3R_e{k}"]) / (d[f"H1R_e{k}"] - d[f"H2R_e{k}"])])
else:
h = heights(d, k, "b")
c = (h.H2 - h.H3) / (h.H1 - h.H2)
c = c.replace([np.inf, -np.inf], np.nan).dropna()
for p7 in P7:
pk = pvs.peak(c, p7)
rows.append({"P1": p1, "P4": p4, "P6": k, "reading": rd, "P7": p7, "n": len(c), "peak": round(pk, 3),
"median": round(c.median(), 3), "share_below_1": round(100 * (c < 1).mean(), 1),
"range": f"{c.min():.3g} to {c.max():.3g}",
"match": within(pk, 0.108, 0.0108) and (c < 1).mean() > 0.5, "partial": (c < 1).mean() > 0.5})
f9f = pd.DataFrame(rows)
print(f9f.to_string(index=False))
print(f"\nNotebook run time {(time.time() - started) / 60:.1f} min")
P1 P4 P6 reading P7 n peak median share_below_1 range match partial a a 1 a a 6612 -278.781 0.171 71.7 -4.05e+03 to 1.47e+05 False True a a 1 a b 6612 0.200 0.171 71.7 -4.05e+03 to 1.47e+05 False True a a 1 b a 3306 -45.436 0.301 67.7 -578 to 6.52e+03 False True a a 1 b b 3306 0.789 0.301 67.7 -578 to 6.52e+03 False True a a 3 a a 6618 34.409 0.090 73.1 -5.89e+03 to 2.04e+04 False True a a 3 a b 6618 0.634 0.090 73.1 -5.89e+03 to 2.04e+04 False True a a 3 b a 3309 180.319 -0.032 70.4 -2.16e+04 to 740 False True a a 3 b b 3309 0.698 -0.032 70.4 -2.16e+04 to 740 False True a a 5 a a 6634 -188.227 0.051 73.1 -1.7e+03 to 1.04e+04 False True a a 5 a b 6634 0.612 0.051 73.1 -1.7e+03 to 1.04e+04 False True a a 5 b a 3317 -177.943 -0.279 72.1 -7.15e+03 to 1.85e+03 False True a a 5 b b 3317 0.490 -0.279 72.1 -7.15e+03 to 1.85e+03 False True a b 1 a a 11556 -278.781 0.402 67.8 -4.05e+03 to 1.47e+05 False True a b 1 a b 11556 1.162 0.402 67.8 -4.05e+03 to 1.47e+05 False True a b 1 b a 5778 -279.093 0.465 64.3 -7.8e+03 to 6.52e+03 False True a b 1 b b 5778 0.936 0.465 64.3 -7.8e+03 to 6.52e+03 False True a b 3 a a 11572 34.409 0.356 69.0 -5.89e+03 to 2.04e+04 False True a b 3 a b 11572 1.088 0.356 69.0 -5.89e+03 to 2.04e+04 False True a b 3 b a 5786 -323.825 0.371 66.2 -2.16e+04 to 4.2e+03 False True a b 3 b b 5786 0.796 0.371 66.2 -2.16e+04 to 4.2e+03 False True a b 5 a a 11592 -188.227 0.346 68.7 -1.7e+03 to 1.04e+04 False True a b 5 a b 11592 1.097 0.346 68.7 -1.7e+03 to 1.04e+04 False True a b 5 b a 5796 -177.943 0.295 67.5 -7.15e+03 to 1.85e+03 False True a b 5 b b 5796 0.711 0.295 67.5 -7.15e+03 to 1.85e+03 False True b a 1 a a 4270 8.060 0.137 71.9 -642 to 1.36e+03 False True b a 1 a b 4270 0.164 0.137 71.9 -642 to 1.36e+03 False True b a 1 b a 2135 -66.323 0.274 68.1 -689 to 7.61e+03 False True b a 1 b b 2135 0.749 0.274 68.1 -689 to 7.61e+03 False True b a 3 a a 4270 25395.203 0.034 73.0 -1.21e+03 to 1.06e+06 False True b a 3 a b 4270 0.604 0.034 73.0 -1.21e+03 to 1.06e+06 False True b a 3 b a 2135 -1016.204 -0.428 71.7 -3.23e+04 to 1.08e+04 False True b a 3 b b 2135 0.550 -0.428 71.7 -3.23e+04 to 1.08e+04 False True b a 5 a a 4266 20.040 -0.005 73.5 -1.05e+03 to 1.47e+03 False True b a 5 a b 4266 0.570 -0.005 73.5 -1.05e+03 to 1.47e+03 False True b a 5 b a 2133 64.629 -0.627 72.8 -621 to 2.42e+03 False True b a 5 b b 2133 0.530 -0.627 72.8 -621 to 2.42e+03 False True b b 1 a a 7432 317.491 0.347 69.0 -4.03e+04 to 1.36e+03 False True b b 1 a b 7432 1.090 0.347 69.0 -4.03e+04 to 1.36e+03 False True b b 1 b a 3716 147.482 0.454 63.5 -9.96e+03 to 7.61e+03 False True b b 1 b b 3716 0.954 0.454 63.5 -9.96e+03 to 7.61e+03 False True b b 3 a a 7458 24102.211 0.319 68.9 -2.54e+03 to 1.06e+06 False True b b 3 a b 7458 1.058 0.319 68.9 -2.54e+03 to 1.06e+06 False True b b 3 b a 3729 -1016.204 0.315 66.0 -3.23e+04 to 1.08e+04 False True b b 3 b b 3729 0.738 0.315 66.0 -3.23e+04 to 1.08e+04 False True b b 5 a a 7464 -11.129 0.294 68.9 -1.64e+03 to 1.47e+03 False True b b 5 a b 7464 1.076 0.294 68.9 -1.64e+03 to 1.47e+03 False True b b 5 b a 3732 47.968 0.259 67.5 -2.58e+03 to 2.42e+03 False True b b 5 b b 3732 0.685 0.259 67.5 -2.58e+03 to 2.42e+03 False True Notebook run time 27.0 min
Result: partially reproduced. More than half the values are below 1 in every combination (63.9 to 73.6%; medians -0.627 to 0.465), so the bottom edge is most often concave. The peak is within 0.0108 of 0.108 in no combination; with outliers out it is 0.164 to 1.162, and over all values the bins are stretched by extreme ratios (down to -4.03e4 and up to 1.06e6) where H1 is close to H2.
Sensitivity to location E¶
The cleaned path of E follows a side branch of the vessel network in its mask (data check), and its PVS lobes are only 29 to 52% of the slice's PVS. E is one of the paper's 14 segments, so every claim keeps it. The cell below computes each claim's table again without E and counts, for each column the outcome rests on, the combinations that meet it with and without E.
In [27]:
FULL = {k: globals()[k].copy() for k in ("cc", "f5a", "f5e", "f5bc", "arl", "f5d", "t2", "f6a", "f6b", "rpl", "f5b",
"eau", "s5", "f7", "vr", "f8", "f9c", "f9e", "lr", "f9a", "f9f")}
COLS = {"cc": ["holds", "counts_in_range"], "f5a": ["match"], "f5e": ["all_within", "median_IQR_within"],
"f5bc": ["reading_a", "reading_b"], "arl": ["match"], "f5d": ["match", "partial"], "t2": ["match", "partial"],
"f6a": ["match", "partial"], "f6b": ["match", "median_within"], "rpl": ["match"], "f5b": ["match", "partial"],
"eau": ["match", "similar"], "s5": ["match"], "f7": ["match"], "vr": ["match"], "f8": ["match", "ratios_within"],
"f9c": ["match"], "f9e": ["match"], "f9a": ["match", "order_holds"], "f9f": ["match", "partial"]}
# the source of each claim's table cell, as run above
SRC = {k: [c for c in In if f"\n{k} = pd.DataFrame(rows)" in c][-1] for k in FULL}
saved = (DATA, META, COUNTS_B, plt.show)
DATA = {p1: d[d.location != "E"] for p1, d in saved[0].items()}
META = {p1: m.drop(index="E") for p1, m in saved[1].items()}
COUNTS_B = {p1: [r for r in rs if r["location"] != "E"] for p1, rs in saved[2].items()}
plt.show = lambda *a, **k: plt.close("all")
NOE = {}
try:
for _tab in FULL: # the claim cells reuse k and other loop names
try:
with contextlib.redirect_stdout(io.StringIO()):
exec(SRC[_tab], globals())
NOE[_tab] = globals()[_tab].copy()
except Exception as err:
print("not recomputed without E:", _tab, repr(err))
finally:
DATA, META, COUNTS_B, plt.show = saved
globals().update({k: v for k, v in FULL.items()})
rows_e = []
for k, cols in COLS.items():
if k not in NOE:
continue
for c in cols:
full, noe = FULL[k][c], NOE[k][c]
if k == "eau":
full, noe = full[FULL[k].bins == "linear"], noe[NOE[k].bins == "linear"]
rows_e.append({"table": k, "column": c, "combinations": len(full),
"with_E": int(full.astype(bool).sum()), "without_E": int(noe.astype(bool).sum())})
print(pd.DataFrame(rows_e).to_string(index=False))
if len(NOE) == len(FULL):
lr_n = NOE["lr"]
print("\nleft-right-no-difference without E: readings finding a difference per combination",
lr_n.differences_found.value_counts().to_dict(), "; p_c", lr_n.p_c.astype(float).min(), "to", lr_n.p_c.astype(float).max())
print("cross-section-counts without E (P2 a, P4 b):", NOE["cc"][(NOE["cc"].P2 == "a") & (NOE["cc"].P4 == "b")].to_dict("records"))
ref = lambda t: t[(t.P1 == "a") & (t.P4 == "b")]
print("\nReference values, P1 a, P4 b, with E and without:")
print(" fig5a slope, R2 (with intercept):", ref(FULL["f5a"]).iloc[0][["slope", "R2"]].tolist(), ref(NOE["f5a"]).iloc[0][["slope", "R2"]].tolist())
q = lambda t: t[(t.P1 == "a") & (t.P4 == "b") & (t.P5 == "a") & (t.P7 == "b")].iloc[0][["peak", "median", "IQR"]].tolist()
print(" fig5e peak, median, IQR (P5 a, outliers out):", q(FULL["f5e"]), q(NOE["f5e"]))
q = lambda t: t[(t.P1 == "a") & (t.P4 == "b") & (t.P7 == "b") & (t.bins == "linear")].iloc[0][["peak", "median"]].tolist()
print(" fig6b peak, median (outliers out, linear):", q(FULL["f6b"]), q(NOE["f6b"]))
q = lambda t: t[(t.P1 == "a") & (t.P4 == "b") & (t.P6 == 3) & (t.units == "a") & (t.P7 == "b")].iloc[0][["R3", "R3/R1", "R2/R1", "R3/R2"]].tolist()
print(" fig8 R3, R3/R1, R2/R1, R3/R2 (erosion 3, outliers out):", q(FULL["f8"]), q(NOE["f8"]))
table column combinations with_E without_E
cc holds 8 0 0
cc counts_in_range 8 1 1
f5a match 8 0 0
f5e all_within 32 0 0
f5e median_IQR_within 32 0 0
f5bc reading_a 4 4 4
f5bc reading_b 4 1 1
arl match 128 0 0
f5d match 96 4 1
f5d partial 96 92 94
t2 match 96 0 0
t2 partial 96 24 18
f6a match 8 0 0
f6a partial 8 2 3
f6b match 16 0 0
f6b median_within 16 0 0
rpl match 24 0 0
f5b match 96 0 0
f5b partial 96 96 96
eau match 48 1 1
eau similar 48 1 1
s5 match 48 15 18
f7 match 24 24 24
vr match 24 13 9
f8 match 72 0 0
f8 ratios_within 72 0 12
f9c match 24 14 14
f9e match 24 0 0
f9a match 48 1 0
f9a order_holds 48 12 11
f9f match 48 0 0
f9f partial 48 48 48
left-right-no-difference without E: readings finding a difference per combination {2: 12} ; p_c 0.0547 to 0.414
cross-section-counts without E (P2 a, P4 b): [{'P1': 'a', 'P2': 'a', 'P4': 'b', 'fewest': 409, 'most': 604, 'segments_with_0': 0, 'mean_spacing_um': 0.877, 'counts_in_range': False, 'spacing_ok': True, 'holds': False}, {'P1': 'b', 'P2': 'a', 'P4': 'b', 'fewest': 264, 'most': 394, 'segments_with_0': 0, 'mean_spacing_um': 1.347, 'counts_in_range': True, 'spacing_ok': False, 'holds': False}]
Reference values, P1 a, P4 b, with E and without:
fig5a slope, R2 (with intercept): [np.float64(0.912), np.float64(0.176)] [np.float64(1.025), np.float64(0.199)]
fig5e peak, median, IQR (P5 a, outliers out): [np.float64(1.568), np.float64(1.565), np.float64(1.287)] [np.float64(1.508), np.float64(1.566), np.float64(1.23)]
fig6b peak, median (outliers out, linear): ['1.68e+15', '5.28e+15'] ['1.63e+15', '5.19e+15']
fig8 R3, R3/R1, R2/R1, R3/R2 (erosion 3, outliers out): [np.float64(27.918), np.float64(1.008), np.float64(3.066), np.float64(0.377)] [np.float64(27.918), np.float64(0.997), np.float64(3.044), np.float64(0.388)]
Location E. Without E most of the counts the outcomes rest on do not change. Those that move: fig8's three ratios are within tolerance in 12 of 72 combinations in place of none, which would meet its partial rule; fig9a loses its one matching combination, while the order H1 above H2 above H3 (its partial rule) holds in 11 of 48; fig5d keeps 1 matching combination of the 4; vessel-radius-peak holds in 9 of 24 in place of 13, spline-resistance-five-times in 18 of 48 (both bin scales) in place of 15, the partial rule of table2 in 18 of 96 in place of 24 and that of fig6a in 3 of 8 in place of 2. Without E the fewest cross sections at 0.648 um (P4 b) is 409 and the most still 604. The readings of left-right-no-difference still split in every combination (p of reading c 0.055 to 0.414). So fig8-ellipse-parameters and fig9a-heights-taper depend on E: without it both would be partially reproduced. E is one of the paper's 14 segments, but its cross sections here are of a side branch, not of the artery the paper measured, so neither keeping nor dropping it gives the paper's E. The outcomes keep E, as the claims file's 14 segments do, and the record states the dependence.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| cross-section-counts | not reproduced (proposed change to reproduced not accepted) |
| segmentation-validation-overlap | blocked: data missing |
| fig5a-pvs-area-scales-with-vessel-area | not reproduced |
| fig5e-area-ratio-distribution | not reproduced |
| fig5bc-area-intersegment-variation | blocked: methods underspecified |
| area-ratio-per-location | not reproduced |
| fig5d-polynomial-closest-area-ratio | reproduced (4 of 96 combinations) |
| table2-fit-area-ratio-peaks | partially reproduced |
| fig6a-resistance-power-law | partially reproduced |
| fig6b-resistance-segmentation-peak | not reproduced |
| resistance-per-location | not reproduced |
| fig5b-polynomial-closest-resistance | partially reproduced |
| ellipse-annulus-underestimate | reproduced (one combination, P4 a) |
| spline-resistance-five-times | reproduced |
| fig7-velocity-profiles | reproduced |
| vessel-radius-peak | reproduced |
| fig8-ellipse-parameters | not reproduced |
| fig9c-pvs-two-lobes | reproduced |
| fig9e-pvs-width-asymmetry | not reproduced |
| left-right-no-difference | blocked: methods underspecified |
| fig9a-heights-taper | reproduced (one combination) |
| fig9f-pvs-concave | partially reproduced |
How the analysis differs from the paper's. It starts from the deposited segmentations. The centerline is skimage's 3D thinning, not MATLAB's bwskel with manual removal of points; the longest path through the skeleton, cleaned of points farther from the polynomial fit than the vessel radius and of points that revisit a stretch already passed, stands in for the manual cleaning. Cross sections are sampled by linear interpolation of the masks, not with obliqueslice, and their PVS is the two largest PVS regions within one vessel radius of the vessel. Planes that the PVS mask fills edge to edge (E and F) are cleared. Pixels in both masks are counted as vessel. The resistance is solved by finite differences on the pixel grid, not by finite elements. Reading (b) of P2 is run for the counts only. Every other reading of the claims file is run, and the tables show each. The ellipse is fit to the PVS and vessel together, holes filled, as the claims file has it since a reader answered the question on fig8-ellipse-parameters. These choices changed five outcomes after the first verification (see the attempt record). After the second verification no outcome changed. The data check now shows how far the cleaned path meanders; counts estimated for a path along the fitted centerline fall inside the criterion's range, and a change of cross-section-counts to reproduced was proposed. The third verification did not agree, since no slices were taken along such a path and its step was not measured, so the outcome stays not reproduced. Location E, whose cleaned path follows a side branch, is kept; the outcomes of fig8-ellipse-parameters and fig9a-heights-taper would both be partially reproduced without it.
What was found about the dataset. The earlier published version lacks the PVS masks of L, M and N. The deposited vessel and PVS masks overlap in 11 of 14 files (up to 344565 voxels at E), though the Methods say overlaps were removed. At E and F the PVS mask fills every plane below a depth (87 and 58 planes). The vessel mask of E holds a network of branching vessels, though the Methods say branches were removed by hand. L, M and N state a depth step of 0.648 um while their session descriptions say 1 um; round vessel cross sections are consistent with the stated 0.648 um. The emission wavelengths put 698 nm on the vessel channel (FITC-dextran), which the Methods give as 525 nm. The 2D time series used to validate the segmentation are not in the dandiset.