Add files via upload

README.md

@@ -1 +1 @@
-# rcal-toolbox
+rcbox : robust calibration toolbox by outlier-aware iterative solvers

build_and_run.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+# build binaries
+cd rcbox
+python setup.py build_ext --inplace
+cd ..
+
+# run unit tests
+py.test rcbox/unit_tests/utest_*
+
+# run script examples & make figures
+python example_rmds.py
+python example_toa.py
+python example_tdoa.py
+

example_rmds.py

@@ -0,0 +1,96 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Aug 20 2018
+
+@author: Charles Vanwynsberghe
+
+Simple simulation example for robust diffuse field calibration by RMDS
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d.art3d import juggle_axes
+from scipy.spatial.distance import squareform
+
+from rcbox.utils import get_D
+from rcbox.rmds import RMDS
+
+np.random.seed(3212)
+
+# %% set up data
+
+X = np.random.randn(16, 3)
+D = get_D(X)
+X_ini = X + 0.1*np.random.randn(*X.shape)
+
+# add outliers to toa
+K = 20
+
+# outlier values follow Rayleigh distribution with mode = D ground truth
+D_noisy_table = np.random.rayleigh(scale=np.triu(D))
+D_noisy_table += D_noisy_table.T
+
+mask = np.zeros((int((D.size - D.shape[0])/2),))
+mask[0:K] = 1
+np.random.shuffle(mask)
+mask = squareform(mask).astype(bool)
+
+D_noisy = D.copy()
+D_noisy[mask] = D_noisy_table[mask]
+o = D_noisy - D
+e = np.triu(0.01*np.random.randn(*D.shape), k=1)
+e = e + e.T
+
+D_noisy = D_noisy + e
+
+# %% find geometry from euclidean distance matrix
+
+# least square case (lambda = infinity)
+solver_ls = RMDS(D_noisy)
+solver_ls.Run(lbda=np.inf, Xinit=X_ini, EpsLim=10**-6, itmax=10000)
+solver_ls.align(X)
+
+# outlier-aware case
+solver_oa = RMDS(D_noisy)
+solver_oa.Run(lbda=0.2, Xinit=X_ini, EpsLim=10**-6, itmax=10000)
+solver_oa.align(X)
+
+# %% plot result
+
+fig = plt.figure(1, figsize=(8, 3.5))
+fig.suptitle("MDS example")
+
+ax_ls = fig.add_subplot(121, projection='3d')
+ax_ls.set_title("least-square (smacof)")
+
+ax_ls.scatter(*solver_ls.x_aligned.T, marker='x', facecolor='RoyalBlue',
+              color='none', s=35, linewidths=1.5, label='r estimated')
+
+ax_ls.scatter(*solver_ls.x_ref_centered.T, marker='o', facecolor='none',
+              color='RoyalBlue', s=35, linewidths=1.5, label='r source')
+
+ax_oa = fig.add_subplot(122, projection='3d')
+ax_oa.set_title("outlier-aware")
+
+ax_oa.scatter(*solver_oa.x_aligned.T, marker='x', facecolor='RoyalBlue',
+              color='none', s=35, linewidths=1.5, label='r estimated')
+
+ax_oa.scatter(*solver_oa.x_ref_centered.T, marker='o', facecolor='none',
+              color='RoyalBlue', s=35, linewidths=1.5, label='r source')
+
+plt.savefig("example_rmds_geometries.pdf")
+
+# %% L-curve choice for lambda
+
+solver_oa.tune_lambda(np.geomspace(1, 0.01, 200), Xinit=X_ini,
+                      EpsLim=10**-6, itmax=10000)
+
+plt.figure()
+plt.plot(solver_oa.lambda_sequence, 2*solver_oa.k_seq)
+plt.axhline(2*K, color='k', linestyle='dashed')
+plt.axhline(solver_oa.Nr*(solver_oa.Nr-1), color='grey', linestyle='dashed')
+
+plt.xlabel(r"$\lambda$")
+plt.ylabel(r"$||\mathbf{O}||_0$")
+
+plt.savefig("example_rmds_Lcurve.pdf")

