From 377bb97d488ca040324280d3df25533fb87c93b5 Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Fri, 21 Mar 2014 18:00:53 +0100
Subject: [PATCH] Blobs: initial experimental code for generating blobs in
grain_%04d.mat files
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
2_difspot/gtReadBlobFromHDF5.m | 15 ++++++++
2_difspot/gtSegDefaultParameters.m | 1 +
2_difspot/gtSegmentationDoubleThreshold.m | 21 ++++++++---
2_difspot/gtWriteBlobToHDF5.m | 23 ++++++++++++
5_reconstruction/gtForwardSimulate_v2.m | 44 +++++++++++++++++++++++
5 files changed, 100 insertions(+), 4 deletions(-)
create mode 100644 2_difspot/gtReadBlobFromHDF5.m
create mode 100644 2_difspot/gtWriteBlobToHDF5.m
diff --git a/2_difspot/gtReadBlobFromHDF5.m b/2_difspot/gtReadBlobFromHDF5.m
new file mode 100644
index 00000000..961af0a8
--- /dev/null
+++ b/2_difspot/gtReadBlobFromHDF5.m
@@ -0,0 +1,15 @@
+function [blob_vol, blob_bb] = gtReadBlobFromHDF5(blob_dir, blob_id)
+ sub_dir = fullfile(blob_dir, sprintf('%06d', blob_id - mod(blob_id, 1e4)));
+ fnameout = fullfile(sub_dir, sprintf('difblob_%06d.hdf5', blob_id) );
+ if (isempty(dir(fnameout)))
+ error('gtReadBlobFromHDF5:file_not_existing', ...
+ 'The file: %s does not exist!', fnameout);
+ end
+
+ blob_vol = h5read(fnameout, '/rawblob');
+ if (nargout == 2)
+ blob_bb = h5readatt(fnameout, '/rawblob', 'bb');
+ blob_bb = [blob_bb(2:3:end), blob_bb(1:3:end), blob_bb(3:3:end)];
+ blob_bb = double([blob_bb(1, :), blob_bb(2, :)]);
+ end
+end
\ No newline at end of file
diff --git a/2_difspot/gtSegDefaultParameters.m b/2_difspot/gtSegDefaultParameters.m
index 39a1f8f1..133d7759 100644
--- a/2_difspot/gtSegDefaultParameters.m
+++ b/2_difspot/gtSegDefaultParameters.m
@@ -22,6 +22,7 @@ par_seg.debug = false; % display messages
par_seg.writeblobs = true; % write difblobs to the table
par_seg.writespots = true; % write difspot metadata to the table
par_seg.writeedfs = true; % save difspots as edf files
+par_seg.writehdf5 = true; % write difblobs to the table
par_seg.wrapping = true; % for 360 degree data, wrap from
% the last image back to the first image
par_seg.segmentation_stack_size = 1000; % how many images in memory? (1000 images approx 32Gb)
diff --git a/2_difspot/gtSegmentationDoubleThreshold.m b/2_difspot/gtSegmentationDoubleThreshold.m
index 44a3284d..74499134 100644
--- a/2_difspot/gtSegmentationDoubleThreshold.m
+++ b/2_difspot/gtSegmentationDoubleThreshold.m
@@ -118,6 +118,11 @@ if isfield(parameters.seg, 'writeedfs')
else
stack.writeedfs = true;
end
+if isfield(parameters.seg, 'writehdf5')
+ stack.writehdf5 = parameters.seg.writehdf5;
+else
+ stack.writehdf5 = true;
+end
% do we want to wrap around? Not for live, but for post mortem we could
if isfield(parameters.seg, 'wrapping')
stack.wrapping = parameters.seg.wrapping;
@@ -131,9 +136,9 @@ else
end
% you must write either blobs or edfs
-if (stack.writeblobs) || (stack.writespots && stack.writeedfs)
- fprintf('Output configuration:\n (writeblobs=%d, writespots=%d, writeedfs=%d)\n',...
- stack.writeblobs,stack.writespots,stack.writeedfs)
+if (stack.writeblobs) || (stack.writehdf5) || (stack.writespots && stack.writeedfs)
+ fprintf('Output configuration:\n (writeblobs=%d, writehdf5=%d, writespots=%d, writeedfs=%d)\n',...
+ stack.writeblobs, stack.writehdf5, stack.writespots, stack.writeedfs)
else
fprintf(['Output configuration:\n (writeblobs=%d, writespots=%d, writeedfs=%d) '...
'\n is not valid!\nYou must save images in some way\n Quitting... \n'], ...
@@ -149,6 +154,7 @@ stack.parameters = parameters;
stack.seg = parameters.seg;
stack.fulldir = fullfile(parameters.acq.dir, '1_preprocessing', 'full');
stack.difspotdir = fullfile(parameters.acq.dir, '2_difspot');
+stack.difblobdir = fullfile(parameters.acq.dir, '2_difblob');
% stack range
if strcmpi(parameters.acq.type, '360degree')
@@ -463,7 +469,14 @@ function sfWriteBlob(blob, stack)
'zIndex', zz, ...
'graylevel', vals);
end
-
+
+ if (stack.writehdf5)
+ blob_bb = blob.bb;
+ % put z into reference of the dataset
+ blob_bb(3) = blob_bb(3) + stack.first - 1;
+ gtWriteBlobToHDF5(stack.difblobdir, blob.difblobID, blob.rawgreyvol, blob_bb);
+ end
+
if (stack.writespots)
gtSegmentationBlob2SpotTable(blob, stack)
end
diff --git a/2_difspot/gtWriteBlobToHDF5.m b/2_difspot/gtWriteBlobToHDF5.m
new file mode 100644
index 00000000..cb0a8c42
--- /dev/null
+++ b/2_difspot/gtWriteBlobToHDF5.m
@@ -0,0 +1,23 @@
+function gtWriteBlobToHDF5(blob_dir, blob_id, blob_vol, blob_bb)
+ sub_dir = fullfile(blob_dir, sprintf('%06d', blob_id - mod(blob_id, 1e4)));
+ if (~exist(sub_dir, 'dir'))
+ [status, message, messageid] = mkdir(sub_dir);
+ if (~status)
+ gtError(messageid, message)
+ elseif (~isempty(messageid))
+ warning(messageid, message)
+ end
+ end
+ fnameout = fullfile(sub_dir, sprintf('difblob_%06d.hdf5', blob_id) );
+
+ if (exist(fnameout, 'file'))
+ delete(fnameout);
+ end
+ vol_size = [size(blob_vol, 1) size(blob_vol, 2) size(blob_vol, 3)];
+ chunk_size = min([16 16 4], vol_size);
+ h5create(fnameout, '/rawblob', vol_size, 'Datatype', 'single', 'Deflate', 1, 'ChunkSize', chunk_size);
+% h5create(fnameout, '/rawblob', size(blob_vol), 'Datatype', 'single');
+
+ h5write(fnameout, '/rawblob', single(blob_vol))
+ h5writeatt(fnameout, '/rawblob', 'bb', uint32(blob_bb))
+end
\ No newline at end of file
diff --git a/5_reconstruction/gtForwardSimulate_v2.m b/5_reconstruction/gtForwardSimulate_v2.m
index 13b6ca9d..82cdfd9c 100644
--- a/5_reconstruction/gtForwardSimulate_v2.m
+++ b/5_reconstruction/gtForwardSimulate_v2.m
@@ -184,6 +184,11 @@ for n = first : last
numspots = numel(onDet); % number of exspected diffraction spots
flag = ones(numspots, 1); % see above
% stack of zero padded diffraction spots (projections used for 3D reconstruction)
+ bl = struct( ...
+ 'intm', cell(numspots, 1), ...
+ 'bbuim', cell(numspots, 1), ...
+ 'bbvim', cell(numspots, 1), ...
+ 'bbwim', cell(numspots, 1));
difstack = zeros(spotsCommProps.stackHSize, numspots, spotsCommProps.stackVSize);
bb = zeros(numspots, 4); % difspot BoundingBox from database
spotid = zeros(numspots, 1); % Mysql difspottable difspotID
@@ -299,6 +304,30 @@ for n = first : last
difstack(:, ii, :) = spot';
end
+ %%% We now build the diffraction blob stack
+ for ii = reshape(matched, [1 numel(matched)])
+ [blob_vol, blob_bb] = buildBlob(spotid(ii), parameters, spotsCommProps);
+% blob_vol = permute(blob_vol, [2 1 3]);
+
+% if (fsim.assemble_figure)
+% full = gtPlaceSubImage2(sum(blob_vol, 3), full, bb(ii, 1), bb(ii, 2), 'summed');
+% end
+
+ % Essentially zero pads the blob to make it fit into ASTRA's
+ % diffraction stack
+ shifts_blob = getBlobStackShifts(spotsCommProps, blob_bb);
+ shifts = [shifts_blob.h + 1, shifts_blob.v + 1, 1];
+
+ bl(ii).intm = gtPlaceSubVolume( ...
+ zeros(spotsCommProps.stackVSize, spotsCommProps.stackHSize, blob_bb(6)), ...
+ blob_vol, shifts);
+ bl(ii).intm = permute(bl(ii).intm, [2 1 3]);
+
+ bl(ii).bbuim = [blob_bb(1), (blob_bb(1)+blob_bb(4)-1)]; % Might be exchanged with v
+ bl(ii).bbvim = [blob_bb(2), (blob_bb(2)+blob_bb(5)-1)]; % Might be exchanged with u
+ bl(ii).bbwim = [blob_bb(3), (blob_bb(3)+blob_bb(6)-1)];
+ end
+
%%% We now Build ASTRA's geometry
for ii = reshape(matched, [1 numel(matched)])
% If available, we want to use the angles as determined from
@@ -400,6 +429,7 @@ for n = first : last
found_reflections = length(find(flag > 1));
out.completeness = found_reflections / numspots;
+ out.bl = bl(included);
out.proj.stack = difstack(:, included, :);
out.proj.geom = proj_geom(included, :);
out.proj.geom_full = proj_geom_full(included, :);
@@ -665,6 +695,15 @@ function spot = buildSpot(candidateID, parameters, info, spotsCommProps, bb)
spot = spot / sum(spot(:)) * volume;
end
+function [blob, bb] = buildBlob(candidateID, parameters, spotsCommProps)
+ base_dir = fullfile(parameters.acq.dir, '2_difblob');
+ [blob, bb] = gtReadBlobFromHDF5(base_dir, candidateID);
+
+ % here we force the intensity to be the same in all spots -
+ volume = spotsCommProps.Xsize * spotsCommProps.Xsize * spotsCommProps.Ysize;
+ blob = blob / sum(blob(:)) * volume;
+end
+
function candidateIDs = findCandidates(spotsProps, reflPos, segmentedSpots)
% at this stage we want to include spots which may overlap with the spot we are looking for.
% the omega difference of candidate spots must be smaller than thr_max_offset,
@@ -700,3 +739,8 @@ function shifts = getDiffStackShifts(spotsCommProps, bb)
shifts.v = round((spotsCommProps.stackVSize - bb(4)) / 2);
end
+function shifts = getBlobStackShifts(spotsCommProps, bb)
+ shifts.h = round((spotsCommProps.stackHSize - bb(4)) / 2);
+ shifts.v = round((spotsCommProps.stackVSize - bb(5)) / 2);
+end
+
--
GitLab