From 5fd1cb0cdbd02186411f34c5a94ab4212ac25324 Mon Sep 17 00:00:00 2001 From: Lars Bittrich Date: Tue, 9 Apr 2019 10:03:58 +0200 Subject: [PATCH] Zeissimporter initial commit with coordinate transfer --- gepard.py | 23 ++++- opticalscan.py | 8 +- sampleview.py | 100 ++++++++++--------- zeissimporter.py | 255 +++++++++++++++++++++++++++++++++++++++++++++++ zeissxml.py | 100 +++++++++++++++++++ 5 files changed, 437 insertions(+), 49 deletions(-) create mode 100644 zeissimporter.py create mode 100644 zeissxml.py diff --git a/gepard.py b/gepard.py index d6c6b03..76d3ac8 100755 --- a/gepard.py +++ b/gepard.py @@ -84,6 +84,16 @@ class MeasureParticleWindow(QtWidgets.QMainWindow): self.fname = str(fileName) self.view.open(self.fname) self.scalingChanged(1.) + + @QtCore.pyqtSlot() + def importProject(self, fileName=False): + if fileName is False: + fileName = QtWidgets.QFileDialog.getOpenFileName(self, "Import Zeiss Zen Project", + defaultPath, "*.xml")[0] + if fileName: + self.fname = str(fileName) + self.view.importProject(self.fname) + self.scalingChanged(1.) @QtCore.pyqtSlot() def new(self, fileName=False): @@ -115,6 +125,10 @@ class MeasureParticleWindow(QtWidgets.QMainWindow): self.openAct.setShortcut("Ctrl+O") self.openAct.triggered.connect(self.open) + self.importAct = QtWidgets.QAction("&Import Project...", self) + self.importAct.setShortcut("Ctrl+I") + self.importAct.triggered.connect(self.importProject) + self.newAct = QtWidgets.QAction("&New Measurement...", self) self.newAct.setShortcut("Ctrl+N") self.newAct.triggered.connect(self.new) @@ -211,12 +225,14 @@ class MeasureParticleWindow(QtWidgets.QMainWindow): def unblockUI(self, connected): self.openAct.setEnabled(True) + self.importAct.setEnabled(True) self.newAct.setEnabled(True) self.updateConnected(connected) self.exitAct.setEnabled(True) def blockUI(self): self.openAct.setEnabled(False) + self.importAct.setEnabled(False) self.newAct.setEnabled(False) self.connectRamanAct.setEnabled(False) self.disconnectRamanAct.setEnabled(False) @@ -237,6 +253,7 @@ class MeasureParticleWindow(QtWidgets.QMainWindow): def createMenus(self): self.fileMenu = QtWidgets.QMenu("&File", self) self.fileMenu.addAction(self.newAct) + self.fileMenu.addAction(self.importAct) self.fileMenu.addAction(self.openAct) self.fileMenu.addSeparator() self.fileMenu.addAction(self.exitAct) @@ -272,6 +289,7 @@ class MeasureParticleWindow(QtWidgets.QMainWindow): self.toolbar.setIconSize(QtCore.QSize(100,50)) self.toolbar.addAction(self.aboutAct) self.toolbar.addAction(self.newAct) + self.toolbar.addAction(self.importAct) self.toolbar.addAction(self.openAct) self.toolbar.addSeparator() self.toolbar.addAction(self.connectRamanAct) @@ -295,6 +313,7 @@ if __name__ == '__main__': logpath = QtCore.QStandardPaths.writableLocation( QtCore.QStandardPaths.AppLocalDataLocation) + fp = None if logpath != "": if not os.path.exists(logpath): os.mkdir(logpath) @@ -309,4 +328,6 @@ if __name__ == '__main__': meas = MeasureParticleWindow(logpath) meas.showMaximized() - ret = app.exec_() \ No newline at end of file + ret = app.exec_() + if fp is not None: + fp.close() \ No newline at end of file diff --git a/opticalscan.py b/opticalscan.py index a1ed149..32f651a 100755 --- a/opticalscan.py +++ b/opticalscan.py @@ -113,11 +113,12 @@ def loadAndPasteImage(srcnames, fullimage, fullzval, width, height, class PointCoordinates(QtWidgets.QGridLayout): readPoint = QtCore.pyqtSignal(float, float, float, name='readPoint') - def __init__(self, N, ramanctrl, parent=None): + def __init__(self, N, ramanctrl, parent=None, names=None): super().__init__(parent) self.dswidgets = [] self.N = 0 self.ramanctrl = ramanctrl + self.names = names self.pimageOnly = QtWidgets.QPushButton("Image") self.addWidget(self.pimageOnly, 0, 6, QtCore.Qt.AlignRight) self.pimageOnly.released.connect(QtCore.pyqtSlot()(lambda : self.read(-1))) @@ -140,7 +141,10 @@ class PointCoordinates(QtWidgets.QGridLayout): self.itemAtPosition(i+1,5).setVisible(True) self.itemAtPosition(i+1,6).setVisible(True) for i in range(self.N, N): - lx = QtWidgets.QLabel(f"{i+1} -> x:") + if self.names is not None: + lx = QtWidgets.QLabel(f"{self.names[i]} -> x:") + else: + lx = QtWidgets.QLabel(f"{i+1} -> x:") ly = QtWidgets.QLabel("y:") lz = QtWidgets.QLabel("z:") wx = QtWidgets.QDoubleSpinBox() diff --git a/sampleview.py b/sampleview.py index 98996e5..c82a073 100644 --- a/sampleview.py +++ b/sampleview.py @@ -27,6 +27,7 @@ from opticalscan import OpticalScan from ramanscanui import RamanScanUI from detectionview import ParticleDetectionView from analysis.analysisview import ParticleAnalysis +from zeissimporter import ZeissImporter from viewitems import FitPosIndicator, Node, Edge, ScanIndicator, RamanScanIndicator, SegmentationContours from helperfunctions import polygoncovering, cv2imread_fix import cv2 @@ -84,7 +85,7 @@ class SampleView(QtWidgets.QGraphicsView): self.oscanwidget.imageUpdate.connect(self.loadPixmap) self.oscanwidget.boundaryUpdate.connect(self.resetBoundary) self.analysiswidget = None - self.setMinimumSize(600,600) + self.setMinimumSize(600, 600) self.darkenPixmap = False self.microscopeMode = None self.coordTestMode = False @@ -237,6 +238,13 @@ class SampleView(QtWidgets.QGraphicsView): self.activateMaxMode(loadnew=True) self.imparent.snapshotAct.setEnabled(True) + def importProject(self, fname): + zimp = ZeissImporter(fname, self.ramanctrl, self) + if zimp.validimport: + zimp.exec() + if zimp.result() == QtWidgets.QDialog.Accepted: + self.open(zimp.gepardname) + def new(self, fname): self.saveDataSet() if self.dataset is not None: @@ -347,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.pixelscale_df if self.microscopeMode == 'df' else self.dataset.pixelscale_bf) if self.dataset is None or pixelscale is None: self.ScalingChanged.emit(-1.0) else: @@ -445,10 +453,10 @@ class SampleView(QtWidgets.QGraphicsView): if len(p2) + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program, see COPYING. +If not, see . +""" +import os +from PyQt5 import QtCore, QtWidgets +from zeissxml import ZeissHandler, make_parser +from opticalscan import PointCoordinates +from helperfunctions import cv2imread_fix, cv2imwrite_fix +from ramancom.ramancontrol import defaultPath +from dataset import DataSet +from scipy.optimize import least_squares +import cv2 +import numpy as np + +class ZeissImporter(QtWidgets.QDialog): + def __init__(self, fname, ramanctrl, parent=None): + super().__init__(parent) + + if not ramanctrl.connect() or not self.readImportData(fname): + msg = QtWidgets.QMessageBox() + msg.setText('Connection failed! Please enable remote control.') + msg.exec() + self.validimport = False + return + else: + self.validimport = True + + self.ramanctrl = ramanctrl + + vbox = QtWidgets.QVBoxLayout() + pointgroup = QtWidgets.QGroupBox('Marker coordinates at Raman spot [µm]', self) + self.points = PointCoordinates(len(self.markers), self.ramanctrl, self, + names = [m.name for m in self.markers]) + self.points.pimageOnly.setVisible(False) + pointgroup.setLayout(self.points) + self.points.readPoint.connect(self.takePoint) + + self.pconvert = QtWidgets.QPushButton('Convert', self) + self.pexit = QtWidgets.QPushButton('Cancel', self) + self.pconvert.released.connect(self.convert) + self.pexit.released.connect(self.reject) + self.pconvert.setEnabled(False) + + btnLayout = QtWidgets.QHBoxLayout() + btnLayout.addStretch() + btnLayout.addWidget(self.pconvert) + btnLayout.addWidget(self.pexit) + + label = QtWidgets.QLabel("Z-Image blur radius", self) + self.blurspinbox = QtWidgets.QSpinBox(self) + self.blurspinbox.setMinimum(3) + self.blurspinbox.setMaximum(99) + self.blurspinbox.setSingleStep(2) + self.blurspinbox.setValue(5) + blurlayout = QtWidgets.QHBoxLayout() + blurlayout.addWidget(label) + blurlayout.addWidget(self.blurspinbox) + blurlayout.addStretch() + + vbox.addWidget(pointgroup) + vbox.addLayout(blurlayout) + vbox.addLayout(btnLayout) + self.setLayout(vbox) + + def readImportData(self, fname): + path = os.path.split(fname)[0] + self.zmapimgname = os.path.join(path, '3D.tif') + self.edfimgname = os.path.join(path, 'EDF.tif') + xmlname = os.path.join(path, '3D.tif_metadata.xml') + + errmsges = [] + if not os.path.exists(self.zmapimgname): + errmsges.append('Depth map image not found: 3D.tif') + if not os.path.exists(self.edfimgname): + errmsges.append('EDF image not found: EDF.tif') + if not os.path.exists(xmlname): + errmsges.append('XML metadata not found: 3D.tif_metadata.xml') + else: + parser = make_parser() + z = ZeissHandler() + parser.setContentHandler(z) + parser.parse(xmlname) + + if len(z.markers)<3: + errmsges.append('Fewer than 3 markers found to adjust coordinates!') + if None in [z.region.centerx, z.region.centery, + z.region.width, z.region.height]: + errmsges.append('Image dimensions incomplete or missing!') + if None in [z.zrange.z0, z.zrange.zn, z.zrange.dz]: + errmsges.append('ZStack information missing or incomplete!') + + if len(errmsges)>0: + QtWidgets.QMessageBox.error(self, 'Error!', + '\n'.join(errmsges), + QtWidgets.QMessageBox.Ok, + QtWidgets.QMessageBox.Ok) + return False + + self.region = z.region + self.zrange = z.zrange + self.markers = z.markers + return True + + @QtCore.pyqtSlot(float, float, float) + def takePoint(self, x, y, z): + points = self.points.getPoints() + if len(points)>=3: + self.pconvert.setEnabled(True) + + @QtCore.pyqtSlot() + def convert(self): + fname = QtWidgets.QFileDialog.getSaveFileName(self, + 'Create New GEPARD Project', defaultPath, '*.pkl')[0] + if fname=='': + return + dataset = DataSet(fname, newProject=True) + T, pc, zpc = self.getTransform() + imgshape, warp_mat = self.convertZimg(dataset, T, pc, zpc) + self.convertImage(dataset, warp_mat) + dataset.save() + self.gepardname = dataset.fname + self.accept() + + def convertImage(self, dataset, warp_mat): + img = cv2imread_fix(self.edfimgname) + img = cv2.warpAffine(img, warp_mat, img.shape[:2][::-1]) + cv2imwrite_fix(dataset.getImageName(), img) + + + def convertZimg(self, dataset, T, pc, zpc): + N = int(round((self.zrange.zn-self.zrange.z0)/self.zrange.dz)) + dataset.zpositions = np.linspace(self.zrange.z0, + self.zrange.zn, N)-zpc[2]+pc[2] + zimg = cv2imread_fix(self.zmapimgname, cv2.IMREAD_GRAYSCALE) + zmdist = zimg.mean() + zm = zmdist/255.*(self.zrange.zn-self.zrange.z0) + self.zrange.z0 + + radius = self.blurspinbox.value() + blur = cv2.GaussianBlur(zimg, (radius, radius), 0) + + pshift = self.ramanctrl.getRamanPositionShift() + dataset.pshift = pshift + pixelscale = self.region.width/zimg.shape[1] + + # use input image as single image aquired in one shot + dataset.imagedim_df = (self.region.width, self.region.height, 0.0) + dataset.pixelscale_df = pixelscale + + dataset.imagedim_bf = (self.region.width, self.region.height, 0.0) + dataset.pixelscale_bf = pixelscale + + # set image center as reference point in data set (transform from Zeiss) + p0 = np.dot((np.array([self.region.centerx, + self.region.centery,zm])-zpc),T)[:2] + pc[:2] + dataset.readin = False + dataset.lastpos = p0 + dataset.maxdim = p0 + p0 + + # pixel triangle for coordinate warping transformation + srcTri = np.array( [[0, 0], [zimg.shape[1] - 1, 0], + [0, zimg.shape[0] - 1]] ).astype(np.float32) + # upper left point (0,0) in Zeiss coordinates: + z0 = np.array([self.region.centerx - self.region.width/2, + self.region.centery + self.region.height/2]) + # transform pixel data to Zeiss coordinates + dstTri = np.array([[p[0]*pixelscale + z0[0], + z0[1] - p[1]*pixelscale, zm] for p in srcTri]).astype(np.double)-zpc + + # transform to Raman coordinates + dstTri = np.dot(dstTri,T) + pc[np.newaxis,:] + + # tilt blur image based on transformend z and adapt zpositions + x = np.linspace(0,1,blur.shape[1]) + y = np.linspace(0,1,blur.shape[0]) + x, y = np.meshgrid(x,y) + zmap = x*(dstTri[1,2]-dstTri[0,2]) + y*(dstTri[2,2]-dstTri[0,2]) + \ + (zimg * ((self.zrange.zn-self.zrange.z0)/255.) - \ + zmdist*((self.zrange.zn-self.zrange.z0)/255.)) + zmin, zmax = zmap.min(), zmap.max() + dataset.zpositions = np.array([zmap.min(), zmap.max()]) + blur = (zmap-zmin)*(255./(zmax-zmin)) + blur[blur>255.] = 255. + blur = np.uint8(blur) + # transform triangle back to pixel + dstTri = np.array([dataset.mapToPixel(p[:2]) for p in dstTri]).astype(np.float32) + + warp_mat = cv2.getAffineTransform(srcTri, dstTri) + + blur = cv2.warpAffine(blur, warp_mat, zimg.shape[::-1]) + zimgname = dataset.getZvalImageName() + cv2imwrite_fix(zimgname, blur) + return zimg.shape, warp_mat + + def getTransform(self): + points = self.points.getPoints() + pshift = self.ramanctrl.getRamanPositionShift() + points[:,0] -= pshift[0] + points[:,1] -= pshift[1] + zpoints = np.array([m.getPos() for m in self.markers], dtype=np.double) + pc = points.mean(axis=0) + zpc = zpoints.mean(axis=0) + + points -= pc[np.newaxis,:] + zpoints -= zpc[np.newaxis,:] + + def getRotMat(angles): + c1, s1 = np.cos(angles[0]), np.sin(angles[0]) + c2, s2 = np.cos(angles[1]), np.sin(angles[1]) + c3, s3 = np.cos(angles[2]), np.sin(angles[2]) + return np.mat([[c1*c3-s1*c2*s3, -c1*s3-s1*c2*c3, s1*s2], + [s1*c3+c1*c2*s3, -s1*s3+c1*c2*c3, -c1*s2], + [s1*s3, s2*c3, c2]]) + + 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:] + T = getRotMat(optangles).T.A + e = (np.dot(zpoints, T)-shift[np.newaxis,:]-points) + print("Transformation angles:", optangles, flush=True) + print("Transformation shift:", shift, flush=True) + print("Transformation err:", e, flush=True) + d = np.linalg.norm(e, axis=1) + if np.any(d>1.): + QtWidgets.QMessageBox.warning(self, 'Warning!', + f'Transformation residuals are large:{d}', + QtWidgets.QMessageBox.Ok, + QtWidgets.QMessageBox.Ok) + return T, pc-shift, zpc + \ No newline at end of file diff --git a/zeissxml.py b/zeissxml.py new file mode 100644 index 0000000..2abb808 --- /dev/null +++ b/zeissxml.py @@ -0,0 +1,100 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +GEPARD - Gepard-Enabled PARticle Detection +Copyright (C) 2018 Lars Bittrich and Josef Brandt, Leibniz-Institut für +Polymerforschung Dresden e. V. + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program, see COPYING. +If not, see . +""" + +from xml.sax import make_parser, handler + +class Marker: + def __init__(self, name, x, y, z): + self.name, self.x, self.y, self.z = name, x, y, z + + def getPos(self): + return self.x, self.y, self.z + + def __repr__(self): + return str(self) + + def __str__(self): + return f'Marker: "{self.name}" at: {self.getPos()} µm' + +class Region: + def __init__(self): + self.centerx, self.centery = None, None + self.width, self.height = None, None + + def __repr__(self): + return str(self) + + def __str__(self): + return f'Region center: {self.centerx, self.centery} µm\n' + \ + f'Region size: {self.width, self.height} µm' + +class ZRange: + def __init__(self): + self.z0, self.zn = None, None + self.dz = None + + def __repr__(self): + return str(self) + + def __str__(self): + return f'Z-range in: {self.z0, self.zn} µm at step size: {self.dz} µm' + +class ZeissHandler(handler.ContentHandler): + def __init__(self): + self.markers = [] + self.region = Region() + self.zrange = ZRange() + self.intag = False + self.subtag = '' + + def characters(self, content): + if self.intag: + self.content += content.strip() + + def startElement(self, name, attrs): + if name == 'Marker': + self.markers.append(Marker(attrs['Id'], attrs['StageXPosition'], + attrs['StageYPosition'], + attrs['FocusPosition'])) + elif name == 'TileRegion' or name == 'ZStackSetup': + self.intag = True + if self.intag: + self.content = '' + if name in ['First','Last','Interval']: + self.subtag = name + + def endElement(self, name): + if name == 'TileRegion' or name == 'ZStackSetup': + self.intag = False + if self.intag and name == 'CenterPosition': + self.region.centerx, self.region.centery = map(float, self.content.split(',')) + elif self.intag and name == 'ContourSize': + self.region.width, self.region.height = map(float, self.content.split(',')) + elif self.intag and name == 'Value' and \ + self.subtag in ['First','Last','Interval']: + attrmap = {'First':'z0','Last':'zn','Interval':'dz'} + # values a given in meter even though µm are set as units? + # -> we convert to µm + self.zrange.__setattr__(attrmap[self.subtag], + float(self.content)*1.0e6) + self.subtag = '' + \ No newline at end of file -- 2.24.2