tranformation fix to separate the previous implicit image orientation in pixel sign values (f1fe95fd) · Commits · GEPARD / GEPARD

analysis/database.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -44,7 +44,7 @@ class DataBaseWindow(QtWidgets.QMainWindow):
		self.path = os.path.join(Path.home(), 'gepard', 'databases')
		self.importPath = self.path
		if not os.path.exists(self.path):
		os.mkdir(self.path)
		os.makedirs(self.path)
		self.activeDatabase = None
		self.activeSpectrum = None
		self.activeSpectrumName = None

dataset.py

+411 −386

Original line number	Diff line number	Diff line
		@@ -124,6 +124,9 @@ class DataSet(object):
		self.zpositions = [] # z-positions for optical scan
		self.heightmap = None
		self.zvalimg = None
		self.coordinatetransform = None # if imported form extern source coordinate system may be rotated
		self.signx = 1.
		self.signy = -1.

		# parameters specifically for raman scan
		self.pshift = None # shift of raman scan position relative to image center
		@@ -260,7 +263,14 @@ class DataSet(object):

		def getZval(self, pixelpos):
		assert self.zvalimg is not None
		zp = self.zvalimg[round(pixelpos[1]), round(pixelpos[0])]
		i, j = int(round(pixelpos[1])), int(round(pixelpos[0]))
		if i>=self.zvalimg.shape[0]:
		print('error in getZval:', self.zvalimg.shape, i, j)
		i = self.zvalimg.shape[0]-1
		if j>=self.zvalimg.shape[1]:
		print('error in getZval:', self.zvalimg.shape, i, j)
		j = self.zvalimg.shape[1]-1
		zp = self.zvalimg[i,j]
		z0, z1 = self.zpositions.min(), self.zpositions.max()
		return zp/255.*(z1-z0) + z0

		@@ -274,44 +284,59 @@ class DataSet(object):
		assert not self.readin
		p0 = copy(self.lastpos)

		if self.coordinatetransform is not None:
		z = 0. if len(p)<3 else p[2]
		T, pc = self.coordinatetransform
		p = (np.dot(np.array([p[0], p[1], z])-pc, T.T))

		if mode == 'df':
		p0[0] -= self.imagedim_df[0]/2
		p0[1] += self.imagedim_df[1]/2
		return (p[0] - p0[0])/self.pixelscale_df, (p0[1] - p[1])/self.pixelscale_df
		p0[0] -= self.signx*self.imagedim_df[0]/2
		p0[1] -= self.signy*self.imagedim_df[1]/2
		x, y = self.signx(p[0] - p0[0])/self.pixelscale_df, self.signy(p[1] - p0[1])/self.pixelscale_df

		elif mode == 'bf':
		p0[0] -= self.imagedim_bf[0]/2
		p0[1] += self.imagedim_bf[1]/2
		return (p[0] - p0[0])/self.pixelscale_bf, (p0[1] - p[1])/self.pixelscale_bf
		p0[0] -= self.signx*self.imagedim_bf[0]/2
		p0[1] += self.signy*self.imagedim_bf[1]/2
		x, y = self.signx(p[0] - p0[0])/self.pixelscale_bf, self.signy(p[1] - p0[1])/self.pixelscale_bf

		else:
		print('mapToPixelMode not understood')
		return
		raise ValueError(f'mapToPixel mode: {mode} not understood')

		return x, y

		def mapToLength(self, pixelpos, mode='df', force=False):
		def mapToLength(self, pixelpos, mode='df', force=False, returnz=False):
		if not force:
		assert not self.readin
		p0 = copy(self.lastpos)
		p0[0] += self.coordOffset[0]
		p0[1] += self.coordOffset[1]

		if mode == 'df':
		p0[0] -= self.imagedim_df[0]/2
		p0[1] += self.imagedim_df[1]/2
		return (pixelpos[0]self.pixelscale_df + p0[0]), (p0[1] - pixelpos[1]self.pixelscale_df)
		p0[0] -= self.signx*self.imagedim_df[0]/2
		p0[1] -= self.signy*self.imagedim_df[1]/2
		x, y = (self.signxpixelpos[0]self.pixelscale_df + p0[0]), (p0[1] + self.signypixelpos[1]self.pixelscale_df)
		elif mode == 'bf':
		p0[0] -= self.imagedim_bf[0]/2
		p0[1] += self.imagedim_bf[1]/2
		return (pixelpos[0]self.pixelscale_bf + p0[0]), (p0[1] - pixelpos[1]self.pixelscale_bf)
		p0[0] -= self.signx*self.imagedim_bf[0]/2
		p0[1] -= self.signy*self.imagedim_bf[1]/2
		x, y = (self.signxpixelpos[0]self.pixelscale_bf + p0[0]), (p0[1] + self.signypixelpos[1]self.pixelscale_bf)
		else:
		raise ValueError(f'mapToLength mode: {mode} not understood')

		def mapToLengthRaman(self, pixelpos, microscopeMode='df', noz=False):
		p0x, p0y = self.mapToLength(pixelpos, mode = microscopeMode)
		x, y = p0x + self.pshift[0], p0y + self.pshift[1]
		z = None
		if not noz:
		if (returnz and self.zvalimg is not None) or self.coordinatetransform is not None:
		z = self.mapHeight(x, y)
		z += self.getZval(pixelpos)

		if self.coordinatetransform is not None:
		T, pc = self.coordinatetransform
		x, y, z = (np.dot(np.array([x,y,z]), T) + pc)

		if returnz:
		return x, y, z
		return x, y

		def mapToLengthRaman(self, pixelpos, microscopeMode='df', noz=False):
		p0x, p0y, z = self.mapToLength(pixelpos, mode=microscopeMode, returnz=True)
		x, y = p0x + self.pshift[0], p0y + self.pshift[1]
		return x, y, z

		def newProject(self, fname):

