vector tomography update

dgursoy · Apr 10, 2018 · bbf5df8 · bbf5df8
1 parent 3a920bd
commit bbf5df8
Show file tree

Hide file tree

Showing 19 changed files with 1,555 additions and 32 deletions.
diff --git a/VERSION b/VERSION
@@ -1 +1 @@
-1.1.2
+1.1.3
diff --git a/gridrec_butterworth.npy b/gridrec_butterworth.npy
diff --git a/gridrec_cosine.npy b/gridrec_cosine.npy
diff --git a/gridrec_hamming.npy b/gridrec_hamming.npy
diff --git a/gridrec_hann.npy b/gridrec_hann.npy
diff --git a/gridrec_none.npy b/gridrec_none.npy
diff --git a/gridrec_parzen.npy b/gridrec_parzen.npy
diff --git a/gridrec_ramlak.npy b/gridrec_ramlak.npy
diff --git a/gridrec_shepp.npy b/gridrec_shepp.npy
diff --git a/meta.yaml b/meta.yaml
@@ -1,6 +1,6 @@
 package:
   name: tomopy
-  version: '1.1.2'
+  version: '1.1.3'
 
 source:
   git_url: https://github.com/tomopy/tomopy.git

diff --git a/setup.py b/setup.py
@@ -80,6 +80,7 @@
         'src/sirt.c',
         'src/morph.c',
         'src/stripe.c',
+        'src/vector.c',
         'src/remove_ring.c'],
     define_macros=[('USE_MKL', None)] if use_mkl else [])
 

diff --git a/src/project.c b/src/project.c
@@ -141,12 +141,113 @@ project(
     free(indi);
 }
 
