codebase / stimlib

Commit 16e854e8d248263dc2e23c57ced58475fa709595

Authored by Pavel Govyadinov
2 parents 2e325f29 e38fc0c3

Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

Showing 2 changed files with 95 additions and 43 deletions   Show diff stats
Inline Side-by-side
stim/envi/agilent_binary.h
  Show/Hide comments   View file @16e854e
... ... @@ -5,6 +5,8 @@
5 5 #include <string>
6 6 #include <fstream>
7 7 #include <complex>
  8 +#include <cstring>
  9 +#include <chrono>
8 10  
9 11 //CUDA
10 12 //#ifdef CUDA_FOUND
... ... @@ -44,6 +46,10 @@ public:
44 46 alloc();
45 47 }
46 48  
  49 + char* data() {
  50 + return (char*)ptr;
  51 + }
  52 +
47 53 size_t dim(size_t i){
48 54 return R[i];
49 55 }
... ... @@ -58,17 +64,20 @@ public:
58 64 /// Default constructor, sets the resolution to zero and the data pointer to NULL
59 65 agilent_binary(){
60 66 memset(R, 0, sizeof(size_t) * 3); //set the resolution to zero
  67 + memset(Z, 0, sizeof(double) * 2);
61 68 ptr = NULL;
62 69 }
63 70  
64 71 /// Constructor with resolution
65 72 agilent_binary(size_t x, size_t y, size_t z){
66 73 alloc(x, y, z);
  74 + memset(Z, 0, sizeof(double) * 2);
67 75 }
68 76  
69 77 /// Constructor with filename
70 78 agilent_binary(std::string filename){
71 79 ptr = NULL;
  80 + memset(Z, 0, sizeof(double) * 2);
72 81 load(filename);
73 82 }
74 83  
... ... @@ -84,6 +93,11 @@ public:
84 93 return *this; //return the result
85 94 }
86 95  
  96 + operator bool() {
  97 + if (R[0] == 0 || R[1] == 0 || R[2] == 0) return false;
  98 + else return true;
  99 + }
  100 +
87 101 ~agilent_binary(){
88 102 free(ptr);
89 103 }
... ... @@ -110,6 +124,8 @@ public:
110 124 ptr = (T*) malloc(bytes()); //allocate space for the data
111 125 infile.read((char*)ptr, bytes()); //read the data
112 126 infile.close();
  127 + Z[0] = 1;
  128 + Z[1] = R[2];
113 129 }
114 130  
115 131 void save(std::string filename){
... ... @@ -189,6 +205,17 @@ public:
189 205 ptr[i] = -log10(ptr[i] / background->ptr[i]);
190 206 }
191 207  
  208 + //crops the image down to a set number of samples
  209 + void crop(size_t n) {
  210 + if (n < R[2]) { //if the requested size is smaller than the image
  211 + R[2] = n; //update the number of bands
  212 + T* old_ptr = ptr; //store the old pointer
  213 + alloc(); //allocate space for the new image
  214 + memcpy(ptr, old_ptr, bytes()); //copy the old data to the new image
  215 + free(old_ptr); //free the old data
  216 + }
  217 + }
  218 +
