Commit 5fd1cb0c authored by Lars Bittrich's avatar Lars Bittrich
Browse files

Zeissimporter initial commit with coordinate transfer

parent eb5a3c56
Loading
Loading
Loading
Loading
+22 −1
Original line number Original line Diff line number Diff line
@@ -85,6 +85,16 @@ class MeasureParticleWindow(QtWidgets.QMainWindow):
            self.view.open(self.fname)  
            self.view.open(self.fname)  
            self.scalingChanged(1.)
            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() 
    @QtCore.pyqtSlot() 
    def new(self, fileName=False):
    def new(self, fileName=False):
        if fileName is False:
        if fileName is False:
@@ -115,6 +125,10 @@ class MeasureParticleWindow(QtWidgets.QMainWindow):
        self.openAct.setShortcut("Ctrl+O")
        self.openAct.setShortcut("Ctrl+O")
        self.openAct.triggered.connect(self.open)
        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 = QtWidgets.QAction("&New Measurement...", self)
        self.newAct.setShortcut("Ctrl+N")
        self.newAct.setShortcut("Ctrl+N")
        self.newAct.triggered.connect(self.new)
        self.newAct.triggered.connect(self.new)
@@ -211,12 +225,14 @@ class MeasureParticleWindow(QtWidgets.QMainWindow):
    
    
    def unblockUI(self, connected):
    def unblockUI(self, connected):
        self.openAct.setEnabled(True)
        self.openAct.setEnabled(True)
        self.importAct.setEnabled(True)
        self.newAct.setEnabled(True)
        self.newAct.setEnabled(True)
        self.updateConnected(connected)
        self.updateConnected(connected)
        self.exitAct.setEnabled(True)
        self.exitAct.setEnabled(True)
    
    
    def blockUI(self):
    def blockUI(self):
        self.openAct.setEnabled(False)
        self.openAct.setEnabled(False)
        self.importAct.setEnabled(False)
        self.newAct.setEnabled(False)
        self.newAct.setEnabled(False)
        self.connectRamanAct.setEnabled(False)
        self.connectRamanAct.setEnabled(False)
        self.disconnectRamanAct.setEnabled(False)
        self.disconnectRamanAct.setEnabled(False)
@@ -237,6 +253,7 @@ class MeasureParticleWindow(QtWidgets.QMainWindow):
    def createMenus(self):
    def createMenus(self):
        self.fileMenu = QtWidgets.QMenu("&File", self)
        self.fileMenu = QtWidgets.QMenu("&File", self)
        self.fileMenu.addAction(self.newAct)
        self.fileMenu.addAction(self.newAct)
        self.fileMenu.addAction(self.importAct)
        self.fileMenu.addAction(self.openAct)
        self.fileMenu.addAction(self.openAct)
        self.fileMenu.addSeparator()
        self.fileMenu.addSeparator()
        self.fileMenu.addAction(self.exitAct)
        self.fileMenu.addAction(self.exitAct)
@@ -272,6 +289,7 @@ class MeasureParticleWindow(QtWidgets.QMainWindow):
        self.toolbar.setIconSize(QtCore.QSize(100,50))
        self.toolbar.setIconSize(QtCore.QSize(100,50))
        self.toolbar.addAction(self.aboutAct)
        self.toolbar.addAction(self.aboutAct)
        self.toolbar.addAction(self.newAct)
        self.toolbar.addAction(self.newAct)
        self.toolbar.addAction(self.importAct)
        self.toolbar.addAction(self.openAct)
        self.toolbar.addAction(self.openAct)
        self.toolbar.addSeparator()
        self.toolbar.addSeparator()
        self.toolbar.addAction(self.connectRamanAct)
        self.toolbar.addAction(self.connectRamanAct)
@@ -295,6 +313,7 @@ if __name__ == '__main__':


    logpath = QtCore.QStandardPaths.writableLocation(
    logpath = QtCore.QStandardPaths.writableLocation(
              QtCore.QStandardPaths.AppLocalDataLocation)
              QtCore.QStandardPaths.AppLocalDataLocation)
    fp = None
    if logpath != "":
    if logpath != "":
        if not os.path.exists(logpath):
        if not os.path.exists(logpath):
            os.mkdir(logpath)
            os.mkdir(logpath)
