It works, png is output. The only problem is that some triangles are missing

Samo Penic

2014-06-13 a61c001cd35ff70a314ef417c4beda9c7e68d3ad

It works, png is output. The only problem is that some triangles are missing

 config.h

@@ -1,9 +1,24 @@
/* config.h.  Generated from config.h.in by configure.  */
/* config.h.in.  Generated from configure.ac by autoheader.  */

/* Define to 1 if you have the <cairo/cairo.h> header file. */
#define HAVE_CAIRO_CAIRO_H 1

/* Define to 1 if you have the <cairo/cairo-png.h> header file. */
/* #undef HAVE_CAIRO_CAIRO_PNG_H */

/* Define to 1 if you have the <cairo.h> header file. */
/* #undef HAVE_CAIRO_H */

/* Define to 1 if you have the <cairo-png.h> header file. */
/* #undef HAVE_CAIRO_PNG_H */

/* Define to 1 if you have the <inttypes.h> header file. */
#define HAVE_INTTYPES_H 1

/* Define to 1 if you have the `cairo' library (-lcairo). */
#define HAVE_LIBCAIRO 1

/* Define to 1 if you have the `confuse' library (-lconfuse). */
#define HAVE_LIBCONFUSE 1


 config.h.in

@@ -1,8 +1,23 @@
/* config.h.in.  Generated from configure.ac by autoheader.  */

/* Define to 1 if you have the <cairo/cairo.h> header file. */
#undef HAVE_CAIRO_CAIRO_H

/* Define to 1 if you have the <cairo/cairo-png.h> header file. */
#undef HAVE_CAIRO_CAIRO_PNG_H

/* Define to 1 if you have the <cairo.h> header file. */
#undef HAVE_CAIRO_H

/* Define to 1 if you have the <cairo-png.h> header file. */
#undef HAVE_CAIRO_PNG_H

/* Define to 1 if you have the <inttypes.h> header file. */
#undef HAVE_INTTYPES_H

/* Define to 1 if you have the `cairo' library (-lcairo). */
#undef HAVE_LIBCAIRO

/* Define to 1 if you have the `confuse' library (-lconfuse). */
#undef HAVE_LIBCONFUSE


 configure.ac

@@ -14,6 +14,17 @@
AC_CHECK_LIB([confuse], [cfg_parse])
# FIXME: Replace `main' with a function in `-lm':
AC_CHECK_LIB([m], [pow])
# Checks for libcairo

## cairographics needs pkg-config
#AC_PATH_PROG(PKGCONF, pkg-config , [], 
#         [$PATH:/usr/local/bin:/ext/bin:/ext:/sw/bin:/opt/bin])
AC_CHECK_HEADERS([cairo/cairo.h cairo/cairo-png.h])
AC_CHECK_HEADERS([cairo.h cairo-png.h])
AC_CHECK_LIB([cairo],[cairo_move_to])
#PKG_CHECK_MODULES(CAIRO, cairo)
#AC_SUBST(CAIRO_CFLAGS)
#AC_SUBST(CAIRO_LIBS)

AC_CHECK_LIB([gslcblas],[cblas_dgemm])
AC_CHECK_LIB([gsl],[gsl_blas_dgemm])

 src/cross-section.c

@@ -1,7 +1,7 @@
#include<general.h>
#include<cross-section.h>
#include<coord.h>

