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

Last change on this file since 3611 was 3611, checked in by ldelgass, 9 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.