contrib/sixs/hkl3d.c

   1 /* This file is part of the hkl library.
   2  *
   3  * The hkl library is free software: you can redistribute it and/or modify
   4  * it under the terms of the GNU General Public License as published by
   5  * the Free Software Foundation, either version 3 of the License, or
   6  * (at your option) any later version.
   7  *
   8  * The hkl library is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11  * GNU General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU General Public License
  14  * along with the hkl library.  If not, see <http://www.gnu.org/licenses/>.
  15  *
  16  * Copyright (C) 2003-2016 Synchrotron SOLEIL
  17  *                         L'Orme des Merisiers Saint-Aubin
  18  *                         BP 48 91192 GIF-sur-YVETTE CEDEX
  19  *
  20  * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
  21  */
  22 #include <stdio.h>
  23
  24 //#undef H5_USE_16_API
  25
  26 #include "hkl.h"
  27 #include <hdf5.h>
  28
  29 #define ROOT "/nfs/ruche-sixs/sixs-soleil/com-sixs/2015/Shutdown4-5/XpadAu111/"
  30 #define FILENAME "align_FLY2_omega_00045.nxs"
  31
  32 #define DATASET_IMAGE "com_113934/scan_data/xpad_image"
  33 #define DATASET_MU "com_113934/scan_data/UHV_MU"
  34 #define DATASET_OMEGA "com_113934/scan_data/UHV_OMEGA"
  35 #define DATASET_DELTA "com_113934/scan_data/UHV_DELTA"
  36 #define DATASET_GAMMA "com_113934/scan_data/UHV_GAMMA"
  37 #define DATASET_UB "com_113934/SIXS/I14-C-CX2__EX__DIFF-UHV__#1/UB"
  38 #define DATASET_WAVELENGTH "com_113934/SIXS/Monochromator/wavelength"
  39 #define DATASET_DIFFRACTOMETER_TYPE "com_113934/SIXS/I14-C-CX2__EX__DIFF-UHV__#1/type"
  40
  41 /* Display all the node informations */
  42
  43 static herr_t attribute_info(hid_t location_id, const char *attr_name, const H5A_info_t *ainfo, void *op_data)
  44 {
  45         printf("    Attribute: %d %s\n", location_id, attr_name);
  46
  47         return 0;
  48 }
  49
  50 static herr_t file_info(hid_t loc_id, const char *name, const H5L_info_t *info, void *opdata)
  51 {
  52         H5O_info_t statbuf;
  53         hsize_t n = 0;
  54
  55         /*
  56          * Get type of the object and display its name and type.
  57          * The name of the object is passed to this function by
  58          * the Library. Some magic :-)
  59          */
  60         H5Oget_info_by_name(loc_id, name, &statbuf, H5P_DEFAULT);
  61         switch (statbuf.type) {
  62         case H5O_TYPE_UNKNOWN:
  63                 printf(" Object with name %s is an unknown type\n", name);
  64                 break;
  65         case H5O_TYPE_GROUP:
  66                 printf(" Object with name %s is a group\n", name);
  67                 break;
  68         case H5O_TYPE_DATASET:
  69                 printf(" Object with name %s is a dataset\n", name);
  70                 break;
  71         case H5O_TYPE_NAMED_DATATYPE:
  72                 printf(" Object with name %s is a named datatype\n", name);
  73                 break;
  74         default:
  75                 printf(" Unable to identify an object ");
  76         }
  77
  78         H5Aiterate_by_name(loc_id,  name, H5_INDEX_NAME, H5_ITER_NATIVE, &n, attribute_info, NULL, H5P_DEFAULT);
  79
  80         return 0;
  81 }
  82
  83 typedef struct _HklDataframeSixsUhv HklDataframeSixsUhv;
  84
  85 struct _HklDataframeSixsUhv {
  86         hid_t file;
  87         hid_t mu;
  88         hid_t omega;
  89         hid_t delta;
  90         hid_t gamma;
  91         hid_t ub;
  92         hid_t wavelength;
  93         hid_t dtype;
  94 };
  95
  96
  97 static HklDataframeSixsUhv hkl_h5_open(const char *filename)
  98 {
  99         /* herr_t status; */
 100         hid_t h5file_id;
 101         HklDataframeSixsUhv dataframe = {0};
 102
 103         /* is it an hdf5 file */
 104         if (H5Fis_hdf5(filename) == 0)
 105                 return dataframe;
 106
 107         h5file_id = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT);
 108         if (h5file_id < 0)
 109                 return dataframe;
 110
 111         /* display all the node informations */
 112         /* status = H5Lvisit(h5file_id, H5_INDEX_CRT_ORDER, H5_ITER_NATIVE, file_info, NULL); */
 113         /* if(status < 0) */
 114         /*      return dataframe; */
 115
 116         dataframe.file = h5file_id;
 117         dataframe.mu = H5Dopen (h5file_id, DATASET_MU, H5P_DEFAULT);
 118         dataframe.omega = H5Dopen (h5file_id, DATASET_OMEGA, H5P_DEFAULT);
 119         dataframe.delta = H5Dopen (h5file_id, DATASET_DELTA, H5P_DEFAULT);
 120         dataframe.gamma = H5Dopen (h5file_id, DATASET_GAMMA, H5P_DEFAULT);
 121
 122         dataframe.ub = H5Dopen (h5file_id, DATASET_UB, H5P_DEFAULT);
 123         dataframe.wavelength = H5Dopen (h5file_id, DATASET_WAVELENGTH, H5P_DEFAULT);
 124         dataframe.dtype = H5Dopen (h5file_id, DATASET_DIFFRACTOMETER_TYPE, H5P_DEFAULT);
 125
 126         return dataframe;
 127 }
 128
 129 static hssize_t check_ndims(hid_t dataset_id, int expected, hid_t *space_id)
 130 {
 131         *space_id = H5Dget_space (dataset_id);
 132         return H5Sget_simple_extent_ndims (*space_id) == expected;
 133 }
 134
 135 static int hkl_h5_len(const HklDataframeSixsUhv *dataframe)
 136 {
 137         hid_t space_id = H5Dget_space (dataframe->mu);
 138         return H5Sget_simple_extent_npoints(space_id);
 139 }
 140
 141
 142 static int hkl_h5_is_valid(const HklDataframeSixsUhv *dataframe)
 143 {
 144         int res = TRUE;
 145         hid_t space_id;
 146         hssize_t n;
 147
 148         /* check the dimensionnality of the axes */
 149         res &= check_ndims(dataframe->mu, 1, &space_id);
 150         n = H5Sget_simple_extent_npoints(space_id);
 151
 152         /* check the dimensionnality of all the axes */
 153         res &= check_ndims(dataframe->omega, 1, &space_id);
 154         res &= n == H5Sget_simple_extent_npoints(space_id);
 155
 156         res &= check_ndims(dataframe->delta, 1, &space_id);
 157         res &= n == H5Sget_simple_extent_npoints(space_id);
 158
 159         res &= check_ndims(dataframe->gamma, 1, &space_id);
 160         res &= n == H5Sget_simple_extent_npoints(space_id);
 161
 162         return res;
 163 }
 164
 165 static herr_t hkl_h5_close(const HklDataframeSixsUhv *dataframe)
 166 {
 167         H5Dclose(dataframe->dtype);
 168         H5Dclose(dataframe->wavelength);
 169         H5Dclose(dataframe->ub);
 170         H5Dclose(dataframe->gamma);
 171         H5Dclose(dataframe->delta);
 172         H5Dclose(dataframe->omega);
 173         H5Dclose(dataframe->mu);
 174
 175         return  H5Fclose (dataframe->file);
 176 }
 177
 178 typedef struct _HklDataframe HklDataframe;
 179
 180 struct _HklDataframe {
 181         int i;
 182         int len;
 183         const HklDataframeSixsUhv *_dataframe;
 184 };
 185
 186 static const HklDataframe hkl_dataframe_first(const HklDataframeSixsUhv *dataframe)
 187 {
 188         HklDataframe frame;
 189
 190         frame.i = 0;
 191         frame.len = hkl_h5_len(dataframe);
 192         frame._dataframe = dataframe;
 193
 194         return frame;
 195 }
 196
 197 static int hkl_dataframe_done(const HklDataframe dataframe)
 198 {
 199         return dataframe.i < dataframe.len;
 200 }
 201
 202 static const HklDataframe hkl_dataframe_next(const HklDataframe dataframe)
 203 {
 204         HklDataframe frame = dataframe;
 205
 206         frame.i++;
 207
 208         return frame;
 209 }
 210
 211
 212 static herr_t get_position(hid_t dataset_id, int idx, double *position)
 213 {
 214         hid_t space_id;
 215         hid_t mem_type_id;
 216         hid_t mem_space_id;
 217         hsize_t count[1];
 218         hsize_t offset[1];
 219
 220         mem_type_id = H5Dget_type(dataset_id);
 221
 222         /* Get the dataspace handle */
 223         if ( (space_id = H5Dget_space( dataset_id )) < 0 )
 224                 goto out;
 225
 226         /* Define a hyperslab in the dataset of the size of the records */
 227         offset[0] = idx;
 228         count[0]  = 1;
 229         if ( H5Sselect_hyperslab(space_id, H5S_SELECT_SET, offset, NULL, count, NULL) < 0 )
 230                 goto out;
 231
 232         /* Create a memory dataspace handle */
 233         if ( (mem_space_id = H5Screate_simple( 1, count, NULL )) < 0 )
 234                 goto out;
 235
 236         if ( H5Dread(dataset_id, mem_type_id, mem_space_id, space_id, H5P_DEFAULT, position ) < 0 )
 237                 goto out;
 238
 239         /* Terminate access to the memory dataspace */
 240         if ( H5Sclose( mem_space_id ) < 0 )
 241                 goto out;
 242
 243         /* Terminate access to the dataspace */
 244         if ( H5Sclose( space_id ) < 0 )
 245                 goto out;
 246
 247         return 0;
 248 out:
 249         H5Tclose(mem_type_id);
 250         return -1;
 251 }
 252
 253 static herr_t hkl_dataframe_geometry_get(const HklDataframe dataframe, HklGeometry **geometry)
 254 {
 255         herr_t status = 0;
 256         hid_t datatype;
 257         double wavelength;
 258         double axes[4];
 259
 260         /* create the HklGeometry */
 261         if((*geometry) == NULL){
 262                 char *name;
 263                 size_t n;
 264                 HklFactory *factory;
 265
 266                 /* read the diffractometer type from the hdf5 file */
 267                 datatype = H5Dget_type(dataframe._dataframe->dtype);
 268                 n = H5Tget_size(datatype);
 269                 name = malloc(n+1);
 270                 status = H5Dread(dataframe._dataframe->dtype,
 271                                  datatype,
 272                                  H5S_ALL, H5S_ALL,
 273                                  H5P_DEFAULT, name);
 274                 if(status >= 0){
 275                         /* remove the last "\n" char */
 276                         name[n-1] = 0;
 277
 278                         factory = hkl_factory_get_by_name(name, NULL);
 279                         *geometry = hkl_factory_create_new_geometry(factory);
 280                 }
 281                 free(name);
 282                 H5Tclose(datatype);
 283         }
 284
 285         /* read the wavelength double */
 286         /* TODO check the right size */
 287         /* TODO how to obtain the unit of the  wavelength */
 288         datatype = H5Dget_type(dataframe._dataframe->wavelength);
 289         status = H5Dread(dataframe._dataframe->wavelength,
 290                          datatype,
 291                          H5S_ALL, H5S_ALL,
 292                          H5P_DEFAULT, &wavelength);
 293         if(status >= 0)
 294                 hkl_geometry_wavelength_set(*geometry, wavelength, HKL_UNIT_USER, NULL);
 295         H5Tclose(datatype);
 296
 297         /* read the axis positions of the ith dataframe */
 298         /* check how to decide about the dataset connection and the hkl axes connection */
 299         /* TODO check the right size */
 300         /* TODO how to obtain the unit of the axes position */
 301         if (get_position(dataframe._dataframe->mu,
 302                          dataframe.i, &axes[0]) < 0)
 303                 goto out;
 304         if (get_position(dataframe._dataframe->omega,
 305                          dataframe.i, &axes[1]) < 0)
 306                 goto out;
 307         if (get_position(dataframe._dataframe->gamma,
 308                          dataframe.i, &axes[2]) < 0)
 309                 goto out;
 310         if (get_position(dataframe._dataframe->delta,
 311                          dataframe.i, &axes[3]) < 0)
 312                 goto out;
 313
 314         hkl_geometry_axis_values_set(*geometry, axes, 4, HKL_UNIT_USER, NULL);
 315         /* hkl_geometry_fprintf(stdout, *geometry); */
 316         /* fprintf(stdout, "\n"); */
 317
 318         return 0;
 319 out:
 320         return -1;
 321
 322 }
 323
 324 int main (int argc, char ** argv)
 325 {
 326         const char *filename =  ROOT FILENAME;
 327         HklDataframeSixsUhv dataframe_h5 = hkl_h5_open(filename);
 328         int res = hkl_h5_is_valid(&dataframe_h5);
 329         /* fprintf(stdout, "h5file is valid : %d\n", res); */
 330         HklGeometry *geometry = NULL;
 331
 332         for(HklDataframe dataframe =  hkl_dataframe_first(&dataframe_h5);
 333             hkl_dataframe_done(dataframe);
 334             dataframe = hkl_dataframe_next(dataframe))
 335         {
 336                 hkl_dataframe_geometry_get(dataframe, &geometry);
 337                 /* fprintf(stdout, " %d", dataframe.i); */
 338         }
 339
 340         hkl_geometry_free(geometry);
 341         hkl_h5_close(&dataframe_h5);
 342
 343         return res;
 344 }