Forum

Author Topic: Script for exporting orthophotos using markers  (Read 70394 times)

_mARCel_

  • Newbie
  • *
  • Posts: 22
  • I'm an archaeologist from Germany
    • View Profile
    • ArchaeoBW GmbH Home
Re: Script for exporting orthophotos using markers
« Reply #60 on: August 28, 2018, 05:13:17 PM »
Hi Dukytony,

here is our script for 1.4.1.

Best regards,
Marcel



Code: [Select]
"""
/********************************************************************************************************************
                                    OPEX - Orthophotoexport (V. 1.10)

                                                    Content:
  Processing code to exports orthophoto in planar projection defined by vektor (from three markers) and Z axes.
  This code is mainly for exporting archaeological profiles; ready for importing in GIS.
  Before using the script you need:
  - A georeferenzed Model
  - An "A" and "B" Profile Nail (Named "A" and "B" in Photoscan)
  - A "Y" Point in Photoscan with the coordinates (X-coord from "A",Y-coord from "A", Z-coord from "A"+ a high number)

  Compatibility - PhotoScan Professional 1.4.1
  Version - 1.10

-------------------
written by         :   Alexey Pasumansky, AgiSoft LLC,
                                              Marcel C. Hagner, ArchaeoBW GmbH
      Florian Tubbesing, ArchaeoBW GmbH
copyright         : 2018 ArchaeoBW GmbH
email                : info@archaeobw.de
 ********************************************************************************************************************/

/*************************************************************************************
 *   This program is free software; you can redistribute it and/or modify it under   *
 *   under the terms of a Creative Commons Attribution-NonCommercial-ShareAlike 4.0  *
 *   International License.                                                          *
 *                                                                                   *
 *   Attribution - You must give a appropriate credit, provide a link to the license,*
 * and indicate of changes were made. You may do so in any reasonable manner, but  *
 * not in any way that suggests the licensor endorses you or your use              *
 *   *
 * NonCommercial - You may not use the material for commercial purposes.           *
 * *
 * ShareAlike - If you remix, transform, or build upon the material, you must      *
 * distribute your contributions under the same license as the original.           *
 *************************************************************************************/
"""


import os
import PhotoScan
from PySide2 import QtCore, QtWidgets

def getMarker(chunk, label):

for marker in chunk.markers:
if marker.label.upper() == label.upper():
return marker
return 0