example_tdoa.py

@@ -0,0 +1,86 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Aug 20 2018
+
+@author: Charles Vanwynsberghe
+
+Simple simulation example for robust TDOA calibration
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+from rcbox.rtdoa import ROno, denoise_tdoa_soft
+from rcbox.utils import make_TDOA
+
+np.random.seed(1234)  # fix random number generator
+
+# %% set up variables
+
+N = 30  # number of sources
+M = 20  # number of microphones
+K = 20  # number of outlier per TDOA matrix (sparsity with K << M*(M-1)/2)
+sig_o = 10  # standard deviation of outlier amplitudes
+reg_l1 = 0.5
+
+# %% generate geometry and TDOA
+
+X = np.random.randn(M, 3)
+S = 2*np.random.randn(N, 3)
+
+tau = make_TDOA(S, X, c0=1.)
+
+# generate initial guess
+r0 = X + 0.5*np.random.randn(*X.shape)
+s0 = S + 0.5*np.random.randn(*S.shape)
+
+# add outliers to tau
+o = np.zeros_like(tau)
+idx = np.triu_indices(X.shape[0], k=1)
+
+for n in range(o.shape[0]):
+    a = np.arange(idx[0].size)
+    np.random.shuffle(a)
+    for k in range(K):
+        idx_k = (idx[0][a[k]], idx[1][a[k]])
+        o[n][idx_k] = sig_o*np.random.randn()
+
+    o[n, ...] = o[n, ...] - o[n, ...].T
+
+tau_noisy = tau + o
+
+# %% compare TDOA solvers
+
+fig = plt.figure(1, figsize=(14, 4))
+fig.suptitle("TDOA example")
+
+ax_0 = fig.add_subplot(131, projection='3d')
+ax_1 = fig.add_subplot(132, projection='3d')
+ax_2 = fig.add_subplot(133, projection='3d')
+
+# Case 1: baseline
+ono = ROno(tau=tau_noisy, c0=1., fast=True)
+ono.Init(r0=r0, s0=s0)
+ono.Run(np.inf, param_converge=1e-6, verbose=0)
+ono.Plot_result(r_gt=X, ax=ax_0)  # , show_init=True, show_iters=True)
+ax_0.set_title("TDOA - baseline\n")
+
+# Case2: outlier-aware denoising + baseline
+tau_denoised = np.zeros_like(tau_noisy)
+for n in range(N):
+    tau_denoised[n, ...], _, _ = denoise_tdoa_soft(tau_noisy[n, ...], reg_l1,
+                                                   eps=1e-6, tmax=200)
+
+dono = ROno(tau=tau_denoised, c0=1., fast=True)
+dono.Init(r0=r0, s0=s0)
+dono.Run(np.inf, param_converge=1e-6, verbose=0)
+dono.Plot_result(r_gt=X, ax=ax_1)  # , show_init=True, show_iters=True)
+ax_1.set_title("TDOA - denoising + baseline\n")
+
+# case 3: outlier-aware R-Ono
+rono = ROno(tau=tau_noisy, c0=1., fast=True)
+rono.Init(r0=r0, s0=s0)
+rono.Run(reg_l1, param_converge=1e-6, verbose=0)
+rono.Plot_result(r_gt=X, ax=ax_2)  # , show_init=True, show_iters=True)
+ax_2.set_title("TDOA - outlier-aware\n")
+
+plt.savefig("example_tdoa_geometries.pdf")

example_toa.py

