source: nanovis/branches/1.2/ZincBlendeReconstructor.cpp @ 5699

Last change on this file since 5699 was 5406, checked in by ldelgass, 9 years ago

backport from trunk: remove unused ctor params

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