Forum

Author Topic: How can I compute the area that a camera is covering?  (Read 21840 times)

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #15 on: March 31, 2015, 05:03:37 PM »
Hello Mark,

Actually Python in PhotoScan is quite slow, so looping across millions of polygons will take a long time.
Best regards,
Alexey Pasumansky,
Agisoft LLC

AndyC

  • Newbie
  • *
  • Posts: 30
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #16 on: April 27, 2015, 01:32:28 PM »
I've a few questions for the latest script (for V1.1.0) (apologies in advance if they're rather basic, I'm rather new to Python).

(1) Is the user input 'camera index' the total number of ‘cameras’ (images) for which footprints should be exported??

(2) When specifying the path for the output file, is is necessary to specify the format? If so which formats are permissible?

(3) Should the 'step' parameter be adjusted for the coordinate system? (e.g., is using a coordinate system in m, with camera footprints of ~30 m, would it be best to reduce the step from 100 down to perhaps 1?)

(4) Will each camera footprint be exported as an individual unit (e.g., multiple polygons in a single shapefile), or will they be amalgamated?

Thanks, Andy

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #17 on: April 27, 2015, 02:20:48 PM »
Hello Andy,

1) Camera index starts from zero, so the last camera in the chunk of N cameras should have (N - 1) index. Using N is incorrect and may crash the script or lead to the unexpected results.

2) It's recommended to specify the format for your own convenience. I suggest to use .txt extension.

3) Step parameter is actually in pixels, going along the image frame perimeter from the top-left corner in clock-wise direction.

4) The output file is for the single camera only (for the one specified by the user) and includes the intersections of the rays going through camera sensor center and pixels on the image frame border.
Best regards,
Alexey Pasumansky,
Agisoft LLC

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #18 on: April 27, 2015, 09:46:28 PM »
Alexei,

it would very interesting to have script reporting area covered on mesh for every camera but giving only the coordinates of 4 corners for each camara frame...
Best Regards,
Paul Pelletier,
Surveyor

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #19 on: April 27, 2015, 11:07:42 PM »
Hello pap1956,

You can modify the script by removing "stepping" and leaving only four courners of the image. For multiple photos it's necessary to add a loop like the folloiwng:
Code: [Select]
for camera in chunk.cameras
Best regards,
Alexey Pasumansky,
Agisoft LLC

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #20 on: April 29, 2015, 07:30:56 PM »
Alexeiu,

i did put in a loop to loop thru all cameras in chunk and write image footprints to file but I got following error:

 File "C:/Users/paul.pelletier/Documents/Photoscan/ImageFootprint_4corners.py", line 26
    for camera in chunk.cameras
                              ^
SyntaxError: invalid syntax

Can you help me?

Best Regards,
Paul Pelletier,
Surveyor

AndyC

  • Newbie
  • *
  • Posts: 30
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #21 on: April 29, 2015, 07:38:38 PM »
Pap1956, probably not the cause of the error, but is your script trying to query the projected coordinates of point 0,0 twice?
« Last Edit: April 29, 2015, 07:40:14 PM by AndyC »

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #22 on: April 29, 2015, 07:49:58 PM »
Pap1956, probably not the cause of the error, but is your script trying to query the projected coordinates of point 0,0 twice?
Yes as I want to output the 4 corners of each camera plus repeated 1st point so to import polygons in Global mapper
Best Regards,
Paul Pelletier,
Surveyor

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #23 on: April 29, 2015, 10:54:29 PM »
Hello pap1956,

Line with "for" statement should end with the semicolon ":".

Also if you want to include the corner export loop for each camera, you need to add indent for all the lines of this loop, to make it included to the upper level loop.
Best regards,
Alexey Pasumansky,
Agisoft LLC

Paulo

  • Hero Member
  • *****
  • Posts: 1303
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #24 on: April 30, 2015, 04:26:07 PM »
Thanks for your help, Alexei...

t works but is very slow due to density of mesh (million faces). As I just want an aproximate view of image coverage on ground, I import a crude mesh from ASTER DEM into PS and the script runs fast and I get a text file with a polygon (5 vertices) for each camera....
Best Regards,
Paul Pelletier,
Surveyor

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #25 on: April 30, 2015, 04:33:51 PM »
Hello pap1956,

Script speed depend on the number of faces in the mesh model.

The fifth vertex of the exported polygon, according to the script, should be the same as the first one, as (0,0) point is used twice.
Best regards,
Alexey Pasumansky,
Agisoft LLC

simon_29

  • Newbie
  • *
  • Posts: 10
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #26 on: June 28, 2016, 11:07:13 AM »
Hello everybody !

I would like to compute the area covered for each of my cameras. I have already tried to use the code given as an example. I added the loop to go through each camera of my chunk. I also tried to select only the four corners. As a result i got ... nothing. It seems like it doesn't find any intersections with the faces of my mesh. Therefore I tried to add some pixels.  I got some results but it is not satisfying as I don't get the vertices positions corresponding to my corners.  I attached the txt file that I have at the end so taht you can see for yourself.

During the process I printed the position of the pixel I am using. It seems like I go through all of my corners [(0,0), (6015,0), (6015,399),(0,3999)] as my sensor.width=6016 and my sensor.height=4000. Nevertheless I don't get any vertices coordinates.

It is always for the same pixels that I manage to have intersections but never the four corners. I don't know what I am doing wrong. I barely changed the code from Alexey. I have tried many things but I don't have any ideas left. This is why I would like to know if there is anybody who met the same problem.