@@ -310,3 +329,5 @@ if __name__ == '__main__':
    meas = MeasureParticleWindow(logpath)
    meas = MeasureParticleWindow(logpath)
    meas.showMaximized()
    meas.showMaximized()
    ret = app.exec_()
    ret = app.exec_()
    if fp is not None:
        fp.close()
 No newline at end of file
+6 −2
Original line number Original line Diff line number Diff line
@@ -113,11 +113,12 @@ def loadAndPasteImage(srcnames, fullimage, fullzval, width, height,
class PointCoordinates(QtWidgets.QGridLayout):
class PointCoordinates(QtWidgets.QGridLayout):
    readPoint = QtCore.pyqtSignal(float, float, float, name='readPoint')
    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)
        super().__init__(parent)
        self.dswidgets = []
        self.dswidgets = []
        self.N = 0
        self.N = 0
        self.ramanctrl = ramanctrl
        self.ramanctrl = ramanctrl
        self.names = names
        self.pimageOnly = QtWidgets.QPushButton("Image")
        self.pimageOnly = QtWidgets.QPushButton("Image")
        self.addWidget(self.pimageOnly, 0, 6, QtCore.Qt.AlignRight)
        self.addWidget(self.pimageOnly, 0, 6, QtCore.Qt.AlignRight)
        self.pimageOnly.released.connect(QtCore.pyqtSlot()(lambda : self.read(-1)))
        self.pimageOnly.released.connect(QtCore.pyqtSlot()(lambda : self.read(-1)))
@@ -140,6 +141,9 @@ class PointCoordinates(QtWidgets.QGridLayout):
            self.itemAtPosition(i+1,5).setVisible(True)
            self.itemAtPosition(i+1,5).setVisible(True)
            self.itemAtPosition(i+1,6).setVisible(True)
            self.itemAtPosition(i+1,6).setVisible(True)
        for i in range(self.N, N):
        for i in range(self.N, N):
            if self.names is not None:
                lx = QtWidgets.QLabel(f"{self.names[i]} -> x:")
            else:
                lx = QtWidgets.QLabel(f"{i+1} -> x:")
                lx = QtWidgets.QLabel(f"{i+1} -> x:")
            ly = QtWidgets.QLabel("y:")
            ly = QtWidgets.QLabel("y:")
            lz = QtWidgets.QLabel("z:")
            lz = QtWidgets.QLabel("z:")
+54 −46
Original line number Original line Diff line number Diff line
@@ -27,6 +27,7 @@ from opticalscan import OpticalScan
from ramanscanui import RamanScanUI
from ramanscanui import RamanScanUI
from detectionview import ParticleDetectionView
from detectionview import ParticleDetectionView
from analysis.analysisview import ParticleAnalysis
from analysis.analysisview import ParticleAnalysis
from zeissimporter import ZeissImporter
from viewitems import FitPosIndicator, Node, Edge, ScanIndicator, RamanScanIndicator, SegmentationContours
from viewitems import FitPosIndicator, Node, Edge, ScanIndicator, RamanScanIndicator, SegmentationContours
from helperfunctions import polygoncovering, cv2imread_fix
from helperfunctions import polygoncovering, cv2imread_fix
import cv2
import cv2
@@ -237,6 +238,13 @@ class SampleView(QtWidgets.QGraphicsView):
        self.activateMaxMode(loadnew=True)
        self.activateMaxMode(loadnew=True)
        self.imparent.snapshotAct.setEnabled(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):
    def new(self, fname):
        self.saveDataSet()
        self.saveDataSet()
        if self.dataset is not None:
        if self.dataset is not None:

zeissimporter.py

0 → 100644
+255 −0
Original line number Original line Diff line number Diff line
# -*- coding: utf-8 -*-
"""
GEPARD - Gepard-Enabled PARticle Detection
Copyright (C) 2018  Lars Bittrich and Josef Brandt, Leibniz-Institut für 
Polymerforschung Dresden e. V. <bittrich-lars@ipfdd.de>    

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 <https://www.gnu.org/licenses/>.
"""
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

zeissxml.py

0 → 100644
+100 −0
Original line number Original line Diff line number Diff line
#!/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. <bittrich-lars@ipfdd.de>    

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 <https://www.gnu.org/licenses/>.
"""

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