def vect(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.normalized()

class ExportOFDlg(QtWidgets.QDialog):

def __init__ (self, parent):
QtWidgets.QDialog.__init__(self, parent)

chunk = doc.chunk
self.supported_types = ["jpg", "tif", "png"]
self.blending = PhotoScan.BlendingMode.AverageBlending


self.setWindowTitle("OPEX")

self.btnQuit = QtWidgets.QPushButton("&Exit")
self.btnQuit.setFixedSize(150,50)

self.btnP1 = QtWidgets.QPushButton("&Export")
self.btnP1.setFixedSize(150,50)

self.resTxt = QtWidgets.QLabel()
self.resTxt.setText("Export resolution (m):")
self.resTxt.setFixedSize(150, 25)

self.vectTxt = QtWidgets.QLabel()
self.vectTxt.setText("Horizontal vector:")
self.vectTxt.setFixedSize(150, 25)

self.vect2Txt = QtWidgets.QLabel()
self.vect2Txt.setText("Vertical vector:")
self.vect2Txt.setFixedSize(150, 25)

self.resEdt = QtWidgets.QLineEdit()
self.resEdt.setText("0.00025")
self.resEdt.setFixedSize(150, 25)

self.radioBtn_v = QtWidgets.QRadioButton("Vertical")
self.radioBtn_h = QtWidgets.QRadioButton("Horizontal")
self.radioBtn_v.setChecked(True)
self.radioBtn_h.setChecked(False)

self.llistV = [None, None]

self.llistV[0] = QtWidgets.QComboBox()
self.llistV[0].resize(150, 30)
self.llistV[1] = QtWidgets.QComboBox()
self.llistV[1].resize(150, 30)

self.llistH = [None, None]

self.llistH[0] = QtWidgets.QComboBox()
self.llistH[0].resize(150, 30)
self.llistH[1] = QtWidgets.QComboBox()
self.llistH[1].resize(150, 30)


self.radio = QtWidgets.QVBoxLayout()
self.radio.addWidget(self.radioBtn_v)
self.radio.addWidget(self.radioBtn_h)

for marker in chunk.markers:
for ilist in self.llistV:
ilist.addItem(marker.label)
for ilist in self.llistH:
ilist.addItem(marker.label)


layout = QtWidgets.QGridLayout()   #creating layout
layout.setSpacing(10)
layout.addWidget(self.resTxt, 0, 0)
layout.addWidget(self.resEdt, 0, 1)

layout.addWidget(self.vectTxt, 1, 0)
layout.addWidget(self.vect2Txt, 2, 0)
layout.addWidget(self.llistH[0], 1, 1)
layout.addWidget(self.llistH[1], 1, 2)
layout.addWidget(self.llistV[0], 2, 1)
layout.addWidget(self.llistV[1], 2, 2)
layout.addLayout(self.radio, 3, 0)

layout.addWidget(self.btnP1, 3, 1)
layout.addWidget(self.btnQuit, 3, 2)
self.setLayout(layout) 

proc_export = lambda : self.procExport()

QtCore.QObject.connect(self.btnP1, QtCore.SIGNAL("clicked()"), proc_export)
QtCore.QObject.connect(self.btnQuit, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("reject()"))

self.exec()

def procExport(self):

chunk = doc.chunk

path = PhotoScan.app.getExistingDirectory("Please choose the folder with your psx files.")
crs = PhotoScan.app.getCoordinateSystem("Please choose the Output Coordinate System.", )
filelist = os.listdir(path)
doclist = list()
for file in filelist:
if file[-3:].upper() == "PSX":
doclist.append(path + "\\" + file)


try:
d_x = d_y = float(self.resEdt.text())
except(ValueError):
PhotoScan.app.messageBox("Incorrect export resolution! Please use point delimiter.\n")
print("Script aborted.")
return 0

if self.llistV[0].currentIndex() == self.llistV[1].currentIndex():
PhotoScan.app.messageBox("Can't use the same marker for vector start and end.\n")
print("Script aborted.")
return 0

if self.llistH[0].currentIndex() == self.llistH[1].currentIndex():
PhotoScan.app.messageBox("Can't use the same marker for vector start and end.\n")
print("Script aborted.")
return 0

if len(chunk.markers) < 2:
PhotoScan.app.messageBox("No markers.\n")
print("Script aborted.")
return 0

print("script is starting.....")

type = "tif"

for doc_path in doclist:
try:
doc_unlock = PhotoScan.Document()
doc_unlock.open(doc_path, read_only=True)
doc_unlock.read_only = False
doc.open(doc_path, read_only=False)
chunk = doc.chunk
chunk.crs=crs
T = chunk.transform.matrix

path = doc_path[:-3] + type


self.btnP1.setDisabled(True)
self.btnQuit.setDisabled(True)

m1 = getMarker(chunk, self.llistH[0].currentText())
m2 = getMarker(chunk, self.llistH[1].currentText())

horizontal = m2.position - m1.position
horizontal = T.mulv(horizontal)

m1 = getMarker(chunk, self.llistV[0].currentText())
m2 = getMarker(chunk, self.llistV[1].currentText())

vertical = m2.position - m1.position
vertical = T.mulv(vertical)

normal = vect(vertical, horizontal)

if self.radioBtn_h.isChecked():
vertical = vect(horizontal, normal)
horizontal = horizontal.normalized()

elif self.radioBtn_v.isChecked():
vertical = vertical.normalized()
horizontal = -vect(vertical, normal)
#horizontal, vertical = -vertical, -horizontal

else:
PhotoScan.app.messageBox("Error! Can't export orthomosaic." + path)
print("Script aborted")
continue

m = -1
for marker in chunk.markers:#Choosing marker by name
m += 1
label = marker.label
if "A" in label:
a = m
R = PhotoScan.Matrix ([horizontal, vertical, -normal])
origin = T.mulp(chunk.markers[a].position)
A = m1.reference.location * PhotoScan.Vector([0,0,1]) #A is m1´s height
X = (-1) * R * origin

horizontal.size = 4
horizontal.w = X.x
vertical.size = 4
vertical.w = X.y + A # A is added
normal.size = 4
normal.w = -X.z

proj = PhotoScan.Matrix ([horizontal, vertical, -normal, PhotoScan.Vector([0,0,0,1])])

if doc.path[-3:].lower() != "psx":  #saving in PSX required to build orthomosaic
doc.save(doc.path[:-3] + "psx")
chunk = doc.chunk

chunk.buildOrthomosaic(surface = PhotoScan.DataSource.ModelData, blending = self.blending, projection = proj, dx = d_x, dy = d_y)
doc.save()
chunk.exportOrthomosaic(path, dx = d_x, dy = d_y)

self.btnP1.setDisabled(False)
self.btnQuit.setDisabled(False)
PhotoScan.app.update()

except:
continue
print("Finished process")
print()
return 1



def main():

global doc
doc = PhotoScan.app.document

app = QtWidgets.QApplication.instance()
parent = app.activeWindow()

dlg = ExportOFDlg(parent)

PhotoScan.app.addMenuItem("ArchaeoBW Toolbox/OPEX", main)

matt07lx

  • Newbie
  • *
  • Posts: 30
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #61 on: July 16, 2019, 01:42:19 AM »
Can anyone confirm if this works in 1.5? I've tried running it and it just freezes upon trying to build the orthomosaic.

_mARCel_

  • Newbie
  • *
  • Posts: 22
  • I'm an archaeologist from Germany
    • View Profile
    • ArchaeoBW GmbH Home
Re: Script for exporting orthophotos using markers
« Reply #62 on: July 16, 2019, 08:15:58 AM »
Hi matt07lx,

we are currently using Metashape 1.5.2 and it works just fine. I suppose you would like to have the current code from us? I ask my colleague if I may insert it.  ;)

