/* * Ray2mesh : software for geophysicists. * Compute various scores attached to the mesh cells, based on geometric   information that rays bring when the traverse the cell. * * Copyright (C) 2003, St閜hane Genaud and Marc Grunberg * * This tool is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU * Library General Public License for more details. * * You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the Free * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *//* * score.c * * Current cooking state for score magic formula : * *            score = hit_score * length_score * where : * * 1. hit_score : number of faces hexadreal cells traversed / 6 *  * 2. length_score :  *    the cell volume is approximated by a bounding box,  *    for which we compute the length d of diagonal. *    We also compute the average length l of all ays traversing the *    cell. lenght_score = l / d  */#include <math.h>#include <ray/raydescartes.h>#include <mesh/mesh.h>#include <mesh/layer.h>#include "ray2mesh.h"#include "raymaths.h"#include "const_ray2mesh.h"/*--------------------------------------------------------------------*//* get_hit_score : given a cell, compute a simple score function of   *//* hits/cell face: the best score is obtained if all faces are        *//* traversed at least once and 0<=score<=1                            *//*--------------------------------------------------------------------*/static float get_hit_score(struct cell_info_t *cell){    int i;    int hits = 0;    for (i = 0; i < NB_HIT_FACES && cell->faces_hit; i++) {	if (cell->faces_hit[i] > 0) {	    hits += 1;	}    }    return (((float) hits) / NB_HIT_FACES);}/*--------------------------------------------------------------------*//* get_center_score :                                                 *//*--------------------------------------------------------------------*/static float get_center_score(int x, int y, int z,	/* to remove later */			      double lat_min,			      double lon_min,			      int zstart,			      int zend,			      float cell_width_lat,			      float cell_width_lon,			      struct coord_geo_t *in,			      struct coord_geo_t *out){    struct coord_geo_t *m,	/* middle of segment [in,out] */     O;				/* cell center */    double ox, oy, oz;		/* cell center coords */    /* find cell center geo coordinates (in rad) */    ox = (lat_min * TO_RAD) + cell_width_lat * (0.5 + x);    oy = (lon_min * TO_RAD) + cell_width_lon * (0.5 + y);    /* O coordinates may have angles > PI ,      * convert it for a better reading*/    /* ox = fmod(ox, PI);       oy = fmod(oy, PI); */    /* layers bounds are expressed with 0 being the surface */    oz = (zend + zstart) / 2;    printf("cell[%d,%d,%d]  center O=(%.2f %.2f %.2f)\n", x, y, z, ox, oy,	   oz);    /* now get the middle of the segment m=[in;out] */    m = geo_middle(in, out);    printf	("middle=(%.2f %.2f %.2f)\nin=(%.2f %.2f %.2f) out=(%.2f %.2f %.2f)\n",	 m->lat, m->lon, m->prof, in->lat, in->lon, in->prof, out->lat,	 out->lon, out->prof);    /* ... and measures distance d=[O,m] */    O.lat = ox;    O.lon = oy;    O.prof = oz;    /* geo_distance takes geo coords in radians and depth=0 at the surface */    return ((float) geo_distance(&O, m));}/*--------------------------------------------------------------------*//* compute_score : given the cell_info array c and mesh parameters m  *//* review each cell in c, and updates the score field of each cell    *//*--------------------------------------------------------------------*/void compute_score(struct cell_info_t ****c, struct mesh_t *m){    int x, y, z;		/* for enumeration of cell ids            */    int i;    int nb_lat_cell, nb_lon_cell;    float diagonal_length, length;    enum {SCORE=0, HIT, LENGTH, DISP, NB };    struct coord_geo_t p1, p2;	/* used to compute the great diagonal */    if (!c) {	fprintf (stderr, "compute_score : no cell traversed !\n");    	return;    }    for (z = 0; z < m->nlayers; z++) {	/* number of cells in lattitude for this layer */	nb_lat_cell = m->layer[z]->nlat;	/* number of cells in longitude for this layer */	nb_lon_cell = m->layer[z]->nlon;	/* compute length of cell edges in kms to compare to ray length */	p1.lat = 0.;	p1.lon = 0.;	p1.prof =  m->layer[z]->zstart;	p2.lat  = m->layer[z]->lat_unit * m->parameter->lat_unit_size * TO_RAD;	p2.lon  = m->layer[z]->lon_unit * m->parameter->lon_unit_size * TO_RAD;	p2.prof = m->layer[z]->zend;        /* geo_distance takes geo coords in radians and depth=0 at the surface */	diagonal_length = geo_distance (&p1, &p2);#ifdef DEBUG	fprintf(stderr,		"compute_score() : cell size (layer=%d): %.3fx%.3fx%.3f (diag lenght=%.3f)\n",		z, cell_lx, cell_ly,		m->layer[z]->zend - m->layer[z]->zstart, diagonal_length);	fflush(stderr);#endif	for (x = 0; x < nb_lat_cell; x++) {	    for (y = 0; y < nb_lon_cell; y++) {		/* length per traversing ray info */		if ( c[z][x][y] && c[z][x][y]->nitems) {		    /* update hit score from face hit information */		    c[z][x][y]->score[HIT] = get_hit_score(c[z][x][y]);		    length = 0;		    /* sum_distOM = 0.; */		    for (i = 0; i < c[z][x][y]->nitems; i++) {			/* sum of P and S ray length */			length += (c[z][x][y]->item[i].P_length + c[z][x][y]->item[i].S_length);		    }		    /* dispersion score based on blocks traversed */		    c[z][x][y]->score[DISP] = (float) c[z][x][y]->nblocks / NB_CELL_BLOCKS3;		    /* average ray length */		    c[z][x][y]->score[LENGTH] = (length / c[z][x][y]->nitems) / diagonal_length;		    /* score */		    c[z][x][y]->score[SCORE] = c[z][x][y]->score[LENGTH] * c[z][x][y]->score[HIT] * 				               c[z][x][y]->score[DISP] * log10f(c[z][x][y]->nitems+1);		    /* average dist from center */		    /* (c[z][x][y]).center_score = sum_distOM / (c[z][x][y]).nitems; */		}	    }	}    }}