192 219 //#ifdef CUDA_FOUND
193 220 /// Perform an FFT and return a binary file with bands in the specified range
194 221 agilent_binary<T> fft(double band_min, double band_max, double ELWN = 15798, int UDR = 2){
... ... @@ -242,18 +269,19 @@ public:
242 269 HANDLE_ERROR(cudaMemcpy(cpu_fft, gpu_fft, R[0] * R[1] * (R[2]/2+1) * sizeof(cufftComplex), cudaMemcpyDeviceToHost)); //copy data from the host to the device
243 270  
244 271 //double int_delta = 0.00012656; //interferogram sample spacing in centimeters
245   - double int_delta = (1.0 / ELWN) * ((double)UDR / 2.0); //calculate the interferogram spacing
246   - double int_length = int_delta * R[2]; //interferogram length in centimeters
247   - double fft_delta = 1/int_length; //spectrum spacing (in inverse centimeters, wavenumber)
248   - double fft_max = fft_delta * R[2]/2; //get the maximum wavenumber value supported by the specified number of interferogram samples
  272 + double int_delta = (1.0 / ELWN) * ((double)UDR / 2.0); //calculate the interferogram spacing
  273 + double int_length = int_delta * R[2]; //interferogram length in centimeters
  274 + double fft_delta = 1/int_length; //spectrum spacing (in inverse centimeters, wavenumber)
  275 + double fft_max = fft_delta * R[2]/2; //get the maximum wavenumber value supported by the specified number of interferogram samples
249 276  
250   - if(band_max > fft_max) band_max = fft_max; //the user gave a band outside of the FFT range, reset the band to the maximum available
  277 + if(band_max > fft_max) band_max = fft_max; //the user gave a band outside of the FFT range, reset the band to the maximum available
  278 + if (band_min < 0) band_min = 0;
251 279  
252 280 size_t start_i = (size_t)std::ceil(band_min / fft_delta); //calculate the first band to store
253 281 size_t size_i = (size_t)std::floor(band_max / fft_delta) - start_i; //calculate the number of bands to store
254   - size_t end_i = start_i + size_i; //last band number
  282 + size_t end_i = start_i + size_i; //last band number
255 283 agilent_binary<T> result(R[0], R[1], size_i);
256   - result.Z[0] = start_i * fft_delta; //set the range for the FFT result
  284 + result.Z[0] = start_i * fft_delta; //set the range for the FFT result
257 285 result.Z[1] = end_i * fft_delta;
258 286  
259 287 for(size_t b = start_i; b < end_i; b++){
... ...
stim/visualization/gl_network.h
  Show/Hide comments   View file @16e854e
... ... @@ -45,14 +45,19 @@ public:
45 45  
46 46 ///render cylinder based on points from the top/bottom hat
47 47 ///@param C1 set of points from one of the hat
48   - void renderCylinder(std::vector< stim::vec3<T> > C1, std::vector< stim::vec3<T> > C2) {
  48 + void renderCylinder(std::vector< stim::vec3<T> > C1, std::vector< stim::vec3<T> > C2, stim::vec3<T> P1, stim::vec3<T> P2) {
49 49 glBegin(GL_QUAD_STRIP);
50   - for (unsigned i = 0; i < C1.size(); i++) { // for every point on the circle
  50 + for (unsigned i = 0; i < C1.size(); i++) { // for every point on the circle
  51 + stim::vec3<T> n1 = C1[i] - P1; // set normal vector for every vertex on the quads
  52 + stim::vec3<T> n2 = C2[i] - P2;
  53 + n1 = n1.norm();
  54 + n2 = n2.norm();
  55 + glNormal3f(n1[0], n1[1], n1[2]);
51 56 glVertex3f(C1[i][0], C1[i][1], C1[i][2]);
52   - glVertex3f(C2[i][0], C2[i][1], C2[i][2]);
  57 + glNormal3f(n2[0], n2[1], n2[2]);
  58 + glVertex3f(C2[i][0], C2[i][1], C2[i][2]);
53 59 }
54 60 glEnd();
55   - //glFlush();
56 61 }
57 62  
58 63 ///render the vertex as sphere based on glut build-in function
... ... @@ -80,12 +85,10 @@ public:
80 85 T angle_xy = 0.0;
81 86  
82 87 glBegin(GL_QUADS);
83   - for (unsigned i = 0; i < slice; i++) // around the z-axis
84   - {
  88 + for (unsigned i = 0; i < slice; i++) { // around the z-axis
85 89 angle_z = i * step_z; // step step_z each time
86 90  
87   - for (unsigned j = 0; j < stack; j++) // along the z-axis
88   - {
  91 + for (unsigned j = 0; j < stack; j++) { // along the z-axis
89 92 angle_xy = j * step_xy; // step step_xy each time, draw floor by floor
90 93  
91 94 xx[0] = radius * std::sin(angle_z) * std::cos(angle_xy); // four vertices
... ... @@ -104,8 +107,7 @@ public:
104 107 yy[3] = radius * std::sin(angle_z) * std::sin(angle_xy + step_xy);
105 108 zz[3] = radius * std::cos(angle_z);
106 109  
107   - for (unsigned k = 0; k < 4; k++)
108   - {
  110 + for (unsigned k = 0; k < 4; k++) {
109 111 glVertex3f(x + xx[k], y + yy[k], z + zz[k]); // draw the floor plane
110 112 }
111 113 }
... ... @@ -162,37 +164,42 @@ public:
162 164 if (!glIsList(dlist)) { // if dlist isn't a display list, create it
163 165 dlist = glGenLists(1); // generate a display list
164 166 glNewList(dlist, GL_COMPILE); // start a new display list
165   - glEnable(GL_COLOR_MATERIAL);
166   - glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
167 167 for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
168 168 for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
169 169 stim::circle<T> C1 = E[e].circ(p - 1);
170 170 stim::circle<T> C2 = E[e].circ(p);
171   - C1.set_R(2*radius); // scale the circle to the same
  171 + C1.set_R(2*radius); // re-scale the circle to the same
172 172 C2.set_R(2*radius);
173   - std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
174   - std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  173 + std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane
  174 + std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20);
175 175 glBegin(GL_QUAD_STRIP);
176   - for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle)
  176 + for (unsigned i = 0; i < Cp1.size(); i++) {
  177 + stim::vec3<T> n1 = Cp1[i] - E[e][p - 1]; // set normal vector for every vertex on the quads
  178 + stim::vec3<T> n2 = Cp2[i] - E[e][p];
  179 + n1 = n1.norm();
  180 + n2 = n2.norm();
  181 +
  182 + glNormal3f(n1[0], n1[1], n1[2]);
177 183 glTexCoord1f(E[e].r(p - 1));
178 184 glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
  185 + glNormal3f(n2[0], n2[1], n2[2]);
179 186 glTexCoord1f(E[e].r(p));
180 187 glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
181 188 }
182 189 glEnd();
183   - } // set the texture coordinate based on the specified magnitude index
  190 + } // set the texture coordinate based on the specified magnitude index
184 191 }
185 192 for (unsigned n = 0; n < V.size(); n++) {
186 193 size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex
187 194 if (num != 0) { // if it has outgoing edge
188   - unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex
  195 + unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex
189 196 glTexCoord1f(E[idx].r(0)); // bind the texture as metric of first point on that edge
190 197 }
191 198 else {
192 199 unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
193 200 glTexCoord1f(E[idx].r(E[idx].size() - 1)); // bind the texture as metric of last point on that edge
194 201 }
195   - renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20, 20);
  202 + renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20);
196 203 }
197 204 glEndList(); // end the display list
198 205 }
... ... @@ -230,7 +237,7 @@ public:
230 237 else {
231 238 unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
232 239 }
233   - renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20, 20);
  240 + renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20);
234 241 }
235 242 glEndList();
236 243 }
... ... @@ -247,51 +254,68 @@ public:
247 254 glDeleteLists(dlist+2, 1);
248 255 if (!glIsList(dlist+2)) { // if dlist isn't a display list, create it
249 256 glNewList(dlist+2, GL_COMPILE); // start a new display list
250   - glEnable(GL_COLOR_MATERIAL);
251   - glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
252 257 for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
253 258 if (map[e] != unsigned(-1)) {
254   - if(I == 0) // if it is to render GT
  259 + if (I == 0) { // if it is to render GT
255 260 glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
256   - else // if it is to render T
  261 + }
  262 + else { // if it is to render T
257 263 glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
  264 + }
258 265  
259 266 for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
260 267 stim::circle<T> C1 = E[e].circ(p - 1);
261 268 stim::circle<T> C2 = E[e].circ(p);
262   - C1.set_R(2*radius); // scale the circle to the same
  269 + C1.set_R(2*radius); // re-scale the circle to the same
263 270 C2.set_R(2*radius);
264 271 std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
265 272 std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
266   - renderCylinder(Cp1, Cp2);
  273 + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
267 274 }
268 275 }
269 276 else {
270   - glColor3f(1.0f, 1.0f, 1.0f); // white color for the un-mapping edges
271   - for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
  277 + glColor3f(0.2, 0.0, 0.0); // red color for the un-mapping edges
  278 + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
272 279 stim::circle<T> C1 = E[e].circ(p - 1);
273 280 stim::circle<T> C2 = E[e].circ(p);
274 281 C1.set_R(2*radius); // scale the circle to the same
275 282 C2.set_R(2*radius);
276 283 std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
277 284 std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
278   - renderCylinder(Cp1, Cp2);
  285 + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
279 286 }
280 287 }
281 288 }
282   - for (unsigned v = 0; v < V.size(); v++) {
283   - size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
  289 + for (unsigned n = 0; n < V.size(); n++) {
  290 + size_t num_edge = V[n].e[0].size() + V[n].e[1].size();
284 291 if (num_edge > 1) { // if it is the joint vertex
285   - glColor3f(0.3, 0.3, 0.3); // gray color
286   - renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20, 20);
  292 + if (V[n].e[0].size() != 0) {
  293 + unsigned idx = V[n].e[0][0]; // find the index of first incoming edge of that vertex
  294 + stim::circle<T> C = E[idx].circ(0);
  295 + stim::vec3<T> normal = C.n();
  296 + glNormal3f(normal[0], normal[1], normal[2]);// for every ball, we only have one normal vector
  297 + }
  298 + else {
  299 + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
  300 + stim::circle<T> C = E[idx].circ(E[idx].size() - 1);
  301 + stim::vec3<T> normal = C.n();
  302 + glNormal3f(normal[0], normal[1], normal[2]);// for every ball, we only have one normal vector
  303 + }
  304 + glColor3f(0.3, 0.3, 0.3); // gray
  305 + renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
287 306 }
288 307 else { // if it is the terminal vertex
  308 + if (V[n].e[0].size() != 0) {
  309 + unsigned idx = V[n].e[0][0]; // find the index of first incoming edge of that vertex
  310 + }
  311 + else {
  312 + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
  313 + }
289 314 glColor3f(0.6, 0.6, 0.6); // more white gray
290   - renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20, 20);
  315 + renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
291 316 }
292 317 }
293 318 glEndList();
294   - glDisable(GL_COLOR_MATERIAL);
295 319 }
296 320 glCallList(dlist+2);
297 321 }
... ... @@ -319,7 +343,7 @@ public:
319 343 glEnd();
320 344 }
321 345 else {
322   - glColor3f(1.0, 1.0, 1.0); // white color
  346 + glColor3f(0.2, 0.0, 0.0); // red color for the un-mapping edges
323 347 glBegin(GL_LINE_STRIP);
324 348 for (unsigned p = 0; p < E[e].size(); p++) {
325 349 glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
... ...