Best regards,
Marcel

Florian T.

  • Newbie
  • *
  • Posts: 10
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #63 on: July 16, 2019, 10:57:56 AM »
Hi matt07lx,

as Marcel mentioned before we use the following script for Metashape 1.5.3. and it works fine.
We hope that it´s usefull for you.

Best regards,

Florian

Code: [Select]
"""
/********************************************************************************************************************
                                    OPEX - Orthophotoexport (V. 2.0)

                                                    Content:
  Processing code to exports orthophoto in planar projection defined by vektor (from three markers) and Z axes.
  This code is mainly for exporting archaeological profiles; ready for importing in GIS.
  Before using the script you need:
  - A georeferenzed Model
  - An "A" and "B" Profile Nail (Named "A" and "B" in Photoscan)
  - A "Y" Point in Photoscan with the coordinates (X-coord from "A",Y-coord from "A", Z-coord from "A"+ a high number)

  Compatibility - Metashape Professional 1.5.3
  Version - 2.0

-------------------
written by         :   Alexey Pasumansky, AgiSoft LLC,
                                              Marcel C. Hagner, ArchaeoBW GmbH
      Florian Tubbesing, ArchaeoBW GmbH
copyright          :  2019 ArchaeoBW GmbH
email                 :   info@archaeobw.de
 ********************************************************************************************************************/

/*************************************************************************************
 *   This program is free software; you can redistribute it and/or modify it under                 *
 *   the terms of a Creative Commons Attribution-NonCommercial-ShareAlike 4.0               *
 *   International License.                                                                                                              *
 *                                                                                                                                                      *
 *   Attribution - You must give a appropriate credit, provide a link to the license,                *
 * and indicate of changes were made. You may do so in any reasonable manner, but   *
 * not in any way that suggests the licensor endorses you or your use                            *
 *     *
 * NonCommercial - You may not use the material for commercial purposes.            *
 *   *
 * ShareAlike - If you remix, transform, or build upon the material, you must       *
 * distribute your contributions under the same license as the original.            *
 *************************************************************************************/
"""
import os
import Metashape
from PySide2 import QtCore, QtWidgets

def getMarker(chunk, label):

for marker in chunk.markers:
if marker.label.upper() == label.upper():
return marker
return 0

def vect(a, b):

result = Metashape.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.normalized()

class ExportOFDlg(QtWidgets.QDialog):

def __init__ (self, parent):
QtWidgets.QDialog.__init__(self, parent)

chunk = doc.chunk
self.supported_types = ["jpg", "tif", "png"]
self.blending = Metashape.BlendingMode.AverageBlending


self.setWindowTitle("OPEX")

self.btnQuit = QtWidgets.QPushButton("&Exit")
self.btnQuit.setFixedSize(150,50)

self.btnP1 = QtWidgets.QPushButton("&Export")
self.btnP1.setFixedSize(150,50)

self.resTxt = QtWidgets.QLabel()
self.resTxt.setText("Export resolution (m):")
self.resTxt.setFixedSize(150, 25)

self.vectTxt = QtWidgets.QLabel()
self.vectTxt.setText("Horizontal vector:")
self.vectTxt.setFixedSize(150, 25)

self.vect2Txt = QtWidgets.QLabel()
self.vect2Txt.setText("Vertical vector:")
self.vect2Txt.setFixedSize(150, 25)

self.resEdt = QtWidgets.QLineEdit()
self.resEdt.setText("0.00025")
self.resEdt.setFixedSize(150, 25)

self.radioBtn_v = QtWidgets.QRadioButton("Vertical")
self.radioBtn_h = QtWidgets.QRadioButton("Horizontal")
self.radioBtn_v.setChecked(True)
self.radioBtn_h.setChecked(False)

self.llistV = [None, None]

self.llistV[0] = QtWidgets.QComboBox()
self.llistV[0].resize(150, 30)
self.llistV[1] = QtWidgets.QComboBox()
self.llistV[1].resize(150, 30)

self.llistH = [None, None]

self.llistH[0] = QtWidgets.QComboBox()
self.llistH[0].resize(150, 30)
self.llistH[1] = QtWidgets.QComboBox()
self.llistH[1].resize(150, 30)


self.radio = QtWidgets.QVBoxLayout()
self.radio.addWidget(self.radioBtn_v)
self.radio.addWidget(self.radioBtn_h)

for marker in chunk.markers:
for ilist in self.llistV:
ilist.addItem(marker.label)
for ilist in self.llistH:
ilist.addItem(marker.label)


layout = QtWidgets.QGridLayout()   #creating layout
layout.setSpacing(10)
layout.addWidget(self.resTxt, 0, 0)
layout.addWidget(self.resEdt, 0, 1)

