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

Last change on this file since 2831 was 2822, checked in by ldelgass, 13 years ago

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