Forum

Author Topic: Setting mapping region  (Read 32065 times)

jmsuomal

  • Newbie
  • *
  • Posts: 17
    • View Profile
Setting mapping region
« on: January 09, 2014, 06:19:32 PM »
I would like to set the mapping region with Pytho scripting, but I seem to have problems because I dont understand the chunk.region object.

Here is what I have managed to do so far:
Code: [Select]
# Define: Region [CenterX, CenterY, SizeX, SizeY, ZRotationDeg]
MapRegionVec = [304331.4, 5740607.3, 173.8, 143.6, -25.0]

# Define: Coordinate system
CoordinateSystemEPSG = "EPSG::32632"

# create chunk...

# Set the coordinate system
CoordinateSystem = PhotoScan.CoordinateSystem()
CoordinateSystem.init(CoordinateSystemEPSG)
chunk.crs = CoordinateSystem
chunk.projection = CoordinateSystem

# (load camera calib...)
# (load photos...)
# (load initial camera orientations...)
# (matchPhotos...)
# (alignPhotos...)

# Set region center:
newregion = PhotoScan.Region()
centerUTM = PhotoScan.Vector([MapRegionVec[0], MapRegionVec[1], chunk.region.center[2]])
centerGEO = chunk.crs.unproject(centerUTM)
centerLocal = chunk.crs.localframe(centerGEO)
# LocalCoordinates are a 4x4 matrix instead of XYZ vector!
# I cant put that to XYZ. I have no idea what to do so I use the original center point for now...
newregion.center = chunk.region.center

# Set region size (this one works fine)
newregion.size = PhotoScan.Vector([MapRegionVec[2], MapRegionVec[3], chunk.region.size[2]])

# Set region rotation
import math
SinRotZ= math.sin(math.radians(MapRegionVec[4]))
CosRotZ= math.cos(math.radians(MapRegionVec[4]))
newregion.rot = PhotoScan.Matrix( [[SinRotZ,-CosRotZ,0], [CosRotZ,SinRotZ,0], [0,0,1]] )
# This works technically but does not do what I want!

# put newregion to chunk
chunk.region = newregion

But on that I cannot get the chunk.region.center or chunk.region.rot configured correctly. Thus I have a couple open questions.

- What are the units/coordinatesystem in the chunk.region.center values?  My chunk coordinate system is WGS84/UTM 32. Here's how the current center numbers look:

Code: [Select]
>>> chunk.region.center
Vector([-4.059450202416892, -0.8626608305524442, 12.09349905182462])
Z could be quite sensible value for altitude, but what are those X and Y units.  They are not UTM coordinates for sure. How do I convert my UTM coordinates to these?

- What is this chunk.region.rot object? Rotation matrix? Using standard z-rotation matrix of:

[sin(RotAngleRad),-cos(RotAngleRad),0],
[cos(RotAngleRad),sin(RotAngleRad),0],
[0,0,1]

does not seem to do what I would assume it does. For example running:

Code: [Select]
region = chunk.region
region.rot = PhotoScan.Matrix([[1,0,0],[0,1,0],[0,0,1]])
chunk.region = region

...does not put the region along the XY axis, but rotates it in an angle (that might be based to the orientation of first photo?) How do I set the region to certain Kappa angle?

Any hints or even an example script how to do this?

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14004
    • View Profile
Re: Setting mapping region
« Reply #1 on: January 09, 2014, 06:32:15 PM »
Hello jmsuomal,

Chunk.region.size and center are in the internal units referred to the internal chunk coordinate system.

So, for example, to convert chunk.region.center vector from this system you need to use chunk.transform matrix:

Code: [Select]
c = chunk.region.center
c.size = 4
c.w = 1
c_temp = chunk.transform * c
c_temp.size = 3
c_geo = chunk.crs.project(c_temp)

and to get vector V from geographic coordinates to the internal chunk coordinate system you need to perform inverted workflow:

Code: [Select]
v_temp = chunk.crs.unproject(v)
v_temp.size = 4
v_temp.w = 1
v_internal = chunk.transform.inv() * v_temp
v_internal.size = 3
« Last Edit: January 09, 2014, 06:34:14 PM by Alexey Pasumansky »
Best regards,
Alexey Pasumansky,
Agisoft LLC

