source: nanovis/branches/1.1/ZincBlendeReconstructor.cpp @ 4904

Last change on this file since 4904 was 4904, checked in by ldelgass, 5 years ago

Merge serveral changes from trunk. Does not include threading, world space
changes, etc.

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