| commit | author | age | ||
| 0af0cf | 1 | #include<general.h> |
| SP | 2 | #include<cross-section.h> |
| 3 | #include<coord.h> | |
| a61c00 | 4 | #include<cairo/cairo.h> |
| 0af0cf | 5 | /** @brief Calculates cross-section of vesicle with plane. |
| SP | 6 | * |
| 7 | * 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: | |
| 8 | * | |
| 9 | * if a distance of point to plane (given by equation $D=\frac{ax_0+by_0+cz_0+d}{\sqrt{a^2+b^2+c^2}}$, where $x_0$, $y_0$ and $z_0$ are coordinates of a given vertex) is less than maximal allowed distance between vertices {\tt sqrt(vesicle->dmax)} than vertex is a candidate for crossection calculation. | |
| 10 | * | |
| 11 | */ | |
| d27c07 | 12 | ts_coord_list *get_crossection_with_plane(ts_vesicle *vesicle,ts_double a,ts_double b,ts_double c, ts_double d){ |
| 0af0cf | 13 | |
| SP | 14 | |
| ad4e70 | 15 | ts_uint i, j,k,l; |
| 0af0cf | 16 | ts_double pp,Dsq; // distance from the plane squared |
| ad4e70 | 17 | ts_double ppn1; // distance from the plane squared of a neighbor |
| 0af0cf | 18 | ts_vertex *vtx; |
| a61c00 | 19 | ts_uint ntria=0; // number triangles |
| SP | 20 | ts_triangle *tria[2]; // list of triangles |
| 0af0cf | 21 | ts_coord_list *pts=init_coord_list(); |
| d27c07 | 22 | for(i=0;i<vesicle->vlist->n;i++){ |
| 0af0cf | 23 | vtx=vesicle->vlist->vtx[i]; |
| SP | 24 | |
| 25 | pp=vtx->x*a+vtx->y*b+vtx->z*c+d; | |
| 26 | Dsq=pp*pp/(a*a+b*b+c*c); | |
| 27 | if(Dsq<vesicle->dmax){ | |
| 28 | for(j=0;j<vtx->neigh_no;j++){ | |
| a61c00 | 29 | ppn1=vtx->neigh[j]->x*a+vtx->neigh[j]->y*b+vtx->neigh[j]->z*c+d; |
| SP | 30 | if(pp*ppn1<=0){ //the combination of vertices are good candidates for a crossection |
| 31 | //find triangle that belongs to the two vertices | |
| 32 | ntria=0; | |
| 33 | for(k=0;k<vtx->tristar_no;k++){ | |
| ad4e70 | 34 | if(vtx->tristar[k]->vertex[0]==vtx && ( vtx->tristar[k]->vertex[1]==vtx->neigh[j] || vtx->tristar[k]->vertex[2]==vtx->neigh[j]) ){ |
| a61c00 | 35 | //triangle found. |
| SP | 36 | tria[ntria]=vtx->tristar[k]; |
| 37 | ntria++; | |
| 38 | } | |
| 39 | } | |
| 40 | // if ntria !=1 there is probably something wrong I would say... | |
| 41 | if(ntria==0) continue; | |
| ad4e70 | 42 | /* if(ntria!=1) { //there should be 2 triangles. of course, all some of them will be doubled. |
| a61c00 | 43 | fprintf(stderr,"ntria=%u\n",ntria); |
| ad4e70 | 44 | fatal ("Error in mesh. 2 triangles not found",123123); |
| SP | 45 | } */ |
| a61c00 | 46 | //find the two intersections (in general) to form a intersection line |
| ad4e70 | 47 | for(l=0;l<ntria;l++){ |
| SP | 48 | //we add intersection line between two points for each of the triangles found above. |
| 49 | add_crosssection_point(pts,a,b,c,d,tria[l]->vertex[0], tria[l]->vertex[1]); | |
| 50 | add_crosssection_point(pts,a,b,c,d,tria[l]->vertex[0], tria[l]->vertex[2]); | |
| 51 | add_crosssection_point(pts,a,b,c,d,tria[l]->vertex[1], tria[l]->vertex[2]); | |
| a61c00 | 52 | |
| SP | 53 | } |
| 0af0cf | 54 | } |
| SP | 55 | } |
| 56 | } | |
| 57 | } | |
| 58 | return pts; | |
| 59 | } | |
| 60 | ||
| 61 | ||
| ad4e70 | 62 | |
| SP | 63 | ts_bool add_crosssection_point(ts_coord_list *pts, ts_double a, ts_double b, ts_double c, ts_double d, ts_vertex *vtx1, ts_vertex *vtx2){ |
| 64 | ts_double pp=vtx1->x*a+vtx1->y*b+vtx1->z*c+d; | |
| 65 | ts_double pp2=vtx2->x*a+vtx1->y*b+vtx2->z*c+d; | |
| 66 | if(pp*pp2<=0){ | |
| 67 | ts_double u=pp/(a*(vtx1->x-vtx2->x)+b*(vtx1->y-vtx2->y)+c*(vtx1->z-vtx2->z)); | |
| 68 | add_coord(pts, vtx1->x+u*(vtx2->x - vtx1->x), | |
| 69 | vtx1->y+u*(vtx2->y - vtx1->y), | |
| 70 | vtx1->z+u*(vtx2->z - vtx1->z), | |
| 71 | TS_COORD_CARTESIAN); | |
| 72 | ||
| 73 | return TS_SUCCESS; | |
| 74 | } else { | |
| 75 | return TS_FAIL; | |
| 76 | } | |
| 77 | } | |
| 78 | ||
| a61c00 | 79 | /** Saves calculated crossection as a png image */ |
| SP | 80 | ts_bool crossection_to_png(ts_coord_list *pts, char *filename){ |
| 81 | ||
| 82 | cairo_surface_t *surface; | |
| 83 | cairo_t *cr; | |
| 84 | ts_uint i; | |
| ad4e70 | 85 | surface = cairo_image_surface_create (CAIRO_FORMAT_RGB24, 1800, 1800); |
| a61c00 | 86 | cr = cairo_create (surface); |
| SP | 87 | cairo_rectangle(cr, 0.0, 0.0, 1800,1800); |
| 88 | cairo_set_source_rgb(cr, 0.3, 0.3, 0.3); | |
| 89 | cairo_fill(cr); | |
| 90 | cairo_set_line_width (cr, 5.0/30.0); | |
| 324ad7 | 91 | cairo_set_line_cap(cr, CAIRO_LINE_CAP_ROUND); |
| SP | 92 | cairo_set_line_join(cr, CAIRO_LINE_JOIN_ROUND); |
| a61c00 | 93 | cairo_translate(cr, 900,900); |
| SP | 94 | cairo_scale (cr, 30, 30); |
| 95 | cairo_set_source_rgb (cr, 1.0, 1.0, 1.0); | |
| 96 | ||
| 97 | for(i=0;i<pts->n;i+=2){ | |
| 98 | cairo_move_to(cr, pts->coord[i]->e1, pts->coord[i]->e2); | |
| 99 | cairo_line_to(cr, pts->coord[i+1]->e1, pts->coord[i+1]->e2); | |
| 100 | } | |
| 101 | cairo_stroke(cr); | |
| 102 | cairo_surface_write_to_png (surface,filename); | |
| 103 | cairo_surface_finish (surface); | |
| 104 | return TS_SUCCESS; | |
| 105 | } | |
| ad4e70 | 106 | |
| SP | 107 | |
| 108 | ts_bool save_crossection_snapshot(ts_coord_list *pts, ts_uint timestepno){ | |
| 109 | char filename[255]; | |
| 110 | sprintf(filename,"timestep_%.6u.png",timestepno); | |
| 111 | crossection_to_png(pts,filename); | |
| 112 | return TS_SUCCESS; | |
| 113 | } | |