jmsuomal

  • Newbie
  • *
  • Posts: 17
    • View Profile
Re: Setting mapping region
« Reply #2 on: January 09, 2014, 06:42:33 PM »
Thanks Alexey,

That solved the center coordinates. :)

Any suggestion how I should handle the rotation?

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14004
    • View Profile
Re: Setting mapping region
« Reply #3 on: January 09, 2014, 06:54:04 PM »
Hello,

Somewhere on forum I have already posted the code that rotates the bounding box in accordance of the chunk coordinate system, so you need to use it and for additional rotations multiply the matrix (m) by rotation 3x3 matrices:

Code: [Select]
T = chunk.transform
v = PhotoScan.Vector( [0,0,0,1] )
v_t = T * v
v_t.size = 3

m = chunk.crs.localframe(v_t)
m = m * T
s = math.sqrt(m[0,0]**2 + m[0,1]**2 + m[0,2]**2) #scale factor
R = PhotoScan.Matrix( [[m[0,0],m[0,1],m[0,2]], [m[1,0],m[1,1],m[1,2]], [m[2,0],m[2,1],m[2,2]]])

R = R * (1. / s)
reg = chunk.region
reg.rot = R.t()
chunk.region = reg

If chunk.transform and chunk.crs are None, you should use identity matrix instead them.
Best regards,
Alexey Pasumansky,
Agisoft LLC

jmsuomal

  • Newbie
  • *
  • Posts: 17
    • View Profile
Re: Setting mapping region
« Reply #4 on: January 10, 2014, 11:37:07 AM »
Thanks Alexey. You are the best! I definitely would not have been able to work that out myself.

For other readers of the forum: here's my final working code:

Code: [Select]
# Define: Map region
#MapRegionVec = [CenterX, CenterY, SizeX0, SizeY0, ZRotationDeg]
MapRegionVec = [304331.4, 5740607.3, 173.8, 143.6, -25.0]

# ALIGN PHOTOS
chunk.alignPhotos()

# SET MAPPING REGION (BOUNDING BOX)
# create new region object
newregion = PhotoScan.Region()
# Set region center:
centerUTM = PhotoScan.Vector([MapRegionVec[0], MapRegionVec[1], 0])
centerGEO = chunk.crs.unproject(centerUTM)
centerGEO.size = 4
centerGEO.w = 1
centerLocal = chunk.transform.inv() * centerGEO
newregion.center = PhotoScan.Vector([centerLocal[0], centerLocal[1], chunk.region.center[2]])
# Set region size
newregion.size = PhotoScan.Vector([MapRegionVec[2], MapRegionVec[3], chunk.region.size[2]])
# Set region rotation
# build rotation matrix in our coordinate system
RotZDeg = -MapRegionVec[4]
SinRotZ= math.sin(math.radians(RotZDeg))
CosRotZ= math.cos(math.radians(RotZDeg))
RotMat = PhotoScan.Matrix([[CosRotZ,-SinRotZ,0,0],[SinRotZ,CosRotZ,0,0],[0,0,1,0],[0,0,0,1]])
#  rotate region bounding box
T = chunk.transform
v = PhotoScan.Vector( [0,0,0,1] )
v_t = T * v
v_t.size = 3
m = chunk.crs.localframe(v_t)
m = m * T
m = RotMat*m
s = math.sqrt(m[0,0]**2 + m[0,1]**2 + m[0,2]**2) #scale factor
R = PhotoScan.Matrix( [[m[0,0],m[0,1],m[0,2]], [m[1,0],m[1,1],m[1,2]], [m[2,0],m[2,1],m[2,2]]])
R = R * (1. / s)
newregion.rot = R.t()
# put newregion to chunk
chunk.region = newregion

ppkong

  • Jr. Member
  • **
  • Posts: 54
    • View Profile
Re: Setting mapping region
« Reply #5 on: April 29, 2014, 11:30:04 AM »
Thanks,Very much!
« Last Edit: May 04, 2014, 11:03:21 AM by ppkong »

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14004
    • View Profile
Re: Setting mapping region
« Reply #6 on: May 04, 2014, 11:23:19 AM »
Hello ppkong,

