Skip to content

Commit

Permalink
vector and project2 are corrected
Browse files Browse the repository at this point in the history
  • Loading branch information
@dgursoy
dgursoy committed Apr 18, 2018
1 parent 7c2dd5f commit 2ea7cec
Show file tree
Hide file tree
Showing 7 changed files with 79 additions and 120 deletions.
26 changes: 15 additions & 11 deletions src/project.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,14 +144,10 @@ project(

void
project2(
const float *obj, int oy, int ox, int oz,
const float *objx, const float *objy, int oy, int ox, int oz,
float *data, int dy, int dt, int dx,
const float *center, const float *theta, int axis)
const float *center, const float *theta)
{
// printf ("oy=%d, ox=%d, oz=%d\n", oy, ox, oz);
// printf ("dy=%d, dt=%d, dx=%d\n", dy, dt, dx);
// printf ("axis=%d\n", axis);

float *gridx = (float *)malloc((ox+1)*sizeof(float));
float *gridy = (float *)malloc((oz+1)*sizeof(float));
float *coordx = (float *)malloc((oz+1)*sizeof(float));
Expand All @@ -174,6 +170,7 @@ project2(
int s, p, d;
int quadrant;
float theta_p, sin_p, cos_p;
float srcx, srcy, detx, dety, dv, vx, vy;
float mov, xi, yi;
int asize, bsize, csize;

Expand All @@ -193,9 +190,20 @@ project2(

for (d=0; d<dx; d++)
{

// Calculate coordinates
xi = -ox-oz;
yi = (1-dx)/2.0+d+mov;

srcx = xi*cos_p-yi*sin_p;
srcy = xi*sin_p+yi*cos_p;
detx = -xi*cos_p-yi*sin_p;
dety = -xi*sin_p+yi*cos_p;

dv = sqrt(pow(srcx-detx, 2)+pow(srcy-dety, 2));
vx = (srcx-detx)/dv;
vy = (srcy-dety)/dv;

calc_coords(
ox, oz, xi, yi, sin_p, cos_p, gridx, gridy,
coordx, coordy);
Expand Down Expand Up @@ -225,7 +233,7 @@ project2(
{
// Calculate simdata
calc_simdata2(s, p, d, ox, oz, dt, dx,
csize, indx, indy, dist, obj, axis,
csize, indx, indy, dist, vx, vy, objx, objy,
data); // Output: simulated data
}
}
Expand Down Expand Up @@ -254,10 +262,6 @@ project3(
float *data, int dy, int dt, int dx,
const float *center, const float *theta, int axis)
{
// printf ("oy=%d, ox=%d, oz=%d\n", oy, ox, oz);
// printf ("dy=%d, dt=%d, dx=%d\n", dy, dt, dx);
// printf ("axis=%d\n", axis);

float *gridx = (float *)malloc((ox+1)*sizeof(float));
float *gridy = (float *)malloc((oz+1)*sizeof(float));
float *coordx = (float *)malloc((oz+1)*sizeof(float));
Expand Down
46 changes: 5 additions & 41 deletions src/utils.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,30 +295,14 @@ calc_simdata2(
int ry, int rz,
int dt, int dx,
int csize, const int *indx, const int *indy, const float *dist,
const float *model, int axis, float *simdata)
float vx, float vy,
const float *modelx, const float *modely, float *simdata)
{
int n;

if (axis==0)
{
for (n=0; n<csize-1; n++)
{
simdata[d + p*dx + s*dt*dx] += model[indy[n] + indx[n]*rz + s*ry*rz]*dist[n];
}
}
else if (axis==1)
{
for (n=0; n<csize-1; n++)
{
simdata[d + p*dx + s*dt*dx] += model[s + indx[n]*rz + indy[n]*ry*rz]*dist[n];
}
}
else if (axis==2)
for (n=0; n<csize-1; n++)
{
for (n=0; n<csize-1; n++)
{
simdata[d + p*dx + s*dt*dx] += model[indx[n] + s*rz + indy[n]*ry*rz]*dist[n];
}
simdata[d + p*dx + s*dt*dx] += (modelx[indy[n] + indx[n]*rz + s*ry*rz] * vx + modely[indy[n] + indx[n]*rz + s*ry*rz] * vy) * dist[n];
}
}

Expand Down Expand Up @@ -355,24 +339,4 @@ calc_simdata3(
simdata[d + p*dx + s*dt*dx] += (modelx[indx[n] + s*rz + indy[n]*ry*rz] * vx + modelz[indx[n] + s*rz + indy[n]*ry*rz] * vy) * dist[n];
}
}
}

// void
// calc_simdata2(
// int s, int p, int d,
// int ry, int rz,
// int dt, int dx,
// int csize, const int *indi, const float *dist,
// float vx, float vy,
// const float *modelx, const float *modely,
// float *simdata)
// {
// int n;

// int index_model = s*ry*rz;
// int index_data = d+p*dx+s*dt*dx;
// for (n=0; n<csize-1; n++)
// {
// simdata[index_data] += (modelx[indi[n]+index_model] * vx + modely[indi[n]+index_model] * vy) * dist[n];
// }
// }
}
31 changes: 9 additions & 22 deletions src/utils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,8 @@ project(

void DLL
project2(
const float *obj,
const float *objx,
const float *objy,
int oy,
int ox,
int oz,
Expand All @@ -89,8 +90,7 @@ project2(
int dt,
int dx,
const float *center,
const float *theta,
int axis);
const float *theta);

void DLL
project3(
Expand Down Expand Up @@ -262,11 +262,11 @@ vector(
int dx,
const float *center,
const float *theta,
float *recon,
float *recon1,
float *recon2,
int ngridx,
int ngridy,
int num_iter,
int axis);
int num_iter);

void DLL
vector2(
Expand Down Expand Up @@ -386,8 +386,9 @@ calc_simdata2(
const int *indx,
const int *indy,
const float *dist,
const float *model,
int axis,
float vx, float vy,
const float *modelx,
const float *modely,
float *simdata);

void DLL
Expand All @@ -406,18 +407,4 @@ calc_simdata3(
int axis,
float *simdata);

// void DLL
// calc_simdata2(
// int s, int p, int d,
// int ngridx, int ngridy,
// int dt, int dx,
// int csize,
// const int *indi,
// const float *dist,
// float vx,
// float vy,
// const float *modelx,
// const float *modely,
// float *simdata);

#endif
54 changes: 28 additions & 26 deletions src/vector.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ void
vector(
const float *data, int dy, int dt, int dx,
const float *center, const float *theta,
float *recon, int ngridx, int ngridy, int num_iter, int axis)
float *recon1, float *recon2, int ngridx, int ngridy, int num_iter)
{
float *gridx = (float *)malloc((ngridx+1)*sizeof(float));
float *gridy = (float *)malloc((ngridy+1)*sizeof(float));
Expand Down Expand Up @@ -76,9 +76,10 @@ vector(
float *simdata;
float upd;
int ind_data;
float srcx, srcy, detx, dety, dv, vx, vy;
float *sum_dist;
float sum_dist2;
float *update;
float *update1, *update2;

for (i=0; i<num_iter; i++)
{
Expand All @@ -91,7 +92,8 @@ vector(
&mov, gridx, gridy); // Outputs: mov, gridx, gridy

sum_dist = (float *)calloc((ngridx*ngridy), sizeof(float));
update = (float *)calloc((ngridx*ngridy), sizeof(float));
update1 = (float *)calloc((ngridx*ngridy), sizeof(float));
update2 = (float *)calloc((ngridx*ngridy), sizeof(float));

// For each projection angle
for (p=0; p<dt; p++)
Expand All @@ -110,36 +112,45 @@ vector(
// Calculate coordinates
xi = -ngridx-ngridy;
yi = (1-dx)/2.0+d+mov;

srcx = xi*cos_p-yi*sin_p;
srcy = xi*sin_p+yi*cos_p;
detx = -xi*cos_p-yi*sin_p;
dety = -xi*sin_p+yi*cos_p;

dv = sqrt(pow(srcx-detx, 2)+pow(srcy-dety, 2));
vx = (srcx-detx)/dv;
vy = (srcy-dety)/dv;

calc_coords(
ngridx, ngridy, xi, yi, sin_p, cos_p, gridx, gridy,
ngridx, ngridy, xi, yi, sin_p, cos_p, gridx, gridy,
coordx, coordy);

// Merge the (coordx, gridy) and (gridx, coordy)
trim_coords(
ngridx, ngridy, coordx, coordy, gridx, gridy,
ngridx, ngridy, coordx, coordy, gridx, gridy,
&asize, ax, ay, &bsize, bx, by);

// Sort the array of intersection points (ax, ay) and
// (bx, by). The new sorted intersection points are
// stored in (coorx, coory). Total number of points
// are csize.
sort_intersections(
quadrant, asize, ax, ay, bsize, bx, by,
quadrant, asize, ax, ay, bsize, bx, by,
&csize, coorx, coory);

// Calculate the distances (dist) between the
// intersection points (coorx, coory). Find the
// Calculate the distances (dist) between the
// intersection points (coorx, coory). Find the
// indices of the pixels on the reconstruction grid.
calc_dist2(
ngridx, ngridy, csize, coorx, coory,
ngridx, ngridy, csize, coorx, coory,
indx, indy, dist);

// Calculate simdata
calc_simdata2(s, p, d, ngridx, ngridy, dt, dx,
csize, indx, indy, dist, recon, axis,
csize, indx, indy, dist, vx, vy, recon1, recon2,
simdata); // Output: simdata


// Calculate dist*dist
sum_dist2 = 0.0;
for (n=0; n<csize-1; n++)
Expand All @@ -151,36 +162,27 @@ vector(
// Update
if (sum_dist2 != 0.0)
{
// ind_data = d + p*dx + s*dt*dx;
ind_data = d + p*dx + s*dt*dx;
upd = (data[ind_data]-simdata[ind_data])/sum_dist2;
for (n=0; n<csize-1; n++)
{
update[indy[n] + indx[n]*ngridy] += upd*dist[n];
update1[indy[n] + indx[n]*ngridy] += upd*dist[n]*vx;
update2[indy[n] + indx[n]*ngridy] += upd*dist[n]*vy;
}
}
}
}

for (m = 0; m < ngridx; m++) {
for (n = 0; n < ngridy; n++) {
if (axis == 0)
{
recon[n + m*ngridy + s*ngridx*ngridy] += update[n + m*ngridy]/sum_dist[n + m*ngridy];
}
else if (axis == 1)
{
recon[s + m*ngridy + n*ngridx*ngridy] += update[n + m*ngridy]/sum_dist[n + m*ngridy];
}
else if (axis == 2)
{
recon[m + s*ngridy + n*ngridx*ngridy] += update[n + m*ngridy]/sum_dist[n + m*ngridy];
}
recon1[n + m*ngridy + s*ngridx*ngridy] += update1[n + m*ngridy]/sum_dist[n + m*ngridy];
recon2[n + m*ngridy + s*ngridx*ngridy] += update1[n + m*ngridy]/sum_dist[n + m*ngridy];
}
}

free(sum_dist);
free(update);
free(update1);
free(update2);
}

free(simdata);
Expand Down
11 changes: 6 additions & 5 deletions tomopy/recon/vector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,21 +73,22 @@
__all__ = ['vector', 'vector2', 'vector3']


def vector(tomo, theta, center=None, num_iter=1, axis=0):
def vector(tomo, theta, center=None, num_iter=1):
tomo = dtype.as_float32(tomo)
theta = dtype.as_float32(theta)

# Initialize tomography data.
tomo = init_tomo(tomo, sinogram_order=False, sharedmem=False)

recon_shape = (tomo.shape[0], tomo.shape[2], tomo.shape[2])
recon = np.zeros(recon_shape, dtype=np.float32)
recon1 = np.zeros(recon_shape, dtype=np.float32)
recon2 = np.zeros(recon_shape, dtype=np.float32)

center_arr = get_center(tomo.shape, center)

extern.c_vector(tomo, center_arr, recon, theta,
num_gridx=tomo.shape[2], num_gridy=tomo.shape[2], num_iter=num_iter, axis=axis)
return recon
extern.c_vector(tomo, center_arr, recon1, recon2, theta,
num_gridx=tomo.shape[2], num_gridy=tomo.shape[2], num_iter=num_iter)
return recon1, recon2


def vector2(tomo1, tomo2, theta1, theta2, center1=None, center2=None, num_iter=1, axis1=1, axis2=2):
Expand Down
9 changes: 5 additions & 4 deletions tomopy/sim/project.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,7 @@ def project(


def project2(
obj, theta, center=None, emission=True, pad=True,
objx, objy, theta, center=None, emission=True, pad=True,
sinogram_order=False, axis=0, ncore=None, nchunk=None):
"""
Project x-rays through a given 3D object.
Expand Down Expand Up @@ -277,11 +277,12 @@ def project2(
ndarray
3D tomographic data.
"""
obj = dtype.as_float32(obj)
objx = dtype.as_float32(objx)
objy = dtype.as_float32(objy)
theta = dtype.as_float32(theta)

# Estimate data dimensions.
oy, ox, oz = obj.shape
oy, ox, oz = objx.shape
dt = theta.size
dy = oy
if pad is True:
Expand All @@ -293,7 +294,7 @@ def project2(
tomo[:] = 0.0
center = get_center(shape, center)

extern.c_project2(obj, center, tomo, theta, axis)
extern.c_project2(objx, objy, center, tomo, theta)

# NOTE: returns sinogram order with emmission=True
if not emission:
Expand Down
Loading

0 comments on commit 2ea7cec

Please sign in to comment.