Commit dd37eeb3 authored by Alessandro Mirone's avatar Alessandro Mirone Committed by Alessandro Mirone

anatomix(h5+nxs)

compila

funziona OK

ripulito
parent f976c069
......@@ -6454,7 +6454,9 @@ void CCspace_read_chunk (CCspace * self , int sn, int ntok , int npbunches,
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
binning,
0, // multiplo
0.0 // threshold
);
} else {
printf(" ERROR : reading from hdf5 not yet adapted to elicoidal scan \n");
......@@ -6501,18 +6503,32 @@ void CCspace_read_chunk (CCspace * self , int sn, int ntok , int npbunches,
self->reading_infos.ff_byteorder,
binning, &(self->filereading_sem)
);
} else if (strcmp(self->reading_infos.proj_reading_type , "h5")){
read_data_from_h5(self->reading_infos.ff_file_list[ 0 ] ,
} else if (strcmp(self->reading_infos.ff_reading_type , "h5")==0){
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,
(FFptr+ posff*size_s_[2*sn_nb +0]*size_s_[2*sn_nb +1]) , // 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],
posff - file_seq_pos ,
(FFptr+ posff*size_s_[2*sn_nb +0]*size_s_[2*sn_nb +1]), // 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
binning,
multiplo,
self->params.REFERENCES_HOTPIXEL_THRESHOLD
);
} else {
fprintf(stderr, "file type MUST be either edf or h5, reading ffs \n");
fprintf(stderr, "file type MUST be either edf or h5, reading ffs it was %s\n", self->reading_infos.ff_reading_type);
exit(1);
}
......@@ -7771,7 +7787,6 @@ void CCD_Filter_Implementation(float * Tptr, float * Rawptr,
int ds;
// char msg[400];
int c0,c1;
int c0_,c1_;
// float temp[9];
......@@ -7785,7 +7800,6 @@ void CCD_Filter_Implementation(float * Tptr, float * Rawptr,
float * a_data;
float * m_data;
a_data = Rawptr ;
m_data = Tptr ;
......@@ -7798,16 +7812,12 @@ void CCD_Filter_Implementation(float * Tptr, float * Rawptr,
x0_ = 0 ;
x1_ = 0 ;
X0_ = x0_ + s0_ ;
X1_ = x1_ + s1_ ;
dim0 = s0_;
dim1 = s1_;
memcpy(m_data ,a_data , dim0*dim1*sizeof(float) );
#define AA(i,j) a_data[(i)*dim1+(j) ]
......@@ -9230,8 +9240,63 @@ void vect_minus_fact_vect(float *A ,float s, float *B, int dim) {
}
}
#define SWAP(x, y) { int temp = x; x = y; y = temp; }
int partition(float *a, int left, int right, int pivotIndex)
{
int pivot = a[pivotIndex];
SWAP(a[pivotIndex], a[right]);
int pIndex = left;
int i;
for (i = left; i < right; i++)
{
if (a[i] <= pivot)
{
SWAP(a[i], a[pIndex]);
pIndex++;
}
}
SWAP(a[pIndex], a[right]);
return pIndex;
}
float quickselect(float * A, int left, int right, int k)
{
while(1) {
if (left == right) {
return A[left];
}
// int pivotIndex = left + rand() % (right - left + 1);
int pivotIndex = left + (right - left + 1)/2;
pivotIndex = partition(A, left, right, pivotIndex);
if (k == pivotIndex) {
return A[k];
} else if (k < pivotIndex) {
right = pivotIndex - 1 ;
} else {
left = pivotIndex + 1 ;
}
}
}
float * medianX( float *A , int dim3 , int dim2 , int dim1 , int nproc ) {
float *res = (float *) malloc(sizeof(float) *dim3*dim2*dim1);
assert(dim1>1);
int i;
#pragma omp parallel for num_threads (nproc)
for ( i=0; i<dim3*dim2; i++) {
res[i] = quickselect(A+i*dim1 , 0, dim1 -(dim1%2), (dim1 -(dim1%2) -1 )/2);
}
return res;
}
/* void passeggiaSLICE(CCspace *self, float **WORK, float **WORKbis, int num_bins, int dim_fft, */
/* float **dumf, fcomplex **dumfC , float * WORK_perproje, float **OVERSAMPLE, */
/* int oversampling, int ncpus , float cpu_offset_x, float cpu_offset_y, */
......
......@@ -39,6 +39,9 @@
#include<cufft.h>
float quickselect(float * A, int left, int right, int k);
typedef float __complex__ fcomplex ;
// typedef fftwf_complex fcomplex ;
......@@ -219,6 +222,7 @@ typedef struct {
/* Cparameters */
typedef struct {
int SUBTRACT_BACKGROUND;
float REFERENCES_HOTPIXEL_THRESHOLD;
int CORRECT_FLATFIELD;
float OFFSETRADIOETFF;
int TAKE_LOGARITHM;
......@@ -1439,3 +1443,8 @@ void passeggiaSLICE(CCspace *self, float **WORK, float **WORKbis, int num_bins,
float denoising_driver(CCspace *self,int dim0,int dim1,float *img,float *result , float weight ) ;
void rotational2zero(CCspace *self,float *SLICE,float *SLICEres) ;
float * medianX( float *A , int dim3 , int dim2 , int dim1 , int nproc );
void CCD_Filter_Implementation(float * Tptr, float * Rawptr,
int Size0, int Size1 ,
float threshold, int ncpus) ;
......@@ -480,18 +480,47 @@ Cspace_new( PyObject *self, PyObject *args)
assert( res->myCCspace->params.dist_ncol = ncol );
}
}
res->Offcorr = cpyutils_get_pydic_Oval(Oarrayspace , "doubleffcorrection" ) ;
Py_INCREF(res->Offcorr);
res->myCCspace->ffcorr = (float *) PyArray_DATA((PyArrayObject *)res->Offcorr) ;
res->myCCspace->params.REFERENCES_HOTPIXEL_THRESHOLD = cpyutils_o2f(cpyutils_getattribute(res->OParameters ,"REFERENCES_HOTPIXEL_THRESHOLD" ) );
res->Obackground = cpyutils_get_pydic_Oval(Oarrayspace , "background" ) ;
assert(PyArray_Check( res->Obackground ));
Py_INCREF(res->Obackground);
res->myCCspace->background = (float *) PyArray_DATA((PyArrayObject *)res->Obackground) ;
{
assert(PyArray_Check( res->Obackground ));
Py_INCREF(res->Obackground);
assert( PyArray_NDIM((PyArrayObject*) res->Obackground ) == 2 || PyArray_NDIM((PyArrayObject*) res->Obackground ) == 3 );
int dim2,dim1;
if( PyArray_NDIM((PyArrayObject*) res->Obackground ) ==3 ){
printf(" PyArray_NDIM((PyArrayObject*) res->Obackground ) ==3\n");
int dimx;
res->myCCspace->background = (float*) cpyutils_PyArray3D_as_array( res->Obackground , &dim2, &dim1, &dimx , NPY_FLOAT32 );
printf(" median npbunches \n" );
int npbunches = cpyutils_o2i(cpyutils_getattribute(res->OParameters ,"NPBUNCHES" ) ); ;
res->myCCspace->background = medianX( res->myCCspace->background , dim2, dim1, dimx , npbunches );
printf(" median applied \n");
} else {
res->myCCspace->background = (float*) cpyutils_PyArray2D_as_array( res->Obackground , &dim2, &dim1, NPY_FLOAT32 );
}
if ( res->myCCspace->params.REFERENCES_HOTPIXEL_THRESHOLD ) {
float * new_back = (float*) malloc( dim2*dim1* sizeof(float) ) ;
CCD_Filter_Implementation( new_back, res->myCCspace->background, dim2 , dim1 ,res->myCCspace->params.REFERENCES_HOTPIXEL_THRESHOLD , 1);
res->myCCspace->background = new_back ;
}
}
if(!PyArray_Check(cpyutils_getattribute(res->OParameters ,"angles" ))) {
res->myCCspace->params.do_custom_angles=0;
} else {
......@@ -516,6 +545,8 @@ Cspace_new( PyObject *self, PyObject *args)
cspace_iproc= res->iproc;
res->myCCspace->params.SUBTRACT_BACKGROUND = cpyutils_o2i(cpyutils_getattribute(res->OParameters ,"SUBTRACT_BACKGROUND" ) );
res->myCCspace->params.CORRECT_FLATFIELD = cpyutils_o2i(cpyutils_getattribute(res->OParameters ,"CORRECT_FLATFIELD" ) );
res->myCCspace->params.OFFSETRADIOETFF = cpyutils_o2f(cpyutils_getattribute(res->OParameters ,"OFFSETRADIOETFF" ) );
......
......@@ -188,6 +188,26 @@ void *cpyutils_PyArray2D_as_array(PyObject *OtmpA, int *dim2, int *dim1,int py
return res;
}
void *cpyutils_PyArray3D_as_array(PyObject *OtmpA,int *dim3, int *dim2, int *dim1,int pyarraytype ) {
PyObject *Otmp;
void *res ;
assert(PyArray_Check(OtmpA ));
assert( PyArray_DESCR((PyArrayObject *)OtmpA)->type_num == pyarraytype );
assert( PyArray_NDIM((PyArrayObject*) OtmpA) ==3 );
Otmp = PyArray_ContiguousFromObject(OtmpA, pyarraytype , 3, 3);
*dim3 = PyArray_DIMS((PyArrayObject *) Otmp)[0];
*dim2 = PyArray_DIMS((PyArrayObject *) Otmp)[1];
*dim1 = PyArray_DIMS((PyArrayObject *) Otmp)[2];
res = (void *) malloc( (*dim2)*(*dim2)*(*dim1) *PyArray_ITEMSIZE((PyArrayObject *) Otmp) ) ;
memcpy( res , PyArray_DATA((PyArrayObject*) Otmp), (*dim2)*(*dim1)*PyArray_ITEMSIZE((PyArrayObject *) Otmp) );
Py_DECREF(Otmp);
return res;
}
void *cpyutils_PyArray5D_as_array(PyObject *OtmpA,int *dim5,
int *dim4, int *dim3,
int *dim2, int *dim1,int pyarraytype ){
......
......@@ -40,6 +40,7 @@ void *cpyutils_PyArray1D_as_array(PyObject *Otmp, int *dim , int pyarraytype );
void *cpyutils_PyArray1D_as_array_TB(PyObject *OtmpA, int *listlenght,int pyarraytype, const char * callFILE, int callLINE );
void *cpyutils_PyArray2D_as_array(PyObject *Otmp, int *dim2, int *dim1,int pyarraytype );
void *cpyutils_PyArray3D_as_array(PyObject *Otmp, int *dim3,int *dim2, int *dim1,int pyarraytype );
void *cpyutils_PyArray5D_as_array(PyObject *Otmp,int *dim5,
int *dim4, int *dim3,
......
......@@ -43,12 +43,12 @@
#include <sys/stat.h>
#include <fcntl.h>
#include<unistd.h>
#include"CCspace.h"
#define H5Dopen_vers 1
///// temporary discativation
// #include <hdf5.h>
#include <hdf5.h>
#undef NDEBUG
#include<assert.h>
......@@ -155,10 +155,6 @@ void prereadEdfHeader( int *sizeImage ,
}
}
}
if(strcmp(c1, "DataType")==0) {
if( strcmp(pos, "UnsignedShort" )==0) {
......@@ -338,6 +334,7 @@ void extended_fread( char *ptr, /* memory to write in */
for(i=0; i<ndims; i++) {
indexes[i]=0;
}
loop=ndims-1;
indexes[ndims-1 ]=-1 ;
pos=-strides[ndims-1];
......@@ -356,7 +353,7 @@ void extended_fread( char *ptr, /* memory to write in */
printf("Error 1/n");
break;
}
res=fread ( ( (char * ) ptr) +(count)*size_of_block ,
size_of_block, 1, stream );
if(res!=1) {
......@@ -365,9 +362,9 @@ void extended_fread( char *ptr, /* memory to write in */
break;
}
count++;
oldpos = pos+ size_of_block;
loop=ndims-1;
/* for(i=0 ; i< ndims; i++) {
printf(" %d ", indexes[i] );
......@@ -391,90 +388,119 @@ void read_data_from_h5( char * filename,
long int pos0 , long int pos1 ,
long int size0 , long int size1,
int rotation_vertical,
int binning
int binning,
int multiplo,
float threshold
) {
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 ;
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
*/
fprintf(stderr, " TEMPORARY DISACTIVATION OF HDF5 in pyhst, do not use hdf5 for the moment ( due to different version on the cluster ) \n");
/* 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 ; */
/* 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); */
/* /\* */
/* * descrive quale parte del dataset leggo */
/* *\/ */
/* offset[0] = nproj; */
/* offset[1] = pos0; */
/* offset[2] = pos1; */
/* count[0] = 1; */
/* count[1] = size0; */
/* count[2] = size1; */
/* status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, */
/* count, NULL); */
/* /\* */
/* * descrive la memoria su cui scrivo */
/* *\/ */
/* dimsm[0] = 1; */
/* dimsm[1] = size0; */
/* dimsm[2] = size1 ; */
/* memspace = H5Screate_simple( 3, dimsm, NULL); */
if(multiplo) {
assert( nproj == 0 ) ;
}
offset[0] = nproj;
offset[1] = pos0;
offset[2] = pos1;
int v_span = 1;
if(multiplo) {
v_span = dims_out[0] ;
}
/* /\* */
/* * descrive quale parte di memoria scrivo */
/* *\/ */
/* offset_out[0] = 0; */
/* offset_out[1] = 0; */
/* offset_out[2] = 0; */
/* count_out[0] = 1; */
/* count_out[1] = size0; */
/* 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. */
/* *\/ */
count[0] = v_span ;
count[1] = size0;
count[2] = size1;
status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL,
count, NULL);
/* status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, */
/* H5P_DEFAULT, target); */
/*
* descrive la memoria su cui scrivo
*/
dimsm[0] = v_span;
dimsm[1] = size0;
dimsm[2] = size1 ;
memspace = H5Screate_simple( 3, dimsm, NULL);
/* status=status; */
/*
* 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;
count_out[2] = size1;
status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, offset_out, NULL,
count_out, NULL);
/* H5Dclose(dataset); */
/* H5Sclose(dataspace); */
/* H5Sclose(memspace); */
/* H5Fclose(file); */
/*
* Read data from hyperslab in the file into the hyperslab in
* memory and display.
*/
float *mytarget;
if(v_span > 1) {
mytarget = (float *) malloc( v_span*size0*size1* sizeof(float)) ;
} else {
mytarget = target;
}
status = H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace,
H5P_DEFAULT, mytarget);
status=status;
if(v_span > 1) {
float line[v_span];
int i,k;
for( i = 0; i< size0*size1; i++) {
for( k = 0 ; k < v_span ; k++ ) {
line[ k ] = mytarget[ i + k * size0*size1 ] ;
target[ i ] = quickselect(line, 0, v_span -(v_span%2), (v_span -(v_span%2) -1 )/2 );
}
}
if(threshold) {
CCD_Filter_Implementation( mytarget, target , size0 , size1 ,threshold , 1);
memcpy( target, mytarget, sizeof(float)*size0*size1 ) ;
}
free( mytarget ) ;
} else {
}
H5Dclose(dataset);
H5Sclose(dataspace);
H5Sclose(memspace);
H5Fclose(file);
}
......
......@@ -89,9 +89,12 @@ void read_data_from_h5( char * filename,
long int pos0 , long int pos1 ,
long int size0 , long int size1,
int rotation_vertical,
int binning
int binning,
int multiplo,
float threshold
);
void read_data_from_edf_eli( char * filename,
float *dataptr , // target_float
int filedatatype, //
......
......@@ -7317,9 +7317,10 @@ __global__ static void gputomo_conicity_kernel(float *d_SLICE, // da allocare
float vz, vx, vy;
float magni, px, pz;
float magnicenter = (source_distance+detector_distance)/ (source_distance ) ;
float axis_position_v = (axis_position-SOURCE_X)/(magnicenter*v2x) ;
// gpu_offset_x -= SOURCE_X;
// gpu_offset_y -= SOURCE_X;
......
......@@ -32,15 +32,10 @@ from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import string
from . import string_six
import traceback
import sys
import numpy
......@@ -115,7 +110,7 @@ def treat_par_file(s):
# if Parameters.VERBOSITY>0 : print( "PYTHONIFYING INPUT FILE")
s=s.replace("!","#")
lines=s.split("\n")
new_lines=""
for line in lines:
......@@ -129,12 +124,12 @@ def treat_par_file(s):
line=add_quotes(line)
elif( line.find("POSTFIX")>=0 and (line.find("POSTFIX") < line.find("=")) and notyetwithquotes(line) ):
line=add_quotes(line)
elif( line.find("DS_NAME")>=0 and (line.find("DS_NAME") < line.find("=")) and notyetwithquotes(line) ):
line=add_quotes(line)
elif( line.find("CURRENT_NAME")>=0 and (line.find("CURRENT_NAME") < line.find("=")) and notyetwithquotes(line) ):
line=add_quotes(line)
elif( line[:3] == "IF_" ):
line=add_quotes(line)
new_lines=new_lines+"\n"+line
# try to execute the input file
......@@ -370,7 +365,7 @@ def get_proj_reading(preanalisi=0):
file_proj_num_list = numpy.array(seqnum_list )* P.FILE_INTERVAL+ P.NUM_FIRST_IMAGE
if P.FILE_PREFIX[-3:]==".h5":
if P.FILE_PREFIX[-3:]==".h5" or P.FILE_PREFIX[-4:]==".nxs":
# the file type is determined by the postfix
projtype = "h5"
file_list = [P.FILE_PREFIX]*len(file_proj_num_list)
......@@ -437,15 +432,20 @@ def get_proj_reading(preanalisi=0):
else: posb=posa
coeffb = i_FF*1.0/FF_FILE_INTERVAL
coeffa = 1.0-coeffb
print(" FILE PREFIX " , P.FILE_PREFIX )
if P.FILE_PREFIX[-3:]==".h5":
if P.FILE_PREFIX[-3:] == ".h5" or P.FILE_PREFIX[-4:] == ".nxs":
fftype="h5"
ff_file_list =[P.FF_PREFIX]
if isinstance( P.FF_PREFIX, tuple) or isinstance( P.FF_PREFIX, list):
ff_file_list = list(P.FF_PREFIX)
else:
ff_file_list = [P.FF_PREFIX]
if posb==posa:
ff_indexes_coefficients.append([posa,-1, 1.0 , 0.0 , (i_pro-i_FF) , -1] )
else:
ff_indexes_coefficients.append([posa,posb, coeffa , coeffb , (i_pro-i_FF) , (i_pro-i_FF) + FF_FILE_INTERVAL ] )
else:
fftype="edf"
edfargs = {"FILE_PREFIX":P.FF_PREFIX, "LENGTH_OF_NUMERICAL_PART":P.FF_LENGTH_OF_NUMERICAL_PART,
......@@ -879,6 +879,16 @@ class Parameters:
"""
REFERENCES_HOTPIXEL_THRESHOLD = 0
"""
if different than zero, then dark and references will be filtered for hotpixel
exceding their 3x3 median by more than such value .
This option is effective for for darks, while for refs it is effective only on hdf5 (or nxs)
files when
they contains multiple slices for the same reference. It is applied after median
"""
CURRENT_NAME = ""
"""
if this variable is set then the current is read in the header of each radiography
......@@ -1674,12 +1684,6 @@ Units are pixels per projection.
A good value is 0.01.
"""
CG_MU = 1e-5
"""
</