layout.addWidget(self.vectTxt, 1, 0)
layout.addWidget(self.vect2Txt, 2, 0)
layout.addWidget(self.llistH[0], 1, 1)
layout.addWidget(self.llistH[1], 1, 2)
layout.addWidget(self.llistV[0], 2, 1)
layout.addWidget(self.llistV[1], 2, 2)
layout.addLayout(self.radio, 3, 0)

layout.addWidget(self.btnP1, 3, 1)
layout.addWidget(self.btnQuit, 3, 2)
self.setLayout(layout) 

proc_export = lambda : self.procExport()

QtCore.QObject.connect(self.btnP1, QtCore.SIGNAL("clicked()"), proc_export)
QtCore.QObject.connect(self.btnQuit, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("reject()"))

self.exec()

def procExport(self):

chunk = doc.chunk

path = Metashape.app.getExistingDirectory("Bitte wählen Sie den Ordner mit den gesammelten Geo-Cut-Modell aus")
crs = Metashape.app.getCoordinateSystem("Wählen Sie das Koordinatensystem für die Ausgabe aus.", )
filelist = os.listdir(path)
doclist = list()
for file in filelist:
if file[-3:].upper() == "PSX":
doclist.append(path + "\\" + file)


try:
d_x = d_y = float(self.resEdt.text())
except(ValueError):
Metashape.app.messageBox("Incorrect export resolution! Please use point delimiter.\n")
print("Script aborted.")
return 0

if self.llistV[0].currentIndex() == self.llistV[1].currentIndex():
Metashape.app.messageBox("Can't use the same marker for vector start and end.\n")
print("Script aborted.")
return 0

if self.llistH[0].currentIndex() == self.llistH[1].currentIndex():
Metashape.app.messageBox("Can't use the same marker for vector start and end.\n")
print("Script aborted.")
return 0

if len(chunk.markers) < 2:
Metashape.app.messageBox("No markers.\n")
print("Script aborted.")
return 0

print("Skript startet...")

type = "_OF.tif"

for doc_path in doclist:
#try:
doc_unlock = Metashape.Document()
doc_unlock.open(doc_path, read_only=True)
doc_unlock.read_only = False
doc.open(doc_path, read_only=False)
chunk = doc.chunk
chunk.crs=crs
T = chunk.transform.matrix

path = doc_path[:-4] + type


self.btnP1.setDisabled(True)
self.btnQuit.setDisabled(True)

m1 = getMarker(chunk, self.llistH[0].currentText())
m2 = getMarker(chunk, self.llistH[1].currentText())

horizontal = m2.position - m1.position
horizontal = T.mulv(horizontal)

m1 = getMarker(chunk, self.llistV[0].currentText())
m2 = getMarker(chunk, self.llistV[1].currentText())

vertical = m2.position - m1.position
vertical = T.mulv(vertical)

normal = vect(vertical, horizontal)

if self.radioBtn_h.isChecked():
vertical = vect(horizontal, normal)
horizontal = horizontal.normalized()

elif self.radioBtn_v.isChecked():
vertical = vertical.normalized()
horizontal = -vect(vertical, normal)
#horizontal, vertical = -vertical, -horizontal

else:
Metashape.app.messageBox("Fehler!Kann das OF nicht exportieren" + path)
print("Skript wird abgebrochen.")
continue

m = -1
for marker in chunk.markers:
m += 1
label = marker.label
if "A" in label:
a = m
R = Metashape.Matrix ([horizontal, vertical, -normal])
origin = T.mulp(chunk.markers[a].position)
A = m1.reference.location * Metashape.Vector([0,0,1]) #A ist der Höhenwert von m1
X = (-1) * R * origin

horizontal.size = 4
horizontal.w = X.x
vertical.size = 4
vertical.w = X.y + A # hier wird A hinzuaddiert
normal.size = 4
normal.w = -X.z

proj = Metashape.Matrix ([horizontal, vertical, -normal, Metashape.Vector([0,0,0,1])])

if doc.path[-3:].lower() != "psx":  #saving in PSX required to build orthomosaic
doc.save(doc.path[:-3] + "psx")
chunk = doc.chunk

chunk.buildOrthomosaic(surface = Metashape.DataSource.ModelData, blending = self.blending, fill_holes=True, projection = proj, dx = d_x, dy = d_y)
doc.save()
chunk.exportOrthomosaic(path, dx = d_x, dy = d_y)

self.btnP1.setDisabled(False)
self.btnQuit.setDisabled(False)
Metashape.app.update()
#print("Orthofoto wird exportiert nach to:\n" + path)
#except:
#continue
print("Skript rum.Viele Grüße MP & FT")
print()
return 1



def main():

global doc
doc = Metashape.app.document

app = QtWidgets.QApplication.instance()
parent = app.activeWindow()

dlg = ExportOFDlg(parent)

Metashape.app.addMenuItem("ArchaeoBW Toolbox/Profile/OPEX/OPEX ab 1.5.3", main)
« Last Edit: July 16, 2019, 11:19:26 AM by Florian T. »

