Camera pose in ECEF, camera matrix <--> View matrix, projection matrix

[rebmaso]

Hello everyone,

Intro

I'm trying to use OsgEarth to generate images taken from specific points of view.
Basically, I want to recreate an aerial image taken from a drone, given camera parameters and the homogenous transform representing the pose of the camera in ECEF ref. system.

Then, I'm also interested in recovering the geographic coordinates from screen coordinates, but that's another thing (see thread #2366)

What I'm struggling with

The point of view I get in OsgEarth is slightly different from what I expect. This tells me that there's an error in the conversion from homogenous transform to view matrix.
I'm sure that the ECEF transform that I'm assigning is correct.

Code

Here's a snippet performing view matrix assignment. Do you see any flaws?
I also attached the code that assigns the projection matrix (which I think is OK).

Init

// ========= Init =========

    osgEarth::initialize();

    osg::ArgumentParser arguments(&argc,argv);

    if ( arguments.read("--help") )
        return usage(argv[0]);

    osgViewer::Viewer viewer(arguments);

    // This is normally called by Viewer::run but we are running our frame loop manually so we need to call it here.
    // viewer.setReleaseContextAtEndOfFrameHint(false);

    // Tell the database pager to not modify the unref settings
    viewer.getDatabasePager()->setUnrefImageDataAfterApplyPolicy( true, false );

    // thread-safe initialization of the OSG wrapper manager. Calling this here
    // prevents the "unsupported wrapper" messages from OSG
    osgDB::Registry::instance()->getObjectWrapperManager()->findWrapper("osg::Image");

    // disable the small-feature culling
    viewer.getCamera()->setSmallFeatureCullingPixelSize(-1.0f);

    // Apply window of custom dimensions
    viewer.apply(new osgViewer::SingleWindow(0, 0, MY_OSG_WIDTH, MY_OSG_HEIGHT));

    viewer.realize();

    // ========= Earth manipulator =========

    // install our default manipulator (do this before calling load)
    EarthManipulator* manip = new EarthManipulator();
    viewer.setCameraManipulator(manip);

    // load an earth file, and support all or our example command-line options
    auto node = MapNodeHelper().load(arguments, &viewer);

    if (!node.valid()) {
        std::cerr << "Invalid OSG node!";
        exit(1);
    }

    viewer.setSceneData( node );

View Matrix

// ========= Manually set osg viewMatrix (VERY CLOSE BUT VIEWPOINT IS A BIT SKEWED)=========

    // REFERENCE: https://stackoverflow.com/questions/44375149/opencv-to-opengl-coordinate-system-transform

    // Homogenous transformation (camera pose in ecef frame)
    osg::Matrixd ecef_T_cam(-0.19242013 ,      -0.88759065   ,     0.41851804 , 4368076.   ,      
                     -0.9395708  ,       0.043574184  ,    -0.33957028 ,  502805.34   ,    
                      0.28316283   ,    -0.45856747    ,   -0.8423388 ,  4605613.  ,       
                      0.       ,         0.        ,        0.        ,        1.);

    osg::Matrixd cam_T_ecef;

    cam_T_ecef.invert_4x4(ecef_T_cam);

    // Rotation
    cam_T_ecef(1,0) *= -1;
    cam_T_ecef(2,0) *= -1;
    cam_T_ecef(1,1) *= -1;
    cam_T_ecef(2,1) *= -1;
    cam_T_ecef(1,2) *= -1;
    cam_T_ecef(2,2) *= -1;

    // Translation
    cam_T_ecef(1,3) *= -1;
    cam_T_ecef(2,3) *= -1;

    osg::Matrixd cam_T_ecef_gl;
    cam_T_ecef_gl.transpose(cam_T_ecef);

    // std::cout << "SETTING W OSG: \n";
    // for(size_t i = 0; i < 4 ; i++) {
    //     for(size_t j = 0; j < 4 ; j++) {
    //         std::cout << cam_T_ecef_gl(i,j) << " ";
    //     }
    //     std::cout << "\n";
    // }
    // std::cout << "\n";

    manip->setByInverseMatrix(cam_T_ecef_gl);

Projection Matrix (THIS PART SHOULD BE OK)

// ========= Manually set osg projectionMatrix (THIS SEEMS TO WORK) =========

    // https://amytabb.com/tips/tutorials/2019/06/28/OpenCV-to-OpenGL-tutorial-essentials/

    double cx = 360.; //principal point x coord
    double cy = 240.; //principal point y coord
    double w = 720.; //image width
    double h = 480.; //image height
    double near = 0.1; //near plane
    double far = 1.0; //far plane
    double f = 480.; // focal length in px (assume fx=fy)

    // Projection matrix
    osg::Matrixd K_gl(2*f/w, 0.0, (w - 2*cx)/w, 0.0,
                        0.0, 2*f/h, (h - 2*cy)/h, 0.0,
                        0.0, 0.0, (-far - near) / (far - near), -2.0*far*near/(far-near),
                        0.0, 0.0, -1.0, 0.0);

    osg::Matrixd projMatrix;
    projMatrix.transpose(K_gl);

    viewer.getCamera()->setProjectionMatrix(projMatrix);

Results

below are what I should see (real), and what I actually see (osgearth). By looking at edges, you can see that images are not aligned. It looks as if the camera is looking from a higher altitude. This tells me that it's not a projection matrix assignment error.

[jasonbeverage]

There are a few things you can try to figure out what is going on:

1. If you are going to use the EarthManipulator, turn off terrain avoidance to make sure it's not moving the camera to avoid going under the terrain. For example do something like this:
   manip ->getSettings()->setTerrainAvoidanceEnabled( false );

2. I'm not sure what the source of the data that you are loading is, but make sure that you are using the correct vertical datum for the elevation. Try adding egm96 to your elevation layer in your earth file and see if things look like they line up better.

3. Your projection matrix is using fixed values for the width and height. Be certain that the window you are rendering to has the same dimensions or else osg will resize the projection matrix and you'll get different results than you are expecting.

[rebmaso]

Looks like #2 did the trick! I was dealing with elevation data from swisstopo. Apparently height is not relative to the wgs84 geoid, but to the egm96 model, as you pointed out!
This was brilliant, thanks a lot!!

Here's part of my earth file incorporating Jason's fix (vdatum = egm96).

    <MBTilesElevation name="elevation">
        <url>../data/swisstopo_dem/output/output.mbtiles</url>
        <min_lod>15</min_lod>
        <vdatum>egm96</vdatum>
    </MBTilesElevation>

[XenonofArcticus]

#3 is likely to be the problem in my experience.

…

On Tue, Oct 31, 2023 at 8:22 AM Jason Beverage ***@***.***> wrote: There are a few things you can try to figure out what is going on: 1. If you are going to use the EarthManipulator, turn off terrain avoidance to make sure it's not moving the camera to avoid going under the terrain. For example do something like this: manip ->getSettings()->setTerrainAvoidanceEnabled( false ); 2. I'm not sure what the source of the data that you are loading is, but make sure that you are using the correct vertical datum for the elevation. Try adding egm96 to your elevation layer in your earth file and see if things look like they line up better. 3. Your projection matrix is using fixed values for the width and height. Be certain that the window you are rendering to has the same dimensions or else osg will resize the projection matrix and you'll get different results than you are expecting. — Reply to this email directly, view it on GitHub <#2369 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAPFE3AEDIU7XUXLG4ZOBHDYCECSJAVCNFSM6AAAAAA6XQQTYWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TIMZVHE2DA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>