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

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