source: trunk/packages/vizservers/nanovis/NvZincBlendeReconstructor.cpp @ 3474

Last change on this file since 3474 was 3452, checked in by ldelgass, 12 years ago

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