<jupyter_start><jupyter_text>Triangulation error vs track length  Build the ETH3D pipeline tool to merge Lidar scans<jupyter_code>%%bash
cd ..
git clone https://github.com/ETH3D/dataset-pipeline/
cd dataset-pipeline/
git clone https://github.com/laurentkneip/opengv
cd opengv
rm -rf build && mkdir build && cd build
cmake .. && make -j
cd ..
rm -rf build && mkdir build/ && cd build
ls ../opengv/build
cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -Dopengv_DIR=$(pwd)/../opengv/build ..
make -j
def compute_accuracy_pcd(tag, tol):
    sfm_dir = Paths().interpolate(scene=scene, outputs=outputs_path, method=keypoints, tag=tag).triangulation
    accuracy_path = sfm_dir / f"reconstruction.tolerance_{tol}.ply"
    cmd = [
        "../multi-view-evaluation/build/ETH3DMultiViewEvaluation",
        "--ground_truth_mlp_path", str(scan_path),
        "--reconstruction_ply_path", str(sfm_dir/"reconstruction.ply"),
        "--tolerances", str(tol),
        "--accuracy_cloud_output_path", str(sfm_dir/"reconstruction"),
    ]
    subprocess.run(cmd, check=True)
    assert accuracy_path.exists()
    return sfm_dir, accuracy_path

def get_error_tl(sfm_path, accuracy_path):
    # extract the accuracy labels from the evaluation outputs
    pcd_acc = open3d.io.read_point_cloud(str(accuracy_path))
    colors = np.asarray(pcd_acc.colors)
    good = np.array([0, 1, 0])
    bad = np.array([1, 0, 0])
    ignore = np.array([0, 0, 1])
    is_ignore = np.all(colors == ignore, -1)
    is_good = np.all(colors == good, -1)
    print(is_good[~is_ignore].mean())

    # compute the SfM track lengths
    rec = pycolmap.Reconstruction(sfm_path)
    pids = np.array(sorted(rec.points3D))
    p3ds = np.array([rec.points3D[i].xyz for i in pids])
    all_tl = np.array([rec.points3D[i].track.length() for i in pids])

    # propagate the accuracy labels from the evaluation pcd to the SfM IDs
    dist_sfm_acc = cdist(p3ds, np.asarray(pcd_acc.points))
    assert np.all(np.min(dist_sfm_acc, 1) < 1e-5)
    nearest = np.argmin(dist_sfm_acc, 1)
    pid2acc = {pids[idx]: i for idx, i in enumerate(nearest)}
    valid = np.array([not is_ignore[pid2acc[i]] for i in pids])

    # compute the 3D error w.r.t the GT Lidar pointcloud
    p3ds_gt = []
    for p3d in tqdm(p3ds):
        _, idx_gt, _ = tree_eval.search_knn_vector_3d(p3d, 1)
        p3d_gt = pcd_eval.points[idx_gt[0]]
        p3ds_gt.append(p3d_gt)
    p3ds_gt = np.array(p3ds_gt)
    all_errs_3d = np.linalg.norm(p3ds - p3ds_gt, axis=-1)

    return all_tl[valid], all_errs_3d[valid]

scan_path = dataset_path / scene / "dslr_scan_eval/scan_alignment.mlp"
scan_merged_path = scan_path.parent / "merged.ply"
!../dataset-pipeline/build/NormalEstimator -i $scan_path -o $scan_merged_path
pcd_eval = open3d.io.read_point_cloud(str(scan_merged_path))
tree_eval = open3d.geometry.KDTreeFlann(pcd_eval)
tol = 0.01

sfm_raw, acc_raw = compute_accuracy_pcd(tag_raw, tol)
tl_raw, errs_3d_raw = get_error_tl(sfm_raw, acc_raw)
print(tag_raw, np.mean(errs_3d_raw < tol))

sfm_ref, acc_ref = compute_accuracy_pcd(tag_refined, tol)
tl_ref, errs_3d_ref = get_error_tl(sfm_ref, acc_ref)
print(tag_refined, np.mean(errs_3d_ref < tol))
max_tl = 8
min_err = 2e-3
fig = go.Figure()

fig.add_shape(
    type='line', x0=0., y0=min_err, x1=1., y1=min_err,
    line=dict(color='Gray', dash='dot'),
    yref='y', xref= 'paper', 
)

fig.add_annotation(x=-.004, y=np.log10(2e-3),
            text="Lidar", xref="paper",
           xanchor="right", yanchor='middle',
            showarrow=False)

fig.add_trace(go.Box(
    y=np.maximum(errs_3d_raw, min_err),
    x=np.minimum(tl_raw, max_tl),
    name='raw',
    boxpoints=False, jitter=0.3,
))

