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

Last change on this file since 3628 was 3611, checked in by ldelgass, 11 years ago

Use nv namespace for classes in nanovis rather than prefixing class names with
Nv (still need to convert shader classes).
  • Property svn:eol-style set to native
File size: 15.7 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(origin.x, origin.y, origin.z,
267                                            width, height, depth, 4,
268                                            fourAnionVolume, fourCationVolume,
269                                            vmin, vmax, nzero_min, cellSize);
270
271    return zincBlendeVolume;
272}
273
274ZincBlendeVolume *
275ZincBlendeReconstructor::buildUp(const Vector3f& origin, const Vector3f& delta,
276                                   int width, int height, int depth,
277                                   int datacount, double emptyvalue, void* data)
278{
279    ZincBlendeVolume *zincBlendeVolume = NULL;
280    float *fourAnionVolume, *fourCationVolume;
281    int size = width * height * depth * 4;
282    fourAnionVolume = new float[size];
283    fourCationVolume = new float[size];
284
285    memset(fourAnionVolume, 0, size * sizeof(float));
286    memset(fourCationVolume, 0, size * sizeof(float));
287
288    _NvAtomInfo *srcPtr = (_NvAtomInfo *) data;
289
290    float *component4A, *component4B;
291    float vmin, vmax, nzero_min;
292    int index;
293    nzero_min = 1e23f;
294    vmin = vmax = srcPtr->atom;
295    for (int i = 0; i < datacount; ++i) {
296        index = srcPtr->getIndex(width, height);
297
298#ifdef _LOADER_DEBUG_
299        TRACE("[%d] index %d (width:%lf height:%lf depth:%lf)", i, index, srcPtr->indexX, srcPtr->indexY, srcPtr->indexZ);
300#endif
301
302        if (index < 0) {
303#ifdef _LOADER_DEBUG_
304            TRACE("There is an invalid data");
305#endif
306            srcPtr +=8;
307            continue;
308        }
309
310        component4A = fourAnionVolume + index;
311        component4B = fourCationVolume + index;
312
313        component4A[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
314        component4A[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
315        component4A[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
316        component4A[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
317     
318        component4B[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
319        component4B[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
320        component4B[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
321        component4B[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
322
323        vmax = _NvMax3(_NvMax4(component4A), _NvMax4(component4B), vmax);
324        vmin = _NvMin3(_NvMin4(component4A), _NvMin4(component4B), vmin);
325
326        if (vmin != 0.0 && vmin < nzero_min) {
327            nzero_min = vmin;   
328        }
329    }
330
331    double dv = vmax - vmin;
332    if (vmax != 0.0f) {
333        for (int i = 0; i < datacount; ++i) {
334            fourAnionVolume[i] = (fourAnionVolume[i] - vmin)/ dv;
335            fourCationVolume[i] = (fourCationVolume[i] - vmin) / dv;
336        }
337    }
338
339    Vector3f cellSize;
340    cellSize.x = 0.25 / width;
341    cellSize.y = 0.25 / height;
342    cellSize.z = 0.25 / depth;
343
344    zincBlendeVolume = new ZincBlendeVolume(origin.x, origin.y, origin.z,
345                                            width, height, depth, 4,
346                                            fourAnionVolume, fourCationVolume,
347                                            vmin, vmax, nzero_min, cellSize);
348
349    return zincBlendeVolume;
350}
351
352void ZincBlendeReconstructor::getLine(std::istream& sin)
353{
354    char ch;
355    int index = 0;
356    do {
357        sin.get(ch);
358        if (ch == '\n') break;
359        buff[index++] = ch;
360        if (ch == '>') {
361            if (buff[1] == '\\')
362                break;
363        }
364    } while (!sin.eof());
365
366    buff[index] = '\0';
367
368#ifdef _LOADER_DEBUG_
369    TRACE("%s", buff);
370#endif
371}
372
373ZincBlendeVolume *
374ZincBlendeReconstructor::loadFromMemory(void *dataBlock)
375{
376    ZincBlendeVolume *volume = NULL;
377    Vector3f origin, delta;
378    int width = 0, height = 0, depth = 0;
379    void *data = NULL;
380    int version = 1;
381
382    unsigned char *stream = (unsigned char *)dataBlock;
383    char str[5][20];
384    do {
385        getLine(stream);
386        if (buff[0] == '#') {
387            continue;
388        } else if (strstr((const char *)buff, "object") != 0) {
389            TRACE("VERSION 1");
390            version = 1;
391            break;
392        } else if (strstr(buff, "record format") != 0) {
393            TRACE("VERSION 2");
394            version = 2;
395            break;
396        }
397    } while (1);
398
399    if (version == 1) {
400        float dummy;
401
402        sscanf(buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4],&width, &height, &depth);
403        getLine(stream);
404        sscanf(buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
405        getLine(stream);
406        sscanf(buff, "%s%f%f%f", str[0], &(delta.x), &dummy, &dummy);
407        getLine(stream);
408        sscanf(buff, "%s%f%f%f", str[0], &dummy, &(delta.y), &dummy);
409        getLine(stream);
410        sscanf(buff, "%s%f%f%f", str[0], &dummy, &dummy, &(delta.z));
411        do {
412            getLine(stream);
413        } while (strcmp(buff, "<\\HDR>") != 0);
414
415        width = width / 4;
416        height = height / 4;
417        depth = depth / 4;
418        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
419        // 8 atom per cell, 4 double (x, y, z, and probability) per atom
420        memcpy(data, stream, width * height * depth * 8 * 4 * sizeof(double));
421
422        volume = buildUp(origin, delta, width, height, depth, data);
423
424        free(data);
425    } else if (version == 2) {
426        const char *pt;
427        int datacount = -1;
428        double emptyvalue;
429        do {
430            getLine(stream);
431            if ((pt = strstr(buff, "delta")) != 0) {   
432                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
433#ifdef _LOADER_DEBUG_
434                TRACE("delta : %f %f %f", delta.x, delta.y, delta.z);
435#endif
436            } else if ((pt = strstr(buff, "datacount")) != 0) {
437                sscanf(pt, "%s%d", str[0], &datacount);
438                TRACE("datacount = %d", datacount);
439            } else if ((pt = strstr(buff, "datatype")) != 0) {
440                sscanf(pt, "%s%s", str[0], str[1]);
441                if (strcmp(str[1], "double64")) {
442                }
443            } else if ((pt = strstr(buff, "count")) != 0) {
444                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
445#ifdef _LOADER_DEBUG_
446                TRACE("width height depth %d %d %d", width, height, depth);
447#endif
448            } else if ((pt = strstr(buff, "emptymark")) != 0) {
449                sscanf(pt, "%s%lf", str[0], &emptyvalue);
450#ifdef _LOADER_DEBUG_
451                TRACE("emptyvalue %lf", emptyvalue);
452#endif
453            } else if ((pt = strstr(buff, "emprymark")) != 0) {
454                sscanf(pt, "%s%lf", str[0], &emptyvalue);
455#ifdef _LOADER_DEBUG_
456                TRACE("emptyvalue %lf", emptyvalue);
457#endif
458            }
459        } while (strcmp(buff, "<\\HDR>") != 0 && strcmp(buff, "</HDR>") != 0);
460
461        if (datacount == -1) datacount = width * height * depth;
462
463        data = malloc(datacount * 8 * 4 * sizeof(double));
464        memset(data, 0, datacount * 8 * 4 * sizeof(double));
465        memcpy(data, stream, datacount * 8 * 4 * sizeof(double));
466
467        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
468
469        free(data);
470    }
471    return volume;
472}
473
474void ZincBlendeReconstructor::getLine(unsigned char*& stream)
475{
476    char ch;
477    int index = 0;
478    do {
479        ch = stream[0];
480        ++stream;
481        if (ch == '\n') break;
482        buff[index++] = ch;
483        if (ch == '>') {
484            if (buff[1] == '\\')
485                break;
486        }
487    } while (1);
488
489    buff[index] = '\0';
490
491#ifdef _LOADER_DEBUG_
492    TRACE("%s", buff);
493#endif
494}
Note: See TracBrowser for help on using the repository browser.