maxind format

[bizartor]

I'm finding it difficult to understand how to specify specific reciprocal planes with just maxind. Say I wished to reconstruct the plane h,k,2±0.1, how would I go about specifying that in meerkat?

[k-eks]

From my point of view, meerkat is a powerful tool but not too user friendly at the moment. Since I'm going to use it a lot, I want to change that a little along the road. For now, this code snippet may help:
f = h5py.File('path/to/reconstructed/file.h5')
rebinned_data = f['rebinned_data']
number_of_pixels_rebinned = f['number_of_pixels_rebinned']
l=450 #see comment below
data = rebinned_data[:,:,l]/number_of_pixels_rebinned[:,:,l]
#now you have to visualise the range data[:,:,l]
in this snippet, l does not directly correspond to Miller indices but rather to the pixel layer. You have to find out which pixel "layer" corresponds to your 2 Miller index.
For example,I made a reconstruction with 801 pixels in all direction from +-16 hkl. My hk2 would correspond to l=450 (801/32=25 --> on multiples of 25 I have integer Miller indices, l=0 --> hk-16 the origin of the reconstructed data)
Hope this helped

[aglie]

Thanks, @k-eks your are absolutely correct.

Now meerkat is focusing on reconstructing all of the reciprocal space as a three dimensional volume. And even if you are only interested in some specific layer the only way to reconstruct everything and then extract that particular layer. So in case you want layers of type hkL+-0.1 for L=-5,-4,-3...5 you first reconstruct everything:

reconstruct_data(
        maxind=[10,10,5], #some numbers for h and k, Lmax = 5 
        number_of_pixels=[1001, 1001, 51] #Nl = 51 means that the step size along z is
                                                                 #step=(Nl-1)/(2*Lmax)=0.2  which corresponds to +-0.1
       )

After the reconstruction is done the layers you want can be found here:

import h5py

l=2
Lmax=5
step=0.2
li =(l+Lmax)*step
f = h5py.File('reconstruction.h5')

layer_hk2 =  f["rebinned_data"][:,:,li]/f["number_of_pixels_rebinned"][:,:,li]

[k-eks]

Oh, and by the way: the calculation
f["rebinned_data"][:,:,li]/f["number_of_pixels_rebinned"][:,:,li]
will raise some NAN "exceptions" (divisions through 0 happen), but these NANs are intentional.
That's what fooled my for some time :)

[aglie]

Well, @k-eks, @bizartor you are the very first users, I can change anything just for the program to suit your needs, so the suggestions are welcome!

[bizartor]

Thank you both for your suggestions and clarifications. I will try them out once I'm back in the lab in a week.

[bizartor]

Thanks to disabling the image mask, meerkat now works and outputs the reconstructed HDF5 file!
However, having got this far I still have some questions remaining:

1. It appears that layer_hk2 = f["rebinned_data"][:,:,li]/f["number_of_pixels_rebinned"][:,:,li] extracts a 2D matrix from a 3D matrix. How do the dimensions of the 3D matrix correspond to the axes a* b* c* and a'* b'* c'*, in crystallographic and orthonormal bases respectively?
2. Why is the division by f["number_of_pixels_rebinned"][:,:,li] needed?
3. Is there a way to save the extracted pixel layer as a separate HDF5 file?
4. What software can I use to view the 3D or 2D HDF5 files once they are generated? I tried adxv, but that didn't recognise the file.
5. Is there a way to sum neighbouring pixel layers? Ideally I would like to view reciprocal planes with fractional indices and arbitrary pixel step size, however, being restricted to integer indices for the full 3D volume bounds makes that unfeasible. Seems like the next best thing would be to reconstruct the reciprocal space with high resolution and then sum neighbouring layers around the required fractional index.

Also, it should probably be noted that include h5py is needed for the above code to work, and I think there might be a mistake in the definition of li and it should be:

li = (l+Lmax)/step + 1

Thank you for your continuing assistance.

[k-eks]

I think, I can answer some of your questions, for the others you have to wait for aglie.

3. This is actually straight forward:

#I assume, you want to save your layer_hk2
outputfile = h5py.File("out_file_name.h5", 'w')
outputfile['data'] = layer_hk2
outputfile.close()

