Commit e002ba91 authored by myron's avatar myron

margin

parent ddfda261
......@@ -2413,8 +2413,14 @@ void CCspace_precalculations( CCspace * self, int ncpus ) {
self->gpu_context->FAN_FACTOR = 0;
}
self->gpu_context->DETECTOR_DISTANCE = self->params.DETECTOR_DISTANCE ;
self->gpu_context->SOURCE_X = self->params.SOURCE_X ;
self->gpu_context->SOURCE_Y = self->params.SOURCE_Y ;
self->gpu_context->SOURCE_X = self->params.SOURCE_X ;
self->gpu_context->SOURCE_Y = self->params.SOURCE_Y ;
self->gpu_context->tilt_psi = self->params.DECT_PSI *M_PI/180.0 ;
self->gpu_context->tilt_tilt = self->params.DECT_TILT *M_PI/180/0 ;
self->gpu_context->CONICITY_MARGIN_UP = self->params.CONICITY_MARGIN_UP ;
self->gpu_context->CONICITY_MARGIN_DOWN = self->params.CONICITY_MARGIN_DOWN ;
......
......@@ -416,7 +416,9 @@ typedef struct {
float DETECTOR_DISTANCE ;
float SOURCE_X ;
float SOURCE_Y ;
float DECT_PSI;
float DECT_TILT;
float DZPERPROJ;
float DXPERPROJ;
......@@ -963,7 +965,7 @@ struct Gpu_Context_struct {
int dimrecx, dimrecy;
float * d_SLICE ;
float * d_cos_s, * d_sin_s, * d_axis_s ;
// set by user
int num_x ;
int num_y ;
......@@ -1017,6 +1019,10 @@ struct Gpu_Context_struct {
float DETECTOR_DISTANCE ;
float SOURCE_X ;
float SOURCE_Y ;
float tilt_psi;
float tilt_tilt;
int * CONICITY_MARGIN_UP ;
int * CONICITY_MARGIN_DOWN ;
......
......@@ -638,8 +638,10 @@ Cspace_new( PyObject *self, PyObject *args)
res->myCCspace->params.SOURCE_DISTANCE = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"SOURCE_DISTANCE" ) );
res->myCCspace->params.DETECTOR_DISTANCE = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"DETECTOR_DISTANCE" ) );
res->myCCspace->params.SOURCE_X = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"SOURCE_X" ) );
// res->myCCspace->params.SOURCE_X = res->myCCspace->params.ROTATION_AXIS_POSITION ;
res->myCCspace->params.SOURCE_Y = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"SOURCE_Y" ) );
res->myCCspace->params.DECT_PSI = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"DECT_PSI" ) );
res->myCCspace->params.DECT_TILT = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"DECT_TILT" ) );
res->myCCspace->params.DZPERPROJ = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"DZPERPROJ" ) );
res->myCCspace->params.DXPERPROJ = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"DXPERPROJ" ) );
......
......@@ -1870,6 +1870,17 @@ typedef struct {
cudaStream_t stream[2];
} Gpu_Streams ;
typedef struct {
float sin_psi;
float cos_psi;
float sin_tilt;
float cos_tilt;
} DectTilt ;
#define USE_R2C_TRANSFORM 1
#define BLOCK_SIZE 16
......@@ -7178,8 +7189,13 @@ __global__ static void multi_pro_gputomo_conicity_kernel(float *d_SINO, // da
float v2z,
float v_size,
float SOURCE_X,
float SOURCE_Z) {
float SOURCE_Z,
DectTilt dectTilt
) {
float cos_angle, sin_angle;
// l'asse e' dato a partire dall' angolo primo pixel, ma i centri dei pixel sono nel mezzo
// float acorr05;
......@@ -7190,6 +7206,23 @@ __global__ static void multi_pro_gputomo_conicity_kernel(float *d_SINO, // da
projection = (blockDim.y * blockIdx.y + threadIdx.y) ; // num_proj
iz = (blockDim.z * blockIdx.z + threadIdx.z )* PROCONO_GPU_MULT + data_start; // nslices_data
float fix, fiz;
if(dectTilt.sin_psi!=0.0f) {
fix = SOURCE_X + dectTilt.cos_psi*(ix-SOURCE_X) - dectTilt.sin_psi * (iz-SOURCE_Z) ;
fiz = SOURCE_Z + dectTilt.sin_psi*(ix-SOURCE_X) + dectTilt.cos_psi * (iz-SOURCE_Z) ;
} else {
fix = ix;
fiz = iz;
}
if(dectTilt.sin_tilt!=0.0f) {
float alpha = ( 1.0e-6 * v_size/v2x) / ( source_distance+detector_distance) ;
fiz = SOURCE_Z + ( -dectTilt.sin_tilt* (fiz-SOURCE_Z)*alpha + dectTilt.cos_tilt ) * (fiz-SOURCE_Z) ;
fix = SOURCE_X + ( -dectTilt.sin_tilt* (fiz-SOURCE_Z)*alpha + 1.0f ) * (fix-SOURCE_Z) ;
} else {
}
// float magni, px, pz;
float magnicenter = (source_distance+detector_distance)/ (source_distance ) ;
......@@ -7235,8 +7268,8 @@ __global__ static void multi_pro_gputomo_conicity_kernel(float *d_SINO, // da
float Area_fz=1.0f ;
{
float L = (source_distance+detector_distance)/(1.0e-6);
float H = (ix-SOURCE_X) * v_size/v2x;
float V = (iz-SOURCE_Z)* v_size/v2x;
float H = (fix-SOURCE_X) * v_size/v2x;
float V = (fiz-SOURCE_Z)* v_size/v2x;
Area_fz = sqrt( (L*L+H*H+V*V)/(L*L+H*H) ) ;
}
......@@ -7248,15 +7281,15 @@ __global__ static void multi_pro_gputomo_conicity_kernel(float *d_SINO, // da
for(ivx=0; ivx<num_x; ivx++) {
Ddeno = Ddeno0 + ivx*sin_angle * v_size *1.0e-6;
d_F_x = cos_angle*v_size*1.0e-6/Dnumex;
F = (sin_angle-(SOURCE_X-ix)*d_F_x);
F = (sin_angle-(SOURCE_X-fix)*d_F_x);
md_fvyF_x = (Ddeno)/Dnumex ;
fvy =( ( Xcorner+ivx*cos_angle ) + (SOURCE_X-ix)*md_fvyF_x );
fvy =( ( Xcorner+ivx*cos_angle ) + (SOURCE_X-fix)*md_fvyF_x );
fvy = fvy / F ;
d_fvy_x = (-md_fvyF_x-fvy*d_F_x) /F ;
Ddeno += fvy*cos_angle * v_size *1.0e-6;
d_z_z = Ddeno / Dnumez;
d_z_x = (iz-SOURCE_Z)*d_fvy_x*cos_angle* v_size *1.0e-6 / Dnumez;
fvz = (iz-SOURCE_Z)*d_z_z -Nfirstslice;
d_z_x = (fiz-SOURCE_Z)*d_fvy_x*cos_angle* v_size *1.0e-6 / Dnumez;
fvz = (fiz-SOURCE_Z)*d_z_z -Nfirstslice;
for(jmult=0; jmult < PROCONO_GPU_MULT; jmult++) {
fvyb = fvy + jmult * d_fvy_x ;
for(imult=0; imult < PROCONO_GPU_MULT; imult++) {
......@@ -7273,15 +7306,15 @@ __global__ static void multi_pro_gputomo_conicity_kernel(float *d_SINO, // da
for(ivy=0; ivy<num_y; ivy++) {
Ddeno = Ddeno0 + ivy*cos_angle * v_size *1.0e-6;
d_F_x = -sin_angle*v_size*1.0e-6/Dnumex;
F = (cos_angle+(SOURCE_X-ix)*(-d_F_x));
F = (cos_angle+(SOURCE_X-fix)*(-d_F_x));
d_fvxF_x = (Ddeno)/Dnumex ;
fvx =( -( Xcorner-ivy*sin_angle ) + (ix-SOURCE_X)*d_fvxF_x );
fvx =( -( Xcorner-ivy*sin_angle ) + (fix-SOURCE_X)*d_fvxF_x );
fvx = fvx / F ;
d_fvx_x = ( d_fvxF_x-fvx*d_F_x) /F ;
Ddeno += fvx*sin_angle * v_size *1.0e-6;
d_z_z = Ddeno / Dnumez;
d_z_x = (iz-SOURCE_Z)*d_fvx_x*sin_angle* v_size *1.0e-6 / Dnumez;
fvz = (iz-SOURCE_Z)*d_z_z -Nfirstslice;
d_z_x = (fiz-SOURCE_Z)*d_fvx_x*sin_angle* v_size *1.0e-6 / Dnumez;
fvz = (fiz-SOURCE_Z)*d_z_z -Nfirstslice;
for(jmult=0; jmult < PROCONO_GPU_MULT; jmult++) {
fvxb = fvx + jmult * d_fvx_x ;
for(imult=0; imult < PROCONO_GPU_MULT; imult++) {
......@@ -7327,7 +7360,8 @@ __global__ static void gputomo_conicity_kernel(float *d_SLICE, // da allocare
float v2z,
float voxel_size,
float SOURCE_X,
float SOURCE_Z) {
float SOURCE_Z,
DectTilt dectTilt) {
......@@ -7420,10 +7454,31 @@ __global__ static void gputomo_conicity_kernel(float *d_SLICE, // da allocare
continue;
// Z = nslices_data-1;
}
res[iy* CONO_GPU_MULT +ix] += tex2D(texProjes,
px + 0.5f + ix*d_px_ix + iy* d_px_iy ,
Z + proje*nslices_data +0.5f
float X = px + ix*d_px_ix + iy* d_px_iy ;
float X1,Z1;
if(dectTilt.sin_tilt!=0.0f) {
float alpha = ( 1.0e-6 * voxel_size/v2x) / ( source_distance+detector_distance) ;
Z1 = SOURCE_Z + ( +dectTilt.sin_tilt* (Z-SOURCE_Z)*alpha + 1.0f/dectTilt.cos_tilt ) * (Z - SOURCE_Z) ;
X1 = SOURCE_X + ( +dectTilt.sin_tilt* (Z-SOURCE_Z)*alpha + 1.0f ) * (X - SOURCE_Z) ;
} else {
X1 = X ;
Z1 = Z;
}
float fix, fiz;
if(dectTilt.sin_psi!=0.0f) {
fix = SOURCE_X + dectTilt.cos_psi*(X1-SOURCE_X) + dectTilt.sin_psi * ( Z1 + data_start - SOURCE_Z ) ;
fiz = SOURCE_Z - dectTilt.sin_psi*(X1-SOURCE_X) + dectTilt.cos_psi * ( Z1 + data_start - SOURCE_Z ) - data_start ;
} else {
fix = X1;
fiz = Z1;
}
res[iy* CONO_GPU_MULT +ix] += tex2D(texProjes,
fix + 0.5f ,
fiz + proje*nslices_data +0.5f
);
}
}
......@@ -7437,7 +7492,26 @@ __global__ static void gputomo_conicity_kernel(float *d_SLICE, // da allocare
continue;
// Z = nslices_data-1;
}
res[0] += tex2D(texProjes, px + 0.5f , Z+ proje*nslices_data + 0.5f ) ; ;
float X = px;
float X1,Z1;
if(dectTilt.sin_tilt!=0.0f) {
float alpha = ( 1.0e-6 * voxel_size/v2x) / ( source_distance+detector_distance) ;
Z1 = SOURCE_Z + ( +dectTilt.sin_tilt* (Z-SOURCE_Z)*alpha + 1.0f/dectTilt.cos_tilt ) * (Z - SOURCE_Z) ;
X1 = SOURCE_X + ( +dectTilt.sin_tilt* (Z-SOURCE_Z)*alpha + 1.0f ) * (X - SOURCE_Z) ;
} else {
X1 = X ;
Z1 = Z;
}
float fix, fiz;
if(dectTilt.sin_psi!=0.0f) {
fix = SOURCE_X + dectTilt.cos_psi*(X1-SOURCE_X) + dectTilt.sin_psi * ( Z1 + data_start - SOURCE_Z ) ;
fiz = SOURCE_Z - dectTilt.sin_psi*(X1-SOURCE_X) + dectTilt.cos_psi * ( Z1 + data_start - SOURCE_Z ) - data_start ;
} else {
fix = X1;
fiz = Z1;
}
res[0] += tex2D(texProjes, fix + 0.5f , fiz + proje*nslices_data + 0.5f ) ; ;
}
}
if(CONO_GPU_MULT>1 ) {
......@@ -7503,6 +7577,8 @@ int gpu_main_conicity(Gpu_Context * self, float * SLICE, float *WORK_perp
texProjes.normalized = false;
CUDA_SAFE_CALL( cudaBindTextureToArray(texProjes,a_Proje_voidptr) );
DectTilt dectTilt = { sin(self->tilt_psi), cos(self->tilt_psi), sin(self->tilt_tilt), cos(self->tilt_tilt) } ;
{
......@@ -7525,7 +7601,7 @@ int gpu_main_conicity(Gpu_Context * self, float * SLICE, float *WORK_perp
data_start, nslices_data,
source_distance, detector_distance,
v2x, v2z, voxel_size,
SOURCE_X, SOURCE_Z);
SOURCE_X, SOURCE_Z, dectTilt);
}
}
......@@ -7613,7 +7689,9 @@ int pro_gpu_main_conicity(Gpu_Context * self, float * SLICE, float *SINO
nx = ( (int ) (self->num_bins*1.0/(pezzo_x * PROCONO_GPU_MULT) +0.999999));
ny = ( (int ) (self->nprojs_span*1.0/pezzo_y +0.999999));
nz = ( (int ) (nslices_data*1.0/(pezzo_z* PROCONO_GPU_MULT) +0.999999));
DectTilt dectTilt = { sin(self->tilt_psi), cos(self->tilt_psi), sin(self->tilt_tilt), cos(self->tilt_tilt) } ;
float *d_SINO;
CUDA_SAFE_CALL(cudaMalloc( &(d_SINO) , sizeof(float) *nslices_data *self->nprojs_span* self->num_bins ));
CUDA_SAFE_CALL(cudaMemset( d_SINO,0, sizeof(float) * nslices_data *self->nprojs_span* self->num_bins ));
......@@ -7634,7 +7712,8 @@ int pro_gpu_main_conicity(Gpu_Context * self, float * SLICE, float *SINO
data_start,
source_distance, detector_distance,
v2x, v2z, voxel_size,
SOURCE_X, SOURCE_Z);
SOURCE_X, SOURCE_Z,
dectTilt);
}
CUDA_SAFE_CALL( cudaMemcpy( SINO, d_SINO, sizeof(float) * nslices_data *self->nprojs_span* self->num_bins , cudaMemcpyDeviceToHost) );
......
......@@ -1885,6 +1885,9 @@ Units are pixels per projection.
"""
DECT_PSI = 0.0
DECT_TILT = 0.0
STEAM_INVERTER = 0
"""
When STEAM INVERTER=1 PyHST goes from volume to projections.
......@@ -2284,7 +2287,7 @@ if not Parameters.DO_V3D_UNSHARP and Parameters.NSLICESATONCE is None and (sys.
if P.CONICITY:
f_tmp = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE)) # ??? /P.IMAGE_PIXEL_SIZE_2
VOXEL_SIZE = P.IMAGE_PIXEL_SIZE_1/f_tmp
Fact = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))*VOXEL_SIZE*1.0/P.IMAGE_PIXEL_SIZE_2
Fact = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))*VOXEL_SIZE*1.0/P.IMAGE_PIXEL_SIZE_2 +(abs(P.DECT_PSI)+abs(P.DECT_PSI) )*3.2/180.0*2
DR1 = (P.END_VOXEL_1 - P.START_VOXEL_1)/2.0
DR2 = (P.END_VOXEL_2 - P.START_VOXEL_2)/2.0
RADIUS = math.sqrt( DR1*DR1+DR2*DR2 )
......@@ -2726,7 +2729,7 @@ if (sys.argv[0][-12:]!="sphinx-build") and not P.DO_V3D_UNSHARP :
f_tmp = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE)) # ??? /P.IMAGE_PIXEL_SIZE_2
P.VOXEL_SIZE = P.IMAGE_PIXEL_SIZE_1/f_tmp
Fact = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))*P.VOXEL_SIZE*1.0/P.IMAGE_PIXEL_SIZE_2
Fact = (P.SOURCE_DISTANCE+P.DETECTOR_DISTANCE)/(float(P.SOURCE_DISTANCE))*P.VOXEL_SIZE*1.0/P.IMAGE_PIXEL_SIZE_2 +(abs(P.DECT_PSI)+abs(P.DECT_PSI) )*3.2/180.0*2
DR1 = (P.END_VOXEL_1 - P.START_VOXEL_1)/2.0
DR2 = (P.END_VOXEL_2 - P.START_VOXEL_2)/2.0
......
......@@ -327,7 +327,7 @@ def do_pyhst():
cython_c_ext = ".pyx"
else:
cython_c_ext = ".cpp"
from distutils.command.build_ext import build_ext
from distutils.command.build_ext import build_ext, Extension
# We subclass the build_ext class in order to handle compiler flags
......
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