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

Last change on this file since 2870 was 2844, checked in by ldelgass, 12 years ago

Cleanups, no functional changes
  • 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\n");
57            version = 1;
58            break;
59        } else if (strstr(buff, "record format") != 0) {
60            TRACE("VERSION 2\n");
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\n");
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\n", 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\n", 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\n", 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\n", 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\n", 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\n", 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, 1, 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)\n", 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\n");
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, 1, 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\n");
380            version = 1;
381            break;
382        } else if (strstr(buff, "record format") != 0) {
383            TRACE("VERSION 2\n");
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\n", 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\n", 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\n", 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\n", 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\n", 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.