source: trunk/packages/vizservers/nanovis/ZincBlendeReconstructor.cpp @ 3884

Last change on this file since 3884 was 3630, checked in by ldelgass, 11 years ago

Nanovis refactoring to fix problems with scaling and multiple results.
Do rendering in world space to properly place and scale multiple data sets.
Also fix flows to reduce resets of animations. More work toward removing
Cg dependency. Fix panning to convert viewport coords to world coords.
  • Property svn:eol-style set to native
File size: 16.0 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * Copyright (c) 2004-2013  HUBzero Foundation, LLC
4 *
5 * Authors:
6 *   Insoo Woo <iwoo@purdue.edu>
7 */
8#include <stdio.h>
9#include <string.h>
10#include <stdlib.h>
11
12#include <vrmath/Vector3f.h>
13
14#include "ZincBlendeReconstructor.h"
15#include "ZincBlendeVolume.h"
16#include "Trace.h"
17
18using namespace nv;
19using namespace vrmath;
20
21ZincBlendeReconstructor *ZincBlendeReconstructor::_instance = NULL;
22
23ZincBlendeReconstructor::ZincBlendeReconstructor()
24{
25}
26
27ZincBlendeReconstructor::~ZincBlendeReconstructor()
28{
29}
30
31ZincBlendeReconstructor *ZincBlendeReconstructor::getInstance()
32{
33    if (_instance == NULL) {
34        return (_instance = new ZincBlendeReconstructor());
35    }
36
37    return _instance;
38}
39
40ZincBlendeVolume *ZincBlendeReconstructor::loadFromFile(const char *fileName)
41{
42    std::ifstream stream;
43    stream.open(fileName, std::ios::binary);
44
45    ZincBlendeVolume *volume = loadFromStream(stream);
46
47    stream.close();
48
49    return volume;
50}
51
52ZincBlendeVolume *ZincBlendeReconstructor::loadFromStream(std::istream& stream)
53{
54    ZincBlendeVolume *volume = NULL;
55    Vector3f origin, delta;
56    int width = 0, height = 0, depth = 0;
57    void *data = NULL;
58    int version = 1;
59
60    char str[5][20];
61    do {
62        getLine(stream);
63        if (buff[0] == '#') {
64            continue;
65        } else if (strstr((const char*) buff, "object") != 0) {
66            TRACE("VERSION 1");
67            version = 1;
68            break;
69        } else if (strstr(buff, "record format") != 0) {
70            TRACE("VERSION 2");
71            version = 2;
72            break;
73        }
74    } while (1);
75
76    if (version == 1) {
77        float dummy;
78
79        sscanf(buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4], &width, &height, &depth);
80        getLine(stream);
81        sscanf(buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
82        getLine(stream);
83        sscanf(buff, "%s%f%f%f", str[0], &(delta.x), &dummy, &dummy);
84        getLine(stream);
85        sscanf(buff, "%s%f%f%f", str[0], &dummy, &(delta.y), &dummy);
86        getLine(stream);
87        sscanf(buff, "%s%f%f%f", str[0], &dummy, &dummy, &(delta.z));
88        do {
89            getLine(stream);
90        } while (strcmp(buff, "<\\HDR>") != 0);
91
92        width = width / 4;
93        height = height / 4;
94        depth = depth / 4;
95        //data = new double[width * height * depth * 8 * 4];
96        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
97        // 8 atom per cell, 4 double (x, y, z, and probability) per atom
98        try {
99            stream.read((char *)data, width * height * depth * 8 * 4 * sizeof(double));
100        } catch (...) {
101            ERROR("Caught exception trying to read stream");
102        }
103
104        volume = buildUp(origin, delta, width, height, depth, data);
105
106        free(data);
107    } else if (version == 2) {
108        const char *pt;
109        int datacount;
110        double emptyvalue;
111        do {
112            getLine(stream);
113            if ((pt = strstr(buff, "delta")) != 0) {
114                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
115#ifdef _LOADER_DEBUG_
116                TRACE("delta : %f %f %f", delta.x, delta.y, delta.z);
117#endif
118            } else if ((pt = strstr(buff, "datacount")) != 0) {
119                sscanf(pt, "%s%d", str[0], &datacount);
120#ifdef _LOADER_DEBUG_
121                TRACE("datacount = %d", datacount);
122#endif
123            } else if ((pt = strstr(buff, "datatype")) != 0) {
124                sscanf(pt, "%s%s", str[0], str[1]);
125                if (strcmp(str[1], "double64")) {
126                }
127            } else if ((pt = strstr(buff, "count")) != 0) {
128                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
129#ifdef _LOADER_DEBUG_
130                TRACE("width height depth %d %d %d", width, height, depth);
131#endif
132            } else if ((pt = strstr(buff, "emptymark")) != 0) {
133                sscanf(pt, "%s%lf", str[0], &emptyvalue);
134#ifdef _LOADER_DEBUG_
135                TRACE("emptyvalue %lf", emptyvalue);
136#endif
137            } else if ((pt = strstr(buff, "emprymark")) != 0) {
138                sscanf(pt, "%s%lf", str[0], &emptyvalue);
139#ifdef _LOADER_DEBUG_
140                TRACE("emptyvalue %lf", emptyvalue);
141#endif
142            }
143        } while (strcmp(buff, "<\\HDR>") != 0 && strcmp(buff, "</HDR>") != 0);
144
145        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
146        memset(data, 0, width * height * depth * 8 * 4 * sizeof(double));
147        stream.read((char *) data, width * height * depth * 8 * 4 * sizeof(double));
148
149        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
150        free(data);
151    }
152    return volume;
153}
154
155struct _NvAtomInfo {
156    double indexX, indexY, indexZ;
157    double atom;
158
159    int getIndex(int width, int height) const
160    {
161        // NOTE
162        // Zinc blende data has different axes from OpenGL
163        // + z -> +x (OpenGL)
164        // + x -> +y (OpenGL)
165        // + y -> +z (OpenGL), But in 3D texture coordinate is the opposite direction of z
166        // The reason why index is multiplied by 4 is that one unit cell has half of eight atoms,
167        // i.e. four atoms are mapped into RGBA component of one texel
168        //return ((int) (indexZ - 1)+ (int) (indexX - 1) * width + (int) (indexY - 1) * width * height) * 4;
169        return ((int)(indexX - 1) + (int)(indexY - 1) * width + (int)(indexZ - 1) * width * height) * 4;
170    }
171};
172
173template<class T>
174inline T _NvMax2(T a, T b)
175{ return ((a >= b)? a : b); }
176
177template<class T>
178inline T _NvMin2(T a, T b)
179{ return ((a >= b)? a : b); }
180
181template<class T>
182inline T _NvMax3(T a, T b, T c)
183{ return ((a >= b)? ((a >= c) ? a : c) : ((b >= c)? b : c)); }
184
185template<class T>
186inline T _NvMin3(T a, T b, T c)
187{ return ((a <= b)? ((a <= c) ? a : c) : ((b <= c)? b : c)); }
188
189template<class T>
190inline T _NvMax9(T* a, T curMax)
191{ return _NvMax3(_NvMax3(a[0], a[1], a[2]), _NvMax3(a[3], a[4], a[5]), _NvMax3(a[6], a[7], curMax)); }
192
193template<class T>
194inline T _NvMin9(T* a, T curMax)
195{ return _NvMin3(_NvMax3(a[0], a[1], a[2]), _NvMin3(a[3], a[4], a[5]), _NvMin3(a[6], a[7], curMax)); }
196
197template<class T>
198inline T _NvMax4(T* a)
199{ return _NvMax2(_NvMax2(a[0], a[1]), _NvMax2(a[2], a[3])); }
200
201template<class T>
202inline T _NvMin4(T* a)
203{ return _NvMin2(_NvMin2(a[0], a[1]), _NvMin2(a[2], a[3])); }
204
205ZincBlendeVolume *
206ZincBlendeReconstructor::buildUp(const Vector3f& origin, const Vector3f& delta,
207                                   int width, int height, int depth, void *data)
208{
209    ZincBlendeVolume *zincBlendeVolume = NULL;
210
211    float *fourAnionVolume, *fourCationVolume;
212    int cellCount = width * height * depth;
213    fourAnionVolume = new float[cellCount * sizeof(float) * 4];
214    fourCationVolume = new float[cellCount * sizeof(float) * 4];
215
216    _NvAtomInfo *srcPtr = (_NvAtomInfo *)data;
217
218    float vmin, vmax, nzero_min;
219    float *component4A, *component4B;
220    int index;
221
222    nzero_min = 0.0f;           /* Suppress compiler warning. */
223    vmin = vmax = srcPtr->atom;
224
225    for (int i = 0; i < cellCount; ++i) {
226        index = srcPtr->getIndex(width, height);
227
228#ifdef _LOADER_DEBUG_
229        TRACE("index %d", index);
230#endif
231
232        component4A = fourAnionVolume + index;
233        component4B = fourCationVolume + index;
234
235        component4A[0] = (float)srcPtr->atom; srcPtr++;
236        component4A[1] = (float)srcPtr->atom; srcPtr++;
237        component4A[2] = (float)srcPtr->atom; srcPtr++;
238        component4A[3] = (float)srcPtr->atom; srcPtr++;
239     
240        component4B[0] = (float)srcPtr->atom; srcPtr++;
241        component4B[1] = (float)srcPtr->atom; srcPtr++;
242        component4B[2] = (float)srcPtr->atom; srcPtr++;
243        component4B[3] = (float)srcPtr->atom; srcPtr++;
244
245        vmax = _NvMax3(_NvMax4(component4A), _NvMax4(component4B), vmax);
246        vmin = _NvMin3(_NvMin4(component4A), _NvMin4(component4B), vmin);
247
248        if (vmin != 0.0 && vmin < nzero_min) {
249            nzero_min = vmin;
250        }
251    }
252
253    double dv = vmax - vmin;
254    if (vmax != 0.0f) {
255        for (int i = 0; i < cellCount; ++i) {
256            fourAnionVolume[i] = (fourAnionVolume[i] - vmin)/ dv;
257            fourCationVolume[i] = (fourCationVolume[i] - vmin) / dv;
258        }
259    }
260
261    Vector3f cellSize;
262    cellSize.x = 0.25 / width;
263    cellSize.y = 0.25 / height;
264    cellSize.z = 0.25 / depth;
265
266    zincBlendeVolume = new ZincBlendeVolume(width, height, depth, 4,
267                                            fourAnionVolume, fourCationVolume,
268                                            vmin, vmax, nzero_min, cellSize);
269
270    zincBlendeVolume->xAxis.setRange(origin.x, origin.x + delta.x * (width-1));
271    zincBlendeVolume->yAxis.setRange(origin.y, origin.y + delta.y * (height-1));
272    zincBlendeVolume->zAxis.setRange(origin.z, origin.z + delta.z * (depth-1));
273
274    return zincBlendeVolume;
275}
276
277ZincBlendeVolume *
278ZincBlendeReconstructor::buildUp(const Vector3f& origin, const Vector3f& delta,
279                                   int width, int height, int depth,
280                                   int datacount, double emptyvalue, void* data)
281{
282    ZincBlendeVolume *zincBlendeVolume = NULL;
283    float *fourAnionVolume, *fourCationVolume;
284    int size = width * height * depth * 4;
285    fourAnionVolume = new float[size];
286    fourCationVolume = new float[size];
287
288    memset(fourAnionVolume, 0, size * sizeof(float));
289    memset(fourCationVolume, 0, size * sizeof(float));
290
291    _NvAtomInfo *srcPtr = (_NvAtomInfo *) data;
292
293    float *component4A, *component4B;
294    float vmin, vmax, nzero_min;
295    int index;
296    nzero_min = 1e23f;
297    vmin = vmax = srcPtr->atom;
298    for (int i = 0; i < datacount; ++i) {
299        index = srcPtr->getIndex(width, height);
300
301#ifdef _LOADER_DEBUG_
302        TRACE("[%d] index %d (width:%lf height:%lf depth:%lf)", i, index, srcPtr->indexX, srcPtr->indexY, srcPtr->indexZ);
303#endif
304
305        if (index < 0) {
306#ifdef _LOADER_DEBUG_
307            TRACE("There is an invalid data");
308#endif
309            srcPtr +=8;
310            continue;
311        }
312
313        component4A = fourAnionVolume + index;
314        component4B = fourCationVolume + index;
315
316        component4A[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
317        component4A[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
318        component4A[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
319        component4A[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
320     
321        component4B[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
322        component4B[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
323        component4B[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
324        component4B[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
325
326        vmax = _NvMax3(_NvMax4(component4A), _NvMax4(component4B), vmax);
327        vmin = _NvMin3(_NvMin4(component4A), _NvMin4(component4B), vmin);
328
329        if (vmin != 0.0 && vmin < nzero_min) {
330            nzero_min = vmin;   
331        }
332    }
333
334    double dv = vmax - vmin;
335    if (vmax != 0.0f) {
336        for (int i = 0; i < datacount; ++i) {
337            fourAnionVolume[i] = (fourAnionVolume[i] - vmin)/ dv;
338            fourCationVolume[i] = (fourCationVolume[i] - vmin) / dv;
339        }
340    }
341
342    Vector3f cellSize;
343    cellSize.x = 0.25 / width;
344    cellSize.y = 0.25 / height;
345    cellSize.z = 0.25 / depth;
346
347    zincBlendeVolume = new ZincBlendeVolume(width, height, depth, 4,
348                                            fourAnionVolume, fourCationVolume,
349                                            vmin, vmax, nzero_min, cellSize);
350
351    zincBlendeVolume->xAxis.setRange(origin.x, origin.x + delta.x * (width-1));
352    zincBlendeVolume->yAxis.setRange(origin.y, origin.y + delta.y * (height-1));
353    zincBlendeVolume->zAxis.setRange(origin.z, origin.z + delta.z * (depth-1));
354
355    return zincBlendeVolume;
356}
357
358void ZincBlendeReconstructor::getLine(std::istream& sin)
359{
360    char ch;
361    int index = 0;
362    do {
363        sin.get(ch);
364        if (ch == '\n') break;
365        buff[index++] = ch;
366        if (ch == '>') {
367            if (buff[1] == '\\')
368                break;
369        }
370    } while (!sin.eof());
371
372    buff[index] = '\0';
373
374#ifdef _LOADER_DEBUG_
375    TRACE("%s", buff);
376#endif
377}
378
379ZincBlendeVolume *
380ZincBlendeReconstructor::loadFromMemory(void *dataBlock)
381{
382    ZincBlendeVolume *volume = NULL;
383    Vector3f origin, delta;
384    int width = 0, height = 0, depth = 0;
385    void *data = NULL;
386    int version = 1;
387
388    unsigned char *stream = (unsigned char *)dataBlock;
389    char str[5][20];
390    do {
391        getLine(stream);
392        if (buff[0] == '#') {
393            continue;
394        } else if (strstr((const char *)buff, "object") != 0) {
395            TRACE("VERSION 1");
396            version = 1;
397            break;
398        } else if (strstr(buff, "record format") != 0) {
399            TRACE("VERSION 2");
400            version = 2;
401            break;
402        }
403    } while (1);
404
405    if (version == 1) {
406        float dummy;
407
408        sscanf(buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4],&width, &height, &depth);
409        getLine(stream);
410        sscanf(buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
411        getLine(stream);
412        sscanf(buff, "%s%f%f%f", str[0], &(delta.x), &dummy, &dummy);
413        getLine(stream);
414        sscanf(buff, "%s%f%f%f", str[0], &dummy, &(delta.y), &dummy);
415        getLine(stream);
416        sscanf(buff, "%s%f%f%f", str[0], &dummy, &dummy, &(delta.z));
417        do {
418            getLine(stream);
419        } while (strcmp(buff, "<\\HDR>") != 0);
420
421        width = width / 4;
422        height = height / 4;
423        depth = depth / 4;
424        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
425        // 8 atom per cell, 4 double (x, y, z, and probability) per atom
426        memcpy(data, stream, width * height * depth * 8 * 4 * sizeof(double));
427
428        volume = buildUp(origin, delta, width, height, depth, data);
429
430        free(data);
431    } else if (version == 2) {
432        const char *pt;
433        int datacount = -1;
434        double emptyvalue;
435        do {
436            getLine(stream);
437            if ((pt = strstr(buff, "delta")) != 0) {   
438                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
439#ifdef _LOADER_DEBUG_
440                TRACE("delta : %f %f %f", delta.x, delta.y, delta.z);
441#endif
442            } else if ((pt = strstr(buff, "datacount")) != 0) {
443                sscanf(pt, "%s%d", str[0], &datacount);
444                TRACE("datacount = %d", datacount);
445            } else if ((pt = strstr(buff, "datatype")) != 0) {
446                sscanf(pt, "%s%s", str[0], str[1]);
447                if (strcmp(str[1], "double64")) {
448                }
449            } else if ((pt = strstr(buff, "count")) != 0) {
450                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
451#ifdef _LOADER_DEBUG_
452                TRACE("width height depth %d %d %d", width, height, depth);
453#endif
454            } else if ((pt = strstr(buff, "emptymark")) != 0) {
455                sscanf(pt, "%s%lf", str[0], &emptyvalue);
456#ifdef _LOADER_DEBUG_
457                TRACE("emptyvalue %lf", emptyvalue);
458#endif
459            } else if ((pt = strstr(buff, "emprymark")) != 0) {
460                sscanf(pt, "%s%lf", str[0], &emptyvalue);
461#ifdef _LOADER_DEBUG_
462                TRACE("emptyvalue %lf", emptyvalue);
463#endif
464            }
465        } while (strcmp(buff, "<\\HDR>") != 0 && strcmp(buff, "</HDR>") != 0);
466
467        if (datacount == -1) datacount = width * height * depth;
468
469        data = malloc(datacount * 8 * 4 * sizeof(double));
470        memset(data, 0, datacount * 8 * 4 * sizeof(double));
471        memcpy(data, stream, datacount * 8 * 4 * sizeof(double));
472
473        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
474
475        free(data);
476    }
477    return volume;
478}
479
480void ZincBlendeReconstructor::getLine(unsigned char*& stream)
481{
482    char ch;
483    int index = 0;
484    do {
485        ch = stream[0];
486        ++stream;
487        if (ch == '\n') break;
488        buff[index++] = ch;
489        if (ch == '>') {
490            if (buff[1] == '\\')
491                break;
492        }
493    } while (1);
494
495    buff[index] = '\0';
496
497#ifdef _LOADER_DEBUG_
498    TRACE("%s", buff);
499#endif
500}
Note: See TracBrowser for help on using the repository browser.