#include<cairo/cairo.h>
/** @brief Calculates cross-section of vesicle with plane.
 *
 * Function returns points of cross-section of vesicle with plane. Plane is described with equation $ax+by+cz+d=0$. Algorithm extracts coordinates of each vertex of a vesicle and then:
@@ -12,12 +12,14 @@
ts_coord_list *get_crossection_with_plane(ts_vesicle *vesicle,ts_double a,ts_double b,ts_double c, ts_double d){


    ts_uint i, j;
    ts_uint i, j,k;
    ts_double pp,Dsq; // distance from the plane squared
    ts_double ppn; // distance from the plane squared of a neighbor
    ts_double ppn1,ppn2; // distance from the plane squared of a neighbor
    ts_double u; //factor to scale vector from first vector to the second to get intersection
    ts_vertex *vtx;

    ts_uint ntria=0; // number triangles
    ts_uint cnt=0;
    ts_triangle *tria[2]; // list of triangles
    ts_coord_list *pts=init_coord_list();
    for(i=0;i<vesicle->vlist->n;i++){
        vtx=vesicle->vlist->vtx[i];
@@ -26,13 +28,67 @@
        Dsq=pp*pp/(a*a+b*b+c*c);
        if(Dsq<vesicle->dmax){
            for(j=0;j<vtx->neigh_no;j++){
                ppn=vtx->neigh[j]->x*a+vtx->neigh[j]->y*b+vtx->neigh[j]->z*c+d;
                if(pp*ppn<0){ //the combination of vertices are good candidates for a crossection
                    u=pp/(a*(vtx->x-vtx->neigh[j]->x)+b*(vtx->y-vtx->neigh[j]->y)+c*(vtx->z-vtx->neigh[j]->z));
                ppn1=vtx->neigh[j]->x*a+vtx->neigh[j]->y*b+vtx->neigh[j]->z*c+d;
                if(pp*ppn1<=0){ //the combination of vertices are good candidates for a crossection
/*                    u=pp/(a*(vtx->x-vtx->neigh[j]->x)+b*(vtx->y-vtx->neigh[j]->y)+c*(vtx->z-vtx->neigh[j]->z));
                    add_coord(pts, vtx->x+u*(vtx->neigh[j]->x - vtx->x),
                            vtx->y+u*(vtx->neigh[j]->y - vtx->y),        
                            vtx->z+u*(vtx->neigh[j]->z - vtx->z),
                            TS_COORD_CARTESIAN);        

*/
                    //find triangle that belongs to the two vertices
                    cnt=0;
                    ntria=0;
                    for(k=0;k<vtx->tristar_no;k++){
                        if(vtx->tristar[k]->vertex[1]==vtx->neigh[j] && vtx->tristar[k]->vertex[0]==vtx){ 
                            //triangle found.
                            tria[ntria]=vtx->tristar[k];
                            ntria++;
                        }
                    }
                    // if ntria !=1 there is probably something wrong I would say...
                    if(ntria==0) continue;
                    if(ntria!=1) {
                        fprintf(stderr,"ntria=%u\n",ntria);
                        fatal ("Error in mesh. 1 triangle not found",123123);
                    }
                    //find the two intersections (in general) to form a intersection line
                    /* we dont know which vertex is which, so we need to recalculate all previously calculated values*/
                    pp=vtx->x*a+vtx->y*b+vtx->z*c+d;
                    ppn1=vtx->neigh[j]->x*a+vtx->neigh[j]->y*b+vtx->neigh[j]->z*c+d;

                    u=pp/(a*(vtx->x-vtx->neigh[j]->x)+b*(vtx->y-vtx->neigh[j]->y)+c*(vtx->z-vtx->neigh[j]->z));
                    add_coord(pts, vtx->x+u*(vtx->neigh[j]->x - vtx->x),
                            vtx->y+u*(vtx->neigh[j]->y - vtx->y),		
                            vtx->z+u*(vtx->neigh[j]->z - vtx->z),
                            TS_COORD_CARTESIAN);
                    ppn2=tria[0]->vertex[2]->x*a+tria[0]->vertex[2]->y*b+tria[0]->vertex[2]->z*c+d;	
                    cnt++;
                    //two more tries to find anothen pair of vertices, one above and one below the intersection plane */	
                    if(pp*ppn2<=0) {
                        u=pp/(a*(vtx->x-tria[0]->vertex[2]->x)+b*(vtx->y-tria[0]->vertex[2]->y)+c*(vtx->z-tria[0]->vertex[2]->z));
                        add_coord(pts, vtx->x+u*(tria[0]->vertex[2]->x - vtx->x),
                            vtx->y+u*(tria[0]->vertex[2]->y - vtx->y),		
                            vtx->z+u*(tria[0]->vertex[2]->z - vtx->z),
                            TS_COORD_CARTESIAN);	
                        cnt++;
                    }

                    if(ppn2*ppn1<=0) {
                        u=ppn1/(a*(vtx->neigh[j]->x-tria[0]->vertex[2]->x)+b*(vtx->neigh[j]->y-tria[0]->vertex[2]->y)+c*(vtx->neigh[j]->z-tria[0]->vertex[2]->z));
                        add_coord(pts, vtx->neigh[j]->x+u*(tria[0]->vertex[2]->x - vtx->neigh[j]->x),
                            vtx->neigh[j]->y+u*(tria[0]->vertex[2]->y - vtx->neigh[j]->y),		
                            vtx->neigh[j]->z+u*(tria[0]->vertex[2]->z - vtx->neigh[j]->z),
                            TS_COORD_CARTESIAN);	
                        cnt++;

                    }


                    if(cnt!=2){
                        fprintf(stderr,"Hey, cnt=%u",cnt);
                    }
                }
            }
        }
@@ -41,3 +97,28 @@
}


/** Saves calculated crossection as a png image */
ts_bool crossection_to_png(ts_coord_list *pts, char *filename){

    cairo_surface_t *surface;
    cairo_t *cr;
    ts_uint i;
    surface = cairo_image_surface_create (CAIRO_FORMAT_ARGB32, 1800, 1800);
    cr = cairo_create (surface);
    cairo_rectangle(cr, 0.0, 0.0, 1800,1800);
    cairo_set_source_rgb(cr, 0.3, 0.3, 0.3);
    cairo_fill(cr);
    cairo_set_line_width (cr, 5.0/30.0);
    cairo_translate(cr, 900,900);
    cairo_scale (cr, 30, 30);
    cairo_set_source_rgb (cr, 1.0, 1.0, 1.0);

    for(i=0;i<pts->n;i+=2){
        cairo_move_to(cr, pts->coord[i]->e1, pts->coord[i]->e2);
        cairo_line_to(cr, pts->coord[i+1]->e1, pts->coord[i+1]->e2);
    }
    cairo_stroke(cr);
    cairo_surface_write_to_png (surface,filename);
    cairo_surface_finish (surface);	
    return TS_SUCCESS;
}

 src/cross-section.h

@@ -1,5 +1,5 @@
#ifndef _H_CROSS_SECTION
#define _H_CROSS_SECTION
ts_coord_list *get_crossection_with_plane(ts_vesicle *vesicle,ts_double a,ts_double b,ts_double c, ts_double d);
ts_bool crossection_to_png(ts_coord_list *pts, char *filename);
#endif


 src/timestep.c

@@ -76,6 +76,7 @@
            for(k=0;k<pts->n;k++){
                fprintf(fdx,"%e, %e, %e\n",pts->coord[k]->e1, pts->coord[k]->e2, pts->coord[k]->e3);
            }
            crossection_to_png(pts, "test.png");
            fclose(fdx);
            free(pts);
            epochtime=get_epoch();