+
 void 
 project2(
-    const float *objx, const float *objy,
+    const float *obj, int oy, int ox, int oz, 
+    float *data, int dy, int dt, int dx,
+    const float *center, const float *theta, int 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));
+    float *coordy = (float *)malloc((ox+1)*sizeof(float));
+    float *ax = (float *)malloc((ox+oz+2)*sizeof(float));
+    float *ay = (float *)malloc((ox+oz+2)*sizeof(float));
+    float *bx = (float *)malloc((ox+oz+2)*sizeof(float));
+    float *by = (float *)malloc((ox+oz+2)*sizeof(float));
+    float *coorx = (float *)malloc((ox+oz+2)*sizeof(float));
+    float *coory = (float *)malloc((ox+oz+2)*sizeof(float));
+    float *dist = (float *)malloc((ox+oz+1)*sizeof(float));
+    int *indx = (int *)malloc((ox+oz+1)*sizeof(int));
+    int *indy = (int *)malloc((ox+oz+1)*sizeof(int));
+    int *indi = (int *)malloc((ox+oz+1)*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;
+    int quadrant;
+    float theta_p, sin_p, cos_p;
+    float mov, xi, yi;
+    int asize, bsize, csize;
+
+    preprocessing(ox, oz, dx, center[0], 
+    &mov, gridx, gridy); // Outputs: mov, gridx, gridy
+
+    // For each projection angle
+    for (p=0; p<dt; p++)
+    {
+        // Calculate the sin and cos values 
+        // of the projection angle and find
+        // at which quadrant on the cartesian grid.
+        theta_p = fmod(theta[p], 2*M_PI);
+        quadrant = calc_quadrant(theta_p);
+        sin_p = sinf(theta_p);
+        cos_p = cosf(theta_p);
+
+        for (d=0; d<dx; d++) 
+        {
+            // Calculate coordinates
+            xi = -ox-oz;
+            yi = (1-dx)/2.0+d+mov;
+            calc_coords(
+                ox, oz, xi, yi, sin_p, cos_p, gridx, gridy, 
+                coordx, coordy);
+
+            // Merge the (coordx, gridy) and (gridx, coordy)
+            trim_coords(
+                ox, oz, 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, 
+                &csize, coorx, coory);
+
+            // Calculate the distances (dist) between the 
+            // intersection points (coorx, coory). Find the 
+            // indices of the pixels on the object grid.
+            calc_dist2(
+                ox, oz, csize, coorx, coory, 
+                indx, indy, dist);
+
+            // For each slice
+            for (s=0; s<dy; s++) 
+            {
+                // Calculate simdata 
+                calc_simdata2(s, p, d, ox, oz, dt, dx,
+                    csize, indx, indy, dist, obj, axis,
+                    data); // Output: simulated data
+            }
+        }
+    }
+
+    free(gridx);
+    free(gridy);
+    free(coordx);
+    free(coordy);
+    free(ax);
+    free(ay);
+    free(bx);
+    free(by);
+    free(coorx);
+    free(coory);
+    free(dist);
+    free(indi);
+}
+
+
+void 
+project3(
+    const float *objx, const float *objy, const float *objz,
     int oy, int ox, int oz, 
     float *data, int dy, int dt, int dx,
-    const float *center, const float *theta)
+    const float *center, const float *theta, int axis)
 {
     float *gridx = (float *)malloc((ox+1)*sizeof(float));
     float *gridy = (float *)malloc((oz+1)*sizeof(float));
@@ -159,6 +260,8 @@ project2(
     float *coorx = (float *)malloc((ox+oz+2)*sizeof(float));
     float *coory = (float *)malloc((ox+oz+2)*sizeof(float));
     float *dist = (float *)malloc((ox+oz+1)*sizeof(float));
+    int *indx = (int *)malloc((ox+oz+1)*sizeof(int));
+    int *indy = (int *)malloc((ox+oz+1)*sizeof(int));
     int *indi = (int *)malloc((ox+oz+1)*sizeof(int));
 
     assert(coordx != NULL && coordy != NULL &&
@@ -171,9 +274,9 @@ project2(
     float srcx, srcy, detx, dety, dv, vx, vy;
     float mov, xi, yi;
     int asize, bsize, csize;
-
+    
     preprocessing(ox, oz, dx, center[0], 
-        &mov, gridx, gridy); // Outputs: mov, gridx, gridy
+    &mov, gridx, gridy); // Outputs: mov, gridx, gridy
 
     // For each projection angle
     for (p=0; p<dt; p++)
@@ -221,16 +324,16 @@ project2(
             // Calculate the distances (dist) between the 
             // intersection points (coorx, coory). Find the 
             // indices of the pixels on the object grid.
-            calc_dist(
+            calc_dist2(
                 ox, oz, csize, coorx, coory, 
-                indi, dist);
+                indx, indy, dist);
 
             // For each slice
             for (s=0; s<dy; s++) 
             {
                 // Calculate simdata 
-                calc_simdata2(s, p, d, ox, oz, dt, dx,
-                    csize, indi, dist, vx, vy, objx, objy,
+                calc_simdata3(s, p, d, ox, oz, dt, dx,
+                    csize, indx, indy, dist, vx, vy, objx, objy, objz, axis,
                     data); // Output: simulated data
             }
         }
@@ -248,4 +351,5 @@ project2(
     free(coory);
     free(dist);
     free(indi);
+
 }
diff --git a/src/utils.c b/src/utils.c
@@ -242,6 +242,33 @@ calc_dist(
 }
 
 
+void 
+calc_dist2(
+    int ry, int rz, 
+    int csize, const float *coorx, const float *coory,
+    int *indx, int *indy, float *dist)
+{
+    int n, i1, i2;
+    float x1, x2;
+    float diffx, diffy, midx, midy;
+
+    for (n=0; n<csize-1; n++) 
+    {
+        diffx = coorx[n+1]-coorx[n];
+        diffy = coory[n+1]-coory[n];
+        dist[n] = sqrt(diffx*diffx+diffy*diffy);
+        midx = (coorx[n+1]+coorx[n])/2;
+        midy = (coory[n+1]+coory[n])/2;
+        x1 = midx+ry/2.;
+        x2 = midy+rz/2.;
+        i1 = (int)(midx+ry/2.);
+        i2 = (int)(midy+rz/2.);
+        indx[n] = i1-(i1>x1);
+        indy[n] = i2-(i2>x2);
+    }
+}
+
+
 void 
 calc_simdata(
     int s, int p, int d,
@@ -266,17 +293,65 @@ calc_simdata2(
     int s, int p, int d,
     int ry, int rz, 
     int dt, int dx,
-    int csize, const int *indi, const float *dist,
+    int csize, const int *indx, const int *indy, const float *dist,
+    const float *model, int axis, 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++) 
+        {
+            simdata[d + p*dx + s*dt*dx] += model[indx[n] + s*rz + indy[n]*ry*rz]*dist[n];
+        }
+    }
+}
+
+
+void 
+calc_simdata3(
+    int s, int p, int d,
+    int ry, int rz, 
+    int dt, int dx,
+    int csize, const int *indx, const int *indy, const float *dist,
     float vx, float vy, 
-    const float *modelx, const float *modely, 
-    float *simdata)
+    const float *modelx, const float *modely, const float *modelz, int axis, 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++) 
+    if (axis==0)
+    {
+        for (n=0; n<csize-1; 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];
+        }
+    }
+    else if (axis==1)
     {
-        simdata[index_data] += (modelx[indi[n]+index_model] * vx + modely[indi[n]+index_model] * vy) * dist[n];
+        for (n=0; n<csize-1; n++) 
+        {
+            simdata[d + p*dx + s*dt*dx] += (modely[s + indx[n]*rz + indy[n]*ry*rz] * vx + modelz[s + indx[n]*rz + indy[n]*ry*rz] * vy) * dist[n];
+        }
+    }
+    else if (axis==2)
+    {
+        for (n=0; n<csize-1; n++) 
+        {
+            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];
+        }
     }
 }