Here is the code:

Code: [Select]
import time
import PhotoScan

def cross(a, b):
result = PhotoScan.Vector([a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y *b.x])
return result

print("Script started")


#cam_index = PhotoScan.app.getInt("Input camera index (starting from zero): ") Manual selection of the camera
save_path = PhotoScan.app.getSaveFileName("Specify output file:")


t0 = time.time()
file = open(save_path, "wt")
file.write('FileName    Pixel x    Pixel y    Vertex Lon    Vertex Lat    Vertex Alt\n') # Header
doc = PhotoScan.app.document
chunk = doc.chunk
model = chunk.model
faces = model.faces
vertices = model.vertices

for camera in chunk.cameras:
    sensor = camera.sensor
    print(camera) #camera label

    step = 1000 #bigger value - faster processing.
    steps = list(zip(list(range(0, sensor.width - 1, step)), [0]*((sensor.width - 1)// step)))
    steps.extend( list(zip([sensor.width - 1]*((sensor.height - 1) // step), list(range(0, sensor.height - 1, step)))) )
    steps.extend( list(zip(list(range((sensor.width - 1), 0, -step)), [sensor.height - 1]*((sensor.width - 1)// step))))
    steps.extend( list(zip([0]*((sensor.height - 1) // step), list(range(sensor.height - 1, 0, -step)))) )
   
# Selection of the four corners:   
    #ltop_corner=PhotoScan.Vector([0, 0]) # left top corner
    #rtop_corner=PhotoScan.Vector([sensor.width - 1, 0]) # right top corner
    #rbottom_corner=PhotoScan.Vector([sensor.width - 1, sensor.height - 1]) # right bottom corner
    #lbottom_corner=PhotoScan.Vector([0, sensor.height - 1]) # left bottom corner
    #
    # List of the four corners
    #steps=[ltop_corner,rtop_corner,rbottom_corner,lbottom_corner]     
   
    print(steps)
 
    for x,y in steps:
        point = PhotoScan.Vector([x, y])
        point = sensor.calibration.unproject(point)
        point = camera.transform.mulv(point)
        vect = point
        p = PhotoScan.Vector(camera.center)
        for face in faces:
            v = face.vertices
            E1 = PhotoScan.Vector(vertices[v[1]].coord - vertices[v[0]].coord)
            E2 = PhotoScan.Vector(vertices[v[2]].coord - vertices[v[0]].coord)
            D = PhotoScan.Vector(vect)
            T = PhotoScan.Vector(p - vertices[v[0]].coord)
            P = cross(D, E2)
            Q = cross(T, E1)
            result = PhotoScan.Vector([Q * E2, P * T, Q * D]) / (P * E1)
            if (0 < result[1]) and (0 < result[2]) and (result[1] + result[2] <= 1):
                t = (1 - result[1] - result[2]) * vertices[v[0]].coord
                u = result[1] * vertices[v[1]].coord
                v_ = result[2] * vertices[v[2]].coord
               
                res = chunk.transform.matrix.mulp(u + v_ + t)
                res = chunk.crs.project(res)
               
                #file.write( "{:>04d}".format(x + 1) + "\t" + "{:04d}".format(y + 1) + "\t" + "{:.8f}".format(res[0]) + "\t" + "{:.8f}".format(res[1]) + "\t" + "{:.4f}".format(res[2]) + "\n")
                file.write("%s""\t""%d""\t""%d""\t""%.08f""\t""%.08f""\t""%.3f \n" % (camera.label, x+1,y+1,res[0],res[1],res[2]))
                break #finish when the first intersection is found

file.close()
t1 = time.time()
t1 -= t0
t1 = float(t1)
print("Script finished in " + "{:.2f}".format(t1) + " seconds.")

I am not sure I have understood everything about the code, especially this step:

Code: [Select]
result = PhotoScan.Vector([Q * E2, P * T, Q * D]) / (P * E1)
 if (0 < result[1]) and (0 < result[2]) and (result[1] + result[2] <= 1):
                t = (1 - result[1] - result[2]) * vertices[v[0]].coord
                u = result[1] * vertices[v[1]].coord
                v_ = result[2] * vertices[v[2]].coord

If someone feels like explaining it to me, I would be very grateful. I am quite new to Photoscan so I don't really know if I should post my question in this topic or somewhere else.

Best wishes,

Simon



Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #27 on: June 28, 2016, 11:10:41 AM »
Hello Simon,

Maybe you can send the PSZ file with the alignment results and model (can be decimated) so that we could reproduce the problem (no need for dense cloud or depth maps in the project)?
Best regards,
Alexey Pasumansky,
Agisoft LLC

simon_29

  • Newbie
  • *
  • Posts: 10
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #28 on: June 28, 2016, 12:19:52 PM »
Hello Alexey, first thank you very much for such a quick answer! I have sent the psz file to support@agisoft.com. You have probably received a download link. Just tell me if I need to sent it back.

Simon

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 14813
    • View Profile
Re: How can I compute the area that a camera is covering?
« Reply #29 on: June 28, 2016, 01:32:32 PM »
Hello Simon,

Thank you for providing the project files.

Are you sure that the rays through the corners of the input images do really intersect with the mesh model? There are no image thumbnails in the project, but I've tried to put the markers in the very corners of the images (using Create Marker option) and haven't got any other automatic projections created.
Best regards,
Alexey Pasumansky,
Agisoft LLC