Laboratoire RDP » Biophysics Team » Vascular Networks in Growing Leaves

#include <string.h>

#include <stdlib.h>

#include <stdio.h>

#include "save_data.h"

#include <math.h>

#include "constants.h"

#include <gsl/gsl_matrix.h>

#include <gsl/gsl_blas.h>

#include <gsl/gsl_eigen.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_cblas.h>

#include <gsl/gsl_linalg.h>

#include "gaussian.h"

/* this function saves the tissue and its chemical and mechanical state */

void save_data(FILE *f, int vertex_number, double *vertices, int edge_number, int *edges_vertices,

               double *lengths, double *rest_lengths, double *edges_width, int *edges_cells, int *periodizer, double width, double height, int cell_number){

fprintf(f,"edge_number= %i;\n",edge_number);

fprintf(f,"vertex_number= %i;\n",vertex_number);

fprintf(f,"cell_number= %i;\n",cell_number);

fprintf(f,"width_before= %f;\n",width);

fprintf(f,"height_before= %f;\n",height);

fprintf(f,"stress_x= %f;\n",stress_x);

fprintf(f,"stress_y= %f;\n",stress_y);

fprintf(f,"sigma= %f;\n",sigma);

fprintf(f,"vertices=[\n");

for (int i=0;i<vertex_number;i++){

    fprintf(f,"%f %f\n",vertices[2*i],vertices[2*i+1]);}

fprintf(f,"];\n");

fprintf(f,"edges_vertices=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%i %i\n",edges_vertices[2*i],edges_vertices[2*i+1]);

}

fprintf(f,"];\n");

fprintf(f,"lengths=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%f\n",lengths[i]);

}

fprintf(f,"];\n");

fprintf(f,"rest_lengths=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%f\n",rest_lengths[i]);

}

fprintf(f,"];\n");

fprintf(f,"edges_width=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%f\n",edges_width[i]);

}

fprintf(f,"];\n");

fprintf(f,"periodizer=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%i %i %i %i\n",periodizer[4*i],periodizer[4*i+1],periodizer[4*i+2],periodizer[4*i+3]);

}

fprintf(f,"];\n");

fprintf(f,"edges_cells=[\n");

for (int i=0;i<edge_number;i++){

    fprintf(f,"%i %i\n",edges_cells[2*i],edges_cells[2*i+1]);

}

fprintf(f,"];\n");

double barycenters[2*cell_number]; for (int i=0; i<2*cell_number;i++){ barycenters[i]=0.;}

double count[cell_number]; for (int i=0; i<cell_number;i++){ count[i]=0.;}

double texture[3*cell_number]; for (int i=0; i<3*cell_number;i++){ texture[i]=0.;}

int box[18]={-1,-1,-1,0,-1,1,0,-1,0,0,0,1,1,-1,1,0,1,1};

double distance=0.; int j_min; double distance_min;

for (int i=0; i<edge_number; i++){

    if (count[edges_cells[2*i]]==0){

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i]]+width*(double)periodizer[4*i])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)periodizer[4*i+1])*lengths[i];

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)periodizer[4*i+2])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)periodizer[4*i+3])*lengths[i];

        count[edges_cells[2*i]]+=2.*lengths[i];

    }

    else {

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i]]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i]+1]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i]]+=1.*lengths[i];

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i]]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i+1]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i]+1]/count[edges_cells[2*i]]-vertices[2*edges_vertices[2*i+1]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i]]+=1.*lengths[i];

    }

    if (count[edges_cells[2*i+1]]==0){

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i]]+width*(double)periodizer[4*i])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)periodizer[4*i+1])*lengths[i];

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)periodizer[4*i+2])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)periodizer[4*i+3])*lengths[i];

        count[edges_cells[2*i+1]]+=2.*lengths[i];

    }

    else {

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i+1]]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i+1]+1]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i+1]]+=1.*lengths[i];

        distance_min=100000;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*edges_cells[2*i+1]]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i+1]]-width*(double)box[2*j],2.)+

                     pow(barycenters[2*edges_cells[2*i+1]+1]/count[edges_cells[2*i+1]]-vertices[2*edges_vertices[2*i+1]+1]-height*(double)box[2*j+1],2.);

            if (distance<distance_min){

                j_min=j;

                distance_min=distance;

            }

        }

        barycenters[2*edges_cells[2*i+1]]+=(vertices[2*edges_vertices[2*i+1]]+width*(double)box[2*j_min])*lengths[i];

        barycenters[2*edges_cells[2*i+1]+1]+=(vertices[2*edges_vertices[2*i+1]+1]+height*(double)box[2*j_min+1])*lengths[i];

        count[edges_cells[2*i+1]]+=1.*lengths[i];

    }