matt07lx

  • Newbie
  • *
  • Posts: 30
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #64 on: July 16, 2019, 12:27:06 PM »
Florian,

Thanks so much for sharing the latest version of the script. When I run it now I get the following error: "UnboundLocalError: local variable 'a' referenced before assignment". Do you have any thoughts as to why I am getting this error?

Florian T.

  • Newbie
  • *
  • Posts: 10
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #65 on: July 19, 2019, 12:15:52 PM »
Hey matt07lx,

in my opinion this error occurs because you don´t have a nail with the name "A" which is needed as the starting point for the profile.
Please ceck the requirements of the script which are written on top of it.

Let me know if it helps.

Best regards,

Florian

#JanS#

  • Newbie
  • *
  • Posts: 1
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #66 on: August 24, 2020, 01:42:00 PM »
Hi all,

we found your script by chance and tried to use it with Metashape 1.6.2 but failed. Is it possible that something changed between Metashape 1.5 and 1.6 that prevents a proper usage?
If so, would it be possible to post an updated version or give us a hint, what we have to change within the script to use it?

Many thanks
Jan

PK_2509

  • Newbie
  • *
  • Posts: 13
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #67 on: December 04, 2024, 12:38:10 PM »
Hello,

is it possible to get an updated version of this script for 2.X?

I have tried to adapt the script for version 2.X myself - but without interface - only hardcode via the console. Unfortunately it doesn't work and I can't validate if the projection matrix fits.

I big difference between the API for 1.5.3 and 2.X is how the projection is used. In 1.5.3 it was still possible to simply pass a projection matrix. Metashape.OrthoProjection is required for the buildOrthomosaic command in 2.X.

So for this it tried this, postet by Alexey (https://www.agisoft.com/forum/index.php?topic=11767.msg52733#msg52733)

Code: [Select]
projection = Metashape.OrthoProjection()
projection.type = Metashape.OrthoProjection.Type.Planar
projection.matrix = proj #the matrix that is calculated based on markers
chunk.buildOrthomosaic(surface_type=Metashape.ModelData, blending = Metashape.MosaicBlending, fill_holes = True, projection = projection)

Overall there was an Ortho which was generated - but looked like turned 90° and also after export it dosent match.

Just to be sure, I'll post my code.

Code: [Select]
T=Metashape.app.document.chunk.transform.matrix
m1 = Metashape.app.document.chunk.markers[4].reference.location
m2 = Metashape.app.document.chunk.markers[5].reference.location
horizontal = m2 - m1
horizontal = T.mulv(horizontal)

m1 = Metashape.app.document.chunk.markers[4].reference.location
m2 = Metashape.app.document.chunk.markers[6].reference.location
vertical = m2 - m1
vertical = T.mulv(vertical)

a = vertical
b = horizontal
result = Metashape.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])
normal = result.normalized()

a = horizontal
b = normal
result = Metashape.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])
vertical = result.normalized()
horizontal = horizontal.normalized()
#horizontal, vertical = -vertical, -horizontal

R = Metashape.Matrix ([horizontal, vertical, -normal])
origin = T.mulp(Metashape.app.document.chunk.markers[4].reference.location)
A = Metashape.app.document.chunk.markers[4].reference.location * Metashape.Vector([0,0,1])
X = (-1) * R * origin

horizontal.size = 4
horizontal.w = X.x
vertical.size = 4
vertical.w = X.y + A
normal.size = 4
normal.w = -X.z


projection = Metashape.OrthoProjection()
projection.type = Metashape.OrthoProjection.Type.Planar
proj=Metashape.Matrix ([horizontal, vertical, -normal, Metashape.Vector([0,0,0,1])])
projection.matrix = proj

Metashape.app.document.chunk.buildOrthomosaic(
          surface_data=Metashape.ModelData,
          blending_mode=Metashape.MosaicBlending,
          projection = projection,
          fill_holes=True,
          resolution_x=0.001,
          resolution_y=0.001
      )

So Marker 4 to 5 are for horizontal and 4 to 6 for vertical.

Or maybe Alexy, can you something about the world_transform in exportRaster? I tried to manipulte the World-File during the export.
The description is 2x3 raster-to-world transformation matrix. But the behavoir is not clear to me - a ([1,0,0],[1,0,0]) results in a world-file ([1,1,0],[0,0.5,0.5]). This was just a test to see the behavoir.
I think there is also a matrix multiplication involved? How can i manipulate just the last entry for the Y value?

Many thanks in advance

PK_2509

  • Newbie
  • *
  • Posts: 13
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #68 on: December 06, 2024, 02:31:52 PM »
I was able to update the code to Version 2.X

The only problem left, is that it only works if i set the CRS to 4978. See Code below.

Code: [Select]
projection = Metashape.OrthoProjection()
projection.type = Metashape.OrthoProjection.Type.Planar
projection.matrix = proj
projection.crs=Metashape.CoordinateSystem("EPSG::4978")