4. As far as I know, there is no single software which does an overall good job. hdfview (https://www.hdfgroup.org/products/java/hdfview/) is quite alright, if you want to look at the "raw" numbers and is in my opinion better suited to have a look on the internal structure of the meerkat data than the python command line. It also does come with a visualisation feature which I do NOT recommend to use since it is very poorly implemented and it will create more confusion than solutions.
   Vor visualisation, there is no suitable program that I am aware of. My personal choice is the python module matplotlib (http://matplotlib.org/index.html), which also integrates perfectly into the ipython notebook (kind of a web browser based command line, if you heard about it). For the command line version:

from matplotlib import pyplot
from matplotlib import colors
<...>
layer_hk2 = f["rebinned_data"][:,:,li]/f["number_of_pixels_rebinned"][:,:,li]
pyplot.imshow(layer_hk2,clim=[0,20],cmap=cm.Greys) #clim sets the contrast
pyplot.show()

I hope, this helps for the moment.

[aglie]

Hi, @bizartor, thank you for the questions, thanks @k-eks for the answers you provided.

1. It appears that layer_hk2 = f["rebinned_data"][:,:,li]/f["number_of_pixels_rebinned"][:,:,li] extracts a 2D matrix from a 3D matrix. How do the dimensions of the 3D matrix correspond to the axes a* b* c* and a'* b'* c'*, in crystallographic and orthonormal bases respectively?

Indeed, the script does extract the 2D sections from 3D dataset. The indices of the 3D dataset go along a_, b_ and c_. If the orthonormal coordinates are selected, then the indices will go aling a_', b_', c_'.

Why is the division by f["number_of_pixels_rebinned"][:,:,li] needed?

The process of reconstruction is done like this: the intensity of each detector pixel gets summed up in the appropriate bin, another array counts how many pixels were reconstructed. Since the scattering in each bin is the average of (corrected) scattering on each detector pixels, the final dataset is calculated by dividing the values from the first dataset, on the values on the second dataset.

In certain cases it is important to know how many pixels were reconstructed in each bin (for example to avoid using the diffuse scattering on the edge of the measured space, or close to detector gaps), I agree that the two datasets are a confusing output. I will change it in next release.

Is there a way to save the extracted pixel layer as a separate HDF5 file?

@k-eks perfectly answered the question

What software can I use to view the 3D or 2D HDF5 files once they are generated? I tried adxv, but that didn't recognise the file.

If you go for the 3D view of the dataset, I would recommend ParaView. It works with hdf5 datasets, but unfortunately it requires preparing an additional .xdmf file.

As for 2D view, I did not find any good gui software either. I too recommend using ipython notebook (or julia, Matlab, R at your convenience). If you are willing to install ipython notebook (which is simple - on windows through anaconda, on linux through pip), I can provide scripts which allow to view sections along a_b_, a_c_ and b_c_ directions.

Is there a way to sum neighbouring pixel layers? Ideally I would like to view reciprocal planes with fractional indices and arbitrary pixel step size, however, being restricted to integer indices for the full 3D volume bounds makes that unfeasible. Seems like the next best thing would be to reconstruct the reciprocal space with high resolution and then sum neighbouring layers around the required fractional index.

Well, actually Yell was intended for reconstructing all of the experimental dataset and then extracting the necessary data as necessary.

By the way, in your case the code which I provided will also work for fractional indexes too (as long as they are compatible with the step size of the array).

import h5py

l=2.2 # This one also works since it is compatible with the step size
Lmax=5
step=0.2
li =(l+Lmax)*step
f = h5py.File('reconstruction.h5')

layer_hk2 =  f["rebinned_data"][:,:,li]/f["number_of_pixels_rebinned"][:,:,li]

Also, it should probably be noted that include h5py is needed for the above code to work,

Thank you, I changed it in the submission. I also edited the step size to 0.2 to match what I put in reconstruction. Sorry for that typos.

and I think there might be a mistake in the definition of li

Nop, the indexing is 0 based, the following is correct:

li = (l+Lmax)/step

[bizartor]

Thank you both @aglie and @k-eks. Being able to work with the HDF5 data in MATLAB helped a lot. Based on some code from Laurent I wrote the following script to do what I needed. Taking the requested plane intercept index and plane thickness it extracts the pixel layers closest to the requested and sums them. The only variables that should need user modification are in the first section and probably caxis_scaling at the end.

% Function to display selected planes of variable thickness from a 
% reconstructed 3D reciprocal space generated by meerkat
% (https://github.com/aglie/meerkat) in HDF5 format.

clear all;    % Clear all variables from workspace

%% USER-DEFINED PARAMETERS
info = struct('Sample',{'1,10-Dichlorodecane/urea'},'Beamsize_x_y_um',{'80,100'},'Diffractometer',{'ALS 8.2.2'},...
    'Temp_K',{'98'},'Distance_mm',{'200'},'Phi_deg',{'0:0.5:360'},'Time_s',{'0.3'},'Energy_keV',{''},'Wavelength_A',{'0.885601'});
% Information for figure title

fname = 'ISF_B_136_09_01.h5';           % Reconstructed reciprocal space HDF5 file name
save_fig = 1;                           % Save-figure flag

plane = 1;          % Axis that the extracted plane intersects: 1=a*, 2=b*, 3=c*
plane_ind = 0;      % Centre index of plane along intersected axis / reciprocal lattice units
plane_thick = 0.04; % Thickness of extracted plane / reciprocal lattice units (set to zero to extract plane 1 px wide)


%% DETERMINE COORDINATES OF REQUESTED PLANE
file_info = hdf5info(fname);    % Extract information about HDF5 file

HKL{1} = -h5read(fname,'/maxind',1,1):h5read(fname,'/step_size',1,1):h5read(fname,'/maxind',1,1); % Array of h values
HKL{2} = -h5read(fname,'/maxind',2,1):h5read(fname,'/step_size',2,1):h5read(fname,'/maxind',2,1); % Array of k values
HKL{3} = -h5read(fname,'/maxind',3,1):h5read(fname,'/step_size',3,1):h5read(fname,'/maxind',3,1); % Array of l values

plane_reverse = 4-plane;% Same as 'plane' but reverse axis order (3=a*, 2=b*, 1=c*) to match order in which h5read imports the 3D volume dimensions
plane_rest = [1 2 3];   % Initialise an array of the remaining 2 axes that are parallel to the plane
plane_rest(plane) = []; % Remove the axis which the plane intersects

[~,low_ind] = min(abs(HKL{plane}-(plane_ind-plane_thick/2)));   % Find index of pixel closes to requested plane_ind-plane_thick/2
[~,high_ind] = min(abs(HKL{plane}-(plane_ind+plane_thick/2)));  % Find index of pixel closes to requested plane_ind+plane_thick/2

start_px_ind = ones(1,3);               % Initialise starting px indices to extract from reciprocal space
start_px_ind(plane_reverse) = low_ind;  % Set low-bound index for extracted plane

width_px_ind = file_info.GroupHierarchy.Datasets(6).Dims;               % Initialise number of px to extract from reciprocal space (i.e. max dimensions)
width_px_ind(plane_reverse) = high_ind-low_ind;                         % Thickness of extracted plane / px
if width_px_ind(plane_reverse)<1; width_px_ind(plane_reverse) = 1; end; % If thickness of extracted plane is less than 1 px, set to 1 px


%% EXTRACT RECIPROCAL PLANE
intens = h5read(fname,'/rebinned_data',start_px_ind,width_px_ind);                  % Extract summed pixel intensity in plane from reciprocal space
num_counts = h5read(fname,'/number_of_pixels_rebinned',start_px_ind,width_px_ind);  % Number of intensity contributions to each voxel in the plane
intens_mean = intens./double(num_counts);                                           % Calculate mean reconstructed intensity for each voxel in plane
intens_mean(isnan(intens_mean)) = 0;                                                % Set NaN or Inf values to zero
intens_sum = squeeze(sum(intens_mean,plane_reverse));                               % Sum voxel intensities along reciprocal axis which plane intercepts


%% PLOTTING PARAMETERS
HKL_plot = HKL(plane_rest);             % Select axis arrays that are parallel to the plane
axes_labels = [{'H'} {'K'} {'L'}];      % Initialise cell array of reciprocal indices for figure axes labels
axes_labels = axes_labels(plane_rest);  % Select the axes that are parallel to the plane

plane_low = HKL{plane}(low_ind);                                        % Reciprocal index of low-bound around selected plane
plane_high = HKL{plane}(high_ind);                                      % Reciprocal index of high-bound around selected plane
plane_ind_title = [{'H'} {'K'} {'L'}];                                  % Initialise cell array of reciprocal indices for figure title
plane_ind_title{plane} = sprintf('%.3f to %.3f',plane_low,plane_high);  % Replace intercepted axis index label by intercept value range
plane_ind_title = strjoin(plane_ind_title,', ');                        % Concatenate cell array into one string

intens_sum_mean = mean(mean(intens_sum(intens_sum>0))); % Calculate mean intensity across reconstructed plane (for positive pixels only)

cell_dims = h5read(fname,'/unit_cell');

fname_print = strrep(fname,'_','-');                                                                    % Filenmae sanitized for figure titles
fig_fname = [fname(1:end-3) sprintf('_plane_(%s)',strrep(strrep(plane_ind_title,', ','_'),' ',''))];    % Construct figure filename

% Construct figure title base
fig_title_base = [info.Sample ', reciprocal plane (H,K,L)=(' plane_ind_title '), '...
    sprintf('cell: (a*,b*,c*)=(%.4f,%.4f,%.4f) \xc5, (\\alpha,\\beta,\\gamma)=(%.4f,%.4f,%.4f)\\circ',cell_dims)...
    char(10) 'Temp=' info.Temp_K 'K, d=' info.Distance_mm 'mm, \phi=' info.Phi_deg '^{\circ}, time=' info.Time_s 's, \lambda=' info.Wavelength_A ...
    char(197) ', beamsize (x,y)=(' info.Beamsize_x_y_um ')\mum on ' info.Diffractometer ' (reconstructed data ' fname_print ')'];


%% PLOT RECIPROCAL PLANE
font_size = 8;                          % Figure font size
caxis_scaling = [0.95,2];                % Scaling of colour axis when viewing plane.  In units of, intens_sum_mean, mean pixel value across entire frame.

handles.h1 = figure;
imagesc(HKL_plot{1},HKL_plot{2},log10(intens_sum),log10(caxis_scaling*intens_sum_mean));    % Plot log10 intensities
% colormap(hot);
colormap(flipud(gray));
set(gca,'Box','on','FontSize',font_size);

% grid on
% set(gca,'XTick',HKL_plot{1}(1):1:HKL_plot{1}(end),'YTick',HKL_plot{2}(1):1:HKL_plot{2}(end),'XColor',[0.2 0.2 0.2],'YColor',[0.2 0.2 0.2])
% axis equal
xlabel(axes_labels{1});
ylabel(axes_labels{2});
title(fig_title_base);

if save_fig==1
   savefig(handles.h1,[fig_fname '.fig']); 
end

The result should look something like this:

I'm not sure why, but for some reason MATLAB imports the HDF5 data in reverse order. So if I specify in meerkat number_of_pixels=[101, 201, 301], then the array in MATLAB will be [301,201,101]. I've worked around it in the script, but I was just wondering if this behaviour is expected?

[aglie]

I'm not sure why, but for some reason MATLAB imports the HDF5 data in reverse order. So if I specify in meerkat number_of_pixels=[101, 201, 301], then the array in MATLAB will be [301,201,101]. I've worked around it in the script, but I was just wondering if this behaviour is expected?

It is the feature of Matlab, which uses the column-major order for array memory layout. Since you have enough memory to load all datasets, the easiest way for you to work around it would be to permute datasets immediately after reading them in:

intens = h5read(fname,'/rebinned_data',start_px_ind,width_px_ind);    
intens = permute(intens,[3,2,1])

Also, a couple of observations:

I noticed from your reconstructed image, that you have visible reconstruction aliasing (the visible lines which follow the frames separated by undefined pixels). Probably your raw dataset was measured with 1 degree increment.

In order to fill in the gaps, I recommend to redo the reconstruction with

microsteps=[1,1,5]

One more thing. From the image, it looks like your crystal has non-orthogonal reciprocal space coordinates. The way it is drawn occurs skewed. In order to draw the image the way it should look like you can either:

• reconstruct the crystal with option reconstruct_in_orthonormal_basis set to True, or
• transfrom the image before drawing. This step can be done with image processing toolbox. For hkX layers this can be done the following way:

metric_tensor = h5read(fname,'/metric_tensor')
[~,normalized_metric_tensor]=cov2corr(metric_tensor);
transfrom_matrix=chol(inv(normalized_metric_tensor))';
Tp = maketform('affine',blkdiag(transfrom_matrix(1:2,1:2),1));

imagesc(imtransform(intens_sum),Tp))

[bizartor]

I noticed from your reconstructed image, that you have visible reconstruction aliasing (the visible lines which follow the frames separated by undefined pixels). Probably your raw dataset was measured with 1 degree increment.

In order to fill in the gaps, I recommend to redo the reconstruction with

microsteps=[1,1,5]

Thanks! The frame steps were 0.5 degrees. The reconstruction looks much better now:

Could you please tell me what the microsteps=[1,1,5] option does and the significance of each value?

Regarding transforming the image in MATLAB, unfortunately I don't have the image processing toolbox, so the functions maketform and imtransform are not available. I guess this route is closed to me at the moment.

[aglie]

Could you please tell me what the microsteps=[1,1,5] option does and the significance of each value?

During reconstruction with default settings each detector pixel gets reconstructed only once. It is like assuming that the pixel has zero size in all three dimensions in reciprocal space. In case when angle increment is relatively large (0.5 for example) this leads to some voxels in reconstructed space not filled by any intensity.

The easy workaround which meerkat uses is to reconstruct pixel more than once. In such case each detector pixel is repeated NxMxL times along detector_x, detector_y and angle increment direction, where N,M,L are defined as microsteps=[N,M,L].

Or, putting it more simply, microsteps=[1,1,5] asks meerkat to reconstruct each frame 5 times, as if it were 5 frames with 0.1 degree angle increment.

Regarding transforming the image in MATLAB, unfortunately I don't have the image processing toolbox, so the functions maketform and imtransform are not available. I guess this route is closed to me at the moment.

In such a case you can use the following function. It is super slow, but does the job (also draws fancy hexagonal pixels):

function hex_imagesc(image,inv_flag)
if nargin<2
    inv_flag='direct';
end

plot(1,1)

[x,y]=ndgrid(1:size(image,1),1:size(image,2));
r=[x(:) y(:)]';

s=1/sqrt(3);
ang=0;
M=[sqrt(3)/2 0;-0.5 1];
if strcmp(inv_flag,'inverse');
    M=[1 0;1/2 sqrt(3)/2]';
    ang=30;
end

hex_x=cosd(linspace(0,360,7)+ang)'*s;
hex_y=sind(linspace(0,360,7)+ang)'*s;

r=M*r;

x=repmat(hex_x,1,size(r,2))+repmat(r(1,:),7,1);
y=repmat(hex_y,1,size(r,2))+repmat(r(2,:),7,1);

patch(x,y,image(:)','EdgeAlpha',0)

to use it try:

hex_imagesc(data)

or

hex_imagesc(data,'inverse')

not sure which one is for direct and which is for reciprocal space.

[bizartor]

Thanks for the hexagonal-pixel code @aglie. It works pretty well.

Regarding the microsteps option, I have compared a selection of reflections from the same plane, with and without microsteps=[1,1,5], and have noticed some slight differences between their positions. The data cursor in the images below is on the same pixel for comparison. When microsteps=[1,1,5] the reflection seems to shift by a few pixels, presumably in the direction of rotation of the detector plane.

Given a frame angle increment of 0.5° and microsteps=[1,1,5], I would have expected each frame to be reconstructed in its original position, and then 4 more times following increments of 0.1°. So the 0° frame would be reconstructed at 0.0, 0.1, 0.2, 0.3 and 0.4°. What the images below seem to imply is that the reconstruction occurs at 0.1, 0.2, 0.3, 0.4 and 0.5° (or something similar), such that there is no-longer intensity in the original frame position of 0.0° from that frame. Is this the case?

microsteps=[1,1,1]:

microsteps=[1,1,5]:

microsteps=[1,1,1]:

microsteps=[1,1,5]:

[aglie]

Wow, thank you, @bizartor. This definitely looks like a shiny bug! I will take a look at it.

[aglie]

Dear @bizartor could you please check the new version of meerkat? I fixed the bug and now the two reconstructions should be internally consistent. Also, the Bragg peak position should be closer to the integer indices, though it is not yet set in stone since I found one possible inconsistency in XDS.

[bizartor]

Thanks @aglie. Unfortunately I am chasing deadlines at the moment, but I shall check the latest version of meerkat as soon as I can, possibly early next week. Is it available via pip install meerkat?

[bizartor]

I just realised I never got back to you about this. I ran out of time to work on that issue, so didn't install the new version of meerkat to see if the bug was fixed. There's not telling when I will next be able to do it, so I hope you have a way of checking. Sorry. :(

[aglie]

Oh, no worries, you see that my ping on Meerkat is also close to a year.

Yes, the new version should be available through pip and should just work. I also found some inconsistencies in XDS which I have not solved yet, but if you will run into a similar problem, just write here I will accelerate the solution.