//    count[edges_cells[2*i]]+=1.;

//    count[edges_cells[2*i+1]]+=1.;

}

for (int i=0; i<cell_number; i++){

    barycenters[2*i]/=count[i];

    barycenters[2*i+1]/=count[i];

}

for (int i=0; i<cell_number;i++){ count[i]=0.;}

for (int i=0; i<edge_number; i++){

    distance_min=100000.;

    for (int j=0; j<9; j++){

        distance=pow(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j]-barycenters[2*edges_cells[2*i+1]],2.)+

                 pow(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j+1]-barycenters[2*edges_cells[2*i+1]+1],2.);

        if (distance<distance_min){

            distance_min=distance;

            j_min=j;

        }

    }

    texture[3*edges_cells[2*i]]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]]);

    texture[3*edges_cells[2*i+1]]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]]);

    texture[3*edges_cells[2*i]+1]+=(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    texture[3*edges_cells[2*i+1]+1]+=(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    texture[3*edges_cells[2*i]+2]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    texture[3*edges_cells[2*i+1]+2]+=(barycenters[2*edges_cells[2*i]]+width*(double)box[2*j_min]-barycenters[2*edges_cells[2*i+1]])*(barycenters[2*edges_cells[2*i]+1]+height*(double)box[2*j_min+1]-barycenters[2*edges_cells[2*i+1]+1]);

    count[edges_cells[2*i]]+=1.;

    count[edges_cells[2*i+1]]+=1.;

}

for (int i=0; i<cell_number; i++){

    texture[3*i]/=count[i];

    texture[3*i+1]/=count[i];

    texture[3*i+2]/=count[i];

}

double average_texture[3*cell_number]; for (int i=0; i<3*cell_number; i++){ average_texture[i]=0.;}

for (int i=0; i<cell_number; i++){ count[i]=0.;}

for (int i=0; i<cell_number; i++){

    for (int k=i; k<cell_number; k++){

        distance_min=100000.;

        for (int j=0; j<9; j++){

            distance=pow(barycenters[2*i]+width*(double)box[2*j]-barycenters[2*k],2.)+

                     pow(barycenters[2*i+1]+height*(double)box[2*j+1]-barycenters[2*k+1],2.);

            if (distance<distance_min){

                distance_min=distance;

            }

        }

        distance_min=sqrt(distance_min);

        average_texture[3*i]+=texture[3*k]*gaussian(distance_min,sigma);

        average_texture[3*i+1]+=texture[3*k+1]*gaussian(distance_min,sigma);

        average_texture[3*i+2]+=texture[3*k+2]*gaussian(distance_min,sigma);

        average_texture[3*k]+=texture[3*i]*gaussian(distance_min,sigma);

        average_texture[3*k+1]+=texture[3*i+1]*gaussian(distance_min,sigma);

        average_texture[3*k+2]+=texture[3*i+2]*gaussian(distance_min,sigma);

        count[i]+=gaussian(distance_min,sigma);

        count[k]+=gaussian(distance_min,sigma);

    }

}

for (int i=0; i<cell_number; i++){

    average_texture[3*i]/=count[i];

    average_texture[3*i+1]/=count[i];

    average_texture[3*i+2]/=count[i];

}

fprintf(f,"barycenters_before=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%f %f\n",barycenters[2*i],barycenters[2*i+1]);

}

fprintf(f,"];\n");

fprintf(f,"texture_before=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%f %f %f\n",texture[3*i],texture[3*i+1],texture[3*i+2]);

}

fprintf(f,"];\n");

fprintf(f,"average_texture_before=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%f %f %f\n",average_texture[3*i],average_texture[3*i+1],average_texture[3*i+2]);

}

fprintf(f,"];\n");

int boundary[cell_number]; for (int i=0; i<cell_number; i++){ boundary[i]=0;}

for (int i=0;i<edge_number;i++){

    if (periodizer[4*i]!=0 || periodizer[4*i+1]!=0 || periodizer[4*i+2]!=0 || periodizer[4*i+3]!=0){

        boundary[edges_cells[2*i]]=1; boundary[edges_cells[2*i+1]]=1;

    }

}

fprintf(f,"boundary=[\n");

for (int i=0;i<cell_number;i++){

    fprintf(f,"%i\n",boundary[i]);

}

fprintf(f,"];\n");

}