source: trunk/packages/vizservers/nanovis/NvZincBlendeReconstructor.cpp @ 3502

Last change on this file since 3502 was 3502, checked in by ldelgass, 6 years ago

Add basic VTK structured points reader to nanovis, update copyright dates.

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