Commit 04e3fad4 authored by Mark Hymers's avatar Mark Hymers

Split apart coil and transform recalculations

Signed-off-by: Mark Hymers's avatarMark Hymers <mark.hymers@ynic.york.ac.uk>
parent fb77dc04
Pipeline #25262 passed with stage
in 7 minutes and 19 seconds
......@@ -22,15 +22,22 @@ def usage():
print_transforms FILENAME
Prints out the transforms from the file.
recalculate_transforms [--debug] [--write] FILENAME CALIBRATIONFILE COHFILE [COILS_TO_USE]
Recalculate the ccs_to_scs_transform for each COH run in the file.
A before and after transform report will be printed.
recalculate_coils [--debug] [--write] FILENAME CALIBRATIONFILE COHFILE [COILS_TO_USE]
Re-fit coil locations and the ccs_to_scs_transform for each COH run in
the file. A before and after transform report will be printed.
If --debug is set, debug output from the fitting routine will be printed.
If --write is set, the ccs_to_scs_transform and fitted_coils will be update
in the given file. Use this with care.
recalculate_transforms [--write] FILENAME [COILS_TO_USE]
Recalculate the ccs_to_scs_transform for each COH run in the file.
A before and after transform report will be printed.
If --write is set, the ccs_to_scs_transform and fitted_coils will be update
in the given file. Use this with care.
copy_transforms [--force] SOURCEFILE TARGETFILE
Copy ccs_to_scs_transform and fitted_coils from one file to another for
all acquisitions.
......@@ -63,7 +70,7 @@ def print_transforms(args):
operations['print_transforms'] = print_transforms
def recalculate_transforms(args):
def recalculate_coils(args):
debug = False
write = False
......@@ -101,10 +108,59 @@ def recalculate_transforms(args):
coils_to_use_inp = None
# Call the relevant routine
fit_data, coilsused, report = proc.recalculate_transforms(filename,
calibration_filename,
coh_filename,
coils_to_use_inp)
fit_data, coilsused, report = proc.recalculate_coils(filename,
calibration_filename,
coh_filename,
coils_to_use_inp)
print(report)
if write:
print("Please read the before and after reports before continuing")
ret = input('Are you sure that you wish to overwrite the existing transforms (y/N)? ')
if ret.lower() != 'y':
print("\nNot saving")
exit(0)
proc.rewrite_coils_and_transforms(filename, fit_data, coilsused)
operations['recalculate_coils'] = recalculate_coils
def recalculate_transforms(args):
write = False
# Options must be first in the list
while True:
if args[0].startswith('--'):
# Parse argument
if args[0] == '--write':
write = True
else:
print("Unknown option {}".format(args[0]))
exit(1)
args.pop(0)
else:
break
# If we don't have at least one argument, just show usage
if len(args) < 1:
usage()
filename = args[0]
if len(args) > 1:
coils_to_use_inp = sorted(args[1].split(','))
else:
coils_to_use_inp = None
# Call the relevant routine
fit_data, coilsused, report = proc.recalculate_transform(filename,
coils_to_use_inp)
print(report)
......@@ -116,6 +172,7 @@ def recalculate_transforms(args):
print("\nNot saving")
exit(0)
# Only replace the transforms
proc.rewrite_transforms(filename, fit_data, coilsused)
......@@ -123,6 +180,7 @@ def recalculate_transforms(args):
operations['recalculate_transforms'] = recalculate_transforms
def copy_transforms(args):
# Options must be first in the list
force = False
......
......@@ -86,3 +86,40 @@ def store_coh_transforms(coh_transforms, datafile):
if acq.attrs.get('acq_type') == 'ACQ' and last_trans is not None:
key = '{acq_name}/ccs_to_scs_transform'.format(acq_name=name)
acq[key] = SoftLink(last_trans)
def update_transforms(coh_transforms, datafile):
"""
Update only the COH transform matrices
"""
index = 0
last_trans = None
for acq in get_ordered_acquisitions(datafile):
name = acq.name
if acq.attrs.get('acq_type') == 'COH':
try:
coh = coh_transforms[index]
except IndexError:
coh = None
logging.warning('COH transform data for index {} does not exist'.format(index))
if coh:
# Create the transform dataset in the hdf5
transform = np.array(coh['transform'])
acq.create_dataset(
name='ccs_to_scs_transform',
data=transform,
dtype='float64',
fletcher32=True,
)
index += 1
# If we have a previously calculated a ccs to scs transform in a
# COH acquisition, create a softlink to it in any subsequent
# 'normal' acquisitions.
if acq.attrs.get('acq_type') == 'ACQ' and last_trans is not None:
key = '{acq_name}/ccs_to_scs_transform'.format(acq_name=name)
acq[key] = SoftLink(last_trans)
......@@ -6,8 +6,9 @@ import logging
import h5py
import numpy as np
from yias_coh.alg03.position import RigidBodyTransform
from yias_coh.coil_positioning import calculate_coil_positions
from yias_coh.data_processing import store_coh_transforms
from yias_coh.data_processing import store_coh_transforms, update_transforms
ARRAY_FORMATTER = {'float_kind': lambda x: "{:8.3f}".format(x)}
......@@ -273,8 +274,158 @@ def extract_transform_report_data(f):
return all_acqs
def recalculate_transforms(filename, calibration_filename, coh_filename,
coils_to_use_inp, debug=False):
def recalculate_transform(filename, coils_to_use_inp, debug=False):
"""
Recalculate ccs_to_scs_transform for a dataset
Parameters
----------
filename
Filename of HDF5 dataset
coils_to_use_inp:
List of coil names to use for fitting (or None if all)
Returns
-------
coils_to_use:
List of coil names which were used
report:
String containing report information
"""
ret = []
# Grab the HDF5 file to read information from
with h5py.File(filename, 'r') as f:
# Load the data into a form suitable for generate_transform_report
# This is used later on
all_acqs = extract_transform_report_data(f)
# Generate the "before" report
ret.append("================================================================")
ret.append(" Before recalculation report ")
ret.append("================================================================")
ret.append(generate_transform_report(all_acqs))
ret.append('')
# Read digitisation coils from HDF5
digitisation_coils = {'sequence': [],
'data': {},
'total_num': 0}
coh_coils = f['geometry']['coils']
coils_to_use = []
for coil_name in sorted(coh_coils.keys()):
full_name = 'coil_{}'.format(coil_name)
dat = {'store': coil_name,
'units': 'mm',
'points': coh_coils[coil_name]['location'][...].tolist(),
'count': coh_coils[coil_name]['location'].shape[0],
'text_point': full_name
}
digitisation_coils['sequence'].append(full_name)
digitisation_coils['data'][full_name] = dat
digitisation_coils['total_num'] += 1
if coils_to_use_inp is not None:
if coil_name in coils_to_use_inp:
coils_to_use.append(full_name)
if coils_to_use_inp is None:
coils_to_use = None
else:
if len(coils_to_use) < 3:
raise Exception("Error: Need at least three coils to use: only have {}".format(coils_to_use))
ret.append("=======================================================")
ret.append("NOTE: Only using coils {} to fit transform".format(','.join([str(x) for x in coils_to_use])))
ret.append("=======================================================")
if coils_to_use is None:
coils_to_use_idx = np.arange(len(digitisation_coils['sequence']))
else:
coils_to_use_idx = []
for cname in coils_to_use:
try:
coils_to_use_idx.append(digitisation_coils['sequence'].index(cname))
except ValueError:
raise ValueError("Coil {} not found in digitisation_coils. Exiting.".format(cname))
coils_to_use_idx = np.array(coils_to_use_idx)
# For each COH acquisition, recalculate the transforms
acqs = f['acquisitions']
ret.append("================================================================")
ret.append(" Running recalculation ")
ret.append("================================================================")
all_fit_data = []
for acq_id, acq_name in enumerate(sorted(acqs.keys())):
# Skip the default symlink
if acq_name == 'default':
continue
acq = acqs[acq_name]
if acq.attrs['acq_type'] == 'COH':
ret.append("Calculating transform for COH acquisition {}".format(acq_name))
# primary is CCS fitted coil positions
fitted_coils = all_acqs[acq_id]['fitted_coils']
num_fitted_coils = len(fitted_coils)
primary = np.zeros((num_fitted_coils, 3))
for idx, name in enumerate(sorted(fitted_coils)):
primary[idx, :] = fitted_coils[name]['location']
# secondary is SCS polhemus positions
pol_coils = all_acqs[acq_id]['pol_coils']
num_pol_coils = len(pol_coils)
secondary = np.zeros((num_pol_coils, 3))
for idx, name in enumerate(sorted(pol_coils)):
secondary[idx, :] = pol_coils[name].mean(axis=0)
if num_fitted_coils != num_pol_coils:
raise Exception("Mismatch between fitted and polhemus coils")
# Extract just the coils to use
primary = primary[coils_to_use_idx, :]
secondary = secondary[coils_to_use_idx, :]
rbt = RigidBodyTransform()
transform, reflection, errors = rbt.calc_transform(primary, secondary)
# Store transform for later use
all_fit_data.append({'transform': transform})
# Update transform so that we can generate a report later
all_acqs[acq_id]['ccs_to_scs_transform'] = np.array(transform)
# Generate the after report
ret.append("================================================================")
ret.append(" After recalculation ")
ret.append("================================================================")
ret.append('')
ret.append(generate_transform_report(all_acqs))
return all_fit_data, coils_to_use, '\n'.join(ret)
def recalculate_coils(filename, calibration_filename, coh_filename,
coils_to_use_inp, debug=False):
"""
Recalculate coil positions and ccs_to_scs_transform for a dataset
......@@ -439,9 +590,9 @@ def recalculate_transforms(filename, calibration_filename, coh_filename,
return all_fit_data, coils_to_use, '\n'.join(ret)
def rewrite_transforms(filename, fit_data, coils_to_use=None):
def rewrite_coils_and_transforms(filename, fit_data, coils_to_use=None):
"""
Rewrite the transform information in an HDF5 file
Rewrite the coil and transform information in an HDF5 file
Parameters
----------
......@@ -449,7 +600,7 @@ def rewrite_transforms(filename, fit_data, coils_to_use=None):
Filename to overwrite transforms in
fit_data
Coil fit data as returned from recalculate_transform
Coil fit data as returned from recalculate_coils
coils_to_use
List of coils which were used or None
......@@ -489,6 +640,53 @@ def rewrite_transforms(filename, fit_data, coils_to_use=None):
store_coh_transforms(fit_data, f)
def rewrite_transforms(filename, fit_data, coils_to_use=None):
"""
Rewrite the transform information in an HDF5 file
Parameters
----------
filename
Filename to overwrite transforms in
fit_data
Transform data as returned from recalculate_transform
coils_to_use
List of coils which were used or None
"""
with h5py.File(filename, 'a') as f:
acqs = f['acquisitions']
# Remove existing datasets and groups
for acq_id, acq_name in enumerate(sorted(acqs.keys())):
# Skip the default symlink
if acq_name == 'default':
continue
acq = acqs[acq_name]
# Remove ccs_to_scs_transform (if any)
if 'ccs_to_scs_transform' in acq:
del acq['ccs_to_scs_transform']
# Get rid of our custom field which tells us which coils
# were used for the transform
if 'transform_coils_used' in acq:
del acq['transform_coils_used']
# We have to put this in manually as store_coh_transforms
# doesn't know about it
if coils_to_use is not None:
acq.create_dataset('transform_coils_used',
shape=(len(coils_to_use), ),
dtype=h5py.special_dtype(vlen=str))
acq['transform_coils_used'][...] = coils_to_use
# Let store_coh_transforms add the datasets back
update_transforms(fit_data, f)
def copy_transforms(src_filename, dest_filename, force=False):
"""
Copy transforms from one file to another
......
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