@@ -0,0 +1,117 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Aug 20 2018
+
+@author: Charles Vanwynsberghe
+
+Simple simulation example for robust TOA calibration by RMDU
+"""
+from __future__ import division
+
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d.art3d import juggle_axes
+
+from rcbox.rmds import RMDU
+from rcbox.utils import make_TOA
+
+np.random.seed(2121)
+
+# %% set up data
+
+Xr = np.random.randn(37, 3)
+Xs = 2*np.random.randn(10, 3)
+X_all = np.concatenate((Xr, Xs), axis=0)
+
+c0 = 1.
+
+toa = make_TOA(Xs, Xr, c0=1.)
+
+# add outliers to toa
+K = 50
+
+# outlier values follow Rayleigh distribution with mode = toa ground truth
+toa_noisy_table = np.random.rayleigh(scale=toa)
+
+mask = np.zeros((toa.size))
+mask[0:K] = 1
+np.random.shuffle(mask)
+mask = mask.reshape(toa.shape).astype(bool)
+
+toa_noisy = toa.copy()
+toa_noisy[mask] = toa_noisy_table[mask]
+o = toa_noisy - toa
+e = 0.001*np.random.randn(*toa.shape)
+
+toa_noisy += e
+
+X_ini = X_all + 0.01*np.random.randn(*X_all.shape)
+
+# %% find geometry from TOA matrix
+
+# least square case (lambda = infinity)
+toa_solver = RMDU(c0*toa_noisy)
+toa_solver.Run(lbda=np.inf, Xinit=X_ini, itmax=15000, EpsLim=1e-6, verbose=0)
+toa_solver.align(X_all)
+
+# outlier-aware case
+rtoa_solver = RMDU(c0*toa_noisy)
+rtoa_solver.Run(lbda=0.15, Xinit=X_ini, itmax=15000, EpsLim=1e-6, verbose=0)
+rtoa_solver.align(X_all)
+
+# %% plot result
+
+fig = plt.figure(1, figsize=(8, 3.5))
+fig.suptitle("MDU example")
+
+ax_ls = fig.add_subplot(121, projection='3d')
+ax_ls.scatter(*toa_solver.x_aligned[0:toa_solver.Nr, :].T,
+              marker='x', facecolor='RoyalBlue',
+              color='none', s=35, linewidths=1.5, label='r estimated')
+
+ax_ls.scatter(*toa_solver.x_aligned[toa_solver.Nr::, :].T,
+              marker='x', facecolor='red',
+              color='none', s=35, linewidths=1.5, label='r estimated')
+
+ax_ls.scatter(*toa_solver.x_ref_centered[0:toa_solver.Nr, :].T, marker='o',
+              facecolor='none',
+              color='RoyalBlue', s=35, linewidths=1.5, label='r source')
+
+ax_ls.scatter(*toa_solver.x_ref_centered[toa_solver.Nr::, :].T, marker='o',
+              facecolor='none',
+              color='red', s=35, linewidths=1.5, label='r source')
+
+ax_oa = fig.add_subplot(122, projection='3d')
+ax_oa.scatter(*rtoa_solver.x_aligned[0:rtoa_solver.Nr, :].T,
+              marker='x', facecolor='RoyalBlue',
+              color='none', s=35, linewidths=1.5, label='r estimated')
+
+ax_oa.scatter(*rtoa_solver.x_aligned[rtoa_solver.Nr::, :].T,
+              marker='x', facecolor='red',
+              color='none', s=35, linewidths=1.5, label='r estimated')
+
+ax_oa.scatter(*rtoa_solver.x_ref_centered[0:rtoa_solver.Nr, :].T, marker='o',
+              facecolor='none',
+              color='RoyalBlue', s=35, linewidths=1.5, label='r source')
+
+ax_oa.scatter(*rtoa_solver.x_ref_centered[rtoa_solver.Nr::, :].T, marker='o',
+              facecolor='none',
+              color='red', s=35, linewidths=1.5, label='r source')
+
+plt.savefig("example_toa_geometries.pdf")
+
+# %% L curve choice of o
+rtoa_solver.tune_lambda(np.geomspace(10, 0.001, 200), Xinit=X_ini,
+                        EpsLim=10**-6, itmax=10000)
+
+plt.figure(2)
+plt.plot(rtoa_solver.lambda_sequence, rtoa_solver.k_seq)
+plt.axhline(K, color='k', linestyle='dashed')
+plt.axhline(rtoa_solver.Nr*rtoa_solver.Ns,
+            color='grey', linestyle='dashed')
+
+plt.xlabel(r"$\lambda$")
+plt.ylabel(r"$||\mathbf{O}||_0$")
+
+plt.savefig("example_toa_Lcurve.pdf")