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

Last change on this file since 3492 was 3492, checked in by ldelgass, 11 years ago

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