Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How to parse the scaling ratio and horizontal rotation angle of 3dtiles data #757

Closed
zysyz opened this issue Sep 23, 2024 · 3 comments
Closed
Labels
question Further information is requested

Comments

@zysyz
Copy link

zysyz commented Sep 23, 2024

Hello, I can currently load 3dtiles data, but I don't know how to overlay the local 3dtiles model with the satellite map. Currently, I have parsed the latitude and longitude of the model from tileset.json, but I don't know how to parse the scaling ratio and horizontal rotation angle of the model. So, I would like to ask for advice. Thank you!

@zysyz zysyz added the question Further information is requested label Sep 23, 2024
@gkjohnson
Copy link
Contributor

I'm not sure what you mean. If your data set is specified at a lat / lon then typically it would be positioned at the correct location on the globe, already. If your method for processing just produced lat / lon to place it at then you'll need to know what coordinate frame orientation relative to the globe ellipsoid it was produced in.

As the question template states - please include a link to a project or data. Otherwise I have to guess at what you're referring to.

@zysyz
Copy link
Author

zysyz commented Sep 24, 2024

I'm not sure what you mean. If your data set is specified at a lat / lon then typically it would be positioned at the correct location on the globe, already. If your method for processing just produced lat / lon to place it at then you'll need to know what coordinate frame orientation relative to the globe ellipsoid it was produced in.我不知道你是什么意思。如果您的数据集指定为纬度/经度,则通常它已经定位在地球仪上的正确位置。如果你的处理方法只是产生了纬度/经度来放置它,那么你需要知道它相对于产生它的地球仪椭球体的坐标系方向。

As the question template states - please include a link to a project or data. Otherwise I have to guess at what you're referring to.如问题模板所述-请包含项目或数据的链接。否则我只能猜测你指的是什么。

I use the geo-tree library to display the map, which is flat. Then, I parse tileset.json to obtain latitude and longitude, and overlay the 3D tiles model on the map. However, the horizontal rotation angle and scaling ratio of the model are incorrect and cannot match the map.
image

My code example is as follows:

var scene = new THREE.Scene();
scene.background = new THREE.Color(0xffffff); 

let needsRender = true;

var camera = new THREE.PerspectiveCamera(
  75,
  window.innerWidth / window.innerHeight,
  0.1,
  100000000
);

camera.position.set(0, 0, 5000); 
scene.add(camera); 

var renderer = new THREE.WebGLRenderer();

const light = new THREE.AmbientLight(0xffffff, 5);
scene.add(light);

var map = new Map.MapView(
  Map.MapView.PLANAR,
  new UrlTileProvider(
    {
      minZoom:0,
      maxZoom:18,
      url:'http://0414.gggis.com/maps/vt?lyrs=s&x={x}&y={y}&z={z}'
    }
  )
);
scene.add(map);
map.updateMatrixWorld(true);

const tilesRenderer = new TilesRenderer( './model/tileset.json' ); 

const dracoLoader = new DRACOLoader(); 
dracoLoader.setDecoderPath( 'https://unpkg.com/[email protected]/examples/jsm/libs/draco/' );
dracoLoader.setDecoderConfig({ type: "js" }); 

const ktx2Loader = new KTX2Loader();
//ktx2Loader.setTranscoderPath( 'https://unpkg.com/[email protected]/examples/jsm/libs/basis/' );
ktx2Loader.detectSupport(renderer);

const loader = new GLTFLoader(tilesRenderer.manager);
loader.setDRACOLoader(dracoLoader);
loader.setKTX2Loader(ktx2Loader);
loader.register(() => new GLTFCesiumRTCExtension());

tilesRenderer.fetchOptions.mode = 'cors';
tilesRenderer.manager.addHandler( /\.(gltf|glb)$/g, loader );
tilesRenderer.setCamera( camera );
tilesRenderer.setResolutionFromRenderer( camera, renderer );

let tilesParent = new THREE.Group();
tilesParent.add(tilesRenderer.group);

let TopParent = new THREE.Group();
TopParent.add(tilesParent);
scene.add(TopParent);

loadJsonData();
let lon,lat,height;

let modelLoaded = false; 

function render() {
  stats.begin();
  requestAnimationFrame(render);
  
  if (tilesRenderer) {
    //if (needsRender) {
      let box = new THREE.Box3();
      if (tilesRenderer.getBoundingBox(box)) {
        let top_position = TopParent.position.clone();
        if (!modelLoaded) {
          camera.position.set(top_position.x+100, top_position.y+300, top_position.z+200);
          controls.target = new THREE.Vector3(top_position.x, top_position.y, top_position.z);
          camera.lookAt( new THREE.Vector3(top_position.x, top_position.y, top_position.z) );
          modelLoaded = true;
        }
      }
    camera.updateMatrixWorld();

    tilesRenderer.update();
    renderer.render(scene,camera) 
    stats.end();
  }
  

  controls.update();
}
render();

