source: nanovis/trunk/ZincBlendeReconstructor.cpp @ 4880

Last change on this file since 4880 was 4063, checked in by ldelgass, 10 years ago

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