mmc_mesh.h

Go to the documentation of this file.

/***************************************************************************/
/***************************************************************************/
#ifndef _MMC_MESH_UNIT_H
#define _MMC_MESH_UNIT_H

#include <stdio.h>
#include <math.h>
#include "mmc_utils.h"

#ifdef MMC_USE_SSE
    #include <smmintrin.h>
#endif

#if !defined(USE_OPENCL) && !defined(__NVCC__)

    #ifdef MMC_SFMT
        #include "mmc_rand_sfmt.c"
    #elif defined MMC_XORSHIFT
        #include "mmc_rand_xorshift128p.c"
    #else
        #include "mmc_rand_posix.c"
    #endif

#elif defined(__NVCC__)
    #include "mmc_rand_xorshift128p.h"
#else
    #include "mmc_rand_xorshift128p.c"
#endif

#define MMC_UNDEFINED (3.40282347e+38F)
#define ID_UNDEFINED  0x7FFFFFFF

#define R_RAND_MAX (1.f/RAND_MAX)
#define TWO_PI     (M_PI*2.0)
#define EPS        1e-6f
#define LOG_MT_MAX 22.1807097779182f
#define R_MIN_MUS  1e9f
#define R_C0       3.335640951981520e-12f  //1/C0 in s/mm
#define DELTA_MUA  1e-4f
#define VERY_BIG   1e30f
#define MIN(a,b)  ((a)<(b)?(a):(b))
#define MAX(a,b)  ((a)>(b)?(a):(b))

#define MESH_ERROR(a)  mesh_error((a),__FILE__,__LINE__)

/***************************************************************************/
typedef struct MMC_mesh {
    int nn;                
    int ne;                
    int nf;                
    int prop;              
    int elemlen;           
    float3* node;          
    int*  elem;            
    int*  elem2;           
    float* edgeroi;        
    float* faceroi;        
    float* noderoi;        
    int*  srcelem;         
    int  srcelemlen;       
    int*  detelem;         
    int  detelemlen;       
    int*  type;            
    int*  facenb;          
    medium* med;           
    float* atte;           
    double* weight;        
    double* dref;          
    float* evol;           
    float* nvol;           
    float4 nmin;           
    float4 nmax;           
} tetmesh;

/***************************************************************************/
typedef struct MMC_raytracer {
    tetmesh* mesh;          
    char method;            
    float3* d;              
    float3* m;              
    float3* n;              
} raytracer;
#ifdef  __cplusplus
extern "C" {
#endif
void mesh_init(tetmesh* mesh);
void mesh_init_from_cfg(tetmesh* mesh, mcconfig* cfg);
void mesh_loadnode(tetmesh* mesh, mcconfig* cfg);
void mesh_loadelem(tetmesh* mesh, mcconfig* cfg);
void mesh_loadfaceneighbor(tetmesh* mesh, mcconfig* cfg);
void mesh_loadmedia(tetmesh* mesh, mcconfig* cfg);
void mesh_loadelemvol(tetmesh* mesh, mcconfig* cfg);
void mesh_loadseedfile(tetmesh* mesh, mcconfig* cfg);

void mesh_clear(tetmesh* mesh);
float mesh_normalize(tetmesh* mesh, mcconfig* cfg, float Eabsorb, float Etotal, int pair);
void mesh_build(tetmesh* mesh);
void mesh_error(const char* msg, const char* file, const int linenum);
void mesh_filenames(const char* format, char* foutput, mcconfig* cfg);
void mesh_saveweight(tetmesh* mesh, mcconfig* cfg, int isref);
void mesh_savedetphoton(float* ppath, void* seeds, int count, int seedbyte, mcconfig* cfg);
void mesh_getdetimage(float* detmap, float* ppath, int count, mcconfig* cfg, tetmesh* mesh);
void mesh_savedetimage(float* detmap, mcconfig* cfg);
float mesh_getdetweight(int photonid, int colcount, float* ppath, mcconfig* cfg);
void mesh_srcdetelem(tetmesh* mesh, mcconfig* cfg);
void mesh_createdualmesh(tetmesh* mesh, mcconfig* cfg);
void mesh_loadroi(tetmesh* mesh, mcconfig* cfg);
double mesh_getreff_approx(double n_in, double n_out);
double mesh_getreff(double n_in, double n_out);

void tracer_init(raytracer* tracer, tetmesh* mesh, char methodid);
void tracer_build(raytracer* tracer);
void tracer_prep(raytracer* tracer, mcconfig* cfg);
void tracer_clear(raytracer* tracer);

float mc_next_scatter(float g, float3* dir, RandType* ran, RandType* ran0, mcconfig* cfg, float* pmom);
#ifdef MCX_CONTAINER
#ifdef __cplusplus
extern "C"
#endif
int mcx_throw_exception(const int id, const char* msg, const char* filename, const int linenum);
#endif

#ifdef  __cplusplus
}
#endif