The script from above (http://www.agisoft.ru/forum/index.php?topic=1886.msg10045#msg10045) orients bounding box in accordance of the coordinate system applied to the chunk.
Best regards,
Alexey Pasumansky,
Agisoft LLC

ppkong

  • Jr. Member
  • **
  • Posts: 54
    • View Profile
Re: Setting mapping region
« Reply #7 on: May 04, 2014, 11:42:13 AM »
Thanks,Very much!
And I have another question. I set the region size with 500*500(meters),and the result region become 475*520,I want to know why?
« Last Edit: May 04, 2014, 02:22:22 PM by ppkong »

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14004
    • View Profile
Re: Setting mapping region
« Reply #8 on: May 05, 2014, 11:36:50 AM »
Hello ppkong,

Could you please send the project file and the script used on support@agisoft.ru for the further investigation?

Please also specify how do you measure the size of the region?
Best regards,
Alexey Pasumansky,
Agisoft LLC

ppkong

  • Jr. Member
  • **
  • Posts: 54
    • View Profile
Re: Setting mapping region
« Reply #9 on: May 05, 2014, 12:06:08 PM »
Hello ppkong,

Could you please send the project file and the script used on support@agisoft.ru for the further investigation?

Please also specify how do you measure the size of the region?
I'm sorry, I found out that I had made a mistake in the calculation process leading to incorrect results. I've solved the problem.
Thank you so much!

theYesMan

  • Newbie
  • *
  • Posts: 8
    • View Profile
Re: Setting mapping region
« Reply #10 on: May 13, 2014, 08:03:53 AM »
Hi Alexey,

>> If chunk.transform and chunk.crs are None, you should use identity matrix instead them.

Would you provide an example for doing this?
My chunk.crs are None.

Thanks.

gatsri

  • Jr. Member
  • **
  • Posts: 77
    • View Profile
Re: Setting mapping region
« Reply #11 on: July 08, 2014, 11:10:56 AM »
how can i implement this code
Quote
T = chunk.transform
v = PhotoScan.Vector( [0,0,0,1] )
v_t = T * v
v_t.size = 3

m = chunk.crs.localframe(v_t)
m = m * T
s = math.sqrt(m[0,0]**2 + m[0,1]**2 + m[0,2]**2) #scale factor
R = PhotoScan.Matrix( [[m[0,0],m[0,1],m[0,2]], [m[1,0],m[1,1],m[1,2]], [m[2,0],m[2,1],m[2,2]]])

R = R * (1. / s)
reg = chunk.region
reg.rot = R.t()
chunk.region = reg

into a MenuItem?! so i can run the script from there?

Thanks for the help!

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14004
    • View Profile
Re: Setting mapping region
« Reply #12 on: July 08, 2014, 11:13:51 AM »
Hello ristag,

Save the file as *.py as plain text (don't forget to add import PhotoScan, math line in the beginning of the code), then use Run Script option in Tools Menu or the corresponding button on the Console pane.
Best regards,
Alexey Pasumansky,
Agisoft LLC

gatsri

  • Jr. Member
  • **
  • Posts: 77
    • View Profile
Re: Setting mapping region
« Reply #13 on: July 08, 2014, 11:26:45 AM »
thanks for the fast answer...

but I would like to run the script with a MenuItem like in the picture:
where i can run the script from there

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14004
    • View Profile
Re: Setting mapping region
« Reply #14 on: July 08, 2014, 11:31:49 AM »
Hello ristag,

Then you should use something similar to the following:

Code: [Select]
import PhotoScan, math

def rotateBb():

doc = PhotoScan.app.document
chunk = doc.activeChunk

T = chunk.transform
v = PhotoScan.Vector( [0,0,0,1] )
v_t = T * v
v_t.size = 3

m = chunk.crs.localframe(v_t)
m = m * T
s = math.sqrt(m[0,0]**2 + m[0,1]**2 + m[0,2]**2) #scale factor
R = PhotoScan.Matrix( [[m[0,0],m[0,1],m[0,2]], [m[1,0],m[1,1],m[1,2]], [m[2,0],m[2,1],m[2,2]]])

R = R * (1. / s)
reg = chunk.region
reg.rot = R.t()
chunk.region = reg

PhotoScan.app.addMenuItem("Custom menu/Bounding box to coordinate system", rotateBb)
Best regards,
Alexey Pasumansky,
Agisoft LLC