Jiaming Guo / flow3

Commit 4bb615eb7ac9970c5056b1306622fece19d8920d

Authored by Jiaming Guo
1 parent 6332d1e7

fixed generating image-stack bugs

Showing 2 changed files with 48 additions and 37 deletions   Show diff stats
Inline Side-by-side
@@ -67,7 +67,7 @@ namespace stim { @@ -67,7 +67,7 @@ namespace stim {
67 #ifdef __CUDACC__ 67 #ifdef __CUDACC__
68 // for sphere 68 // for sphere
69 template <typename T> 69 template <typename T>
70 - __global__ void inside_sphere(const stim::sphere<T> *V, size_t num, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) { 70 + __global__ void inside_sphere(const stim::sphere<T> *V, size_t num, T Z,size_t *R, T *S, unsigned char *ptr, int x, int y, int z) {
71 71
72 unsigned ix = blockDim.x * blockIdx.x + threadIdx.x; 72 unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
73 unsigned iy = blockDim.y * blockIdx.y + threadIdx.y; 73 unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
@@ -78,7 +78,7 @@ namespace stim { @@ -78,7 +78,7 @@ namespace stim {
78 stim::vec3<T> world_pixel; 78 stim::vec3<T> world_pixel;
79 world_pixel[0] = (T)ix * S[1] - x; // translate origin to center of the network 79 world_pixel[0] = (T)ix * S[1] - x; // translate origin to center of the network
80 world_pixel[1] = (T)iy * S[2] - y; 80 world_pixel[1] = (T)iy * S[2] - y;
81 - world_pixel[2] = ((T)z - R[3] / 2) * S[3]; // ???center of box minus half width 81 + world_pixel[2] = ((T)z - Z / 2.0f) * S[3]; // ???center of box minus half width
82 82
83 float distance = FLT_MAX; 83 float distance = FLT_MAX;
84 float tmp_distance; 84 float tmp_distance;
@@ -97,7 +97,7 @@ namespace stim { @@ -97,7 +97,7 @@ namespace stim {
97 97
98 // for cone 98 // for cone
99 template <typename T> 99 template <typename T>
100 - __global__ void inside_cone(const stim::cone<T> *E, size_t num, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) { 100 + __global__ void inside_cone(const stim::cone<T> *E, size_t num, T Z, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) {
101 101
102 unsigned ix = blockDim.x * blockIdx.x + threadIdx.x; 102 unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
103 unsigned iy = blockDim.y * blockIdx.y + threadIdx.y; 103 unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
@@ -107,7 +107,7 @@ namespace stim { @@ -107,7 +107,7 @@ namespace stim {
107 stim::vec3<T> world_pixel; 107 stim::vec3<T> world_pixel;
108 world_pixel[0] = (T)ix * S[1] - x; 108 world_pixel[0] = (T)ix * S[1] - x;
109 world_pixel[1] = (T)iy * S[2] - y; 109 world_pixel[1] = (T)iy * S[2] - y;
110 - world_pixel[2] = ((T)z - R[3] / 2) * S[3]; 110 + world_pixel[2] = ((T)z - Z / 2.0f) * S[3];
111 111
112 float distance = FLT_MAX; 112 float distance = FLT_MAX;
113 float tmp_distance; 113 float tmp_distance;
@@ -133,7 +133,7 @@ namespace stim { @@ -133,7 +133,7 @@ namespace stim {
133 133
134 // for source bus 134 // for source bus
135 template <typename T> 135 template <typename T>
136 - __global__ void inside_cuboid(const stim::cuboid<T> *B, size_t num, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) { 136 + __global__ void inside_cuboid(const stim::cuboid<T> *B, size_t num, T Z, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) {
137 137
138 unsigned ix = blockDim.x * blockIdx.x + threadIdx.x; 138 unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
139 unsigned iy = blockDim.y * blockIdx.y + threadIdx.y; 139 unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
@@ -143,7 +143,7 @@ namespace stim { @@ -143,7 +143,7 @@ namespace stim {
143 stim::vec3<T> world_pixel; 143 stim::vec3<T> world_pixel;
144 world_pixel[0] = (T)ix * S[1] - x; 144 world_pixel[0] = (T)ix * S[1] - x;
145 world_pixel[1] = (T)iy * S[2] - y; 145 world_pixel[1] = (T)iy * S[2] - y;
146 - world_pixel[2] = ((T)z - R[3] / 2) * S[3]; 146 + world_pixel[2] = ((T)z - Z / 2.0f) * S[3];
147 147
148 for (unsigned i = 0; i < num; i++) { 148 for (unsigned i = 0; i < num; i++) {
149 bool left_outside = false; // flag indicates point is outside the left bound 149 bool left_outside = false; // flag indicates point is outside the left bound
@@ -1042,8 +1042,9 @@ namespace stim { @@ -1042,8 +1042,9 @@ namespace stim {
1042 void glSolidCuboid(GLint subdivision, bool manufacture = false, T length = 40.0f, T height = 10.0f) { 1042 void glSolidCuboid(GLint subdivision, bool manufacture = false, T length = 40.0f, T height = 10.0f) {
1043 1043
1044 T width; 1044 T width;
1045 - stim::vec3<T> L = bb.A; // get the bottom left corner  
1046 - stim::vec3<T> U = bb.B; // get the top right corner 1045 + stim::gl_aaboundingbox<T> BB = (*this).boundingbox();
  1046 + stim::vec3<T> L = BB.A; // get the bottom left corner
  1047 + stim::vec3<T> U = BB.B; // get the top right corner
1047 width = U[2] - L[2] + 10.0f; 1048 width = U[2] - L[2] + 10.0f;
1048 1049
1049 if (manufacture) 1050 if (manufacture)
@@ -1887,8 +1888,8 @@ namespace stim { @@ -1887,8 +1888,8 @@ namespace stim {
1887 1888
1888 // there are cases that the fragment can not satisfy the requirement for width 1889 // there are cases that the fragment can not satisfy the requirement for width
1889 if (width < times * radius || n == 0) { // check feasibility 1890 if (width < times * radius || n == 0) { // check feasibility
1890 - ratio = 0.0f; // load original lengths  
1891 - desire_l = tmp_d; 1891 + ratio = 1.0f; // load original lengths
  1892 + desire_l = tmp_d; // loading back-ups
1892 origin_l = tmp_l; 1893 origin_l = tmp_l;
1893 1894
1894 std::cout << "Warning: current ratio is not feasible, use full original line." << std::endl; 1895 std::cout << "Warning: current ratio is not feasible, use full original line." << std::endl;
@@ -1923,11 +1924,21 @@ namespace stim { @@ -1923,11 +1924,21 @@ namespace stim {
1923 height = (T)((desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2)); 1924 height = (T)((desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2));
1924 } 1925 }
1925 1926
1926 - // degenerated case, compromise 1927 + /// degenerated case, use original extending method
1927 if (n == 0) { 1928 if (n == 0) {
1928 - n = 1;  
1929 - width = (T)(origin_l) / (2 * n);  
1930 - height = (T)((desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2)); 1929 + height = (T)((desire_l - origin_l) / 2);
  1930 + // "up"
  1931 + tmp_v = tmp_v + stim::vec3<T>(0, coef_up * height, 0);
  1932 + if (feeder == 1)
  1933 + inlet[i].V.push_back(tmp_v);
  1934 + else if (feeder == 0)
  1935 + outlet[i].V.push_back(tmp_v);
  1936 + // "left"
  1937 + tmp_v = tmp_v + stim::vec3<T>(-coef_left * origin_l, 0, 0);
  1938 + if (feeder == 1)
  1939 + inlet[i].V.push_back(tmp_v);
  1940 + else if (feeder == 0)
  1941 + outlet[i].V.push_back(tmp_v);
1931 } 1942 }
1932 1943
1933 // cube-like structure construction 1944 // cube-like structure construction
@@ -2628,7 +2639,7 @@ namespace stim { @@ -2628,7 +2639,7 @@ namespace stim {
2628 2639
2629 /// make full-synthetic binary image stack 2640 /// make full-synthetic binary image stack
2630 // prepare for image stack 2641 // prepare for image stack
2631 - void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, bool prototype = false, T length = 40.0f, T height = 10.0f, T radius = 5.0f) { 2642 + void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, bool prototype = false, T length = 40.0f, T height = 10.0f, T radius = 5.0f, T scale = 1.0f) {
2632 2643
2633 T max_radius = 0.0f; 2644 T max_radius = 0.0f;
2634 T top = FLT_MIN; 2645 T top = FLT_MIN;
@@ -2658,13 +2669,13 @@ namespace stim { @@ -2658,13 +2669,13 @@ namespace stim {
2658 for (unsigned i = 0; i < num_edge; i++) { 2669 for (unsigned i = 0; i < num_edge; i++) {
2659 for (unsigned j = 0; j < E[i].size(); j++) { 2670 for (unsigned j = 0; j < E[i].size(); j++) {
2660 new_sphere.c = E[i][j]; 2671 new_sphere.c = E[i][j];
2661 - new_sphere.r = E[i].r(j); 2672 + new_sphere.r = E[i].r(j) * scale;
2662 A.push_back(new_sphere); 2673 A.push_back(new_sphere);
2663 if (j != E[i].size() - 1) { 2674 if (j != E[i].size() - 1) {
2664 new_cone.c1 = E[i][j]; 2675 new_cone.c1 = E[i][j];
2665 new_cone.c2 = E[i][j + 1]; 2676 new_cone.c2 = E[i][j + 1];
2666 - new_cone.r1 = E[i].r(j);  
2667 - new_cone.r2 = E[i].r(j + 1); 2677 + new_cone.r1 = E[i].r(j) * scale;
  2678 + new_cone.r2 = E[i].r(j + 1) * scale;
2668 B.push_back(new_cone); 2679 B.push_back(new_cone);
2669 } 2680 }
2670 } 2681 }
@@ -2675,14 +2686,14 @@ namespace stim { @@ -2675,14 +2686,14 @@ namespace stim {
2675 for (unsigned i = 0; i < inlet.size(); i++) { 2686 for (unsigned i = 0; i < inlet.size(); i++) {
2676 for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) { 2687 for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) {
2677 new_sphere.c = inlet[i].V[j]; 2688 new_sphere.c = inlet[i].V[j];
2678 - new_sphere.r = inlet[i].r; 2689 + new_sphere.r = inlet[i].r * scale;
2679 A.push_back(new_sphere); 2690 A.push_back(new_sphere);
2680 } 2691 }
2681 } 2692 }
2682 for (unsigned i = 0; i < outlet.size(); i++) { 2693 for (unsigned i = 0; i < outlet.size(); i++) {
2683 for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) { 2694 for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) {
2684 new_sphere.c = outlet[i].V[j]; 2695 new_sphere.c = outlet[i].V[j];
2685 - new_sphere.r = outlet[i].r; 2696 + new_sphere.r = outlet[i].r * scale;
2686 A.push_back(new_sphere); 2697 A.push_back(new_sphere);
2687 } 2698 }
2688 } 2699 }
@@ -2691,8 +2702,8 @@ namespace stim { @@ -2691,8 +2702,8 @@ namespace stim {
2691 for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) { 2702 for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) {
2692 new_cone.c1 = inlet[i].V[j]; 2703 new_cone.c1 = inlet[i].V[j];
2693 new_cone.c2 = inlet[i].V[j + 1]; 2704 new_cone.c2 = inlet[i].V[j + 1];
2694 - new_cone.r1 = inlet[i].r;  
2695 - new_cone.r2 = inlet[i].r; 2705 + new_cone.r1 = inlet[i].r * scale;
  2706 + new_cone.r2 = inlet[i].r * scale;
2696 B.push_back(new_cone); 2707 B.push_back(new_cone);
2697 } 2708 }
2698 } 2709 }
@@ -2700,15 +2711,15 @@ namespace stim { @@ -2700,15 +2711,15 @@ namespace stim {
2700 for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) { 2711 for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) {
2701 new_cone.c1 = outlet[i].V[j]; 2712 new_cone.c1 = outlet[i].V[j];
2702 new_cone.c2 = outlet[i].V[j + 1]; 2713 new_cone.c2 = outlet[i].V[j + 1];
2703 - new_cone.r1 = outlet[i].r;  
2704 - new_cone.r2 = outlet[i].r; 2714 + new_cone.r1 = outlet[i].r * scale;
  2715 + new_cone.r2 = outlet[i].r * scale;
2705 B.push_back(new_cone); 2716 B.push_back(new_cone);
2706 } 2717 }
2707 } 2718 }
2708 2719
2709 // find out the image stack size 2720 // find out the image stack size
2710 - Xl = main_feeder[0][0] - length / 2; // left bound x coordinate  
2711 - Xr = main_feeder[1][0] + length / 2; // right bound x coordinate 2721 + Xl = main_feeder[0][0] - length / 2 - 2 * radius; // left bound x coordinate
  2722 + Xr = main_feeder[1][0] + length / 2 + 2 * radius; // right bound x coordinate
2712 2723
2713 for (unsigned i = 0; i < A.size(); i++) { 2724 for (unsigned i = 0; i < A.size(); i++) {
2714 if (A[i].c[1] > top) 2725 if (A[i].c[1] > top)
@@ -2720,21 +2731,21 @@ namespace stim { @@ -2720,21 +2731,21 @@ namespace stim {
2720 bb.B[2] = A[i].c[2]; 2731 bb.B[2] = A[i].c[2];
2721 if (A[i].c[2] < bb.A[2]) 2732 if (A[i].c[2] < bb.A[2])
2722 bb.A[2] = A[i].c[2]; 2733 bb.A[2] = A[i].c[2];
2723 - if (A[i].r > max_radius)  
2724 - max_radius = A[i].r; 2734 + if (A[i].r * scale > max_radius)
  2735 + max_radius = A[i].r * scale;
2725 } 2736 }
2726 2737
2727 Yt = top + 2 * radius; // top bound y coordinate 2738 Yt = top + 2 * radius; // top bound y coordinate
2728 Yb = bottom - 2 * radius; // bottom bound y coordinate 2739 Yb = bottom - 2 * radius; // bottom bound y coordinate
2729 - Z = (bb.B[2] - bb.A[2] + 2 * max_radius); // bounding box width(along z-axis) 2740 + Z = (bb.B[2] - bb.A[2] + 4 * radius); // bounding box width(along z-axis)
2730 } 2741 }
2731 2742
2732 /// making image stack main function 2743 /// making image stack main function
2733 - void make_image_stack(stim::image_stack<unsigned char, T> &I, T dx, T dy, T dz, std::string stackdir, bool prototype = false, T radius = 5.0f) { 2744 + void make_image_stack(stim::image_stack<unsigned char, T> &I, T dx, T dy, T dz, std::string stackdir, bool prototype = false, T radius = 5.0f, T scale = 1.0f) {
2734 2745
2735 /// preparation for making image stack 2746 /// preparation for making image stack
2736 T X, Xl, Xr, Y, Yt, Yb, Z; 2747 T X, Xl, Xr, Y, Yt, Yb, Z;
2737 - preparation(Xl, Xr, Yt, Yb, Z, prototype); 2748 + preparation(Xl, Xr, Yt, Yb, Z, prototype, 40.0f, 10.0f, radius, scale);
2738 X = Xr - Xl; // bounding box length(along x-axis) 2749 X = Xr - Xl; // bounding box length(along x-axis)
2739 Y = Yt - Yb; // bounding box height(along y-axis) 2750 Y = Yt - Yb; // bounding box height(along y-axis)
2740 stim::vec3<T> center = bb.center(); // get the center of bounding box 2751 stim::vec3<T> center = bb.center(); // get the center of bounding box
@@ -2744,7 +2755,7 @@ namespace stim { @@ -2744,7 +2755,7 @@ namespace stim {
2744 /// make 2755 /// make
2745 size_x = (int)(X / dx + 1); // set the size of image 2756 size_x = (int)(X / dx + 1); // set the size of image
2746 size_y = (int)(Y / dy + 1); 2757 size_y = (int)(Y / dy + 1);
2747 - size_z = (int)(Z / dz + 1); 2758 + size_z = (int)(Z / dz + 1); // +3 in order to deal with reminder
2748 /// initialize image stack object 2759 /// initialize image stack object
2749 I.init(1, size_x, size_y, size_z); 2760 I.init(1, size_x, size_y, size_z);
2750 I.set_dim(dx, dy, dz); 2761 I.set_dim(dx, dy, dz);
@@ -2811,11 +2822,11 @@ namespace stim { @@ -2811,11 +2822,11 @@ namespace stim {
2811 2822
2812 dim3 block((unsigned)(size_x / max_thread + 1), (unsigned)(size_y / max_thread + 1)); 2823 dim3 block((unsigned)(size_x / max_thread + 1), (unsigned)(size_y / max_thread + 1));
2813 dim3 thread((unsigned)max_thread, (unsigned)max_thread); 2824 dim3 thread((unsigned)max_thread, (unsigned)max_thread);
2814 - inside_sphere << <block, thread >> > (d_V, A.size(), d_R, d_S, d_ptr, x, y, z); 2825 + inside_sphere << <block, thread >> > (d_V, A.size(), Z, d_R, d_S, d_ptr, x, y, z);
2815 cudaDeviceSynchronize(); 2826 cudaDeviceSynchronize();
2816 - inside_cone << <block, thread >> > (d_E, B.size(), d_R, d_S, d_ptr, x, y, z); 2827 + inside_cone << <block, thread >> > (d_E, B.size(), Z, d_R, d_S, d_ptr, x, y, z);
2817 cudaDeviceSynchronize(); 2828 cudaDeviceSynchronize();
2818 - inside_cuboid << <block, thread >> > (d_B, CU.size(), d_R, d_S, d_ptr, x, y, z); 2829 + inside_cuboid << <block, thread >> > (d_B, CU.size(), Z, d_R, d_S, d_ptr, x, y, z);
2819 2830
2820 HANDLE_ERROR(cudaMemcpy(ptr, d_ptr, num * sizeof(unsigned char), cudaMemcpyDeviceToHost)); 2831 HANDLE_ERROR(cudaMemcpy(ptr, d_ptr, num * sizeof(unsigned char), cudaMemcpyDeviceToHost));
2821 2832
@@ -833,7 +833,7 @@ void glut_wheel(int wheel, int direction, int x, int y) { @@ -833,7 +833,7 @@ void glut_wheel(int wheel, int direction, int x, int y) {
833 GLdouble posX, posY, posZ; 833 GLdouble posX, posY, posZ;
834 window_to_world(posX, posY, posZ); // get the world coordinates 834 window_to_world(posX, posY, posZ); // get the world coordinates
835 835
836 - if (!to_select_pressure && (simulation || build_inlet_outlet)) { // check current pixel position only in simualtion and build_inlet_outlet modes 836 + if (!to_select_pressure && (simulation || build_inlet_outlet || manufacture)) { // check current pixel position only in simualtion and build_inlet_outlet modes
837 bool flag = flow.epsilon_vertex((float)posX, (float)posY, (float)posZ, eps, scale, pressure_index); 837 bool flag = flow.epsilon_vertex((float)posX, (float)posY, (float)posZ, eps, scale, pressure_index);
838 if (flag && simulation && !glyph_mode) { 838 if (flag && simulation && !glyph_mode) {
839 float tmp_r; 839 float tmp_r;
@@ -947,7 +947,7 @@ void glut_keyboard(unsigned char key, int x, int y) { @@ -947,7 +947,7 @@ void glut_keyboard(unsigned char key, int x, int y) {
947 case 'm': 947 case 'm':
948 if (manufacture) { 948 if (manufacture) {
949 #ifdef __CUDACC__ 949 #ifdef __CUDACC__
950 - flow.make_image_stack(S, dx, dy, dz, stackdir, image_stack); 950 + flow.make_image_stack(S, dx, dy, dz, stackdir, image_stack, default_radius, scale);
951 951
952 #else 952 #else
953 std::cout << "You need to have a gpu to make image stack, sorry." << std::endl; 953 std::cout << "You need to have a gpu to make image stack, sorry." << std::endl;