Skip to content

Commit

Permalink
MNT: structures removed
Browse files Browse the repository at this point in the history
  • Loading branch information
@dgursoy
dgursoy committed Apr 13, 2015
1 parent c41a060 commit 6591c14
Show file tree
Hide file tree
Showing 16 changed files with 1,214 additions and 946 deletions.
637 changes: 476 additions & 161 deletions tomopy/recon.py
View file

Large diffs are not rendered by default.

60 changes: 26 additions & 34 deletions tomopy/src/art.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,48 +46,40 @@

void
art(
float *data, data_pars *dpars,
float *recon, recon_pars *rpars)
float *data, int dx, int dy, int dz, float center, float *theta,
float *recon, int ngridx, int ngridy, int num_iter)
{
int dx, dy, dz, ry, rz;

dx = dpars->dx;
dy = dpars->dy;
dz = dpars->dz;
ry = rpars->ry;
rz = rpars->rz;

float *gridx = (float *)malloc((ry+1)*sizeof(float));
float *gridy = (float *)malloc((rz+1)*sizeof(float));
float *coordx = (float *)malloc((rz+1)*sizeof(float));
float *coordy = (float *)malloc((ry+1)*sizeof(float));
float *ax = (float *)malloc((ry+rz)*sizeof(float));
float *ay = (float *)malloc((ry+rz)*sizeof(float));
float *bx = (float *)malloc((ry+rz)*sizeof(float));
float *by = (float *)malloc((ry+rz)*sizeof(float));
float *coorx = (float *)malloc((ry+rz)*sizeof(float));
float *coory = (float *)malloc((ry+rz)*sizeof(float));
float *dist = (float *)malloc((ry+rz)*sizeof(float));
int *indi = (int *)malloc((ry+rz)*sizeof(int));
float *gridx = (float *)malloc((ngridx+1)*sizeof(float));
float *gridy = (float *)malloc((ngridy+1)*sizeof(float));
float *coordx = (float *)malloc((ngridy+1)*sizeof(float));
float *coordy = (float *)malloc((ngridx+1)*sizeof(float));
float *ax = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *ay = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *bx = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *by = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *coorx = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *coory = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *dist = (float *)malloc((ngridx+ngridy)*sizeof(float));
int *indi = (int *)malloc((ngridx+ngridy)*sizeof(int));

assert(coordx != NULL && coordy != NULL &&
ax != NULL && ay != NULL && by != NULL && bx != NULL &&
coorx != NULL && coory != NULL && dist != NULL && indi != NULL);

int s, p, d, i, n;
int quadrant;
float proj_angle, sin_p, cos_p;
float theta_p, sin_p, cos_p;
float mov, xi, yi;
int asize, bsize, csize;
float* simdata;
float upd;
int ind_data, ind_recon;


preprocessing(ry, rz, dz, dpars->center,
preprocessing(ngridx, ngridy, dz, center,
&mov, gridx, gridy); // Outputs: mov, gridx, gridy

for (i=0; i<rpars->num_iter; i++)
for (i=0; i<num_iter; i++)
{
printf("ART iteration : %i\n", i+1);

Expand All @@ -102,10 +94,10 @@ art(
// Calculate the sin and cos values
// of the projection angle and find
// at which quadrant on the cartesian grid.
proj_angle = fmod(dpars->proj_angle[p], 2*M_PI);
quadrant = calc_quadrant(proj_angle);
sin_p = sinf(proj_angle);
cos_p = cosf(proj_angle);
theta_p = fmod(theta[p], 2*M_PI);
quadrant = calc_quadrant(theta_p);
sin_p = sinf(theta_p);
cos_p = cosf(theta_p);

// For each detector pixel
for (d=0; d<dz; d++)
Expand All @@ -114,12 +106,12 @@ art(
xi = -1e6;
yi = -(dz-1)/2.0+d+mov;
calc_coords(
ry, rz, 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(
ry, rz, 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
Expand All @@ -134,11 +126,11 @@ art(
// intersection points (coorx, coory). Find the
// indices of the pixels on the reconstruction grid.
calc_dist(
ry, rz, csize, coorx, coory,
ngridx, ngridy, csize, coorx, coory,
indi, dist);

// Calculate simdata
calc_simdata(p, s, d, ry, rz, dy, dz,
calc_simdata(p, s, d, ngridx, ngridy, dy, dz,
csize, indi, dist, recon,
simdata); // Output: simdata

Expand All @@ -151,7 +143,7 @@ art(
}

// Update
ind_recon = s*ry*rz;
ind_recon = s*ngridx*ngridy;
ind_data = d+s*dz+p*dy*dz;
if (sum_dist2 != 0.0)
{
Expand Down
77 changes: 35 additions & 42 deletions tomopy/src/bart.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,37 +46,30 @@

void
bart(
float *data, data_pars *dpars,
float *recon, recon_pars *rpars)
float *data, int dx, int dy, int dz, float center, float *theta,
float *recon, int ngridx, int ngridy, int num_iter,
int num_block, float *ind_block)
{
int dx, dy, dz, ry, rz;

dx = dpars->dx;
dy = dpars->dy;
dz = dpars->dz;
ry = rpars->ry;
rz = rpars->rz;

float *gridx = (float *)malloc((ry+1)*sizeof(float));
float *gridy = (float *)malloc((rz+1)*sizeof(float));
float *coordx = (float *)malloc((rz+1)*sizeof(float));
float *coordy = (float *)malloc((ry+1)*sizeof(float));
float *ax = (float *)malloc((ry+rz)*sizeof(float));
float *ay = (float *)malloc((ry+rz)*sizeof(float));
float *bx = (float *)malloc((ry+rz)*sizeof(float));
float *by = (float *)malloc((ry+rz)*sizeof(float));
float *coorx = (float *)malloc((ry+rz)*sizeof(float));
float *coory = (float *)malloc((ry+rz)*sizeof(float));
float *dist = (float *)malloc((ry+rz)*sizeof(float));
int *indi = (int *)malloc((ry+rz)*sizeof(int));
float *gridx = (float *)malloc((ngridx+1)*sizeof(float));
float *gridy = (float *)malloc((ngridy+1)*sizeof(float));
float *coordx = (float *)malloc((ngridy+1)*sizeof(float));
float *coordy = (float *)malloc((ngridx+1)*sizeof(float));
float *ax = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *ay = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *bx = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *by = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *coorx = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *coory = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *dist = (float *)malloc((ngridx+ngridy)*sizeof(float));
int *indi = (int *)malloc((ngridx+ngridy)*sizeof(int));

assert(coordx != NULL && coordy != NULL &&
ax != NULL && ay != NULL && by != NULL && bx != NULL &&
coorx != NULL && coory != NULL && dist != NULL && indi != NULL);

int s, q, p, d, i, m, n, os;
int quadrant;
float proj_angle, sin_p, cos_p;
float theta_p, sin_p, cos_p;
float mov, xi, yi;
int asize, bsize, csize;
float *simdata;
Expand All @@ -87,14 +80,14 @@ bart(
float *update;
int subset_ind1, subset_ind2;

preprocessing(ry, rz, dz, dpars->center,
preprocessing(ngridx, ngridy, dz, center,
&mov, gridx, gridy); // Outputs: mov, gridx, gridy

for (i=0; i<rpars->num_iter; i++)
for (i=0; i<num_iter; i++)
{
printf("BART iteration : %i\n", i+1);

subset_ind1 = dx/rpars->num_block;
subset_ind1 = dx/num_block;
subset_ind2 = subset_ind1;

simdata = (float *)calloc((dx*dy*dz), sizeof(float));
Expand All @@ -103,28 +96,28 @@ bart(
for (s=0; s<dy; s++)
{
// For each ordered-subset num_subset
for (os=0; os<rpars->num_block+1; os++)
for (os=0; os<num_block+1; os++)
{
if (os == rpars->num_block)
if (os == num_block)
{
subset_ind2 = dx%rpars->num_block;
subset_ind2 = dx%num_block;
}

sum_dist = (float *)calloc((ry*rz), sizeof(float));
update = (float *)calloc((ry*rz), sizeof(float));
sum_dist = (float *)calloc((ngridx*ngridy), sizeof(float));
update = (float *)calloc((ngridx*ngridy), sizeof(float));

// For each projection angle
for (q=0; q<subset_ind2; q++)
{
p = rpars->ind_block[q+os*subset_ind1];
p = ind_block[q+os*subset_ind1];

// Calculate the sin and cos values
// of the projection angle and find
// at which quadrant on the cartesian grid.
proj_angle = fmod(dpars->proj_angle[p], 2*M_PI);
quadrant = calc_quadrant(proj_angle);
sin_p = sinf(proj_angle);
cos_p = cosf(proj_angle);
theta_p = fmod(theta[p], 2*M_PI);
quadrant = calc_quadrant(theta_p);
sin_p = sinf(theta_p);
cos_p = cosf(theta_p);

// For each detector pixel
for (d=0; d<dz; d++)
Expand All @@ -133,12 +126,12 @@ bart(
xi = -1e6;
yi = -(dz-1)/2.0+d+mov;
calc_coords(
ry, rz, 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(
ry, rz, 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
Expand All @@ -153,11 +146,11 @@ bart(
// intersection points (coorx, coory). Find the
// indices of the pixels on the reconstruction grid.
calc_dist(
ry, rz, csize, coorx, coory,
ngridx, ngridy, csize, coorx, coory,
indi, dist);

// Calculate simdata
calc_simdata(p, s, d, ry, rz, dy, dz,
calc_simdata(p, s, d, ngridx, ngridy, dy, dz,
csize, indi, dist, recon,
simdata); // Output: simdata

Expand All @@ -184,9 +177,9 @@ bart(
}

m = 0;
for (n = 0; n < ry*rz; n++) {
for (n = 0; n < ngridx*ngridy; n++) {
if (sum_dist[n] != 0.0) {
ind_recon = s*ry*rz;
ind_recon = s*ngridx*ngridy;
recon[m+ind_recon] += update[m]/sum_dist[n];
}
m++;
Expand Down
58 changes: 25 additions & 33 deletions tomopy/src/fbp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,44 +46,36 @@

void
fbp(
float *data, data_pars *dpars,
float *recon, recon_pars *rpars)
float *data, int dx, int dy, int dz, float center, float *theta,
float *recon, int ngridx, int ngridy)
{
int dx, dy, dz, ry, rz;

dx = dpars->dx;
dy = dpars->dy;
dz = dpars->dz;
ry = rpars->ry;
rz = rpars->rz;

float *gridx = (float *)malloc((ry+1)*sizeof(float));
float *gridy = (float *)malloc((rz+1)*sizeof(float));
float *coordx = (float *)malloc((rz+1)*sizeof(float));
float *coordy = (float *)malloc((ry+1)*sizeof(float));
float *ax = (float *)malloc((ry+rz)*sizeof(float));
float *ay = (float *)malloc((ry+rz)*sizeof(float));
float *bx = (float *)malloc((ry+rz)*sizeof(float));
float *by = (float *)malloc((ry+rz)*sizeof(float));
float *coorx = (float *)malloc((ry+rz)*sizeof(float));
float *coory = (float *)malloc((ry+rz)*sizeof(float));
float *dist = (float *)malloc((ry+rz)*sizeof(float));
int *indi = (int *)malloc((ry+rz)*sizeof(int));
float *gridx = (float *)malloc((ngridx+1)*sizeof(float));
float *gridy = (float *)malloc((ngridy+1)*sizeof(float));
float *coordx = (float *)malloc((ngridy+1)*sizeof(float));
float *coordy = (float *)malloc((ngridx+1)*sizeof(float));
float *ax = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *ay = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *bx = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *by = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *coorx = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *coory = (float *)malloc((ngridx+ngridy)*sizeof(float));
float *dist = (float *)malloc((ngridx+ngridy)*sizeof(float));
int *indi = (int *)malloc((ngridx+ngridy)*sizeof(int));

assert(coordx != NULL && coordy != NULL &&
ax != NULL && ay != NULL && by != NULL && bx != NULL &&
coorx != NULL && coory != NULL && dist != NULL && indi != NULL);

int s, p, d, n;
int quadrant;
float proj_angle, sin_p, cos_p;
float theta_p, sin_p, cos_p;
float mov, xi, yi;
int asize, bsize, csize;
float* simdata;
int ind_data, ind_recon;


preprocessing(ry, rz, dz, dpars->center,
preprocessing(ngridx, ngridy, dz, center,
&mov, gridx, gridy); // Outputs: mov, gridx, gridy

simdata = (float *)calloc((dx*dy*dz), sizeof(float));
Expand All @@ -97,10 +89,10 @@ fbp(
// Calculate the sin and cos values
// of the projection angle and find
// at which quadrant on the cartesian grid.
proj_angle = fmod(dpars->proj_angle[p], 2*M_PI);
quadrant = calc_quadrant(proj_angle);
sin_p = sinf(proj_angle);
cos_p = cosf(proj_angle);
theta_p = fmod(theta[p], 2*M_PI);
quadrant = calc_quadrant(theta_p);
sin_p = sinf(theta_p);
cos_p = cosf(theta_p);

// For each detector pixel
for (d=0; d<dz; d++)
Expand All @@ -109,12 +101,12 @@ fbp(
xi = -1e6;
yi = -(dz-1)/2.0+d+mov;
calc_coords(
ry, rz, 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(
ry, rz, 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
Expand All @@ -129,11 +121,11 @@ fbp(
// intersection points (coorx, coory). Find the
// indices of the pixels on the reconstruction grid.
calc_dist(
ry, rz, csize, coorx, coory,
ngridx, ngridy, csize, coorx, coory,
indi, dist);

// Calculate simdata
calc_simdata(p, s, d, ry, rz, dy, dz,
calc_simdata(p, s, d, ngridx, ngridy, dy, dz,
csize, indi, dist, recon,
simdata); // Output: simdata

Expand All @@ -148,7 +140,7 @@ fbp(
// Update
if (sum_dist2 != 0.0)
{
ind_recon = s*ry*rz;
ind_recon = s*ngridx*ngridy;
ind_data = d+s*dz+p*dy*dz;
for (n=0; n<csize-1; n++)
{
Expand Down

0 comments on commit 6591c14

Please sign in to comment.