#5 works without any odd errors

jamesavery · May 5, 2022 · 5c3b831 · 5c3b831
1 parent cd66f4d
commit 5c3b831
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 14 deletions.
diff --git a/src/generate_gauss_c.py b/src/generate_gauss_c.py
@@ -56,12 +56,12 @@ def tobyt(arr):
 if __name__ == '__main__':
     sample, scale = commandline_args({"sample":"<required>","scale":1})
     outpath = 'dummy'
-    display = 500
+    display = 10
     sigma = 13
-    reps = 20
+    reps =20
 
     vf = h5py.File(f'{hdf5_root_fast}/processed/implant-edt/{scale}x/{sample}.h5', 'r')
-    voxels = vf['voxels'][:,:,:]
+    voxels = vf['voxels'][512:576,:,:]
     vf.close()
 
     Image.fromarray(tobyt(voxels[display,:,:])).save(f"{outpath}/original.png")
@@ -81,8 +81,12 @@ def tobyt(arr):
     histograms.gauss_filter_par_cpu(implant_mask, implant_mask.shape, kernel, reps, result)
     print (f'Parallel C edition: {timeit.default_timer() - start}')
     np.save(f'{outpath}/mine', result)
-
-    Image.fromarray(tobyt(result[display,:,:])).save(f'{outpath}/gauss1.png')
+
+    nz,ny,nx = result.shape    
+    Image.fromarray(tobyt(result[nz//2,:,:])).save(f'{outpath}/gauss-xy.png')
+    Image.fromarray(tobyt(result[:,ny//2,:])).save(f'{outpath}/gauss-xz.png')
+    Image.fromarray(tobyt(result[:,:,nx//2])).save(f'{outpath}/gauss-yz.png')
+    Image.fromarray(tobyt((np.max(np.abs(result),axis=0)!=0).astype(float))).save(f'{outpath}/gauss-xy-nonzero.png')
 
     control = implant_mask.copy()
     start = timeit.default_timer()

diff --git a/src/pybind_kernels/histograms.cc b/src/pybind_kernels/histograms.cc
@@ -52,8 +52,8 @@ void gauss_filter_par_cpu(const py::array_t<float> np_voxels,
         Rx = result_info.shape[2];
 
     float 
-        *tmp0 = (float *) malloc(sizeof(float) * Rz*Ry*Rx),
-        *tmp1 = (float *) malloc(sizeof(float) * Rz*Ry*Rx);
+        *tmp0 = (float *) calloc(Rz*Ry*Rx, sizeof(float)),
+        *tmp1 = (float *) calloc(Rz*Ry*Rx, sizeof(float));
 
     memcpy(tmp0, voxels, Rz*Ry*Rx * sizeof(float));
     uint64_t iters = 3 * reps; // 1 pass for each dimension
@@ -80,14 +80,19 @@ void gauss_filter_par_cpu(const py::array_t<float> np_voxels,
                     const int64_t X[3] = {z,y,x};
                     const int64_t N[3] = {Pz,Py,Px};
 
-                    const int64_t stride = strides[dim], i_start = -min(padding,X[dim]), i_end = min(padding,N[dim]-X[dim]);
-
-                    #pragma omp simd reduction(+:sum)
-                    for (int64_t i = i_start; i <= i_end; i++) {
-                        uint64_t voxel_index = output_index + stride*i;
-                        sum += tin[voxel_index] * kernel[i+padding];
+                    const int64_t stride = strides[dim], i_start = -min(padding,X[dim]-1), i_end = min(padding,N[dim]-X[dim]-1);
+
+                    auto mask_value = voxels[output_index];
+                    if(mask_value>0) {
+                        tout[output_index] = 1;
+                    } else {
+                        #pragma omp simd reduction(+:sum)
+                        for (int64_t i = i_start; i <= i_end; i++) {
+                            uint64_t voxel_index = output_index + stride*i;
+                            sum += tin[voxel_index] * kernel[i+padding];
+                        }
+                        tout[output_index] = sum;
                     }
-                    tout[output_index] = sum;
                 }
             }
         }