Skip to content

GitLab

P
pyhst2
Commit 40ef8641 authored by myron's avatar myron
Browse files
Options

scritto h5eli. Da compilare

parent 605281c0
Hide whitespace changes
Inline Side-by-side
Showing with 366 additions and 65 deletions
PyHST/Cspace/CCspace.c
View file @ 40ef8641
g
#/*##########################################################################
# Copyright (C) 2001-2013 European Synchrotron Radiation Facility
#
......@@ -5155,13 +5155,6 @@ void CCspace_preprocess_chunk(CCspace * self , int sn, int ntok, int ntokt ,
}
int num_myprojs_reduced_forAllThreads = num_myprojs_reduced_forthisproc/npbunches ;
// @@@@@@@
Size0 = size_s_[2*sn_nb +0];
Size1 = size_s_[2*sn_nb +1];
......@@ -6395,24 +6388,8 @@ void CCspace_read_chunk (CCspace * self , int sn, int ntok , int npbunches,
);
}
}
// printf(" DEBUG iread %d npj %d self->reading_infos.headerSizes[iread] %ld \n", iread, npj, self->reading_infos.headerSizes[iread] );
/* { */
/* FILE *f=fopen("del.txt","a"); */
/* fprintf(f,"%p\n", target ); */
/* fclose(f); */
/* } */
/* printf(" MEMSET %p \n", target ); */
memset(target, 0, size_s_[2*sn_nb +0]*size_s_[2*sn_nb +1]*sizeof(float) ) ;
// printf(" LEGGO ELI %s %d %d %d %d \n", proj_file_list[iread], pos_s_[2*sn_nb], pos_s_[2*sn_nb+1], size_s_[2*sn_nb +0], size_s_[2*sn_nb +1] );
int npj_absolute =self->reading_infos.tot_proj_num_list[npj+num_previous_projs] + self->reading_infos. proj_num_offset_list[sn_nb] ;
read_data_from_edf_eli( proj_file_list[iread],
target , // target_float
self->reading_infos.datatype, //
......@@ -6436,24 +6413,10 @@ void CCspace_read_chunk (CCspace * self , int sn, int ntok , int npbunches,
self->reading_infos.currents,
self->reading_infos.currents[npj_absolute] , &(self->filereading_sem)
);
/* FILE *f=fopen("pippo.dat", "a"); */
/* // printf("DEBUG iread %d %d %d %s \n", iread, size_s_[2*sn_nb +0], size_s_[2*sn_nb +1],proj_file_list[iread] ); */
/* fwrite( target + size_s_[2*sn_nb +1] , size_s_[2*sn_nb +1]*sizeof(float) , 1, f ) ; */
/* fclose(f); */
/* { */
/* float SUM=0.0; */
/* int k; */
/* for(k=0; k< size_s_[2*sn_nb +1]*size_s_[2*sn_nb +1]; k++) { */
/* SUM += target[k]; */
/* } */
/* // printf("DEBUG SUMA iread %d %e %s\n",iread, SUM, proj_file_list[iread]); */
/* } */
}
} else if (strcmp(self->reading_infos.proj_reading_type , "h5")==0){
if(self->params.DZPERPROJ==0) {
read_data_from_h5( proj_file_list[ 0 ],
self->reading_infos.proj_h5_dsname,
......@@ -6467,14 +6430,69 @@ void CCspace_read_chunk (CCspace * self , int sn, int ntok , int npbunches,
0.0 // threshold
);
} else {
printf(" ERROR : reading from hdf5 not yet adapted to elicoidal scan \n");
exit(1);
}
float *Bptr = NULL;
if(self->params.SUBTRACT_BACKGROUND ) { Bptr = self->background ; }
int iff;
long int posff;
for(iff=0; iff<2; iff++) {
posff = (long int) self->reading_infos.ff_indexes_coefficients[ 6*iread+iff ];
if(posff!=-1 && posff!=last_posFF[iff]) {
last_posFF[iff]=posff;
if( FF[iff]==NULL) {
FF[iff]= (float *) malloc( self->reading_infos.Dim_2*self->reading_infos.Dim_1 * sizeof(float));
}
int file_seq_pos = min(self->reading_infos.ff_file_list_lenght-1, posff ) ;
int multiplo = 1 ;
if( self->reading_infos.ff_file_list_lenght == 1 ) {
multiplo = 0 ;
}
read_data_from_h5(self->reading_infos.ff_file_list[file_seq_pos] ,
self->reading_infos.ff_h5_dsname,
posff - file_seq_pos ,
FF[iff] , // target_float
0 ,
0 ,
self->reading_infos.Dim_2 ,
self->reading_infos.Dim_1 ,
rotation_vertical,
binning,
multiplo,
self->params.REFERENCES_HOTPIXEL_THRESHOLD
);
}
}
memset(target, 0, size_s_[2*sn_nb +0]*size_s_[2*sn_nb +1]*sizeof(float) ) ;
read_data_from_h5_eli( proj_file_list[ 0 ],
self->reading_infos.proj_h5_dsname,
iread ,
target , // target_float
pos_s_[2*sn_nb], pos_s_[2*sn_nb+1],
size_s_[2*sn_nb +0], size_s_[2*sn_nb +1],
rotation_vertical,
binning,
0, // multiplo
0.0, // threshold
self->params.DZPERPROJ,
Bptr,
self->params.CORRECT_FLATFIELD,
self->reading_infos.ff_indexes_coefficients + 6*iread,
FF,
self->params.DOUBLEFFCORRECTION!=NULL ,
self->params.TAKE_LOGARITHM,
self->ffcorr,
&(self->filereading_sem)
);
}// fine conico con H5
} else {
fprintf(stderr, "file type can be either edf or h5 reading projs but it was %s \n",self->reading_infos.proj_reading_type );
exit(1);
}
} else {
// printf(" metto a NAN\n");
int k;
......
PyHST/Cspace/edftools.c
View file @ 40ef8641
......@@ -405,6 +405,9 @@ void read_data_from_h5( char * filename,
herr_t status;
int rank ;
assert(binning==1);
printf(" reading %s dataset %s at nproj %d \n", filename, dsname, nproj ) ;
file = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
......@@ -505,6 +508,196 @@ void read_data_from_h5( char * filename,
}
void read_data_from_h5_eli( char * filename,
char * dsname,
int nproj,
float *target , // target_float
long int pos0 , long int pos1 ,
long int size0 , long int size1,
int rotation_vertical,
int binning,
float DZPERPROJ,
float *Ddark,
int correct_ff,
float *ff_infos,
float **FFs,
int doubleffcorrection,
int takelog,
float *ffcorr,
sem_t * semaforo_ptr
) {
assert(binning==1);
for(i=0; i < size0; i++){
for(j=0; j < size1; j++){
target[i*size1+j]= NAN ;
}
}
hid_t file; /* handles */
hid_t dataset;
hid_t memspace;
hid_t dataspace;
hsize_t dimsm[3]; /* dataset and chunk dimensions*/
hsize_t count[3];
hsize_t offset[3];
hsize_t count_out[3];
hsize_t offset_out[3];
herr_t status;
int rank ;
float shift = iread * DZPERPROJ;
int ishift = ((int) shift);
shift = shift-ishift; // si assume dzperproj sempre positivo
pos0 = pos0 -ishift ;
int pos0_ = pos0-1;
int pos0_clipped=pos0;
if(pos0<0) pos0_clipped=0;
int pos0_clipped_ = pos0_;
if(pos0_<0) pos0_clipped_=0;
int size0_clipped = size0;
if((size0_clipped + pos0) > size1 ) {
size0_clipped = size1 - pos0;
}
int size0_clipped_ = size0_clipped+ 1 ;
if( pos0_ + size0_clipped_ >= pos0_clipped_ ) {
sem_wait( semaforo_ptr );
printf(" reading %s dataset %s at nproj %d \n", filename, dsname, nproj ) ;
file = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
dataset = H5Dopen(file, dsname);
dataspace = H5Dget_space(dataset);
rank = H5Sget_simple_extent_ndims(dataspace);
assert(rank==3);
hsize_t dims_out[rank];
H5Sget_simple_extent_dims( dataspace , dims_out , NULL );
/*
* descrive quale parte del dataset leggo
*/
offset[0] = nproj;
offset[1] = pos0_clipped_;
offset[2] = pos1;
int v_span = 1;
count[0] = v_span ;
count[1] = size0_clipped_-(pos0_clipped_-pos0_) ;
count[2] = size1;
status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL,
count, NULL);
/*
* descrive la memoria su cui scrivo
*/
dimsm[0] = v_span;
dimsm[1] = size0_clipped_-(pos0_clipped_-pos0_);
dimsm[2] = size1 ;
memspace = H5Screate_simple( 3, dimsm, NULL);
/*
* descrive quale parte di memoria scrivo
*/
offset_out[0] = 0;
offset_out[1] = 0;
offset_out[2] = 0;
count_out[0] = v_span;
count_out[1] = size0_clipped_-(pos0_clipped_-pos0_);
count_out[2] = size1;
status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, offset_out, NULL,
count_out, NULL);
/*
* Read data from hyperslab in the file into the hyperslab in
* memory and display.
*/
float *mytarget;
mytarget = (float *) malloc( v_span*( size0_clipped_ )*size1* sizeof(float)) ;
status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace,
H5P_DEFAULT, mytarget +(pos0_clipped_-pos0_)*size1 );
H5Dclose(dataset);
H5Sclose(dataspace);
H5Sclose(memspace);
H5Fclose(file);
sem_post( semaforo_ptr );
if(pos0_clipped_-pos0_) {
memset( mytarget , 0, (pos0_clipped_-pos0_)* sizeof(float) *((long int) size1) );
}
for(j=0; j< size0_clipped_ ; j++){
for(i=0; i< size1 ; i++){
if((j+pos0_)>=0) {
if(Ddark) {
mytarget [ j*size1+i ] = (mytarget[j*size1+i ] - Ddark[(j+pos0_)*size1+i]) ;
}
if(correct_ff) {
float deno=1;
if(ff_infos[0]!=-1.0f) {
deno = FFs[0][ (j+pos0_)*size1+i] ;
}
if(ff_infos[1]!=-1.0f) {
deno = FFs[1][ (j+pos0_)*size1+i]*ff_infos[3] + ff_infos[2]*deno ;
}
mytarget[j*size1+i ] /= deno;
}
}
}
}
if(doubleffcorrection) {
if(takelog ) {
for(j=0; j<size0; j++) {
for(i=0; i< size1; i++) {
if((j+pos0_)>=0) mytarget[j*size1+i ] *= ffcorr[(j+pos0_)*size1+i] ;
}
}
} else {
for(j=0; j<size0; j++) {
for(i=0; i< size1; i++) {
if((j+pos0_)>=0) mytarget[j*size1+i ] -= ffcorr[(j+pos0_)*size1+i] ;
}
}
}
}
if( rotation_vertical==0 ) {
printf("MISSING FEATURE: helicoidal scan with horizontal axis not yet implemented\n");// da rivedere accuratamente
fprintf(stderr, "MISSING FEATURE: helicoidal scan with horizontal axis not yet implemented \n");
}
for(j=0; j< (pos0_clipped_-pos0_); j++){
for(i=0; i< size1; i++){
mytarget[j*size1+i]= NAN;
}
}
for(i=pos0; i < pos0 + size0_clipped ; i++){
for(j=0; j < size1; j++){
target[(i-pos0)*size1+j]= mytarget[ (i-pos0_) * size1 + j ] *(1-shift) + shift* mytarget[ (i-pos0_-1) * size1 + j ] ;
}
}
} else {
}
sem_post( semaforo_ptr );
}
void read_data_from_edf( char * filename,
......@@ -834,7 +1027,7 @@ void read_data_from_edf_eli( char * filename,
}
}
if(pos0_clipped-pos0) {
if(pos0_clipped_-pos0_) {
memset( buffer , 0, (pos0_clipped_-pos0_)* sizeofdatatype *((long int) size1) * binning*binning );
}
......@@ -869,27 +1062,7 @@ void read_data_from_edf_eli( char * filename,
/* memcpy( target_float, bufferF,size0*size1*sizeof(float) ); */
/* return; */
for(j=0; j<n1*binning; j++){
for(i=0; i< n2*binning; i++){
if((j+pos0_)>=0) {
if(Ddark) {
bufferF[j*n2*binning+i ] = (bufferF[j*n2*binning+i ] - Ddark[(j+pos0_)*Dim_1+i])/curr ;
} else {
bufferF[j*n2*binning+i ] = (bufferF[j*n2*binning+i ] )/ curr ;
}
if(correct_ff) {
float deno=1;
if(ff_infos[0]!=-1.0f) {
deno = FFs[0][ (j+pos0_)*Dim_1+i] ;
}
if(ff_infos[1]!=-1.0f) {
deno = FFs[1][ (j+pos0_)*Dim_1+i]*ff_infos[3] + ff_infos[2]*deno ;
}
bufferF[j*n2*binning+i ] /= deno;
}
}
}
}
if(doubleffcorrection) {
if(takelog ) {
......
PyHST/Cspace/edftools.h
View file @ 40ef8641
......@@ -121,6 +121,29 @@ void read_data_from_edf_eli( char * filename,
);
void read_data_from_h5_eli( char * filename,
char * dsname,
int nproj,
float *target , // target_float
long int pos0 , long int pos1 ,
long int size0 , long int size1,
int rotation_vertical,
int binning,
float DZPERPROJ,
float *Ddark,
int correct_ff,
float *ff_infos,
float **FFs,
int doubleffcorrection,
int takelog,
float *ffcorr,
sem_t * semaforo_ptr
) ;
void write_data_to_edf( float *SLICE, int num_y,int num_x, const char * nomeout ) ;
......
PyHST/conversion.py 0 → 100644
View file @ 40ef8641
import h5py
import sys
import numpy as np
import os
def create_virtual_file( vsource, pro_indxs, datasource_shape, mydtype, filename ):
layout = h5py.VirtualLayout(shape= tuple([ len( pro_indxs ) ]+ list(datasource_shape)[1:]), dtype=mydtype)
for n in range( len( pro_indxs ) ):
layout[n] = vsource[ pro_indxs[n] ]
with h5py.File(filename, "w", libver="latest") as f:
f.create_virtual_dataset("data", layout=layout, fillvalue=-5)
fn = sys.argv[1]
dn = sys.argv[2]
rep = sys.argv[3]
pyhst2_input = """
"""
entry = h5py.File(fn,"r")[dn]
pro_indxs = []
ff_indxs = []
ff_poss = []
bck_indxs = []
keys = entry["instrument/detector"]["image_key"]
old_ff_pos = -1
ffblock = None
ipro = 0
for i,k in enumerate(list(keys)):
if k== 2:
bck_indxs.append(i)
if k== 0:
pro_indxs.append(i)
ipro += 1
if k==1:
if ipro== old_ff_pos:
ffblock.append(i)
else:
ffblock = []
ff_indxs.append( ffblock )
ff_poss.append(ipro )
old_ff_pos = ipro
datasource_shape = entry["instrument/detector/data"] .shape
mydtype = entry["instrument/detector/data"].dtype
vsource = h5py.VirtualSource(fn, dn+"/instrument/detector/data", shape= datasource_shape )
filename = os.path.join(rep,"pyhst2_projs.h5")
create_virtual_file( vsource, pro_indxs, datasource_shape, mydtype, filename )
pyhst2_input += """
FILE_PREFIX={}
PROJ_DS_NAME = "data"
""".format( os.path.join(rep,"pyhst2_projs.h5") )
if len(bck_indxs):
filename = os.path.join(rep,"pyhst2_background.h5")
create_virtual_file( vsource, bck_indxs, datasource_shape, mydtype, filename )
pyhst2_input += """
BACKGROUND_FILE={}
BACKGROUND_DS_NAME = "data"
""".format( os.path.join(rep,"pyhst2_projs.h5") )
s_files="""["""
s_intervals="""["""
if len(ff_indxs):
for ffblock, pos in zip(ff_indxs, ff_pos ) :
s_intervals += """ {} ,""".format(pos)
fname = "pyhst2_ff%05d.h5"%ff_pos
fname = os.path.join(rep, fname)
s_files+=""" "{}",""".format(fname=
create_virtual_file( vsource, ffblock, datasource_shape, mydtype, fname )
pyhst2_input += """
FF_PREFIX={ffprefix}]
FF_DS_NAME = "data"
FF_INTERVALS={intervals}]
""".format( ffprefix=s_files[:-1] , intervals=s_intervals[:-1] )
print (pyhst2_input)
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