fig.add_trace(go.Box(
    y=np.maximum(errs_3d_ref, min_err),
    x=np.minimum(tl_ref, max_tl),
    name='refined',
    boxpoints=False, jitter=0.3,
))

bins = np.unique(np.minimum(tl_raw, max_tl))
ax = dict(linewidth=1, linecolor='black', mirror=True, gridcolor='rgb(210, 210, 210)')
fig.update_layout(
    yaxis=dict(
        type="log",
        tickvals=[1e-3, 1e-2, 1e-1, 1],
        ticktext=['1mm', '1cm', '10cm', '1m'],
        range=[-3.1, 0.9],
        **ax,
    ),
    xaxis = dict(
        tickmode = 'array', tickvals = bins,
        ticktext = list(map(str, bins[:-1])) + [str(bins[-1])+'+'],
        **ax,
    )
)

fig.update_layout(
    xaxis_title='track length',
    boxmode='group',
    paper_bgcolor='rgba(0,0,0,0)',
    plot_bgcolor='rgba(0,0,0,0)',
    legend=dict(
        orientation="v",
        yanchor="top",
        y=0.98,
        xanchor="right",
        x=0.93,
    ),
    font=dict(
        size=18,
        color="Black"
    ),
    title=dict(
        text='3D triangulation error',
        y=0.99,
        x=0.54,
        xanchor='center',
        yanchor='top',
    ),
    width=700,
    height=600,
    margin=go.layout.Margin(
        l=0, #left margin
        r=0, #right margin
        b=0, #bottom margin
        t=35, #top margin
    ),
)

fig.show()
# fig.write_image("plots/boxplot_error_vs_trancklength.pdf")<jupyter_output><empty_output><jupyter_text>Distribution of point displacements<jupyter_code>rec = pycolmap.Reconstruction(sfm_ref)
rec_raw = pycolmap.Reconstruction(sfm_raw)

def gather_projections(images, points3D):
    im2xy = {}
    for imid, image in rec.images.items():
        cam = rec.cameras[image.camera_id]
        for idx, p2d in enumerate(image.points2D):
            if not p2d.has_point3D():
                continue
            p3d = rec.points3D[p2d.point3D_id]
            xy, = image.project([p3d])
            xy = cam.world_to_image(xy)
            im2xy[(imid, idx)] = xy
    return im2xy

def gather_detections(images):
    im2xy = {}
    for imid, image in rec.images.items():
        for idx, p2d in enumerate(image.points2D):
            if not p2d.has_point3D():
                continue
            im2xy[(imid, idx)] = p2d.xy
    return im2xy

im2xy_ref = gather_projections(images, points3D)
im2xy_raw = gather_detections(images_raw)
# im2xy_raw = gather_projections(images_raw, points3D_raw)  # uncomment for movement vs geometric refinement
# im2xy_raw = gather_detections(images)  # uncomment for movement of FBA only

keys = list(im2xy_ref.keys() & im2xy_raw.keys())
err = np.array([np.linalg.norm(im2xy_ref[k] - im2xy_raw[k]) for k in keys])
v = np.sort(err)
cum = np.arange(len(v)) / (len(v)-1) * 100
# TODO: subsample based on density

fig = go.Figure(data=go.Scatter(
    x=v, y=cum,
    fill='tozeroy',
    fillcolor='rgba(99, 110, 250, .2)',
))

ax = dict(linewidth=1, linecolor='black', mirror=True, gridcolor='rgb(210, 210, 210)')
fig.update_layout(
    yaxis=dict(
        range=[0, 100],
        **ax,
    ),
    xaxis = dict(
        type='log',
        range=np.log10([0.2, 12]).tolist(),
        **ax,
    )
)

fig.update_layout(
    xaxis_title='point movement [px]',
    yaxis_title='cumulative distr. [%]',
    paper_bgcolor='rgba(0,0,0,0)',
    plot_bgcolor='rgba(0,0,0,0)',
    font=dict(
        size=18,
        color="Black"
    ),
    title=dict(
        text='Distribution of 2D point movement',
        y=0.99,
        x=0.54,
        xanchor='center',
        yanchor='top',
    ),
    width=720,
    height=400,
    margin=go.layout.Margin(
        l=0, #left margin
        r=0, #right margin
        b=0, #bottom margin
        t=35, #top margin
    ),
)

fig.show()
# fig.write_image("plots/cumulative_movement.pdf")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14778856526429343, 'fraction_numerical': 0.03087918015102481, 'mean_word_length': 3.1412618648799553, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 9, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '654348', 'license': 'Apache-2.0', 'n_tokens_mistral': 3011, 'n_tokens_neox': 2888, 'n_words': 523, 'path': 'notebooks/plot_eth3d_triangulation.ipynb', 'repo_name': 'Tangshengku/pixel-perfect-sfm', 'size': '7416'}