I'm also having the same issue. Is anyone planning on implementing this, or are there any other suggestions for extracting an uncertainty of the IMU-Camera calibration?

Here I list my solution to get the covariance of estimated result.

First of all, the kalibr don't calculate the covariance of the IMU-Camera calibration. If we modify the code to calculate the covariance, it seems much effort and cannot make sure the modification is right. Therefore, I bypass the direct calculation method by calculating the inverse the Hessian block.

In the linear system solver, kalibr will build a problem like H * x = b, where the H is the Hessian block, and the last 8 parameters in x are the quaternion (N=4), translation(N=3) and time offset(N=1). We could add log to print the last 8 parameters and the calculated IMU-Camera transformation to check whether they are the same value. I checked they are the same. Then we extract the corresponding Hessian block by the index. Because Hessian is the sparse matrix, the inverse of the block is the covariance what we needed.

Although I obtain the covariance of T_BC(transformation from camera to IMU), but this value is about 1e-12, which don't seem any sense for application.

Thanks for the insight! If you can upload the changes to made to extract the relevant block from the Hessian I can look through it and see if I can figure out why the values are so extreme.

Extract Hessian in below file,

add below code before the final return

if (_H._M.rows() > 100) {
    printInfo("BlockCholeskyLinearSystemSolver::solveSystem()");
    cout << "Hessian size = " << _H._M.rows() << " X " << _H._M.cols() << ", block size = " << _H._M.bRows()
         << " X " << _H._M.bCols();
    if (_H._M.rows() > 3000 || _H._M.cols() > 3000) {
        cout << endl;
    } else {
        cout << "|H| = " << _H._M.toDense().norm() << endl;
    }
    cout << "rhs size = " << _rhs.size() << ", |rhs| = " << _rhs.norm() << endl;
    cout << "dx size = " << outDx.size() << ", |dx| = " << outDx.norm() << endl;
    // last 3-2 block
    int startBlockRow = _H._M.bRows() - 3;
    int startBlockCol = _H._M.bCols() - 3;
    int r0 = _H._M.rowBaseOfBlock(startBlockRow);
    int r1 = _H._M.rowBaseOfBlock(startBlockRow   2);
    int c0 = _H._M.colBaseOfBlock(startBlockCol);
    int c1 = _H._M.colBaseOfBlock(startBlockCol   2);
    Eigen::MatrixXd hBlock =
        _H._M.slice(startBlockRow, startBlockRow   2, startBlockCol, startBlockCol   2)->toDense();
    cout << "H[" << r0 << ":" << r1 << ", " << c0 << ":" << c1 << "] = " << endl << hBlock << endl;
    cout << "H[" << r0 << ":" << r1 << ", " << c0 << ":" << c1 << "]^-1 = " << endl << hBlock.inverse() << endl;
}

Here list my repository I modified for reference.
kalibr-JefferyModification

Please enable some macro definition to show the final result.

Have you looked into this bug yet? It would be really helpful to be able to extract the variance of the time offset in particular.

I tried to recover covariances of the estimated parameters here a long while ago. You have get some info therein.