score.c

用于2维的射线追踪

/*
 * 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; */

		}
	    }
	}
    }
}