Than I can import the Ortho to QGIS and set CRS to what i want over there - e.g. 32632.
But this means antoher not automated step :)

My Question, can we use an Input CRS during buildOrthomosaic?. Like it used to work before.
I tried to set it, but this screws up the projecion.
If anyone can help, I will implement it and post the updated version.

Paulo

  • Hero Member
  • *****
  • Posts: 1392
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #69 on: December 06, 2024, 08:58:37 PM »

Just to be sure, I'll post my code.

Code: [Select]
T=Metashape.app.document.chunk.transform.matrix
m1 = Metashape.app.document.chunk.markers[4].reference.location
m2 = Metashape.app.document.chunk.markers[5].reference.location
horizontal = m2 - m1
horizontal = T.mulv(horizontal)

m1 = Metashape.app.document.chunk.markers[4].reference.location
m2 = Metashape.app.document.chunk.markers[6].reference.location
vertical = m2 - m1
vertical = T.mulv(vertical)

a = vertical
b = horizontal
result = Metashape.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])
normal = result.normalized()

a = horizontal
b = normal
result = Metashape.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])
vertical = result.normalized()
horizontal = horizontal.normalized()
#horizontal, vertical = -vertical, -horizontal

R = Metashape.Matrix ([horizontal, vertical, -normal])
origin = T.mulp(Metashape.app.document.chunk.markers[4].reference.location)
A = Metashape.app.document.chunk.markers[4].reference.location * Metashape.Vector([0,0,1])
X = (-1) * R * origin

horizontal.size = 4
horizontal.w = X.x
vertical.size = 4
vertical.w = X.y + A
normal.size = 4
normal.w = -X.z


projection = Metashape.OrthoProjection()
projection.type = Metashape.OrthoProjection.Type.Planar
proj=Metashape.Matrix ([horizontal, vertical, -normal, Metashape.Vector([0,0,0,1])])
projection.matrix = proj

Metashape.app.document.chunk.buildOrthomosaic(
          surface_data=Metashape.ModelData,
          blending_mode=Metashape.MosaicBlending,
          projection = projection,
          fill_holes=True,
          resolution_x=0.001,
          resolution_y=0.001
      )


Hello PK_2509,

to get a correct Planar OrthoProjection, think you should change each instance of reference.location by position in above code.

Also to correctly the projection crs add following to code:
Code: [Select]
if Metashape.app.document.chunk.crs.geoccs:
    projection.crs = Metashape.app.document.chunk.crs.geoccs # case of georeferenced project
else:
    projection.crs = Metashape.app.document.chunk.crs

Hope this helps
« Last Edit: December 07, 2024, 05:01:54 AM by Paulo »
Best Regards,
Paul Pelletier,
Surveyor

PK_2509

  • Newbie
  • *
  • Posts: 13
    • View Profile
Re: Script for exporting orthophotos using markers
« Reply #70 on: December 09, 2024, 03:41:27 PM »
Hello Paulo,

thank you for your reply and mentioning that referene.location only means the input in the reference section and position means the actuel postion of the marker in the project.

I have check crs.geoccs - but this does actullay nothing in an already georeferenced project. The output of the orthomoasic is still EPSG: 4978

But here is the updated code for OPEX in Version 3.0 tested with Metashape 2.1.3

If anyone have an idee how we can assign the coordinate systeme in Metashape - please let us know.
A workaroud at the moment would be to load the orthomosaic to GIS an apply or overwrite the coordinate system over there.

Code: [Select]
"""
/********************************************************************************************************************
                                    OPEX - Orthophotoexport (V. 3.0)

                                                    Content:
  Processing code to exports orthophoto in planar projection defined by vektor (from three markers) and Z axes.
  This code is mainly for exporting archaeological profiles; ready for importing in GIS.
  Before using the script you need:
  - A georeferenzed Model
  - An "A" and "B" Profile Nail (Named "A" and "B" in Photoscan)
  - A "Y" Point in Photoscan with the coordinates (X-coord from "A",Y-coord from "A", Z-coord from "A"+ a high number)

  Compatibility - Metashape Professional 2.1.3
  Version - 3.0

        -----------------------------------------------------------------------------------------------
        written by          :   Alexey Pasumansky, AgiSoft LLC,
                                Marcel C. Hagner, ArchaeoBW GmbH
                                Florian Tubbesing, ArchaeoBW GmbH
        updated by          :   Philippe Pathé, Bonn Center for Digital Humanities - 2024
        copyright           :   2019 ArchaeoBW GmbH
        email               :   info@archaeobw.de
 ********************************************************************************************************************/

/***************************************************************************************
 *  This program is free software; you can redistribute it and/or modify it under      *
 *  the terms of a Creative Commons Attribution-NonCommercial-ShareAlike 4.0           *
 *  International License.                                                             *
 *                                                                                     *
 *  Attribution - You must give a appropriate credit, provide a link to the license,   *
 *  and indicate of changes were made. You may do so in any reasonable manner, but     *
 *  not in any way that suggests the licensor endorses you or your use                 *
 *                                                                                     *
 *  NonCommercial - You may not use the material for commercial purposes.              *
 *                                                                                     *
 *  ShareAlike - If you remix, transform, or build upon the material, you must         *
 *  distribute your contributions under the same license as the original.              *
 **************************************************************************************/
"""
import os
import Metashape
from PySide2 import QtCore, QtWidgets