helperfunctions.py

+11 −0

Original line number	Diff line number	Diff line
		@@ -23,7 +23,16 @@ import numpy as np
		import cv2
		import os

		try:
		from skimage.io import imread as skimread
		from skimage.io import imsave as skimsave
		except ImportError:
		skimread = None
		skimsave = None

		def cv2imread_fix(fname, flags=cv2.IMREAD_COLOR):
		if skimread is not None:
		return skimread(fname, as_gray=(flags==cv2.IMREAD_GRAYSCALE))
		with open(fname, "rb") as fp:
		cont = fp.read()
		img = cv2.imdecode(np.fromstring(cont, dtype=np.uint8), flags)
		@@ -31,6 +40,8 @@ def cv2imread_fix(fname, flags=cv2.IMREAD_COLOR):
		return None

		def cv2imwrite_fix(fname, img, params=None):
		if skimsave is not None:
		skimsave(fname, img)
		pathname, ext = os.path.splitext(fname)
		if params is None:
		ret, data = cv2.imencode(ext, img)

zeissimporter.py

+258 −265

Original line number	Diff line number	Diff line
		@@ -59,6 +59,10 @@ class ZeissImporter(QtWidgets.QDialog):
		self.pexit.released.connect(self.reject)
		self.pconvert.setEnabled(False)

		self.xinvert = QtWidgets.QCheckBox('Invert x-axis')
		self.yinvert = QtWidgets.QCheckBox('Invert y-axis')
		self.zinvert = QtWidgets.QCheckBox('Invert z-axis')

		btnLayout = QtWidgets.QHBoxLayout()
		btnLayout.addStretch()
		btnLayout.addWidget(self.pconvert)
		@@ -77,22 +81,30 @@ class ZeissImporter(QtWidgets.QDialog):

		vbox.addWidget(pointgroup)
		vbox.addLayout(blurlayout)
		vbox.addWidget(self.xinvert)
		vbox.addWidget(self.yinvert)
		vbox.addWidget(self.zinvert)
		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')
		self.edfimgname, self.zmapimgname, xmlname = '', '', ''
		for name in os.listdir(path):
		if name.lower().endswith('_meta.xml'):
		xmlname = os.path.join(path, name)
		elif name.lower().endswith('_c1.tif'):
		self.edfimgname = os.path.join(path, name)
		elif name.lower().endswith('_c2.tif'):
		self.zmapimgname = os.path.join(path, name)

		errmsges = []
		if not os.path.exists(self.zmapimgname):
		errmsges.append('Depth map image not found: 3D.tif')
		errmsges.append('Depth map image not found: NAME_c2.tif')
		if not os.path.exists(self.edfimgname):
		errmsges.append('EDF image not found: EDF.tif')
		errmsges.append('EDF image not found: NAME_c1.tif')
		if not os.path.exists(xmlname):
		errmsges.append('XML metadata not found: 3D.tif_metadata.xml')
		errmsges.append('XML metadata not found: NAME_meta.xml')
		else:
		parser = make_parser()
		z = ZeissHandler()
		@@ -108,7 +120,7 @@ class ZeissImporter(QtWidgets.QDialog):
		errmsges.append('ZStack information missing or incomplete!')

		if len(errmsges)>0:
		QtWidgets.QMessageBox.error(self, 'Error!',
		QtWidgets.QMessageBox.critical(self, 'Error!',
		'\n'.join(errmsges),
		QtWidgets.QMessageBox.Ok,
		QtWidgets.QMessageBox.Ok)
		@@ -117,6 +129,9 @@ class ZeissImporter(QtWidgets.QDialog):
		self.region = z.region
		self.zrange = z.zrange
		self.markers = z.markers
		print(self.region)
		print(self.zrange)
		print(self.markers, flush=True)
		return True

		@QtCore.pyqtSlot(float, float, float)
		@@ -127,32 +142,33 @@ class ZeissImporter(QtWidgets.QDialog):

		@QtCore.pyqtSlot()
		def convert(self):
		T, pc, zpc, accept = self.getTransform()
		if accept:
		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)
		self.convertZimg(dataset, T, pc, zpc)
		self.convertImage(dataset)
		dataset.save()
		self.gepardname = dataset.fname
		self.accept()

		def convertImage(self, dataset, warp_mat):
		def convertImage(self, dataset):
		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]
		N = int(round(abs(self.zrange.zn-self.zrange.z0)/self.zrange.dz))
		z0, zn = self.zrange.z0, self.zrange.zn
		if zn<z0:
		zn, z0 = z0, zn
		dataset.zpositions = np.linspace(z0, zn, N)-zpc[2]
		dataset.heightmap = np.zeros(3)
		dataset.signy = 1.
		zimg = cv2imread_fix(self.zmapimgname, cv2.IMREAD_GRAYSCALE)
		zmdist = zimg.mean()
		zm = zmdist/255.*(self.zrange.zn-self.zrange.z0) + self.zrange.z0

		print("zimg shape:", zimg.shape, flush=True)
		radius = self.blurspinbox.value()
		blur = cv2.GaussianBlur(zimg, (radius, radius), 0)

		@@ -167,53 +183,32 @@ class ZeissImporter(QtWidgets.QDialog):
		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.coordinatetransform = T, pc

		# set image center as reference point in data set (use Zeiss coordinates)
		p0 = np.array([self.region.centerx, self.region.centery]) - zpc[: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))
		zmin, zmax = dataset.zpositions.min(), dataset.zpositions.max()
		blur = (blur)*(255.)
		blur[blur>255.] = 255.
		blur[np.isnan(blur)] = 0.
		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
		dataset.zvalimg = "saved"

		def getTransform(self):
		points = self.points.getPoints()
		pshift = self.ramanctrl.getRamanPositionShift()
		points[:,0] -= pshift[0]
		points[:,1] -= pshift[1]
		Parity = np.mat(np.diag([-1. if self.xinvert.isChecked() else 1.,
		-1. if self.yinvert.isChecked() else 1.,
		-1. if self.zinvert.isChecked() else 1.]))
		zpoints = np.array([m.getPos() for m in self.markers], dtype=np.double)
		pc = points.mean(axis=0)
		zpc = zpoints.mean(axis=0)
		@@ -233,11 +228,10 @@ class ZeissImporter(QtWidgets.QDialog):
		# [-45°,45°] for all permutation of markers
		permbest = None
		pointsbest = None
		for perm in permutations(range(points.shape[0])):
		ppoints = points[perm,:]
		ppoints = points[:,:].copy()

		def err(angles_shift):
		T = getRotMat(angles_shift[:3]).T.A
		T = (getRotMat(angles_shift[:3]).T*Parity).A
		return (np.dot(zpoints, T) - angles_shift[np.newaxis,3:] \
		- ppoints).ravel()

		@@ -246,25 +240,24 @@ class ZeissImporter(QtWidgets.QDialog):
		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
		T = (getRotMat(optangles).T*Parity).A
		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)
		d = np.linalg.norm(e, axis=1)
		accept = True
		if np.any(d>1.):
		QtWidgets.QMessageBox.warning(self, 'Warning!',
		ret = QtWidgets.QMessageBox.warning(self, 'Warning!',
		f'Transformation residuals are large:{d}',
		QtWidgets.QMessageBox.Ok,
		QtWidgets.QMessageBox.Ok\|QtWidgets.QMessageBox.Cancel,
		QtWidgets.QMessageBox.Ok)
		return T, pc-shift, zpc
		if ret==QtWidgets.QMessageBox.Cancel:
		accept = False
		return T, pc-shift, zpc, accept