Commit 4da0119f authored by kulvait's avatar kulvait

Improved writing functions to den in Python

parent 2897b915
......@@ -13,6 +13,7 @@ import numpy as np
import os
import cv2
def timeToSeconds(naivetime):
strtime = str(naivetime).split('.', 1)[0]
secondfraction = float("0.%s" % str(naivetime).split('.', 1)[1])
......@@ -64,19 +65,18 @@ def storeNdarrayAsDEN(fileName, dataFrame, force=False):
raise ValueError(
'Dimension of dataFrame should be 2 or 3 but is %d.' % len(dataFrame.shape))
shape = dataFrame.shape # Now len is for sure 3
writeEmptyDEN(fileName, dimx=shape[0], dimy=shape[1], dimz=shape[2],
rows = shape[0]
columns = shape[1]
writeEmptyDEN(fileName, dimx=shape[1], dimy=shape[0], dimz=shape[2],
force=force) # No effect
header = np.fromfile(fileName, np.dtype('<i2'), 3)
rows = np.uint32(header[0])
columns = np.uint32(header[1])
f = open(file, "wb")
header.tofile(f)
for frame in range(i):
newdata = np.array(data[:, :, frame], np.dtype('<f4'))
toWrite = dataFrame.astype(np.float32)
f = open(fileName, "r+b")
for frame in range(shape[2]):
newdata = np.array(toWrite[:, :, frame], np.dtype('<f4'))
newdata = newdata.reshape((rows * columns, ))
# put header in front of image data
f.seek(6 + rows * columns * frame * 4, os.SEEK_SET)
newdata.tofile(f)
f.write(newdata.tobytes())
f.close()
return True
......@@ -104,8 +104,8 @@ def writeFrame(fileName, k, data, force=False):
def writeEmptyDEN(fileName, dimx, dimy, dimz, force=False):
if not force and os.path.exists(fileName):
raise IOError('File %s already exists, no header written' % fileName)
outfile = open(file, "w")
header = np.array([dimx, dimy, dimz])
outfile = open(fileName, "w")
header = np.array([dimy, dimx, dimz])
header = np.array(header, dtype='<i2')
header.tofile(outfile)
fileSize = dimx * dimy * dimz * 4 + 6
......
......@@ -14,32 +14,43 @@ import glob
import os
import pydicom
def getDicomFiles(d, suffixes=["IMA"]):
dicomFiles = []
for s in suffixes:
dicomFiles.extend(glob.glob(os.path.join(d, "*.%s"%s)))
dicomFiles.extend(glob.glob(os.path.join(d, "*.%s" % s)))
dicomFiles = glob.glob(os.path.join(d, "*.IMA"))
dicomFiles.sort()
return dicomFiles
#Get list of pydicom objects for each file with given suffix in the directory
# Get list of pydicom objects for each file with given suffix in the directory
def getDicoms(d, suffixes=["IMA"]):
dicomFiles = getDicomFiles(d, suffixes);
dicomFiles = getDicomFiles(d, suffixes)
dicoms = [pydicom.read_file(x) for x in dicomFiles]
return dicoms
def getPtsIds(dicoms):
patientIDs = list(set([x.PatientID for x in dicoms if "PatientID" in x.dir("PatientID")]))
patientIDs = list(
set([x.PatientID for x in dicoms if "PatientID" in x.dir("PatientID")]))
return patientIDs
def getStudyDates(dicoms):
patientIDs = list(set([x.PatientID for x in dicoms if "StudyDate" in x.dir("StudyDate")]))
patientIDs = list(
set([x.PatientID for x in dicoms if "StudyDate" in x.dir("StudyDate")]))
return patientIDs
def getSeriesDescriptions(dicoms):
IDs = list(set([x.SeriesDescription for x in dicoms if "SeriesDescription" in x.dir("SeriesDescription")]))
IDs = list(
set([x.SeriesDescription for x in dicoms if "SeriesDescription" in x.dir("SeriesDescription")]))
return IDs
def getInstanceUIDs(dicoms):
IDs = list(set([x.SeriesInstanceUID for x in dicoms if "SeriesDescription" in x.dir("SeriesDescription")]))
IDs = list(
set([x.SeriesInstanceUID for x in dicoms if "SeriesDescription" in x.dir("SeriesDescription")]))
return IDs
......@@ -13,67 +13,86 @@ import os
# Estimate rigid transformation
# Most naive approach
def getRigidTransformation(frame1, frame2):
minvalue = min(np.percentile(frame1, 10), np.percentile(frame2, 10))
maxvalue = min(np.percentile(frame1, 90), np.percentile(frame2, 90))
frame1 = (frame1 - minvalue) / (maxvalue - minvalue)
frame2 = (frame2 - minvalue) / (maxvalue - minvalue)
frame1 = np.maximum(0, np.minimum(1, frame1))
frame2 = np.maximum(0, np.minimum(1, frame2))
warp_matrix = cv2.estimateRigidTransform((frame1 * 255 + 0.5).astype('uint8'), (frame2 * 255 + 0.5).astype('uint8'), False)
return warp_matrix
minvalue = min(np.percentile(frame1, 10), np.percentile(frame2, 10))
maxvalue = min(np.percentile(frame1, 90), np.percentile(frame2, 90))
frame1 = (frame1 - minvalue) / (maxvalue - minvalue)
frame2 = (frame2 - minvalue) / (maxvalue - minvalue)
frame1 = np.maximum(0, np.minimum(1, frame1))
frame2 = np.maximum(0, np.minimum(1, frame2))
warp_matrix = cv2.estimateRigidTransform(
(frame1 * 255 + 0.5).astype('uint8'), (frame2 * 255 + 0.5).astype('uint8'), False)
return warp_matrix
# Estimate of corelation coefficient
# Gets result after first step
# There is no function for direct computation
def estimateCC(frame1, frame2):
warp_mode = cv2.MOTION_EUCLIDEAN
warp_matrix = np.eye(2, 3, dtype=np.float32)
max_iterations = 1
termination_eps = 1e-15
criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, max_iterations, termination_eps)
(cc, warp_matrix) = cv2.findTransformECC(frame1, frame2, warp_matrix, warp_mode, criteria)
return (cc, warp_matrix)
warp_mode = cv2.MOTION_EUCLIDEAN
warp_matrix = np.eye(2, 3, dtype=np.float32)
max_iterations = 1
termination_eps = 1e-15
criteria = (cv2.TERM_CRITERIA_EPS |
cv2.TERM_CRITERIA_COUNT, max_iterations, termination_eps)
(cc, warp_matrix) = cv2.findTransformECC(
frame1, frame2, warp_matrix, warp_mode, criteria)
return (cc, warp_matrix)
# Get optimal transformation
def getTransformECC(frame1, frame2, warp_mode=None, max_iterations=None, termination_eps=None, init_warp_matrix=None):
if warp_mode is None:
warp_mode = cv2.MOTION_EUCLIDEAN
if max_iterations is None:
max_iterations = 5000 # Default of findTransformECC is 50
if termination_eps is None:
termination_eps = 1e-10
# Default of findTransformECC is 0.001 if init_warp_matrix is None
if warp_mode == cv2.MOTION_HOMOGRAPHY:
init_warp_matrix = np.eye(3, 3, dtype=np.float32)
else:
init_warp_matrix = np.eye(2, 3, dtype=np.float32)
criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, max_iterations, termination_eps)
(cc, warp_matrix) = cv2.findTransformECC(frame1, frame2, init_warp_matrix, warp_mode, criteria)
return (cc, warp_matrix)
if warp_mode is None:
warp_mode = cv2.MOTION_EUCLIDEAN
if max_iterations is None:
max_iterations = 5000 # Default of findTransformECC is 50
if termination_eps is None:
termination_eps = 1e-10
# Default of findTransformECC is 0.001 if init_warp_matrix is None
if warp_mode == cv2.MOTION_HOMOGRAPHY:
init_warp_matrix = np.eye(3, 3, dtype=np.float32)
else:
init_warp_matrix = np.eye(2, 3, dtype=np.float32)
criteria = (cv2.TERM_CRITERIA_EPS |
cv2.TERM_CRITERIA_COUNT, max_iterations, termination_eps)
(cc, warp_matrix) = cv2.findTransformECC(frame1,
frame2, init_warp_matrix, warp_mode, criteria)
return (cc, warp_matrix)
# Calculate the x and y gradients using Sobel operator
# Then using Euclidean norm of both
def getGradient(frame):
grad_x = cv2.Sobel(frame, cv2.CV_32F, 1, 0, ksize=3)
grad_y = cv2.Sobel(frame, cv2.CV_32F, 0, 1, ksize=3)
return np.sqrt(np.add(np.square(grad_x), np.square(grad_y)))
grad_x = cv2.Sobel(frame, cv2.CV_32F, 1, 0, ksize=3)
grad_y = cv2.Sobel(frame, cv2.CV_32F, 0, 1, ksize=3)
return np.sqrt(np.add(np.square(grad_x), np.square(grad_y)))
# Calculate Laplacian
def getLaplacian(frame):
return (cv2.Laplacian(frame, cv2.CV_32F, ksize=23))
return (cv2.Laplacian(frame, cv2.CV_32F, ksize=23))
# Apply the same algorithm on gradient data
def gradientTransformECC(frame1, frame2, warp_mode=None, max_iterations=None, termination_eps=None, init_warp_matrix=None):
g1 = getGradient(frame1)
g2 = getGradient(frame2)
return getTransformECC(g1, g2, warp_mode, max_iterations, termination_eps,init_warp_matrix)
g1 = getGradient(frame1)
g2 = getGradient(frame2)
return getTransformECC(g1, g2, warp_mode, max_iterations, termination_eps, init_warp_matrix)
# Apply the same algorithm on laplacian data
def laplacianTransformECC(frame1, frame2, warp_mode=None, max_iterations=None, termination_eps=None, init_warp_matrix=None):
l1 = getLaplacian(frame1)
l2 = getLaplacian(frame2)
return getTransformECC(l1, l2, warp_mode, max_iterations, termination_eps, init_warp_matrix)
l1 = getLaplacian(frame1)
l2 = getLaplacian(frame2)
return getTransformECC(l1, l2, warp_mode, max_iterations, termination_eps, init_warp_matrix)
# For solving the problem that we have warp_matrix in the coordinates(j, i)
# here(x, y) = (col, row)
......@@ -82,24 +101,32 @@ def laplacianTransformECC(frame1, frame2, warp_mode=None, max_iterations=None, t
# We have that(Ax) ' = Ax' - Aoffset + offset
# Here offset is(offset_row, offset_col, 0)
# So the operator that acts on(i ',j') space has modified third column.
def shiftWarpMatrix(warp_matrix, offset_row, offset_col):
output_matrix = np.copy(warp_matrix)
offsetVector = np.array([offset_col, offset_row, 0], dtype=np.float32)
thirdRowAdd = -np.matmul(warp_matrix, offsetVector) + offsetVector[0:warp_matrix.shape[0]]
output_matrix[:, 2] = output_matrix[:, 2] + thirdRowAdd
return output_matrix
output_matrix = np.copy(warp_matrix)
offsetVector = np.array([offset_col, offset_row, 0], dtype=np.float32)
thirdRowAdd = - \
np.matmul(warp_matrix, offsetVector) + \
offsetVector[0:warp_matrix.shape[0]]
output_matrix[:, 2] = output_matrix[:, 2] + thirdRowAdd
return output_matrix
# Perform the 'backtransformation' of the frame to the template space
# It seems that this function takes(x, y) indexing(cols, rows)
# Interpolation mode see
# Interpolation mode see
#
#https://docs.opencv.org/3.1.0/da/d54/group__imgproc__transform.html#gga5bb5a1fea74ea38e1a5445ca803ff121aa5521d8e080972c762467c45f3b70e6c
# https://docs.opencv.org/3.1.0/da/d54/group__imgproc__transform.html#gga5bb5a1fea74ea38e1a5445ca803ff121aa5521d8e080972c762467c45f3b70e6c
def warpFrame(frame, warp_matrix, warp_mode, interpolation_mode=None):
if interpolation_mode is None:
interpolation_mode = cv2.INTER_LANCZOS4
(rows, cols) = frame.shape
if warp_mode == cv2.MOTION_HOMOGRAPHY:
frameAligned = cv2.warpPerspective(frame, warp_matrix, (cols, rows), flags=interpolation_mode + cv2.WARP_INVERSE_MAP)
else:
frameAligned = cv2.warpAffine(frame, warp_matrix, (cols, rows), flags=interpolation_mode + cv2.WARP_INVERSE_MAP)
return frameAligned
if interpolation_mode is None:
interpolation_mode = cv2.INTER_LANCZOS4
(rows, cols) = frame.shape
if warp_mode == cv2.MOTION_HOMOGRAPHY:
frameAligned = cv2.warpPerspective(frame, warp_matrix, (
cols, rows), flags=interpolation_mode + cv2.WARP_INVERSE_MAP)
else:
frameAligned = cv2.warpAffine(frame, warp_matrix, (
cols, rows), flags=interpolation_mode + cv2.WARP_INVERSE_MAP)
return frameAligned
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment