...
 
Commits (4)
......@@ -220,9 +220,11 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
self.layout.addLayout(self.menuLayout)
self.layout.addLayout(viewLayout)
if self.datastats.config['minHQI'] is not None:
self.hqiSpinBox.setValue(self.datastats.config['minHQI'])
self.compHqiSpinBox.setValue(self.datastats.config['compHQI'])
minHQI = self.datastats.dataset.resultParams['minHQI']
compHQI = self.datastats.dataset.resultParams['compHQI']
if minHQI is not None:
self.hqiSpinBox.setValue(minHQI)
self.compHqiSpinBox.setValue(compHQI)
self.createActions()
self.createMenus()
......@@ -399,7 +401,8 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
self.resultCheckBoxesLayout.addWidget(self.showTotalSelector)
#generate new checkboxes
self.polymerCheckBoxes = []
for index, polymer in enumerate(self.datastats.uniquePolymers):
uniquePolymers = self.datastats.getUniquePolymers()
for index, polymer in enumerate(uniquePolymers):
self.polymerCheckBoxes.append(QtWidgets.QCheckBox(self))
self.polymerCheckBoxes[index].setText(polymer)
self.resultCheckBoxesLayout.addWidget(self.polymerCheckBoxes[index])
......@@ -427,7 +430,7 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
self.navigationGroup.setEnabled(True)
self.polymerComboBox.currentIndexChanged.disconnect()
self.polymerComboBox.clear()
self.polymerComboBox.addItems(self.datastats.uniquePolymers)
self.polymerComboBox.addItems(uniquePolymers)
self.polymerComboBox.currentIndexChanged.connect(self.displayNewPolymerType)
self.polymerIndex = self.polymerComboBox.currentIndex()
......@@ -453,6 +456,7 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
#draw Sample Spectrum
specIndex = self.currentSpectrumIndex
spectra = self.datastats.spectra
particlestats = self.datastats.getParticleStats()
self.spec_ax.axis("on")
self.spec_ax.clear()
self.spec_ax.plot(spectra[:, 0], spectra[:, specIndex+1])
......@@ -460,7 +464,7 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
self.spec_ax.set_xlabel('Wavenumber (cm-1)', fontsize = 15)
self.spec_ax.set_ylabel('Counts', fontsize = 15)
self.spec_ax.set_title('ScanPoint Number {}, Size = {} µm'.format(specIndex+1,
np.round(self.datastats.particlestats[self.currentParticleIndex][2], 1)))
np.round(particlestats[self.currentParticleIndex][2], 1)))
self.spec_ax.set_xbound(100, (3400 if spectra[-1, 0] > 3400 else spectra[-1, 0]))
wavenumber_diff = list(spectra[:, 0]-100)
y_start = wavenumber_diff.index(min(wavenumber_diff))
......@@ -490,6 +494,48 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
self.parent.highLightRamanIndex(specIndex)
self.lastSpectrumInFocus = specIndex
def selectContour(self, index, centerOn=True):
uniquePolymers = self.datastats.getUniquePolymers()
if uniquePolymers is not None:
#the index is the contour index, find particle index:
specIndex = self.datastats.particles2spectra[index][0] #select first spectrum of partoicle
self.currentParticleIndex = index
self.currentSpectrumIndex = specIndex
selectedPolymer = self.datastats.currentPolymers[specIndex]
try:
self.polymerIndex = uniquePolymers.index(selectedPolymer)
except:
print(selectedPolymer)
raise
#subparticleIndex
partIndicesOfThatPolymer = self.datastats.indices[self.polymerIndex]
subPartInd = partIndicesOfThatPolymer.index(index)
#disconnect analysis widgets:
self.particleSelector.valueChanged.disconnect()
self.spectrumSelector.valueChanged.disconnect()
self.polymerComboBox.currentIndexChanged.disconnect()
#set widgets...
self.particleSelector.setValue(subPartInd+1)
self.particleSelector.setMaximum(len(partIndicesOfThatPolymer))
self.spectrumSelector.setValue(1)
self.spectrumSelector.setMaximum(len(self.datastats.particles2spectra[index]))
selectedPolymer = self.datastats.currentPolymers[specIndex]
self.polymerIndex = uniquePolymers.index(selectedPolymer)
self.polymerComboBox.setCurrentIndex(self.polymerIndex)
#reconnect all widgets:
self.particleSelector.valueChanged.connect(self.selectParticle)
self.spectrumSelector.valueChanged.connect(self.selectSpectrum)
self.polymerComboBox.currentIndexChanged.connect(self.displayNewPolymerType)
self.updateSpecPlot(centerOn=centerOn)
def displayNewPolymerType(self, resetCurrentIndex=True):
self.polymerIndex = self.polymerComboBox.currentIndex()
self.particleSelector.setMaximum(len(self.datastats.indices[self.polymerIndex]))
......@@ -604,7 +650,7 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
#general size histogram
self.bins = np.logspace(0.1, 3, 20)
self.sizes = [i[0] if np.isnan(i[2]) else i[2] for i in self.datastats.particlestats] #extract long size (if ellipse fit is nan -> box fit)
self.sizes = [i[0] if np.isnan(i[2]) else i[2] for i in self.datastats.getParticleStats()] #extract long size (if ellipse fit is nan -> box fit)
sizehist = np.histogram(self.sizes, self.bins)
self.totalhistx = []
for i in range(19):
......@@ -676,6 +722,7 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
return color
def createPolymerOverlay(self):
uniquePolymers = self.datastats.getUniquePolymers()
if not self.noOverlayAct.isChecked() and self.datastats.indices is not None:
if len(self.datastats.indices) > 0:
......@@ -686,9 +733,9 @@ class ParticleAnalysis(QtWidgets.QMainWindow):
for index, indexList in enumerate(self.datastats.indices):
if self.fullOverlayAct.isChecked() or (self.selOverlayAct.isChecked() and self.polymerCheckBoxes[index].isChecked()):
color = self.getColorFromName(self.datastats.uniquePolymers[index], base255=True)
color = self.getColorFromName(uniquePolymers[index], base255=True)
color = QtGui.QColor(color[0], color[1], color[2], alpha=alpha)
legendItems.append((self.datastats.uniquePolymers[index], color))
legendItems.append((uniquePolymers[index], color))
for i in indexList:
colorList[i] = color
......
......@@ -36,7 +36,7 @@ class ExpExcelDialog(QtWidgets.QDialog):
self.setGeometry(200,200, 300, 300)
self.datastats = datastats
self.particles = self.datastats.particlestats
self.particles = self.datastats.getParticleStats()
self.polymers = self.datastats.particleResults
self.additives = self.datastats.currentAdditives
self.hqis = self.datastats.hqis
......
......@@ -21,20 +21,30 @@ If not, see <https://www.gnu.org/licenses/>.
import os
import numpy as np
import operator
from dataset import loadData, recursiveDictCompare
def readDataStats(fname):
ds = loadData(fname)
datastats = DataStats(ds)
datastats.update()
datastats.loadParticleData()
minHQI = datastats.dataset.resultParams['minHQI']
compHQI = datastats.dataset.resultParams['compHQI']
datastats.formatResults(minHQI, compHQI)
datastats.createHistogramData()
return datastats
class DataStats(object):
def __init__(self, dataset):
self.dataset = dataset
self.config = dataset.resultParams
self.spectraResults = None #entire List of all spectra assignments
self.additiveResults = None #entire List of all additives
self.particlestats = None
self.particleResults = None #final assignment for each particle
self.currentPolymers = None #list of polymers after setting entries with low hqi to unknown
self.currentAdditives = None #same thing for the additives
self.uniquePolymers = None #list of present polymer types
self.spectra = None #acquired spectra
self.indices = None #assignment of what spectra-indices belong to what substance
......@@ -102,15 +112,6 @@ class DataStats(object):
return specs
def loadParticleData(self):
self.particlestats = np.array(self.dataset.particlestats)
pixelscale = (self.dataset.pixelscale_df if self.dataset.imagescanMode == 'df' else self.dataset.pixelscale_bf)
#convert to mikrometer scale
for index in range(len(self.particlestats)):
for subindex in range(5):
self.particlestats[index][subindex] = self.particlestats[index][subindex] * pixelscale #multiply by pixelscale
if subindex == 4:
self.particlestats[index][subindex] = self.particlestats[index][subindex] * pixelscale #again for the area...
self.particles2spectra = self.dataset.particles2spectra
sortindices = self.dataset.ramanscansortindex
......@@ -145,12 +146,26 @@ class DataStats(object):
if self.currentAdditives is not None:
self.currentAdditives[self.addhqis < compHqi] = 'unknown'
def getUniquePolymers(self):
if self.currentPolymers is None:
return None
return self.uniquePolymers
def getParticleStats(self):
particlestats = np.array(self.dataset.particlestats)
pixelscale = self.dataset.getPixelScale()
#convert to mikrometer scale
particlestats[:,:5] *= pixelscale
particlestats[:,4] *= pixelscale #again for the area...
return particlestats
def createHistogramData(self):
particlestats = self.getParticleStats()
self.uniquePolymers = np.unique(self.currentPolymers)
self.particleResults = [None]*len(self.particlestats)
self.particleResults = [None]*len(particlestats)
self.typehistogram = {i: 0 for i in self.uniquePolymers}
if len(self.particles2spectra) != len(self.particlestats):
if len(self.particles2spectra) != len(particlestats):
return False
for particleID, specList in enumerate(self.particles2spectra):
......@@ -192,4 +207,12 @@ class DataStats(object):
self.dataset.resultParams = {'minHQI': minHQI,
'compHQI': compHQI}
self.dataset.save()
print('saved dataset')
\ No newline at end of file
testresult = self.testRead()
print('saved dataset; Valid:', testresult)
return testresult
def testRead(self):
statsread = readDataStats(self.dataset.fname)
return recursiveDictCompare(self.__dict__, statsread.__dict__)
\ No newline at end of file
......@@ -261,7 +261,12 @@ class LoadWITecResults(QtWidgets.QDialog):
self.parent.updateBtn.clicked.connect(self.parent.formatResults)
self.parent.formatResults()
self.parent.show_hide_labels()
self.parent.saveAnalysisResults()
minHQI = self.parent.hqiSpinBox.value()
compHQI = self.parent.compHqiSpinBox.value()
if not self.parent.datastats.saveAnalysisResults(minHQI, compHQI):
QtWidgets.QMessageBox.warning(self.parent, 'Error!',
'Data inconsistency after saving!',
QtWidgets.QMessageBox.Ok, QtWidgets.QMessageBox.Ok)
self.parent.setEnabled(True)
event.accept()
......
......@@ -24,6 +24,7 @@ You should have received a copy of the GNU General Public License
along with this program, see COPYING.
If not, see <https://www.gnu.org/licenses/>.
"""
import numpy as np
import cv2
from PyQt5 import QtWidgets
......@@ -56,8 +57,6 @@ class ParticleEditor(object):
return
contourIndices = sorted(contourIndices) #we want to keep the contour with lowest index
print('merging contours:', contourIndices)
self.createSafetyBackup()
#get contours:
contours = [self.datastats.dataset.particlecontours[i] for i in contourIndices]
cnt = np.vstack(tuple(contours)) #combine contous
......@@ -72,7 +71,7 @@ class ParticleEditor(object):
img = np.zeros((rangey, rangex))
for i in contourIndices:
curCnt = self.datastats.dataset.particlecontours[i]
curCnt = self.datastats.dataset.particlecontours[i].copy()
for i in range(len(curCnt)):
curCnt[i][0][0] -= xmin-padding
curCnt[i][0][1] -= ymin-padding
......@@ -87,6 +86,11 @@ class ParticleEditor(object):
else:
temp, contours, hierarchy = cv2.findContours(img, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE)
if len(contours)>1:
QtWidgets.QMessageBox.critical(self.parent, 'ERROR!',
'Particle contours are not connected and cannot be combined!')
return
newContour = contours[0]
stats = self.characterizeParticle(newContour)
......@@ -94,6 +98,8 @@ class ParticleEditor(object):
newContour[i][0][0] += xmin-padding
newContour[i][0][1] += ymin-padding
print('merging contours:', contourIndices)
self.createSafetyBackup()
#check, if dataset contains (already modified) particle2spectra, otherwise create new.
if self.datastats.dataset.particles2spectra is None: #create default assignment
......@@ -133,12 +139,17 @@ class ParticleEditor(object):
print('removing index from particles2spectra:', index)
del self.datastats.dataset.particles2spectra[index]
#save data
self.datastats.saveAnalysisResults()
#update contours in sampleview
self.parent.parent.contouritem.resetContours(self.datastats.dataset.particlecontours)
self.parent.loadParticleData()
#save data
minHQI = self.parent.hqiSpinBox.value()
compHQI = self.parent.compHqiSpinBox.value()
if not self.datastats.saveAnalysisResults(minHQI, compHQI):
QtWidgets.QMessageBox.warning(self.parent, 'Error!',
'Data inconsistency after saving!', QtWidgets.QMessageBox.Ok,
QtWidgets.QMessageBox.Ok)
def reassignParticles(self, contourindices, new_assignment):
......@@ -155,10 +166,16 @@ class ParticleEditor(object):
self.datastats.spectraResults[specIndex] = new_assignment
self.datastats.hqis[specIndex] = 100
#save data
self.datastats.saveAnalysisResults()
#update contours in sampleview
self.parent.parent.contouritem.resetContours(self.datastats.dataset.particlecontours)
self.parent.loadParticleData()
#save data
minHQI = self.parent.hqiSpinBox.value()
compHQI = self.parent.compHqiSpinBox.value()
if not self.datastats.saveAnalysisResults(minHQI, compHQI):
QtWidgets.QMessageBox.warning(self.parent, 'Error!',
'Data inconsistency after saving!',
QtWidgets.QMessageBox.Ok, QtWidgets.QMessageBox.Ok)
def deleteParticles(self):
......
......@@ -21,16 +21,16 @@ class SQLExport(QtWidgets.QDialog):
self.datastats = datastats
self.polymerList = self.datastats.particleResults
self.longSizes = np.round(np.array([i[0] if np.isnan(i[2]) else i[2] for i in self.datastats.particlestats]), 1)
self.shortSize = np.round(np.array([i[1] if np.isnan(i[3]) else i[3] for i in self.datastats.particlestats]), 1)
particlestats = self.datastats.getParticleStats()
self.longSizes = np.round(np.array([i[0] if np.isnan(i[2]) else i[2] for i in particlestats]), 1)
self.shortSize = np.round(np.array([i[1] if np.isnan(i[3]) else i[3] for i in particlestats]), 1)
#spectra can be quite some data size, they are not copied here but referenced later on...
self.particleImages = None
self.log = []
configfilename = os.path.join(os.path.dirname(os.path.split(__file__)[0]),
'database_config.txt')
configfilename = os.path.join(os.path.split(__file__)[0], 'database_config.txt')
if not os.path.exists(configfilename):
QtWidgets.QMessageBox.warning(self, 'Warning!',
......@@ -338,6 +338,9 @@ class SQLExport(QtWidgets.QDialog):
self.cnx = mysql.connector.connect(**self.config) #port: 3306
except mysql.connector.Error as err:
print(err)
QtWidgets.QMessageBox.warning(self, 'Error!',
str(err),
QtWidgets.QMessageBox.Ok, QtWidgets.QMessageBox.Ok)
self.cnx = None
def getEntireTable(self, tablename):
......
......@@ -53,6 +53,58 @@ def saveData(dataset, fname):
pickle.dump(dataset, fp, protocol=-1)
dataset.zvalimg = zvalimg
def arrayCompare(a1, a2):
if a1.shape!=a2.shape:
return False
if a1.dtype!=np.float32 and a1.dtype!=np.float64:
return np.all(a1==a2)
ind = np.isnan(a1)
if not np.any(ind):
return np.all(a1==a2)
return np.all(a1[~ind]==a2[~ind])
def listCompare(l1, l2):
if len(l1)!=len(l2):
return False
for l1i, l2i in zip(l1, l2):
if isinstance(l1i, np.ndarray):
if not isinstance(l2i, np.ndarray) or not arrayCompare(l1i, l2i):
return False
elif isinstance(l1i, (list, tuple)):
if not isinstance(l2i, (list, tuple)) or not listCompare(l1i, l2i):
return False
elif l1i!=l2i and ((~np.isnan(l1i)) or (~np.isnan(l2i))):
return False
return True
def recursiveDictCompare(d1, d2):
for key in d1:
if not key in d2:
print("key missing in d2:", key, flush=True)
return False
a = d1[key]
b = d2[key]
print(key, type(a), type(b), flush=True)
if isinstance(a, np.ndarray):
if not isinstance(b, np.ndarray) or not arrayCompare(a, b):
print("data is different!", a, b)
return False
elif isinstance(a, dict):
if not isinstance(b, dict):
print("data is different!", a, b)
return False
if not recursiveDictCompare(a, b):
return False
elif isinstance(a, (list, tuple)):
if not isinstance(b, (list, tuple)) or not listCompare(a, b):
print("data is different!", a, b)
return False
elif a != b:
if (a is not None) and (b is not None):
print("data is different!", a, b)
return False
return True
class DataSet(object):
def __init__(self, fname, newProject=False):
self.fname = fname
......@@ -115,6 +167,14 @@ class DataSet(object):
self.fname = self.newProject(fname)
self.updatePath()
def __eq__(self, other):
return recursiveDictCompare(self.__dict__, other.__dict__)
def getPixelScale(self, mode=None):
if mode is None:
mode = self.imagescanMode
return (self.pixelscale_df if mode == 'df' else self.pixelscale_bf)
def saveZvalImg(self):
if self.zvalimg is not None:
cv2imwrite_fix(self.getZvalImageName(), self.zvalimg)
......@@ -243,8 +303,7 @@ class DataSet(object):
p0[1] += self.imagedim_bf[1]/2
return (pixelpos[0]*self.pixelscale_bf + p0[0]), (p0[1] - pixelpos[1]*self.pixelscale_bf)
else:
print('mapToRamanMode not understood')
return
raise ValueError(f'mapToLength mode: {mode} not understood')
def mapToLengthRaman(self, pixelpos, microscopeMode='df', noz=False):
p0x, p0y = self.mapToLength(pixelpos, mode = microscopeMode)
......@@ -315,9 +374,7 @@ class DataSet(object):
saveData(self, self.fname)
def saveBackup(self):
# backupNameNotFound = True
inc = 0
# while backupNameNotFound:
while True:
directory = os.path.dirname(self.fname)
filename = self.name + '_backup_' + str(inc) + '.pkl'
......@@ -327,5 +384,4 @@ class DataSet(object):
else:
saveData(self, path)
return filename
# backupNameNotFound = False
......@@ -355,7 +355,7 @@ class SampleView(QtWidgets.QGraphicsView):
self.announceScaling()
def announceScaling(self):
pixelscale = (self.dataset.pixelscale_df if self.microscopeMode == 'df' else self.dataset.pixelscale_bf)
pixelscale = self.dataset.getPixelScale(self.microscopeMode)
if self.dataset is None or pixelscale is None:
self.ScalingChanged.emit(-1.0)
else:
......@@ -453,10 +453,8 @@ class SampleView(QtWidgets.QGraphicsView):
if len(p2)<len(p1):
p1 = [[pi[1],pi[0]] for pi in p2]
self.dataset.grid = p1
if self.microscopeMode == 'df':
wxs, wys = width/self.dataset.pixelscale_df, height/self.dataset.pixelscale_df
else:
wxs, wys = width/self.dataset.pixelscale_bf, height/self.dataset.pixelscale_bf
pixelscale = self.dataset.getPixelScale(self.microscopeMode)
wxs, wys = width/pixelscale, height/pixelscale
self.scanitems = []
for i, p in enumerate(p1):
......@@ -489,46 +487,7 @@ class SampleView(QtWidgets.QGraphicsView):
@QtCore.pyqtSlot(int, bool)
def selectContour(self, index, centerOn=True):
if self.analysiswidget is not None:
if self.analysiswidget.uniquePolymers is not None:
#the index is the contour index, find particle index:
specIndex = self.analysiswidget.particles2spectra[index][0] #select first spectrum of partoicle
self.analysiswidget.currentParticleIndex = index
self.analysiswidget.currentSpectrumIndex = specIndex
selectedPolymer = self.analysiswidget.currentPolymers[specIndex]
try:
self.analysiswidget.polymerIndex = self.analysiswidget.uniquePolymers.index(selectedPolymer)
except:
print(selectedPolymer)
raise
#subparticleIndex
partIndicesOfThatPolymer = self.analysiswidget.indices[self.analysiswidget.polymerIndex]
subPartInd = partIndicesOfThatPolymer.index(index)
#disconnect analysis widgets:
self.analysiswidget.particleSelector.valueChanged.disconnect()
self.analysiswidget.spectrumSelector.valueChanged.disconnect()
self.analysiswidget.polymerComboBox.currentIndexChanged.disconnect()
#set widgets...
self.analysiswidget.particleSelector.setValue(subPartInd+1)
self.analysiswidget.particleSelector.setMaximum(len(partIndicesOfThatPolymer))
self.analysiswidget.spectrumSelector.setValue(1)
self.analysiswidget.spectrumSelector.setMaximum(len(self.analysiswidget.particles2spectra[index]))
selectedPolymer = self.analysiswidget.currentPolymers[specIndex]
self.analysiswidget.polymerIndex = self.analysiswidget.uniquePolymers.index(selectedPolymer)
self.analysiswidget.polymerComboBox.setCurrentIndex(self.analysiswidget.polymerIndex)
#reconnect all widgets:
self.analysiswidget.particleSelector.valueChanged.connect(self.analysiswidget.selectParticle)
self.analysiswidget.spectrumSelector.valueChanged.connect(self.analysiswidget.selectSpectrum)
self.analysiswidget.polymerComboBox.currentIndexChanged.connect(self.analysiswidget.displayNewPolymerType)
self.analysiswidget.updateSpecPlot(centerOn=centerOn)
self.analysiswidget.selectContour(index, centerOn)
def prepareAnalysis(self):
self.clearItems()
......
......@@ -22,7 +22,7 @@ import numpy as np
from PyQt5 import QtCore, QtWidgets, QtGui
class SegmentationContours(QtWidgets.QGraphicsItem):
def __init__(self, parent=None, contours=[], pos=(0,0)):
def __init__(self, parent, contours=[], pos=(0,0)):
super().__init__()
self.parent = parent
self.setPos(pos[0], pos[1])
......@@ -122,14 +122,14 @@ class SegmentationContours(QtWidgets.QGraphicsItem):
for index in self.selectedContours:
# partIndex = int(np.where(self.parent.dataset.ramanscansortindex == index)[0])
partIndex = index
assignments.append(self.parent.analysiswidget.particleResults[partIndex])
assignments.append(self.parent.analysiswidget.datastats.particleResults[partIndex])
assignments.append("other")
for assignment in np.unique(np.array(assignments)):
combineActs.append(combineMenu.addAction(assignment))
reassignActs = []
reassignMenu = QtWidgets.QMenu("Reassign particle(s) into")
for polymer in self.parent.analysiswidget.uniquePolymers:
for polymer in self.parent.analysiswidget.datastats.getUniquePolymers():
reassignActs.append(reassignMenu.addAction(polymer))
reassignActs.append(reassignMenu.addAction("other"))
......@@ -145,25 +145,16 @@ class SegmentationContours(QtWidgets.QGraphicsItem):
elif numParticles == 1:
combineMenu.setDisabled(True)
action = contextMenu.exec_(event.screenPos())
if action == deleteAct:
print('deleting')
elif action in combineActs:
newAssignment = action.text()
# if newAssignment == "other":
# QtWidgets.QMessageBox.about(self.parent, "Not yet implemented", "we are getting there...")
# return
self.parent.analysiswidget.editor.combineParticles(self.selectedContours, newAssignment)
elif action in reassignActs:
newAssignment = action.text()
# if newAssignment == "other":
# QtWidgets.QMessageBox.about(self.parent, "Not yet implemented", "we are getting there...")
# return
self.parent.analysiswidget.editor.reassignParticles(self.selectedContours, newAssignment)
......
......@@ -26,6 +26,7 @@ from helperfunctions import cv2imread_fix, cv2imwrite_fix
from ramancom.ramancontrol import defaultPath
from dataset import DataSet
from scipy.optimize import least_squares
from itertools import permutations
import cv2
import numpy as np
......@@ -228,20 +229,34 @@ class ZeissImporter(QtWidgets.QDialog):
[s1*c3+c1*c2*s3, -s1*s3+c1*c2*c3, -c1*s2],
[s1*s3, s2*c3, c2]])
# find the transformation matrix with best fit for small angles in
# [-45°,45°] for all permutation of markers
permbest = None
pointsbest = None
for perm in permutations(range(points.shape[0])):
ppoints = points[perm,:]
def err(angles_shift):
T = getRotMat(angles_shift[:3]).T.A
return (np.dot(zpoints, T)-angles_shift[np.newaxis,3:]-points).ravel()
best = None
for i in range(100):
angle = np.pi*np.random.rand(3)-np.pi/2.
opt = least_squares(err, np.concatenate((angle, np.zeros(3))), method='lm')
if best is None or best.cost>opt.cost:
best = opt
optangles = best.x[:3]
shift = best.x[3:]
return (np.dot(zpoints, T) - angles_shift[np.newaxis,3:] \
- ppoints).ravel()
angle = np.zeros(3)
opt = least_squares(err, np.concatenate((angle, np.zeros(3))),
bounds=(np.array([-np.pi/4]*3+[-np.inf]*3),
np.array([np.pi/4]*3+[np.inf]*3)),
method='dogbox')
if permbest is None or \
permbest.cost>opt.cost:
print("Current best permutation:", perm, flush=True)
permbest = opt
pointsbest = ppoints
optangles = permbest.x[:3]
shift = permbest.x[3:]
T = getRotMat(optangles).T.A
e = (np.dot(zpoints, T)-shift[np.newaxis,:]-points)
e = (np.dot(zpoints, T)-shift[np.newaxis,:]-pointsbest)
print("Transformation angles:", optangles, flush=True)
print("Transformation shift:", shift, flush=True)
print("Transformation err:", e, flush=True)
......