static inline void vec_add(float3* a, float3* b, float3* res) {
    res->x = a->x + b->x;
    res->y = a->y + b->y;
    res->z = a->z + b->z;
}

static inline void vec_diff(float3* a, float3* b, float3* res) {
    res->x = b->x - a->x;
    res->y = b->y - a->y;
    res->z = b->z - a->z;
}

static inline void vec_mult(float3* a, float sa, float3* res) {
    res->x = sa * a->x;
    res->y = sa * a->y;
    res->z = sa * a->z;
}

static inline void vec_mult_add(float3* a, float3* b, float sa, float sb, float3* res) {
    res->x = sb * b->x + sa * a->x;
    res->y = sb * b->y + sa * a->y;
    res->z = sb * b->z + sa * a->z;
}

static inline void vec_cross(float3* a, float3* b, float3* res) {
    res->x = a->y * b->z - a->z * b->y;
    res->y = a->z * b->x - a->x * b->z;
    res->z = a->x * b->y - a->y * b->x;
}

static inline void mmc_sincosf(float x, float* sine, float* cosine) {
#if defined(__GNUC__) && defined(__linux__) && !defined(__clang__)
    __builtin_sincosf(x, sine, cosine);
#else
    *sine = sinf(x);
    *cosine = cosf(x);
#endif
}

//#ifndef MMC_USE_SSE
static inline float vec_dot(float3* a, float3* b) {
    return a->x * b->x + a->y * b->y + a->z * b->z;
}/*
#else

#ifndef __SSE4_1__
static inline float vec_dot(float3 *a,float3 *b){
        float dot;
        __m128 na,nb,res;
        na=_mm_load_ps(&a->x);
        nb=_mm_load_ps(&b->x);
        res=_mm_mul_ps(na,nb);
        res=_mm_hadd_ps(res,res);
        res=_mm_hadd_ps(res,res);
        _mm_store_ss(&dot,res);
        return dot;
}
#else
static inline float vec_dot(float3 *a,float3 *b){
        float dot;
        __m128 na,nb,res;
        na=_mm_load_ps(&a->x);
        nb=_mm_load_ps(&b->x);
        res=_mm_dp_ps(na,nb,0x7f);
        _mm_store_ss(&dot,res);
        return dot;
}
#endif
#endif
*/

static inline float pinner(float3* Pd, float3* Pm, float3* Ad, float3* Am) {
    return vec_dot(Pd, Am) + vec_dot(Pm, Ad);
}


static inline float dist2(float3* p0, float3* p1) {
    return (p1->x - p0->x) * (p1->x - p0->x) + (p1->y - p0->y) * (p1->y - p0->y) + (p1->z - p0->z) * (p1->z - p0->z);
}

static inline float dist(float3* p0, float3* p1) {
    return sqrtf(dist2(p0, p1));
}

static inline float dist2d2(float* p0, float* p1) {
    return (p1[0] - p0[0]) * (p1[0] - p0[0]) + (p1[1] - p0[1]) * (p1[1] - p0[1]);
}

static inline float mmc_rsqrtf(float a) {
#ifdef MMC_USE_SSE
    return _mm_cvtss_f32( _mm_rsqrt_ss( _mm_set_ss( a ) ) );
#else
    return 1.f / sqrtf(a);
#endif
}

#endif