source: nanovis/trunk/ZincBlendeReconstructor.cpp @ 5399

Last change on this file since 5399 was 4894, checked in by ldelgass, 9 years ago

sync with release branch

  • Property svn:eol-style set to native
File size: 16.1 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            } else if ((pt = strstr(_buff, "emprymark")) != 0) {
139                sscanf(pt, "%s%lf", str[0], &emptyvalue);
140#ifdef _LOADER_DEBUG_
141                TRACE("emptyvalue %lf", emptyvalue);
142#endif
143            }
144        } while (strcmp(_buff, "<\\HDR>") != 0 && strcmp(_buff, "</HDR>") != 0);
145
146        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
147        memset(data, 0, width * height * depth * 8 * 4 * sizeof(double));
148        stream.read((char *) data, width * height * depth * 8 * 4 * sizeof(double));
149
150        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
151        free(data);
152    }
153    return volume;
154}
155
156struct _NvAtomInfo {
157    double indexX, indexY, indexZ;
158    double atom;
159
160    int getIndex(int width, int height) const
161    {
162        // NOTE
163        // Zinc blende data has different axes from OpenGL
164        // + z -> +x (OpenGL)
165        // + x -> +y (OpenGL)
166        // + y -> +z (OpenGL), But in 3D texture coordinate is the opposite direction of z
167        // The reason why index is multiplied by 4 is that one unit cell has half of eight atoms,
168        // i.e. four atoms are mapped into RGBA component of one texel
169        //return ((int) (indexZ - 1)+ (int) (indexX - 1) * width + (int) (indexY - 1) * width * height) * 4;
170        return ((int)(indexX - 1) + (int)(indexY - 1) * width + (int)(indexZ - 1) * width * height) * 4;
171    }
172};
173
174template<class T>
175inline T _NvMax2(T a, T b)
176{ return ((a >= b)? a : b); }
177
178template<class T>
179inline T _NvMin2(T a, T b)
180{ return ((a >= b)? a : b); }
181
182template<class T>
183inline T _NvMax3(T a, T b, T c)
184{ return ((a >= b)? ((a >= c) ? a : c) : ((b >= c)? b : c)); }
185
186template<class T>
187inline T _NvMin3(T a, T b, T c)
188{ return ((a <= b)? ((a <= c) ? a : c) : ((b <= c)? b : c)); }
189
190template<class T>
191inline T _NvMax9(T* a, T curMax)
192{ return _NvMax3(_NvMax3(a[0], a[1], a[2]), _NvMax3(a[3], a[4], a[5]), _NvMax3(a[6], a[7], curMax)); }
193
194template<class T>
195inline T _NvMin9(T* a, T curMax)
196{ return _NvMin3(_NvMax3(a[0], a[1], a[2]), _NvMin3(a[3], a[4], a[5]), _NvMin3(a[6], a[7], curMax)); }
197
198template<class T>
199inline T _NvMax4(T* a)
200{ return _NvMax2(_NvMax2(a[0], a[1]), _NvMax2(a[2], a[3])); }
201
202template<class T>
203inline T _NvMin4(T* a)
204{ return _NvMin2(_NvMin2(a[0], a[1]), _NvMin2(a[2], a[3])); }
205
206ZincBlendeVolume *
207ZincBlendeReconstructor::buildUp(const Vector3f& origin, const Vector3f& delta,
208                                   int width, int height, int depth, void *data)
209{
210    ZincBlendeVolume *zincBlendeVolume = NULL;
211
212    float *fourAnionVolume, *fourCationVolume;
213    int cellCount = width * height * depth;
214    fourAnionVolume = new float[cellCount * sizeof(float) * 4];
215    fourCationVolume = new float[cellCount * sizeof(float) * 4];
216
217    _NvAtomInfo *srcPtr = (_NvAtomInfo *)data;
218
219    float vmin, vmax, nzero_min;
220    float *component4A, *component4B;
221    int index;
222
223    nzero_min = 0.0f;           /* Suppress compiler warning. */
224    vmin = vmax = srcPtr->atom;
225
226    for (int i = 0; i < cellCount; ++i) {
227        index = srcPtr->getIndex(width, height);
228
229#ifdef _LOADER_DEBUG_
230        TRACE("index %d", index);
231#endif
232
233        component4A = fourAnionVolume + index;
234        component4B = fourCationVolume + index;
235
236        component4A[0] = (float)srcPtr->atom; srcPtr++;
237        component4A[1] = (float)srcPtr->atom; srcPtr++;
238        component4A[2] = (float)srcPtr->atom; srcPtr++;
239        component4A[3] = (float)srcPtr->atom; srcPtr++;
240     
241        component4B[0] = (float)srcPtr->atom; srcPtr++;
242        component4B[1] = (float)srcPtr->atom; srcPtr++;
243        component4B[2] = (float)srcPtr->atom; srcPtr++;
244        component4B[3] = (float)srcPtr->atom; srcPtr++;
245
246        vmax = _NvMax3(_NvMax4(component4A), _NvMax4(component4B), vmax);
247        vmin = _NvMin3(_NvMin4(component4A), _NvMin4(component4B), vmin);
248
249        if (vmin != 0.0 && vmin < nzero_min) {
250            nzero_min = vmin;
251        }
252    }
253
254    double dv = vmax - vmin;
255    if (vmax != 0.0f) {
256        for (int i = 0; i < cellCount; ++i) {
257            fourAnionVolume[i] = (fourAnionVolume[i] - vmin)/ dv;
258            fourCationVolume[i] = (fourCationVolume[i] - vmin) / dv;
259        }
260    }
261
262    Vector3f cellSize;
263    cellSize.x = 0.25 / width;
264    cellSize.y = 0.25 / height;
265    cellSize.z = 0.25 / depth;
266
267    zincBlendeVolume = new ZincBlendeVolume(width, height, depth, 4,
268                                            fourAnionVolume, fourCationVolume,
269                                            vmin, vmax, nzero_min, cellSize);
270
271    zincBlendeVolume->xAxis.setRange(origin.x, origin.x + delta.x * (width-1));
272    zincBlendeVolume->yAxis.setRange(origin.y, origin.y + delta.y * (height-1));
273    zincBlendeVolume->zAxis.setRange(origin.z, origin.z + delta.z * (depth-1));
274
275    return zincBlendeVolume;
276}
277
278ZincBlendeVolume *
279ZincBlendeReconstructor::buildUp(const Vector3f& origin, const Vector3f& delta,
280                                   int width, int height, int depth,
281                                   int datacount, double emptyvalue, void* data)
282{
283    ZincBlendeVolume *zincBlendeVolume = NULL;
284    float *fourAnionVolume, *fourCationVolume;
285    int size = width * height * depth * 4;
286    fourAnionVolume = new float[size];
287    fourCationVolume = new float[size];
288
289    memset(fourAnionVolume, 0, size * sizeof(float));
290    memset(fourCationVolume, 0, size * sizeof(float));
291
292    _NvAtomInfo *srcPtr = (_NvAtomInfo *) data;
293
294    float *component4A, *component4B;
295    float vmin, vmax, nzero_min;
296    int index;
297    nzero_min = 1e23f;
298    vmin = vmax = srcPtr->atom;
299    for (int i = 0; i < datacount; ++i) {
300        index = srcPtr->getIndex(width, height);
301
302#ifdef _LOADER_DEBUG_
303        TRACE("[%d] index %d (width:%lf height:%lf depth:%lf)", i, index, srcPtr->indexX, srcPtr->indexY, srcPtr->indexZ);
304#endif
305
306        if (index < 0) {
307#ifdef _LOADER_DEBUG_
308            TRACE("There is an invalid data");
309#endif
310            srcPtr +=8;
311            continue;
312        }
313
314        component4A = fourAnionVolume + index;
315        component4B = fourCationVolume + index;
316
317        component4A[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
318        component4A[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
319        component4A[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
320        component4A[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
321     
322        component4B[0] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
323        component4B[1] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
324        component4B[2] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
325        component4B[3] = (srcPtr->atom != emptyvalue)? (float) srcPtr->atom : 0.0f; srcPtr++;
326
327        vmax = _NvMax3(_NvMax4(component4A), _NvMax4(component4B), vmax);
328        vmin = _NvMin3(_NvMin4(component4A), _NvMin4(component4B), vmin);
329
330        if (vmin != 0.0 && vmin < nzero_min) {
331            nzero_min = vmin;   
332        }
333    }
334
335    double dv = vmax - vmin;
336    if (vmax != 0.0f) {
337        for (int i = 0; i < datacount; ++i) {
338            fourAnionVolume[i] = (fourAnionVolume[i] - vmin)/ dv;
339            fourCationVolume[i] = (fourCationVolume[i] - vmin) / dv;
340        }
341    }
342
343    Vector3f cellSize;
344    cellSize.x = 0.25 / width;
345    cellSize.y = 0.25 / height;
346    cellSize.z = 0.25 / depth;
347
348    zincBlendeVolume = new ZincBlendeVolume(width, height, depth, 4,
349                                            fourAnionVolume, fourCationVolume,
350                                            vmin, vmax, nzero_min, cellSize);
351
352    zincBlendeVolume->xAxis.setRange(origin.x, origin.x + delta.x * (width-1));
353    zincBlendeVolume->yAxis.setRange(origin.y, origin.y + delta.y * (height-1));
354    zincBlendeVolume->zAxis.setRange(origin.z, origin.z + delta.z * (depth-1));
355
356    return zincBlendeVolume;
357}
358
359void ZincBlendeReconstructor::getLine(std::istream& sin)
360{
361    char ch;
362    int index = 0;
363    do {
364        sin.get(ch);
365        if (ch == '\n') break;
366        _buff[index++] = ch;
367        if (ch == '>') {
368            if (_buff[1] == '\\')
369                break;
370        }
371    } while (!sin.eof());
372
373    _buff[index] = '\0';
374
375#ifdef _LOADER_DEBUG_
376    TRACE("%s", _buff);
377#endif
378}
379
380ZincBlendeVolume *
381ZincBlendeReconstructor::loadFromMemory(const void *dataBlock)
382{
383    ZincBlendeVolume *volume = NULL;
384    Vector3f origin, delta;
385    int width = 0, height = 0, depth = 0;
386    void *data = NULL;
387    int version = 1;
388
389    const unsigned char *stream = (const unsigned char *)dataBlock;
390    char str[5][20];
391    do {
392        getLine(stream);
393        if (_buff[0] == '#') {
394            continue;
395        } else if (strstr((const char *)_buff, "object") != 0) {
396            TRACE("VERSION 1");
397            version = 1;
398            break;
399        } else if (strstr(_buff, "record format") != 0) {
400            TRACE("VERSION 2");
401            version = 2;
402            break;
403        }
404    } while (1);
405
406    if (version == 1) {
407        float dummy;
408
409        sscanf(_buff, "%s%s%s%s%s%d%d%d", str[0], str[1], str[2], str[3], str[4],&width, &height, &depth);
410        getLine(stream);
411        sscanf(_buff, "%s%f%f%f", str[0], &(origin.x), &(origin.y), &(origin.z));
412        getLine(stream);
413        sscanf(_buff, "%s%f%f%f", str[0], &(delta.x), &dummy, &dummy);
414        getLine(stream);
415        sscanf(_buff, "%s%f%f%f", str[0], &dummy, &(delta.y), &dummy);
416        getLine(stream);
417        sscanf(_buff, "%s%f%f%f", str[0], &dummy, &dummy, &(delta.z));
418        do {
419            getLine(stream);
420        } while (strcmp(_buff, "<\\HDR>") != 0);
421
422        width = width / 4;
423        height = height / 4;
424        depth = depth / 4;
425        data = malloc(width * height * depth * 8 * 4 * sizeof(double));
426        // 8 atom per cell, 4 double (x, y, z, and probability) per atom
427        memcpy(data, stream, width * height * depth * 8 * 4 * sizeof(double));
428
429        volume = buildUp(origin, delta, width, height, depth, data);
430
431        free(data);
432    } else if (version == 2) {
433        const char *pt;
434        int datacount = -1;
435        double emptyvalue;
436        do {
437            getLine(stream);
438            if ((pt = strstr(_buff, "delta")) != 0) {   
439                sscanf(pt, "%s%f%f%f", str[0], &(delta.x), &(delta.y), &(delta.z));
440#ifdef _LOADER_DEBUG_
441                TRACE("delta : %f %f %f", delta.x, delta.y, delta.z);
442#endif
443            } else if ((pt = strstr(_buff, "datacount")) != 0) {
444                sscanf(pt, "%s%d", str[0], &datacount);
445                TRACE("datacount = %d", datacount);
446            } else if ((pt = strstr(_buff, "datatype")) != 0) {
447                sscanf(pt, "%s%s", str[0], str[1]);
448                if (strcmp(str[1], "double64")) {
449                }
450            } else if ((pt = strstr(_buff, "count")) != 0) {
451                sscanf(pt, "%s%d%d%d", str[0], &width, &height, &depth);
452#ifdef _LOADER_DEBUG_
453                TRACE("width height depth %d %d %d", width, height, depth);
454#endif
455            } else if ((pt = strstr(_buff, "emptymark")) != 0) {
456                sscanf(pt, "%s%lf", str[0], &emptyvalue);
457#ifdef _LOADER_DEBUG_
458                TRACE("emptyvalue %lf", emptyvalue);
459#endif
460            } else if ((pt = strstr(_buff, "emprymark")) != 0) {
461                sscanf(pt, "%s%lf", str[0], &emptyvalue);
462#ifdef _LOADER_DEBUG_
463                TRACE("emptyvalue %lf", emptyvalue);
464#endif
465            }
466        } while (strcmp(_buff, "<\\HDR>") != 0 && strcmp(_buff, "</HDR>") != 0);
467
468        if (datacount == -1) datacount = width * height * depth;
469
470        data = malloc(datacount * 8 * 4 * sizeof(double));
471        memset(data, 0, datacount * 8 * 4 * sizeof(double));
472        memcpy(data, stream, datacount * 8 * 4 * sizeof(double));
473
474        volume =  buildUp(origin, delta, width, height, depth, datacount, emptyvalue, data);
475
476        free(data);
477    }
478    return volume;
479}
480
481void ZincBlendeReconstructor::getLine(const unsigned char*& stream)
482{
483    char ch;
484    int index = 0;
485    do {
486        ch = stream[0];
487        ++stream;
488        if (ch == '\n') break;
489        _buff[index++] = ch;
490        if (ch == '>') {
491            if (_buff[1] == '\\')
492                break;
493        }
494    } while (1);
495
496    _buff[index] = '\0';
497
498#ifdef _LOADER_DEBUG_
499    TRACE("%s", _buff);
500#endif
501}
Note: See TracBrowser for help on using the repository browser.