From c0eb87dd3b6609f0ddbfbaee6ba4ff33cd05354a Mon Sep 17 00:00:00 2001 From: Adrien Thob Date: Tue, 1 Jan 2019 16:56:15 +0000 Subject: [PATCH 1/7] Update rendermodule.c Simplified the code for dealing with the remainder particles --- sphviewer/extensions/rendermodule.c | 22 +--------------------- 1 file changed, 1 insertion(+), 21 deletions(-) diff --git a/sphviewer/extensions/rendermodule.c b/sphviewer/extensions/rendermodule.c index f1012ed..05bb5df 100644 --- a/sphviewer/extensions/rendermodule.c +++ b/sphviewer/extensions/rendermodule.c @@ -96,7 +96,7 @@ void c_render(float *x, float *y, float *t, float *mass, } } - + if((r-thread_id) > 0) ppt+=1; //to include the remainder particle // Let's compute the local image // for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){ for(l=0;l 0){ - i = nth*ppt+thread_id; - xx = (int) x[i]; - yy = (int) y[i]; - tt = (int) t[i]; - mm = mass[i]; - - if(tt < 1) tt = 1; - if(tt > size_lim) tt = size_lim; - - for(j=-tt; j= 0) && ( (xx+j) < xsize) && ( (yy+k) >=0) && ( (yy+k) < ysize)){ - local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt); - } - } - } - } // Let's merge the local images From 5324ee7c9e1fa10d979a9702e030e21ff47f066f Mon Sep 17 00:00:00 2001 From: Adrien Thob Date: Tue, 1 Jan 2019 18:23:28 +0000 Subject: [PATCH 2/7] Update scenemodule.c Simplified the kview computation to get the inclusion of particles with a kernel that intersects the field of view, despite the particle being outside of it --- sphviewer/extensions/scenemodule.c | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/sphviewer/extensions/scenemodule.c b/sphviewer/extensions/scenemodule.c index ae6e488..0bd538f 100644 --- a/sphviewer/extensions/scenemodule.c +++ b/sphviewer/extensions/scenemodule.c @@ -142,16 +142,6 @@ long int compute_scene(float *x, float *y, float *z, float *hsml, ymin = extent[2]; ymax = extent[3]; - - for(i=0;i= xmin) & (x[i] <= xmax) & - (y[i] >= ymin) & (y[i] <= ymax) ) { - - kview[idx] = i; - idx += 1; - } - } - #pragma omp parallel for firstprivate(n,xmin,xmax,ymin,ymax,xsize,ysize) for(i=0;i 0) & - (fabsf(x[i]) <= fabsf(zpart)*xmax/zoom) & - (fabsf(y[i]) <= fabsf(zpart)*ymax/zoom) ) { - - kview[idx] = i; - idx += 1; - } - } - #pragma omp parallel for firstprivate(n,xmin,xmax,ymin,ymax,xsize,ysize,lbin,zoom,r,zpart) for(i=0;i= -hsml[i]) & (x[i] <= xsize+hsml[i]) & + (y[i] >= -hsml[i]) & (y[i] <= ysize+hsml[i]) ) { //this considers edge-particles with intersecting kernels + + kview[idx] = i; + idx += 1; + } + } return idx; } From f08467461a0dec17b4c603cbedf311607d12691e Mon Sep 17 00:00:00 2001 From: Adrien Thob Date: Wed, 2 Jan 2019 12:26:47 +0000 Subject: [PATCH 3/7] Update scenemodule.c Correct kview computation to remove the inclusion of particles with hsml[i] < 0 (equivalent to the old zpart > 0 condition) --- sphviewer/extensions/scenemodule.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sphviewer/extensions/scenemodule.c b/sphviewer/extensions/scenemodule.c index 0bd538f..63e72eb 100644 --- a/sphviewer/extensions/scenemodule.c +++ b/sphviewer/extensions/scenemodule.c @@ -181,7 +181,8 @@ long int compute_scene(float *x, float *y, float *z, float *hsml, } for(i=0;i= -hsml[i]) & (x[i] <= xsize+hsml[i]) & + if( (hsml[i] > 0) & + (x[i] >= -hsml[i]) & (x[i] <= xsize+hsml[i]) & (y[i] >= -hsml[i]) & (y[i] <= ysize+hsml[i]) ) { //this considers edge-particles with intersecting kernels kview[idx] = i; From c04b46042beb3b373bb938710baa2e26ef89ad83 Mon Sep 17 00:00:00 2001 From: Adrien Thob Date: Thu, 3 Jan 2019 08:08:49 +0000 Subject: [PATCH 4/7] Update rendermodule.c Detangle tt between the kernel pixel size (the tt variable) with the kernel projected size (the new tt_f float variable), allowing kernels with sizes under the resolution to be computed properly --- sphviewer/extensions/rendermodule.c | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/sphviewer/extensions/rendermodule.c b/sphviewer/extensions/rendermodule.c index 05bb5df..d9a385f 100644 --- a/sphviewer/extensions/rendermodule.c +++ b/sphviewer/extensions/rendermodule.c @@ -77,7 +77,7 @@ void c_render(float *x, float *y, float *t, float *mass, float *local_image; int i,j,k,l; int xx, yy, tt; - float mm; + float tt_f, mm; int r, nth, ppt, thread_id; nth = omp_get_num_threads(); // Get number of threads @@ -103,18 +103,19 @@ void c_render(float *x, float *y, float *t, float *mass, i = thread_id+nth*l; xx = (int) x[i]; yy = (int) y[i]; - tt = (int) t[i]; + tt_f = t[i]; + tt = (int) ceil(tt_f); mm = mass[i]; - + if(tt < 1) tt = 1; if(tt > size_lim) tt = size_lim; - + // Let's compute the convolution with the Kernel for(j=-tt; j= 0) && ( (xx+j) < xsize) && ( (yy+k) >=0) && ( (yy+k) < ysize)){ - local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt); - } + if( ( (xx+j) >= 0) && ( (xx+j) < xsize) && ( (yy+k) >=0) && ( (yy+k) < ysize)){ + local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt_f); + } } } } From 824db23a628b7320e5a06aedc0a5b29fa079e5a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C2=96=C2=96=C2=96Adrien=20Thob?= Date: Fri, 4 Jan 2019 23:12:37 +0000 Subject: [PATCH 5/7] Introducing the projections feature with the fish-eye camera case --- sphviewer/Camera.py | 6 +- sphviewer/Render.py | 10 ++- sphviewer/Scene.py | 32 +++++++- sphviewer/extensions/rendermodule.c | 115 +++++++++++++++++++++------- sphviewer/extensions/scenemodule.c | 27 ++++++- 5 files changed, 154 insertions(+), 36 deletions(-) diff --git a/sphviewer/Camera.py b/sphviewer/Camera.py index c21fd46..a6d7cc4 100644 --- a/sphviewer/Camera.py +++ b/sphviewer/Camera.py @@ -27,11 +27,11 @@ class Camera(object): def __init__(self, x = None, y = None, z = None, r = None, t = None, p = None, zoom = None, roll = None, - xsize = None, ysize = None, extent = None): + xsize = None, ysize = None, extent = None, projection = None): self._name = 'CAMERA' self.__params = {'x':x,'y':y,'z':z,'r':r, 't':t,'p':p,'zoom':zoom, 'roll':roll, - 'xsize':xsize, 'ysize':ysize, 'extent':extent} + 'xsize':xsize, 'ysize':ysize, 'extent':extent, 'projection': projection} def get_params(self): return self.__params @@ -115,4 +115,4 @@ def set_autocamera(self,Particles, mode='minmax'): self.__params = {'x':xmean,'y':ymean,'z':zmean,'r':r, 't':0,'p':0,'zoom':1, 'roll':0, - 'xsize':500, 'ysize':500, 'extent':None} + 'xsize':500, 'ysize':500, 'extent':None, 'projection': None} diff --git a/sphviewer/Render.py b/sphviewer/Render.py index ef8d0ac..debecb1 100644 --- a/sphviewer/Render.py +++ b/sphviewer/Render.py @@ -86,9 +86,15 @@ def __init__(self,Scene): def __make_render(self,x,y,t,kview,xsize,ysize): from .extensions import render - + + if(self.Scene.Camera.get_params()['r'] == 'infinity'): + projection = 0 + elif(self.Scene.Camera.get_params()['projection'] == 'fisheye'): + projection = 2 + else: + projection = 1 image = render.render(self.Scene._x, self.Scene._y, self.Scene._hsml, - self.Scene._m,np.int32(xsize),np.int32(ysize)) + self.Scene._m,np.int32(xsize),np.int32(ysize),np.int32(projection),np.float32(self.Scene.get_extent())) return np.reshape(image,[ysize,xsize]) diff --git a/sphviewer/Scene.py b/sphviewer/Scene.py index cddead8..a01f9da 100644 --- a/sphviewer/Scene.py +++ b/sphviewer/Scene.py @@ -188,9 +188,39 @@ def __compute_scene(self): else: self.__extent = np.float32(self._camera_params['extent']) else: - projection = 1 rcam = np.float32(self._camera_params['r']) self.__extent = np.array([0,0,0,0],dtype=np.float32) + if(isinstance(self._camera_params['projection'],str)): + if('fisheye' in self._camera_params['projection']): + projection = 2 + try: + shift = float(self._camera_params['projection'][7:]) + theta_rad=np.pi*np.array(theta)/180. + phi_rad=np.pi*np.array(phi)/180. + roll_rad=np.pi*np.array(roll)/180. + kaxis = np.array([np.cos(phi_rad)*np.cos(roll_rad), + np.cos(theta_rad)*np.sin(roll_rad)+np.sin(theta_rad)*np.sin(phi_rad)*np.cos(roll_rad), + np.sin(theta_rad)*np.sin(roll_rad)-np.cos(theta_rad)*np.sin(phi_rad)*np.cos(roll_rad)]).T + Kchange = np.cross(kaxis,np.eye(3)) + Rchange = np.eye(3)-np.sin(shift*np.pi/180)*Kchange+(1-np.cos(shift*np.pi/180))*np.dot(Kchange,Kchange) + zaxis = np.array([np.sin(phi_rad),-np.sin(theta_rad)*np.cos(phi_rad),np.cos(theta_rad)*np.cos(phi_rad)]).T + jaxis = np.matmul(Rchange,zaxis) + xcopy,ycopy,zcopy = np.array([xcam,ycam,zcam])-rcam*(zaxis-jaxis) + iaxis = np.cross(jaxis,kaxis) + tcopy_rad = np.arctan2(-jaxis[1],jaxis[2]); tcopy = tcopy_rad/np.pi*180 + cossintcopy=np.array([np.cos(tcopy_rad),np.sin(tcopy_rad)]) + rollcopy_rad = np.arctan2(-iaxis[0],kaxis[0]); rollcopy = rollcopy_rad/np.pi*180 + cossinrollcopy=np.array([np.cos(rollcopy_rad),np.sin(rollcopy_rad)]) + pcopy_rad=np.arctan2(jaxis[0],np.nanmean([([kaxis[0],-iaxis[0]]/cossinrollcopy)[np.argmax(np.abs(cossintcopy))], + ([jaxis[2],-jaxis[1]]/cossintcopy)[np.argmax(np.abs(cossintcopy))]])) + pcopy=pcopy_rad/np.pi*180 + xcam,ycam,zcam,theta,phi,roll = xcopy,ycopy,zcopy,tcopy,pcopy,rollcopy + except: + None + else: + raise ValueError("'"+self._camera_params['projection']+"' isn't a valid projection type") + else: + projection = 1 xx,yy,hh,kk = scene.scene(pos[:,0],pos[:,1], pos[:,2],hsml, diff --git a/sphviewer/extensions/rendermodule.c b/sphviewer/extensions/rendermodule.c index d9a385f..be23e60 100644 --- a/sphviewer/extensions/rendermodule.c +++ b/sphviewer/extensions/rendermodule.c @@ -58,9 +58,19 @@ float cubic_kernel(float r, float h){ return sigma*func/(h*h*h); } +float sinc(float x){ + float value; + + if(x==0) + value = 1.0; + else + value = sin(x)/x; + return value; +} + void c_render(float *x, float *y, float *t, float *mass, - int xsize, int ysize, int n, float *image){ + int xsize, int ysize, int n, int projection, float *extent, float *image){ // C function calculating the image of the particles convolved with our kernel int size_lim; @@ -95,27 +105,71 @@ void c_render(float *x, float *y, float *t, float *mass, local_image[k*xsize+j] = 0.0; } } - - if((r-thread_id) > 0) ppt+=1; //to include the remainder particle - // Let's compute the local image - // for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){ - for(l=0;l size_lim) tt = size_lim; + // Dome projection/Fish-eye camera : Azimuth Equidistant projection distortions + if(projection==2){ + float xxrad, yyrad, ttrad, rho, sincrho, cosrho; + if((r-thread_id) > 0) ppt+=1; //to include the remainder particle + // Let's compute the local image + // for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){ + for(l=0;l size_lim) tt = size_lim; - // Let's compute the convolution with the Kernel - for(j=-tt; j=0) && (jxx=0) && (kyy=1) cosd = 1; else if(cosd<=-1) cosd = -1; + local_image[kyy*xsize+jxx] += mm*cubic_kernel(acos(cosd), ttrad); //I should normalize here by the surface area + } + } + } + } + } + // + else{ + if((r-thread_id) > 0) ppt+=1; //to include the remainder particle + // Let's compute the local image + // for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){ + for(l=0;l size_lim) tt = size_lim; + + // Let's compute the convolution with the Kernel + for(j=-tt; j= 0) && ( (xx+j) < xsize) && ( (yy+k) >=0) && ( (yy+k) < ysize)){ local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt_f); } + } } } } @@ -138,10 +192,10 @@ void c_render(float *x, float *y, float *t, float *mass, void test_C(){ // This function if for testing purposes only. It writes a file called image_test.bin - int *x, *y, *t; - float *mass; - int xsize, ysize, n; - float *image; + int *x, *y; + float *t, *mass; + int xsize, ysize, n, projection; + float *extent, *image; int i; xsize = 1000; @@ -150,8 +204,9 @@ void test_C(){ x = (int *)malloc( n * sizeof(int) ); y = (int *)malloc( n * sizeof(int) ); - t = (int *)malloc( n * sizeof(int) ); + t = (float *)malloc( n * sizeof(int) ); mass = (float *)malloc( n * sizeof(float) ); + image = (float *)malloc( 4 * sizeof(float) ); image = (float *)malloc( xsize * ysize * sizeof(float) ); srand( time(NULL) ); @@ -163,7 +218,7 @@ void test_C(){ mass[i] = rand() % 499; } - c_render(x,y,t,mass,xsize,ysize,n,image); + c_render(x,y,t,mass,xsize,ysize,n,projection,extent,image); FILE *output; @@ -176,15 +231,15 @@ void test_C(){ static PyObject *rendermodule(PyObject *self, PyObject *args){ PyArrayObject *x_obj, *y_obj, *t_obj; - PyArrayObject *m_obj; + PyArrayObject *m_obj, *extent_obj; float *x, *y, *t; float *mass; int xsize, ysize; - int n; - float *image; + int n, projection; + float *extent, *image; int DOUBLE = 0; - if(!PyArg_ParseTuple(args, "OOOOii",&x_obj, &y_obj, &t_obj, &m_obj, &xsize, &ysize)) + if(!PyArg_ParseTuple(args, "OOOOiiiO",&x_obj, &y_obj, &t_obj, &m_obj, &xsize, &ysize, &projection, &extent_obj)) return NULL; // Let's check the size of the 1-dimensions arrays. @@ -209,6 +264,8 @@ static PyObject *rendermodule(PyObject *self, PyObject *args){ return NULL; } + extent = (float *)extent_obj->data; + image = (float *)malloc( xsize * ysize * sizeof(float) ); int i, j; @@ -219,7 +276,7 @@ static PyObject *rendermodule(PyObject *self, PyObject *args){ } // Here we do the work - c_render(x,y,t,mass,xsize,ysize,n,image); + c_render(x,y,t,mass,xsize,ysize,n,projection,extent,image); if(DOUBLE) free(mass); diff --git a/sphviewer/extensions/scenemodule.c b/sphviewer/extensions/scenemodule.c index 63e72eb..339dfa4 100644 --- a/sphviewer/extensions/scenemodule.c +++ b/sphviewer/extensions/scenemodule.c @@ -151,7 +151,7 @@ long int compute_scene(float *x, float *y, float *z, float *hsml, } // If the camera is not at the infinity, let's put it at a certain // distance from the object. - else + else if(projection == 1) { float FOV = 2.0*fabsf(atanf(1.0/zoom)); xmax = 1.0; @@ -179,6 +179,31 @@ long int compute_scene(float *x, float *y, float *z, float *hsml, hsml[i] = (hsml[i]/zpart)/(xmax-xmin) * xsize; } } + // Dome projection/Fish-eye camera + else if(projection == 2) + { + float xfovmax = 2.0*fabsf(atanf(1.0/zoom))/3.141592; + float xfovmin = -xfovmax; + float yfovmax = 0.5*(xfovmax-xfovmin)*ysize/xsize; + float yfovmin = -yfovmax; + + extent[0] = xfovmin; + extent[1] = xfovmax; + extent[2] = yfovmin; + extent[3] = yfovmax; + + + float zpart, radius, polar; +#pragma omp parallel for firstprivate(n,xfovmin,xfovmax,yfovmin,yfovmax,xsize,ysize,lbin,zoom,r,zpart,radius,polar) + for(i=0;i 0) & From 5e0da149f81af61c3b44358b4d2d19a141fcc6da Mon Sep 17 00:00:00 2001 From: Alejandro Benitez-Llambay Date: Wed, 5 Jun 2019 14:11:06 +0100 Subject: [PATCH 6/7] Update setup.py Adding "-ffast-math" compiler flag by default. This speeds up the code by more than a factor of 3. --- setup.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/setup.py b/setup.py index ef83800..1cb1646 100644 --- a/setup.py +++ b/setup.py @@ -7,15 +7,15 @@ import numpy as np module_scene = Extension('sphviewer/extensions/scene', sources = ['sphviewer/extensions/scenemodule.c'], - extra_compile_args=['-fopenmp','-w', '-std=c99'], + extra_compile_args=['-fopenmp','-w', '-std=c99', '-ffast-math'], extra_link_args=['-lgomp']) module_render = Extension('sphviewer/extensions/render', sources = ['sphviewer/extensions/rendermodule.c'], - extra_compile_args=['-fopenmp','-w', '-std=c99'], + extra_compile_args=['-fopenmp','-w', '-std=c99', '-ffast-math'], extra_link_args=['-lgomp']) module_makehsv = Extension('sphviewer/tools/makehsv', sources = ['sphviewer/tools/makehsvmodule.c'], - extra_compile_args=['-fopenmp','-w', '-std=c99'], + extra_compile_args=['-fopenmp','-w', '-std=c99', '-ffast-math'], extra_link_args=['-lgomp']) From 5a1c157261aa37ba1eab7d9ae2f24428816e8a05 Mon Sep 17 00:00:00 2001 From: Adrien Thob Date: Tue, 3 Dec 2019 01:53:42 +0000 Subject: [PATCH 7/7] Corrected cause of 'segmentation fault (core dumped)' errors --- .gitignore | 2 + sphviewer/extensions/rendermodule.c | 99 +++++++++++++++-------------- 2 files changed, 54 insertions(+), 47 deletions(-) diff --git a/.gitignore b/.gitignore index 3ff90a6..deba437 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,8 @@ PKG-INFO requires.txt SOURCES.txt top_level.txt +build/ +dist/ # Python bytecode *.pyc diff --git a/sphviewer/extensions/rendermodule.c b/sphviewer/extensions/rendermodule.c index 66e7c88..093126c 100644 --- a/sphviewer/extensions/rendermodule.c +++ b/sphviewer/extensions/rendermodule.c @@ -126,57 +126,62 @@ void c_render(float *x, float *y, float *t, float *mass, mm = mass[i]; if(tt <= 1) { - local_image[yy*xsize+xx] += mm; - continue; - } - if(tt > size_lim) tt = size_lim; - - int jxx,kyy; - float jxxrad, kyyrad, rhopix, cosd; - // Let's compute the convolution with the Kernel - for(j=-tt; j=0) && (jxx=0) && (kyy=1) cosd = 1; else if(cosd<=-1) cosd = -1; - local_image[kyy*xsize+jxx] += mm*cubic_kernel(acos(cosd), ttrad); //I should normalize here by the surface area - } + if( (rho <= 3.141592) && ( xx >= 0) && ( xx < xsize) && ( yy >= 0) && ( yy < ysize)){ + local_image[yy*xsize+xx] += mm; } - } + } + else{ + if(tt > size_lim) tt = size_lim; + + int jxx,kyy; + float jxxrad, kyyrad, rhopix, cosd; + // Let's compute the convolution with the Kernel + for(j=-tt; j=0) && (jxx=0) && (kyy=1) cosd = 1; else if(cosd<=-1) cosd = -1; + local_image[kyy*xsize+jxx] += mm*cubic_kernel(acos(cosd), ttrad); //I should normalize here by the surface area } + } + } + } + } + } + // + else{ + if((r-thread_id) > 0) ppt+=1; //to include the remainder particle + // Let's compute the local image + // for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){ + for(l=0;l= 0) && ( xx < xsize) && ( yy >= 0) && ( yy < ysize)){ + local_image[yy*xsize+xx] += mm; + } + } + else{ + if(tt > size_lim) tt = size_lim; + + // Let's compute the convolution with the Kernel + for(j=-tt; j= 0) && ( (xx+j) < xsize) && ( (yy+k) >= 0) && ( (yy+k) < ysize)){ + local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt_f); + } } - // - else{ - if((r-thread_id) > 0) ppt+=1; //to include the remainder particle - // Let's compute the local image - // for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){ - for(l=0;l size_lim) tt = size_lim; - - // Let's compute the convolution with the Kernel - for(j=-tt; j= 0) && ( (xx+j) < xsize) && ( (yy+k) >=0) && ( (yy+k) < ysize)){ - local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt_f); - } - } + } } } }