def getMarker(chunk, label):
    for marker in chunk.markers:
        if marker.label.upper() == label.upper():
            return marker
    return 0

def vect(a, b):
    result = Metashape.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.normalized()

class ExportOFDlg(QtWidgets.QDialog):

    def __init__(self, parent):
        QtWidgets.QDialog.__init__(self, parent)

        chunk = doc.chunk
        self.supported_types = ["jpg", "tif", "png"]
        self.blending = Metashape.BlendingMode.AverageBlending

        self.setWindowTitle("OPEX")
        # GUI initialization remains unchanged...
        self.btnQuit = QtWidgets.QPushButton("&Exit")
        self.btnQuit.setFixedSize(150,50)
       
        self.btnP1 = QtWidgets.QPushButton("&Export")
        self.btnP1.setFixedSize(150,50)
       
        self.resTxt = QtWidgets.QLabel()
        self.resTxt.setText("Export resolution (m):")
        self.resTxt.setFixedSize(150, 25)
       
        self.vectTxt = QtWidgets.QLabel()
        self.vectTxt.setText("Horizontal vector:")
        self.vectTxt.setFixedSize(150, 25)
       
        self.vect2Txt = QtWidgets.QLabel()
        self.vect2Txt.setText("Vertical vector:")
        self.vect2Txt.setFixedSize(150, 25)
       
        self.resEdt = QtWidgets.QLineEdit()
        self.resEdt.setText("0.00025")
        self.resEdt.setFixedSize(150, 25)
       
        self.radioBtn_v = QtWidgets.QRadioButton("Vertical")
        self.radioBtn_h = QtWidgets.QRadioButton("Horizontal")
        self.radioBtn_v.setChecked(True)
        self.radioBtn_h.setChecked(False)
       
        self.llistV = [None, None]
       
        self.llistV[0] = QtWidgets.QComboBox()
        self.llistV[0].resize(150, 30)
        self.llistV[1] = QtWidgets.QComboBox()
        self.llistV[1].resize(150, 30)
       
        self.llistH = [None, None]
       
        self.llistH[0] = QtWidgets.QComboBox()
        self.llistH[0].resize(150, 30)
        self.llistH[1] = QtWidgets.QComboBox()
        self.llistH[1].resize(150, 30)
       
       
        self.radio = QtWidgets.QVBoxLayout()
        self.radio.addWidget(self.radioBtn_v)
        self.radio.addWidget(self.radioBtn_h)
       
        for marker in chunk.markers:
        for ilist in self.llistV:
        ilist.addItem(marker.label)
        for ilist in self.llistH:
        ilist.addItem(marker.label)
       
       
        layout = QtWidgets.QGridLayout()   #creating layout
        layout.setSpacing(10)
        layout.addWidget(self.resTxt, 0, 0)
        layout.addWidget(self.resEdt, 0, 1)
       
        layout.addWidget(self.vectTxt, 1, 0)
        layout.addWidget(self.vect2Txt, 2, 0)
        layout.addWidget(self.llistH[0], 1, 1)
        layout.addWidget(self.llistH[1], 1, 2)
        layout.addWidget(self.llistV[0], 2, 1)
        layout.addWidget(self.llistV[1], 2, 2)
        layout.addLayout(self.radio, 3, 0)
       
        layout.addWidget(self.btnP1, 3, 1)
        layout.addWidget(self.btnQuit, 3, 2)
        self.setLayout(layout)
       
        proc_export = lambda : self.procExport()
       
        QtCore.QObject.connect(self.btnP1, QtCore.SIGNAL("clicked()"), proc_export)
        QtCore.QObject.connect(self.btnQuit, QtCore.SIGNAL("clicked()"), self, QtCore.SLOT("reject()"))
       
        self.exec()

    def procExport(self):
        chunk = doc.chunk
        print("Starting export process...")

        path = Metashape.app.getExistingDirectory("Bitte wählen Sie den Ordner mit den gesammelten Geo-Cut-Modell aus")
        print(f"Selected path: {path}")

        #crs = Metashape.app.getCoordinateSystem("Wählen Sie das Koordinatensystem für die Ausgabe aus.")
        #print(f"Selected Coordinate System (CRS): {crs}")

        filelist = os.listdir(path)
        doclist = [os.path.join(path, file) for file in filelist if file.endswith(".psx")]
        print(f"Identified Metashape projects: {doclist}")

        try:
            d_x = d_y = float(self.resEdt.text())
            print(f"Set export resolution: d_x = {d_x}, d_y = {d_y}")
        except ValueError:
            Metashape.app.messageBox("Incorrect export resolution! Please use point delimiter.\n")
            print("Script aborted due to invalid resolution input.")
            return 0

        if self.llistV[0].currentIndex() == self.llistV[1].currentIndex():
            Metashape.app.messageBox("Can't use the same marker for vector start and end.\n")
            print("Script aborted due to invalid vertical vector markers.")
            return 0

        if self.llistH[0].currentIndex() == self.llistH[1].currentIndex():
            Metashape.app.messageBox("Can't use the same marker for vector start and end.\n")
            print("Script aborted due to invalid horizontal vector markers.")
            return 0

        if len(chunk.markers) < 2:
            Metashape.app.messageBox("No markers.\n")
            print("Script aborted due to insufficient markers in the chunk.")
            return 0

        print("Script prerequisites verified. Starting processing...")

        type = "_OF.tif"
        for doc_path in doclist:
            print(f"Processing project: {doc_path}")
            doc_unlock = Metashape.Document()
            doc_unlock.open(doc_path, read_only=True)
            doc_unlock.read_only = False
            doc.open(doc_path, read_only=False)
            chunk = doc.chunk
            #chunk.crs = crs
            T = chunk.transform.matrix

            path = doc_path[:-4] + type
            print(f"Export path for orthomosaic: {path}")

            self.btnP1.setDisabled(True)
            self.btnQuit.setDisabled(True)

            m1 = getMarker(chunk, self.llistH[0].currentText())
            m2 = getMarker(chunk, self.llistH[1].currentText())
            horizontal = m2.position - m1.position
            horizontal = T.mulv(horizontal)
            print(f"Horizontal vector calculated: {horizontal}")

            m1 = getMarker(chunk, self.llistV[0].currentText())
            m2 = getMarker(chunk, self.llistV[1].currentText())
            vertical = m2.position - m1.position
            vertical = T.mulv(vertical)
            print(f"Vertical vector calculated: {vertical}")

            normal = vect(vertical, horizontal)
            print(f"Normal vector calculated: {normal}")

            if self.radioBtn_h.isChecked():
                vertical = vect(horizontal, normal)
                horizontal = horizontal.normalized()
                print("Horizontal orientation selected.")
            elif self.radioBtn_v.isChecked():
                vertical = vertical.normalized()
                horizontal = -vect(vertical, normal)
                print("Vertical orientation selected.")
            else:
                Metashape.app.messageBox("Fehler! Kann das OF nicht exportieren" + path)
                print("Script aborted due to invalid orientation selection.")
                continue

            R = Metashape.Matrix([horizontal, vertical, -normal])
            print(f"Rotation matrix R: {R}")

            m = -1
            for marker in chunk.markers:
                m += 1
                label = marker.label
                if "A" in label:
                    a = m

            origin = T.mulp(chunk.markers[a].position)
            print(f"Origin point: {origin}")

            A = m1.reference.location * Metashape.Vector([0, 0, 1])
            X = (-1) * R * origin

            horizontal.size = 4
            horizontal.w = X.x
            vertical.size = 4
            vertical.w = X.y + A
            normal.size = 4
            normal.w = -X.z

            proj = Metashape.Matrix([horizontal, vertical, -normal, Metashape.Vector([0, 0, 0, 1])])
            print(f"Projection matrix: {proj}")

            projection = Metashape.OrthoProjection()
            projection.type = Metashape.OrthoProjection.Type.Planar
            projection.matrix = proj
            if Metashape.app.document.chunk.crs.geoccs:
                projection.crs = Metashape.app.document.chunk.crs.geoccs # case of georeferenced project
            else:
                projection.crs = Metashape.app.document.chunk.crs
            print(f"Projection object created with matrix: {projection.matrix}")
            print(f"CRS for project: {projection.crs}")

            chunk.buildOrthomosaic(
                surface_data=Metashape.DataSource.ModelData,
                blending_mode=self.blending,
                fill_holes=True,
                projection=projection,
                resolution_x=d_x,
                resolution_y=d_y,
                ghosting_filter=False,
                cull_faces=False,
                refine_seamlines=False
            )
            print("Orthomosaic built successfully.")

            chunk.exportRaster(
                path=path,
                source_data=Metashape.DataSource.OrthomosaicData,
                resolution_x=d_x,
                resolution_y=d_y,
                save_alpha=True,
                image_compression=Metashape.ImageCompression(),
                white_background=True
            )
            print(f"Orthomosaic exported to {path}")

            self.btnP1.setDisabled(False)
            self.btnQuit.setDisabled(False)
            Metashape.app.update()

        print("Script completed successfully. Viele Grüße MP & FT")
        return 1

def main():
    global doc
    doc = Metashape.app.document
    app = QtWidgets.QApplication.instance()
    parent = app.activeWindow()
    dlg = ExportOFDlg(parent)

Metashape.app.addMenuItem("ArchaeoBW Toolbox/Profile/OPEX/OPEX ab 2.1.3", main)