window.addEventListener("resize", () => {
  camera.aspect = window.innerWidth / window.innerHeight;
  camera.updateProjectionMatrix();
  renderer.setSize(window.innerWidth, window.innerHeight);
  renderer.setPixelRatio(window.devicePixelRatio);
});

tilesRenderer.onLoadTileSet = () => {
  const sphere = new THREE.Sphere();
  tilesRenderer.getBoundingSphere(sphere);
  tilesRenderer.group.position.copy( sphere.center ).multiplyScalar( - 1 );

  const position = sphere.center.clone();
  const distanceToEllipsoidCenter = position.length();

  const surfaceDirection = position.normalize();
  const up = new THREE.Vector3(0, 1, 0);
  const rotationToNorthPole = rotationBetweenDirections(surfaceDirection, up);
  tilesParent.quaternion.copy(rotationToNorthPole);
  
  let coordConvert = new CoordConvert();
  let change_potision = coordConvert.lnglat_to_mercator(lon,lat); // Convert latitude and longitude to Mercator coordinates
  TopParent.position.x = change_potision.x;
  TopParent.position.z = -change_potision.y;
  TopParent.position.y = height;  // Move the position of the model along the Y-axis so that it is on the surface of the Earth

  // Adjust the horizontal rotation angle of the model
  TopParent.rotation.set(0,Math.PI/1-0.153,0);

};

// Calculate the rotation quaternion between two directional vectors
function rotationBetweenDirections( dir1, dir2 ) {

	const rotation = new THREE.Quaternion();
	const a = new THREE.Vector3().crossVectors( dir1, dir2 );
	rotation.x = a.x;
	rotation.y = a.y;
	rotation.z = a.z;
	rotation.w = 1 + dir1.clone().dot( dir2 );
	rotation.normalize();

	return rotation;

}

//Analyze tileset.json to obtain the latitude, longitude, and height of the model
function loadJsonData() {
  fetch(tilesRenderer.rootURL)
    .then(response => response.json())
    .then(data => {
      
      if (!data || !data.root || !data.root.boundingVolume) {
        throw new Error('Invalid JSON data structure');
      }

      let root = data.root;
      let boundingVolume = root.boundingVolume;
      let modelCenter;
      if (boundingVolume.hasOwnProperty("box")) {
          let x = boundingVolume.box[0];
          let y = boundingVolume.box[1];
          let z = boundingVolume.box[2];
          let w = 1;
          modelCenter = new THREE.Vector4(x, y, z, w);
      }
      else if (boundingVolume.hasOwnProperty("sphere")) {
        let x = boundingVolume.sphere[0];
        let y = boundingVolume.sphere[1];
        let z = boundingVolume.sphere[2];
        let w = 1;
        modelCenter = new THREE.Vector4(x, y, z, w);
      }

      if (root.hasOwnProperty("transform")) {
        let transform = root.transform;
        let matrix4 = new THREE.Matrix4();
        matrix4.set(
            transform[0], transform[4], transform[8], transform[12],
            transform[1], transform[5], transform[9], transform[13],
            transform[2], transform[6], transform[10], transform[14],
            transform[3], transform[7], transform[11], transform[15]
        );

        let wgs84Vector4 = modelCenter.applyMatrix4(matrix4);
        let wgs84Vector3 = new THREE.Vector3(wgs84Vector4.x, wgs84Vector4.y, wgs84Vector4.z);
    
        let wgs84Cartographic = fromCartesianToWGS84(wgs84Vector3);
        lon = THREE.MathUtils.radToDeg(wgs84Cartographic.longitude);
        lat = THREE.MathUtils.radToDeg(wgs84Cartographic.latitude);
        height = wgs84Cartographic.height;

      }
      else{
        let wgs84Vector3 = new THREE.Vector3(modelCenter.x, modelCenter.y, modelCenter.z);
        let wgs84Cartographic = fromCartesianToWGS84(wgs84Vector3);
        lon = THREE.MathUtils.radToDeg(wgs84Cartographic.longitude);
        lat = THREE.MathUtils.radToDeg(wgs84Cartographic.latitude);
        height = wgs84Cartographic.height;
      }
    })
    .catch(error => console.error('Error loading JSON:', error));
}

// Convert Cartesian coordinates to WGS84 coordinates
function fromCartesianToWGS84(cartesian) {
    const ellipsoid = {
        oneOverRadii: new THREE.Vector3(1.0 / 6378137.0, 1.0 / 6378137.0, 1.0 / 6356752.3142451793),
        oneOverRadiiSquared: new THREE.Vector3(1.0 / (6378137.0 * 6378137.0), 1.0 / (6378137.0 * 6378137.0), 1.0 / (6356752.3142451793 * 6356752.3142451793)),
        _centerToleranceSquared: 1e-10
    };

    const p = scaleToGeodeticSurface(cartesian, ellipsoid);

    if (!p) {
        return undefined;
    }

    const n = p.clone().multiply(ellipsoid.oneOverRadiiSquared).normalize();
    const h = cartesian.clone().sub(p);

    const longitude = Math.atan2(n.y, n.x);
    const latitude = Math.asin(n.z);
    const height = Math.sign(h.dot(cartesian)) * h.length();

    return {
        longitude: longitude,
        latitude: latitude,
        height: height
    };
}

