src/Tools/filter/mrc2hdf/src/mrc2hdf.c

   1 /*
   2 # mrc2hdf : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : mrc2hdf
   6 # Attention
   7 #   $Loccker$
   8 #       $State$
   9 #
  10 */
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <math.h>
  15 #include <float.h>
  16 #define GLOBAL_DECLARATION
  17 #include "../inc/config.h"
  18
  19 #define DEBUG
  20 #include "genUtil.h"
  21 #include "Memory.h"
  22 #include "mrcImage.h"
  23 #include "lmrcImageProjection.h"
  24 #include <hdf5.h>
  25
  26 typedef struct lmrc2hdfInfo {
  27         int   numFile;
  28     int   nResolution;
  29     int flagIn;
  30     mrcImage* in;
  31     float inColor[4];
  32     int flagIn2;
  33     mrcImage* in2;
  34     float in2Color[4];
  35 } lmrc2hdfInfo;
  36
  37 typedef enum lmrc2hdfMode {
  38         a=0,
  39         b=1
  40 } lmrc2hdfMode;
  41
  42 extern int lmrc2hdfimaris(hid_t* out, mrcImage* in, lmrc2hdfInfo* linfo, int mode);
  43 extern int lmrc2hdf(hid_t* out, mrcImage* in, lmrc2hdfInfo* linfo, int mode);
  44
  45 int
  46 main(int argc, char* argv[])
  47 {
  48         mrc2hdfInfo info;
  49     lmrc2hdfInfo linfo;
  50
  51     hid_t   out;
  52     herr_t  status;
  53
  54     int i;
  55
  56         init0(&info);
  57     argCheck(&info, argc, argv);
  58     init1(&info);
  59
  60         DEBUGPRINT("Program Start\n");
  61     if(info.flagIn) {
  62         linfo.numFile = 1;
  63         linfo.in = (mrcImage*)memoryAllocate(sizeof(mrcImage), "in main");
  64         mrcFileRead(linfo.in, info.In, "in main", 0);
  65         linfo.flagIn = 1;
  66         linfo.flagIn2 = 0;
  67     } else if(info.flagInList) {
  68         linfo.numFile = info.flagInList;
  69         linfo.flagIn = info.flagInList;
  70         linfo.in = (mrcImage*)memoryAllocate(sizeof(mrcImage)*linfo.numFile, "in main");
  71         for(i=0; i<linfo.numFile; i++) {
  72             mrcFileRead(&(linfo.in[i]), info.InList[i], "in main", 0);
  73         }
  74         linfo.inColor[0] = info.IR;
  75         linfo.inColor[1] = info.IG;
  76         linfo.inColor[2] = info.IB;
  77         linfo.inColor[3] = info.IA;
  78         if(info.flagInList2) {
  79             if(linfo.numFile != info.flagInList2) {
  80                 fprintf(stderr, "Different Number between -I and -I2: %ld %ld\n", info.flagInList, info.flagInList2);
  81                 exit(EXIT_FAILURE);
  82             }
  83             linfo.flagIn2 = info.flagInList2;
  84             linfo.in2 = (mrcImage*)memoryAllocate(sizeof(mrcImage)*linfo.numFile, "in main");
  85             for(i=0; i<linfo.numFile; i++) {
  86                 mrcFileRead(&(linfo.in2[i]), info.InList2[i], "in main", 0);
  87             }
  88             linfo.in2Color[0] = info.I2R;
  89             linfo.in2Color[1] = info.I2G;
  90             linfo.in2Color[2] = info.I2B;
  91             linfo.in2Color[3] = info.I2A;
  92         } else {
  93             linfo.flagIn2 = 0;
  94         }
  95     } else {
  96         fprintf(stderr, "-i or -I is necessary\n");
  97         usage(argv[0]);
  98     }
  99     linfo.nResolution = info.nResolution;
 100     //out = H5Fopen(info.Out, H5F_ACC_TRUNC, H5P_DEFAULT);
 101     out = H5Fcreate(info.Out, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
 102     if(out<0) {
 103         fprintf(stderr, "Not openend:%s\n", info.Out);
 104         exit(EXIT_FAILURE);
 105     }
 106
 107     switch(info.hdf5mode) {
 108         case 0: {
 109             lmrc2hdf(&out, linfo.in, &linfo, info.mode);
 110             break;
 111         }
 112         case 1: {
 113             lmrc2hdfimaris(&out, linfo.in, &linfo, info.mode);
 114             break;
 115         }
 116         default: {
 117             fprintf(stderr, "not supported hdf5mode: %ld\n", info.hdf5mode);
 118             usage(argv[0]);
 119             exit(EXIT_FAILURE);
 120         }
 121     }
 122     status = H5Fclose(out);
 123
 124         exit(EXIT_SUCCESS);
 125 }
 126
 127 void
 128 additionalUsage()
 129 {
 130         fprintf(stderr, "----- Additional Usage -----\n");
 131 }
 132
 133
 134 int
 135 lhdfAttributeString(hid_t gIDChannel, char tag[], char* stmp)
 136 {
 137     hsize_t sdim;
 138     hid_t attrspaceID, attrID;
 139     hid_t stringType;
 140     herr_t status;
 141
 142     sdim = strlen(stmp);
 143     attrspaceID = H5Screate_simple(1, &sdim, NULL);
 144     stringType = H5Tcopy(H5T_C_S1);
 145     status = H5Tset_size(stringType, 1);
 146     status = H5Tset_strpad(stringType, H5T_STR_NULLTERM);
 147     DEBUGPRINT3("space: %lld type: %lld status: %d \n", attrspaceID, stringType, status);
 148     attrID = H5Acreate(gIDChannel, tag, stringType, attrspaceID, H5P_DEFAULT, H5P_DEFAULT);
 149     DEBUGPRINT1("attrID: %lld\n", attrID);
 150     status = H5Awrite(attrID, stringType, stmp);
 151     status = H5Aclose(attrID);
 152     status = H5Sclose(attrspaceID);
 153
 154     return status;
 155 }
 156
 157 int
 158 lhdfAttributeUInt32(hid_t gIDChannel, char tag[], uint32_t* i, hsize_t num)
 159 {
 160     hid_t attrspaceID, attrID;
 161     herr_t status;
 162
 163     attrspaceID = H5Screate_simple(1, &num, NULL);
 164     attrID = H5Acreate(gIDChannel, tag, H5T_STD_U32LE, attrspaceID, H5P_DEFAULT, H5P_DEFAULT);
 165     status = H5Awrite(attrID, H5T_STD_U32LE, i);
 166     status = H5Aclose(attrID);
 167     status = H5Sclose(attrspaceID);
 168
 169     return status;
 170 }
 171
 172 int
 173 lhdfAttributeInt2String(hid_t gIDChannel, char tag[], hsize_t i)
 174 {
 175     char stmp[1024];
 176     hsize_t sdim;
 177     hid_t attrspaceID, attrID;
 178     hid_t stringType;
 179     herr_t status;
 180
 181     sprintf(stmp, "%lld", i);
 182     sdim = strlen(stmp);
 183     attrspaceID = H5Screate_simple(1, &sdim, NULL);
 184     stringType = H5Tcopy(H5T_C_S1);
 185     status = H5Tset_size(stringType, 1);
 186     status = H5Tset_strpad(stringType, H5T_STR_NULLTERM);
 187     DEBUGPRINT3("space: %lld type: %lld status: %d \n", attrspaceID, stringType, status);
 188     attrID = H5Acreate(gIDChannel, tag, stringType, attrspaceID, H5P_DEFAULT, H5P_DEFAULT);
 189     DEBUGPRINT1("attrID: %lld\n", attrID);
 190     status = H5Awrite(attrID, stringType, stmp);
 191     status = H5Aclose(attrID);
 192     status = H5Sclose(attrspaceID);
 193
 194     return status;
 195 }
 196
 197 int
 198 lhdfAttributeReal2String(hid_t gIDChannel, char tag[], double f)
 199 {
 200     char stmp[1024];
 201     hsize_t sdim;
 202     hid_t attrspaceID, attrID;
 203     hid_t stringType;
 204     herr_t status;
 205
 206     sprintf(stmp, "%lf", f);
 207     sdim = strlen(stmp);
 208     attrspaceID = H5Screate_simple(1, &sdim, NULL);
 209     stringType = H5Tcopy(H5T_C_S1);
 210     status = H5Tset_size(stringType, 1);
 211     status = H5Tset_strpad(stringType, H5T_STR_NULLTERM);
 212     DEBUGPRINT3("space: %lld type: %lld status: %d \n", attrspaceID, stringType, status);
 213     attrID = H5Acreate(gIDChannel, tag, stringType, attrspaceID, H5P_DEFAULT, H5P_DEFAULT);
 214     DEBUGPRINT1("attrID: %lld\n", attrID);
 215     status = H5Awrite(attrID, stringType, stmp);
 216     status = H5Aclose(attrID);
 217     status = H5Sclose(attrspaceID);
 218
 219     return status;
 220 }
 221
 222 int
 223 lmrc2hdfimaris(hid_t* out, mrcImage* in, lmrc2hdfInfo* linfo, int mode)
 224 {
 225     // String for tags
 226     static hsize_t sdim=1024;
 227     char tmp[sdim];
 228     char stmp[sdim];
 229     // Histogram
 230     static hsize_t nHist1024=1024;
 231     static hsize_t nHist=256;
 232     hsize_t iLevel;
 233     uint64_t hist1024[nHist1024];
 234     uint64_t hist[nHist];
 235
 236     // image
 237     double data;
 238     float* image;
 239
 240     hsize_t dims[4];
 241     hsize_t chunk[4];
 242
 243     // for hdf5
 244     herr_t status;
 245     hid_t datasetID;
 246     hid_t datasetID2;
 247     hid_t dataspaceID;
 248     hid_t dataspaceID2;
 249     hid_t datasetPropertyID;
 250     hid_t dset;
 251     static uint32_t numOfDataSets = 1;
 252
 253     int nGroup=4;
 254     hid_t groupID[4];
 255     char* groupName[] = {"/DataSet", "/DataSetInfo", "/DataSetTimes", "/Thumbnail"};
 256     int val;
 257     int i, ires, itime, ichannel;
 258     uint64_t i64;
 259     int nres, ntime, nchannel;
 260     hid_t dsetID, dspaceID, gIDRes, gIDTime, gIDChannel, attrID, attrspaceID, stringType;
 261     mrcImageParaTypeReal x, y, z;
 262     mrcImageParaTypeReal xorg, yorg, zorg;
 263     int irange, numrange;
 264     unsigned char* thumbImage;
 265     double sum, min, max, mean;
 266     mrcImage tmpImage;
 267     mrcImage* tmpIn;
 268     lmrcImageProjectionInfo lpro;
 269
 270     DEBUGPRINT("Start lmrc2hdfimaris\n");
 271     // Top Level Attribute
 272     lhdfAttributeString(*out, "DataSetDirectoryName",       "DataSet");
 273     lhdfAttributeString(*out, "DataSetInfoDirectoryName",   "DataSetInfo");
 274     lhdfAttributeString(*out, "ImarisDataSet",              "ImarisDataSet");
 275     lhdfAttributeString(*out, "ImarisVersion",              "5.5.0");
 276     lhdfAttributeUInt32(*out, "NumberOfDataSets",           &numOfDataSets, 1);
 277     lhdfAttributeString(*out, "ThumbnailDirectoryName",     "Thumbnail");
 278
 279     // Top Level Groups
 280     for(i=0; i<nGroup; i++) {
 281         DEBUGPRINT2("Group: %s %d\n", groupName[i], i);
 282         groupID[i] = H5Gcreate2(*out, groupName[i], H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 283         if(groupID[i]<0) {
 284             fprintf(stderr, "Not create Group: %s\n", groupName[i]);
 285             exit(EXIT_FAILURE);
 286         } else {
 287             DEBUGPRINT1("Success: %d\n", i);
 288         }
 289     }
 290
 291     // DataSet/.../Data
 292     dims[0] = in[0].HeaderN.z;
 293     dims[1] = in[0].HeaderN.y;
 294     dims[2] = in[0].HeaderN.x;
 295
 296     /*
 297     chunk[0] = dims[0];
 298     chunk[1] = dims[1];
 299     chunk[2] = dims[2];
 300     */
 301     chunk[0] = 32;
 302     chunk[1] = 64;
 303     chunk[2] = 64;
 304
 305     datasetPropertyID = H5Pcreate(H5P_DATASET_CREATE);
 306     status = H5Pset_chunk(datasetPropertyID, 3, chunk);
 307
 308     val = 3;
 309     status = H5Pset_fill_value(datasetPropertyID, H5T_NATIVE_INT, &val);
 310     status = H5Pset_alloc_time(datasetPropertyID, H5D_ALLOC_TIME_EARLY);
 311     status = H5Pclose(datasetPropertyID);
 312
 313     // DataSet
 314     nres = linfo->nResolution;
 315     ntime = linfo->numFile;
 316     if(0<linfo->flagIn) {
 317         nchannel = 1;
 318     }
 319     if(0<linfo->flagIn2) {
 320         nchannel = 2;
 321         if(in[0].HeaderN.x != linfo->in2[0].HeaderN.x
 322          ||in[0].HeaderN.y != linfo->in2[0].HeaderN.y
 323          ||in[0].HeaderN.z != linfo->in2[0].HeaderN.z) {
 324             fprintf(stderr, "Different Size : In (%d %d %d) and I2 (%d, %d, %d)\n",
 325                 linfo->in[0].HeaderN.x,  linfo->in[0].HeaderN.y,  linfo->in[0].HeaderN.z,
 326                 linfo->in2[0].HeaderN.x, linfo->in2[0].HeaderN.y, linfo->in2[0].HeaderN.z);
 327             exit(EXIT_FAILURE);
 328         }
 329     }
 330
 331     // ResolutionLevel
 332     for(ires=0; ires<nres; ires++) {
 333         sprintf(tmp, "%s/%s %d", groupName[0], "ResolutionLevel", ires);
 334         gIDRes = H5Gcreate2(groupID[0], tmp, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 335         DEBUGPRINT1("%s\n", tmp);
 336     for(itime=0; itime<ntime; itime++) {
 337         sprintf(tmp, "%s/%s %d/%s %d", groupName[0], "ResolutionLevel", ires, "TimePoint", itime);
 338         gIDTime = H5Gcreate2(gIDRes, tmp, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 339         DEBUGPRINT1("%s\n", tmp);
 340     for(ichannel=0; ichannel<nchannel; ichannel++) {
 341         switch(ichannel) {
 342             case 0:
 343                 tmpIn = linfo->in;
 344                 break;
 345             case 1:
 346                 tmpIn = linfo->in2;
 347                 break;
 348             default:
 349                 fprintf(stderr, "Not supported ichannel: %d\n", ichannel);
 350                 exit(EXIT_FAILURE);
 351                 break;
 352         }
 353         sprintf(tmp, "%s/%s %d/%s %d/%s %d", groupName[0], "ResolutionLevel", ires, "TimePoint", itime, "Channel", ichannel);
 354         gIDChannel = H5Gcreate2(gIDTime, tmp, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 355         DEBUGPRINT1("%s\n", tmp);
 356         // Channel Attribute
 357         lhdfAttributeInt2String(gIDChannel, "ImageSizeX", dims[2]);
 358         lhdfAttributeInt2String(gIDChannel, "ImageSizeY", dims[1]);
 359         lhdfAttributeInt2String(gIDChannel, "ImageSizeZ", dims[0]);
 360
 361         // TimePoint
 362         dspaceID = H5Screate_simple(3, dims, NULL);
 363         sprintf(tmp, "%s/%s %d/%s %d/%s %d/%s", groupName[0], "ResolutionLevel", ires, "TimePoint", itime, "Channel", ichannel, "Data");
 364         DEBUGPRINT1("%s\n", tmp);
 365         DEBUGPRINT1("mode: %d\n", tmpIn[0].HeaderMode);
 366         image = (float*)memoryAllocate(sizeof(float)*dims[0]*dims[1]*dims[2], "in lmrc2hdfimars");
 367
 368         for(i=0; i<256; i++) hist[i] = 0;
 369         for(i=0; i<1024; i++) hist1024[i] = 0;
 370         // image
 371         DEBUGPRINT("image\n");
 372         if(tmpIn[itime].HeaderAMax<=tmpIn[itime].HeaderAMin) {
 373             fprintf(stderr, "Image data is flat\n");
 374             for(i=0; i<dims[0]*dims[1]*dims[2]; i++) {
 375                 image[i] = tmpIn[itime].HeaderAMin;
 376                 hist[0]++;
 377                 hist1024[0]++;
 378             }
 379         } else {
 380             irange=1<<ires;
 381             DEBUGPRINT2("res: %d irange %d\n", ires, irange);
 382             for(z=0; z<dims[0]; z++) {
 383                 if(0==(((int)(z/(double)(dims[0]-1)*100))%10)) {
 384                     DEBUGPRINT2("%d/%llu\n", (int)z, dims[0]);
 385                 }
 386             for(y=0; y<dims[1]; y++) {
 387             for(x=0; x<dims[2]; x++) {
 388                 i64 = (double)x + (double)y*(double)dims[2] + (double)z*(double)dims[2]*(double)dims[1];
 389                 if(dims[0]*dims[1]*dims[2]<=i64) {
 390                     DEBUGPRINT1("%f \n", x + (double)y*(double)dims[2] + (double)z*(double)dims[2]*(double)dims[1]);
 391                     DEBUGPRINT4("%lld %f %f %f\n", i64, x, y, z);
 392                     DEBUGPRINT4("%lld %lld %lld %lld\n", dims[0]*dims[1]*dims[2], dims[2], dims[1], dims[0]);
 393                 }
 394                 sum = numrange = 0;
 395                 for(zorg=z*irange;zorg<(z+1)*irange; zorg++) {
 396                 for(yorg=y*irange;yorg<(y+1)*irange; yorg++) {
 397                 for(xorg=x*irange;xorg<(x+1)*irange; xorg++) {
 398                     numrange++;
 399                     mrcPixelDataGet(&(tmpIn[itime]), xorg, yorg, zorg, &data, mrcPixelRePart, mrcPixelHowNearest);
 400                     sum+=data;
 401                 }
 402                 }
 403                 }
 404                 if(0<numrange) {
 405                     data = sum/numrange;
 406 #undef DEBUG2
 407 #ifdef DEBUG2
 408                     if(0<data) {
 409                         DEBUGPRINT4("%llu %f = %f / %d\n", i64, data, sum, numrange);
 410                     }
 411 #endif
 412 #undef DEBUG2
 413                 } else {
 414                     //DEBUGPRINT("numrange is zero: Something wrong\n");
 415                     data = 0;
 416                 }
 417                 image[i64] = (float)data;
 418             } // x
 419             } // y
 420             } // z
 421
 422             DEBUGPRINT("MIN/MAX\n");
 423             min= max = image[0];
 424             for(i64=1; i64<dims[2]*dims[1]*dims[0]; i64++) {
 425                 if(image[i64]<min) min = image[i64];
 426                 if(max<image[i64]) max = image[i64];
 427             }
 428             DEBUGPRINT2("%f/%f\n", min, max);
 429             DEBUGPRINT("HISTGRAM\n");
 430             for(i64=1; i64<dims[2]*dims[1]*dims[0]; i64++) {
 431                 data = image[i64];
 432                 iLevel = MAX(0,MIN(1023,(int)((data-min)/(max-min)*1023+0.5)));
 433                 hist1024[iLevel]++;
 434                 iLevel = MAX(0,MIN( 255,(int)((data-min)/(max-min)* 255+0.5)));
 435                 hist[iLevel]++;
 436             }
 437         }
 438         // Image
 439         DEBUGPRINT("H5Dcreate");
 440         dsetID = H5Dcreate(gIDChannel, tmp, H5T_IEEE_F32LE, dspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 441         status = H5Dwrite(dsetID, H5T_IEEE_F32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, image);
 442         DEBUGPRINT3("%s End: dataset %lld status %d\n", tmp, dsetID, status);
 443         memoryFree(image);
 444         status = H5Dclose(dsetID);
 445         status = H5Sclose(dspaceID);
 446         // Histogram
 447         lhdfAttributeInt2String(gIDChannel, "HistogramMax", max);
 448         lhdfAttributeInt2String(gIDChannel, "HistogramMin", min);
 449         lhdfAttributeInt2String(gIDChannel, "HistogramMax1024", max);
 450         lhdfAttributeInt2String(gIDChannel, "HistogramMin1024", min);
 451
 452         // Histogram
 453         dspaceID = H5Screate_simple(1, &nHist, NULL);
 454         sprintf(tmp, "%s/%s %d/%s %d/%s %d/%s", groupName[0], "ResolutionLevel", ires, "TimePoint", itime, "Channel", ichannel, "Histogram");
 455         dsetID = H5Dcreate(gIDChannel, tmp, H5T_STD_U64LE, dspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 456         status = H5Dwrite(dsetID, H5T_STD_U64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, hist);
 457         status = H5Dclose(dsetID);
 458         status = H5Sclose(dspaceID);
 459
 460         // Histogram1024
 461         dspaceID = H5Screate_simple(1, &nHist1024, NULL);
 462         sprintf(tmp, "%s/%s %d/%s %d/%s %d/%s", groupName[0], "ResolutionLevel", ires, "TimePoint", itime, "Channel", ichannel, "Histogram1024");
 463         dsetID = H5Dcreate(gIDChannel, tmp, H5T_STD_U64LE, dspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 464         status = H5Dwrite(dsetID, H5T_STD_U64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, hist1024);
 465         status = H5Dclose(dsetID);
 466         status = H5Sclose(dspaceID);
 467
 468         // Group Close
 469         status = H5Gclose(gIDChannel);
 470     }
 471         status = H5Gclose(gIDTime);
 472     }
 473         dims[0]/=2;
 474         dims[1]/=2;
 475         dims[2]/=2;
 476         status = H5Gclose(gIDRes);
 477     }
 478     // DataSetInfo
 479     DEBUGPRINT("DataSetInfo\n");
 480     {
 481         char* dsiGroupName[] = {"Channel 0", "Image", "Imaris", "ImarisDataSet", "Log", "TImeInfo"};
 482         int dsiNGroup = 6;
 483         hid_t dsiGroupID[dsiNGroup];
 484
 485         for(i=0; i<dsiNGroup; i++) {
 486             sprintf(stmp, "%s/%s", groupName[1], dsiGroupName[i]);
 487             dsiGroupID[i] = H5Gcreate2(*out, stmp, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 488         }
 489         // Channel 0
 490         lhdfAttributeString(dsiGroupID[0], "Name", "No");
 491         lhdfAttributeString(dsiGroupID[0], "Color", "1.000 0.000 0.000");
 492         lhdfAttributeString(dsiGroupID[0], "ColorMode", "BaseColor");
 493         lhdfAttributeReal2String(dsiGroupID[0], "ColorOpacity", 1.000);
 494         lhdfAttributeReal2String(dsiGroupID[0], "Min", in[0].HeaderAMax);
 495         lhdfAttributeReal2String(dsiGroupID[0], "Max", in[0].HeaderAMin);
 496         // Image
 497         lhdfAttributeString(dsiGroupID[1], "Description", "No");
 498         lhdfAttributeString(dsiGroupID[1], "Name", "No");
 499         lhdfAttributeInt2String(dsiGroupID[1], "Noc", 1);
 500         lhdfAttributeReal2String(dsiGroupID[1], "ExtMin0", in[0].HeaderStartN.x*in[0].HeaderLength.x);
 501         lhdfAttributeReal2String(dsiGroupID[1], "ExtMin1", in[0].HeaderStartN.y*in[0].HeaderLength.y);
 502         lhdfAttributeReal2String(dsiGroupID[1], "ExtMin2", in[0].HeaderStartN.z*in[0].HeaderLength.z);
 503         lhdfAttributeReal2String(dsiGroupID[1], "ExtMax0", (in[0].HeaderStartN.x+dims[2]-1)*in[0].HeaderLength.x);
 504         lhdfAttributeReal2String(dsiGroupID[1], "ExtMax1", (in[0].HeaderStartN.y+dims[1]-1)*in[0].HeaderLength.y);
 505         lhdfAttributeReal2String(dsiGroupID[1], "ExtMax2", (in[0].HeaderStartN.z+dims[0]-1)*in[0].HeaderLength.z);
 506         if(1<nres) {
 507             lhdfAttributeString(dsiGroupID[1], "ResampleDimensionX", "true");
 508             lhdfAttributeString(dsiGroupID[1], "ResampleDimensionY", "true");
 509             lhdfAttributeString(dsiGroupID[1], "ResampleDimensionZ", "true");
 510         }
 511         lhdfAttributeString(dsiGroupID[1], "Unit", "um");
 512         lhdfAttributeInt2String(dsiGroupID[1], "X", dims[2]);
 513         lhdfAttributeInt2String(dsiGroupID[1], "Y", dims[1]);
 514         lhdfAttributeInt2String(dsiGroupID[1], "Z", dims[0]);
 515         // Imaris
 516         lhdfAttributeString(dsiGroupID[2], "Version", "5,5");
 517         // ImarisDataSet
 518         lhdfAttributeString(dsiGroupID[3], "Creator", "lmrc2hdf5imaris");
 519         lhdfAttributeString(dsiGroupID[3], "NumberOfImages", "1");
 520         lhdfAttributeString(dsiGroupID[3], "Version", "1.0");
 521
 522         for(i=dsiNGroup-1; i>=0; i--) {
 523             status = H5Gclose(dsiGroupID[i]);
 524         }
 525
 526     }
 527     DEBUGPRINT("DataSetTimes\n");
 528     DEBUGPRINT("Thumbnail\n");
 529     // Thumbnail
 530     {
 531       hsize_t dimsT[4];
 532       double zoomx, zoomy;
 533       hsize_t dimsX, dimsY;
 534       float orgx, orgy, orgz;
 535       zoomx = zoomy = 1;
 536       if(in[0].HeaderN.x<in[0].HeaderN.y) {
 537           zoomx *= in[0].HeaderN.x/(double)in[0].HeaderN.y;
 538       } else {
 539           zoomy *= in[0].HeaderN.y/(double)in[0].HeaderN.x;
 540       }
 541       dimsT[0] = dimsX = ((int)(256*zoomx));
 542       dimsT[1] = ((int)(256*zoomy))*4; // RGBA
 543       dimsY = dimsT[1]/4;
 544       dimsT[2]=dimsT[3]=1;
 545         DEBUGPRINT2("zoom: %f %f\n", zoomx, zoomy);
 546         DEBUGPRINT4("dimsT: %llu %llu - %llu %llu\n", dimsT[0], dimsT[1], dimsX, dimsY);
 547         sprintf(tmp, "%s/%s", groupName[3], "Data");
 548         thumbImage = (unsigned char*) memoryAllocate(sizeof(unsigned char)*dimsT[0]*dimsT[1], "in imaris");
 549         lpro.mode = 0;
 550         lmrcImageProjectionMIP(&tmpImage, &in[0], &lpro);
 551         {
 552 #ifdef DEBUG2
 553             mrcFileWrite(&tmpImage, "/tmp/tmpImage.mrc", "in tet", 0);
 554 #endif
 555         }
 556
 557         min = FLT_MAX;
 558         max = FLT_MIN;
 559         sum = 0;
 560         for(x=0; x<in[0].HeaderN.x; x++) {
 561         for(y=0; y<in[0].HeaderN.y; y++) {
 562             mrcPixelDataGet(&tmpImage, x, y, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
 563             sum += data;
 564             if(data<min) min = data;
 565             if(max<data) max = data;
 566         }
 567         }
 568         mean = sum/in[0].HeaderN.x/in[0].HeaderN.y;
 569
 570         DEBUGPRINT1("min:  %f\n", min);
 571         DEBUGPRINT1("max:  %f\n", max);
 572         DEBUGPRINT1("mean: %f\n", mean);
 573         for(x=0; x<dimsX; x++) {
 574         for(y=0; y<dimsY; y++) {
 575             orgx =  x/(dimsX - 1)*(in[0].HeaderN.x-1);
 576             orgy =  y/(dimsY - 1)*(in[0].HeaderN.y-1);
 577             mrcPixelDataGet(&tmpImage, orgx, orgy, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
 578             if(max==min) {
 579                 thumbImage[(int)(y*4+x*dimsT[1])] = 0;
 580             } else {
 581                 data = (int)(255.0*(data-(mean - (mean - min)/3))/((mean + (max-mean)/3)  - (mean - (mean-min)/3))+0.5);
 582                 if(data<0) data = 0;
 583                 if(255<data) data = 255;
 584                 thumbImage[(int)(y*4+x*dimsT[1])] = data;
 585             }
 586             thumbImage[(int)(y*4+1+x*dimsT[1])] = 0;
 587             thumbImage[(int)(y*4+2+x*dimsT[1])] = 0;
 588             thumbImage[(int)(y*4+3+x*dimsT[1])] = 255;
 589         }
 590         }
 591
 592         dspaceID = H5Screate_simple(2, &(dimsT[0]), NULL);  // H5S dimension
 593         sprintf(tmp, "%s/%s", groupName[3], "Data");        // Group Name
 594         dsetID = H5Dcreate2(groupID[3], tmp, H5T_STD_U8LE, dspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 595         status = H5Dwrite(dsetID, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL, H5P_DEFAULT, thumbImage);  // Write
 596         mrcImageFree(&tmpImage, "in lmrc2hdfimaris");
 597
 598         status = H5Dclose(dsetID);
 599         status = H5Sclose(dspaceID);
 600     }
 601     for(i=nGroup-1; i>=0; i--) {
 602         status = H5Gclose(groupID[i]);
 603     }
 604
 605     return 0;
 606 }
 607 int
 608 lmrc2hdf(hid_t* out, mrcImage* in, lmrc2hdfInfo* linfo, int mode)
 609 {
 610     char tmp[1024];
 611     hsize_t dims[4];
 612     hsize_t chunk[4];
 613     herr_t status;
 614     hid_t datasetID;
 615     hid_t datasetID2;
 616     hid_t dataspaceID;
 617     hid_t dataspaceID2;
 618     hid_t datasetPropertyID;
 619     hid_t dset;
 620     hid_t nGroup=4;
 621     hid_t groupID[4];
 622     char* groupName[] = {"/DataSet", "/DataSet/ResolutionLevel 0", "/DataSet/ResolutionLevel 0/TimePoint 0", "/DataSet/ResolutionLevel 0/TimePoint 0/Channel 0"};
 623     //char* groupName[] = {"/DataSet"};
 624     int val;
 625     int i;
 626
 627     nGroup = 1;
 628     for(i=0; i<nGroup; i++) {
 629         DEBUGPRINT1("Group: %s\n", groupName[i]);
 630         groupID[i] = H5Gcreate2(*out, groupName[i], H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 631         if(groupID[i]<0) {
 632             fprintf(stderr, "Not create Group: %s\n", groupName[i]);
 633             exit(EXIT_FAILURE);
 634         } else {
 635             DEBUGPRINT("Success\n");
 636         }
 637     }
 638
 639     dims[0] = in->HeaderN.x;
 640     dims[1] = in->HeaderN.y;
 641     dims[2] = in->HeaderN.z;
 642     dims[3] = linfo->numFile;
 643
 644     chunk[0] = dims[0];
 645     chunk[1] = dims[1];
 646     chunk[2] = dims[2];
 647     chunk[3] = dims[3];
 648
 649     datasetPropertyID = H5Pcreate(H5P_DATASET_CREATE);
 650     status = H5Pset_chunk(datasetPropertyID, 4, chunk);
 651
 652     val = 4;
 653     status = H5Pset_fill_value(datasetPropertyID, H5T_NATIVE_INT, &val);
 654     status = H5Pset_alloc_time(datasetPropertyID, H5D_ALLOC_TIME_EARLY);
 655
 656     dataspaceID2 = H5Screate_simple(4, dims, NULL);
 657     //datasetID2 = H5Dcreate(*out, "/prop", H5T_STD_I32LE, dataspaceID2, H5P_DEFAULT, datasetPropertyID, H5P_DEFAULT);
 658     i=0;
 659     sprintf(tmp, "%s/%s", groupName[i], "Data");
 660     datasetID2 = H5Dcreate2(groupID[i], tmp, H5T_STD_I32LE, dataspaceID2, H5P_DEFAULT, datasetPropertyID, H5P_DEFAULT);
 661     status = H5Pclose(datasetPropertyID);
 662     status = H5Dclose(datasetID2);
 663     status = H5Sclose(dataspaceID2);
 664
 665
 666     dataspaceID = H5Screate_simple(4, dims, NULL);
 667     switch(in->HeaderMode) {
 668         case mrcCharImage:
 669             datasetID = H5Dcreate(*out, "/mrc", H5T_STD_U8LE, dataspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 670             status = H5Dwrite(datasetID, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, in->Image);
 671             break;
 672         case mrcShortImage:
 673             datasetID = H5Dcreate(*out, "/mrc", H5T_STD_I16LE, dataspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 674             status = H5Dwrite(datasetID, H5T_NATIVE_SHORT, H5S_ALL, H5S_ALL, H5P_DEFAULT, in->Image);
 675             break;
 676         case mrcUShortImage:
 677             datasetID = H5Dcreate(*out, "/mrc", H5T_STD_U16LE, dataspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 678             status = H5Dwrite(datasetID, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL, H5P_DEFAULT, in->Image);
 679             break;
 680         case mrcFloatImage:
 681             datasetID = H5Dcreate(*out, "/mrc", H5T_IEEE_F32LE, dataspaceID, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 682             status = H5Dwrite(datasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, in->Image);
 683             break;
 684         default:
 685             fprintf(stderr, "Not supported: %d\n", in->HeaderMode);
 686             exit(EXIT_FAILURE);
 687             break;
 688     }
 689
 690     status = H5Dclose(datasetID);
 691     status = H5Sclose(dataspaceID);
 692
 693     for(i=nGroup-1; i>=0; i--) {
 694         status = H5Gclose(groupID[i]);
 695     }
 696
 697     return 0;
 698 }