//This function imitates the implementation of Cesium, which converts Cartesian coordinates into projection points on the surface of the Earth
function scaleToGeodeticSurface(cartesian, ellipsoid) {
  const oneOverRadii = ellipsoid.oneOverRadii;
  const oneOverRadiiSquared = ellipsoid.oneOverRadiiSquared;
  const centerToleranceSquared = ellipsoid._centerToleranceSquared;

  const positionX = cartesian.x;
  const positionY = cartesian.y;
  const positionZ = cartesian.z;

  const oneOverRadiiX = oneOverRadii.x;
  const oneOverRadiiY = oneOverRadii.y;
  const oneOverRadiiZ = oneOverRadii.z;

  const x2 = positionX * positionX * oneOverRadiiX * oneOverRadiiX;
  const y2 = positionY * positionY * oneOverRadiiY * oneOverRadiiY;
  const z2 = positionZ * positionZ * oneOverRadiiZ * oneOverRadiiZ;

  const squaredNorm = x2 + y2 + z2;
  const ratio = Math.sqrt(1.0 / squaredNorm);

  const intersection = new THREE.Vector3(
      positionX * ratio,
      positionY * ratio,
      positionZ * ratio
  );

  if (squaredNorm < centerToleranceSquared) {
      return !isFinite(ratio) ? undefined : intersection;
  }

  const oneOverRadiiSquaredX = oneOverRadiiSquared.x;
  const oneOverRadiiSquaredY = oneOverRadiiSquared.y;
  const oneOverRadiiSquaredZ = oneOverRadiiSquared.z;

  const gradient = new THREE.Vector3();
  gradient.x = intersection.x * oneOverRadiiSquaredX * 2.0;
  gradient.y = intersection.y * oneOverRadiiSquaredY * 2.0;
  gradient.z = intersection.z * oneOverRadiiSquaredZ * 2.0;

  let lambda = ((1.0 - ratio) * cartesian.length()) / (0.5 * gradient.length());
  let correction = 0.0;

  let func;
  let denominator;
  let xMultiplier;
  let yMultiplier;
  let zMultiplier;
  let xMultiplier2;
  let yMultiplier2;
  let zMultiplier2;
  let xMultiplier3;
  let yMultiplier3;
  let zMultiplier3;

  do {
      lambda -= correction;

      xMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredX);
      yMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredY);
      zMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquaredZ);

      xMultiplier2 = xMultiplier * xMultiplier;
      yMultiplier2 = yMultiplier * yMultiplier;
      zMultiplier2 = zMultiplier * zMultiplier;

      xMultiplier3 = xMultiplier2 * xMultiplier;
      yMultiplier3 = yMultiplier2 * yMultiplier;
      zMultiplier3 = zMultiplier2 * zMultiplier;

      func = x2 * xMultiplier2 + y2 * yMultiplier2 + z2 * zMultiplier2 - 1.0;

      denominator = x2 * xMultiplier3 * oneOverRadiiSquaredX +
                    y2 * yMultiplier3 * oneOverRadiiSquaredY +
                    z2 * zMultiplier3 * oneOverRadiiSquaredZ;

      const derivative = -2.0 * denominator;

      correction = func / derivative;
  } while (Math.abs(func) > 1e-12);

  return new THREE.Vector3(
      positionX * xMultiplier,
      positionY * yMultiplier,
      positionZ * zMultiplier
  );
}

export {renderer};

@gkjohnson
Copy link
Contributor

I'm happy to help but please make things easy to understand. Pasting 300 lines of code with no imports or data, as I asked for, is not something I can to spend time digging in to to understand or figure out how to run.

If you're just asking how to extract a coordinate frame from a lat / lon position you can use the Ellipsoid.getEastNorthUpAxes or getEastNorthUpFrame function to get a matrix with orthogonal vectors for the eastern and northern directions in X and Y respectively. These functions aren't currently documented but a WGS84_ELLIPSOID object is exported with these functions and more used for the globe functionality.

Regarding the scale differences this is going to depend on the 2d projection of your map. The scale will change depending on the location of the tile set due to the projection but that's outside the scope of this project.

Otherwise if you'd like help please make a simple, live code example using something like codesandbox or jsfiddle including all the data necessary and the specific portion of code you have a question about made clear.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
question Further information is requested
Projects
None yet
Development

No branches or pull requests

2 participants