C++ MIN函数代码示例

void TCOD_map_compute_fov_restrictive_shadowcasting_quadrant (map_t *m, int player_x, int player_y, int max_radius, bool light_walls, int maxObstacles, int dx, int dy) {
    static double *startAngle=NULL, *endAngle=NULL;
    static int angleArraySize=0;

    if ( angleArraySize > 0 && angleArraySize < maxObstacles ) {
        free(startAngle);
        startAngle=NULL;
    }
    if ( startAngle == NULL ) {
        angleArraySize = maxObstacles;
        startAngle = (double *)malloc(sizeof(double) * 2 * maxObstacles);
        endAngle = &startAngle[maxObstacles];
    }
    //octant: vertical edge
    {
        int iteration = 1; //iteration of the algo for this octant
        bool done = false;
        int totalObstacles = 0;
        int obstaclesInLastLine = 0;
		double minAngle = 0.0f;
		int x,y;

        //do while there are unblocked slopes left and the algo is within the map's boundaries
        //scan progressive lines/columns from the PC outwards
        y = player_y+dy; //the outer slope's coordinates (first processed line)
        if (y < 0 || y >= m->height) done = true;
		while(!done) {
            //process cells in the line
			double slopesPerCell = 1.0f/(double)(iteration+1);
			double halfSlopes = slopesPerCell*0.5f;
			int processedCell = (int)(minAngle / slopesPerCell);
            int minx = MAX(0,player_x-iteration), maxx = MIN(m->width-1,player_x+iteration);
            done = true;
            for (x = player_x + (processedCell * dx); x >= minx && x <= maxx; x+=dx) {
                int c = x + (y * m->width);
                //calculate slopes per cell
                bool visible = true;
                double startSlope = (double)processedCell*slopesPerCell;
                double centreSlope = startSlope+halfSlopes;
                double endSlope = startSlope+slopesPerCell;
                if (obstaclesInLastLine > 0 && m->cells[c].fov == 0) {
                    int idx = 0;
                    while(visible && idx < obstaclesInLastLine) {
                        if (m->cells[c].transparent == true) {
                            if (centreSlope > startAngle[idx] && centreSlope < endAngle[idx])
                                visible = false;
                            }
                        else {
                            if (startSlope >= startAngle[idx] && endSlope <= endAngle[idx])
                                visible = false;
                        }
                        if (visible && (m->cells[c-(m->width*dy)].fov == 0 || !m->cells[c-(m->width*dy)].transparent) && (x-dx >= 0 && x-dx < m->width && (m->cells[c-(m->width*dy)-dx].fov == 0 || !m->cells[c-(m->width*dy)-dx].transparent))) visible = false;
                        idx++;
                    }
                }
                if (visible) {
                    m->cells[c].fov = 1;
                    done = false;
                    //if the cell is opaque, block the adjacent slopes
                    if (!m->cells[c].transparent) {
                        if (minAngle >= startSlope) minAngle = endSlope;
                        else {
                        	startAngle[totalObstacles] = startSlope;
                        	endAngle[totalObstacles++] = endSlope;
                        }
                        if (!light_walls) m->cells[c].fov = 0;
                    }
                }
                processedCell++;
            }
            if (iteration == max_radius) done = true;
            iteration++;
            obstaclesInLastLine = totalObstacles;
            y += dy;
            if (y < 0 || y >= m->height) done = true;
			if ( minAngle == 1.0f ) done=true;
        }
    }
    //octant: horizontal edge
    {
        int iteration = 1; //iteration of the algo for this octant
        bool done = false;
        int totalObstacles = 0;
        int obstaclesInLastLine = 0;
		double minAngle = 0.0f;
		int x,y;

        //do while there are unblocked slopes left and the algo is within the map's boundaries
        //scan progressive lines/columns from the PC outwards
        x = player_x+dx; //the outer slope's coordinates (first processed line)
        if (x < 0 || x >= m->width) done = true;
		while(!done) {
            //process cells in the line
			double slopesPerCell = 1.0f/(double)(iteration+1);
			double halfSlopes = slopesPerCell*0.5f;
			int processedCell = (int)(minAngle / slopesPerCell);
            int miny = MAX(0,player_y-iteration), maxy = MIN(m->height-1,player_y+iteration);
            done = true;
            for (y = player_y + (processedCell * dy); y >= miny && y <= maxy; y+=dy) {
                int c = x + (y * m->width);
//.........这里部分代码省略.........

int computeArea(int A, int B, int C, int D, int E, int F, int G, int H) {
    int i=MAX(A,E),j=MAX(B,F),m=MIN(C,G),n=MIN(D,H),area=(D-B)*(C-A)+(H-F)*(G-E);
	if(i<m && j<n)
		area-=(n-j)*(m-i);
	return area;
}

/*
 * Set "p_ptr->food", notice observable changes
 *
 * The "p_ptr->food" variable can get as large as 20000, allowing the
 * addition of the most "filling" item, Elvish Waybread, which adds
 * 7500 food units, without overflowing the 32767 maximum limit.
 *
 * Perhaps we should disturb the player with various messages,
 * especially messages about hunger status changes.  XXX XXX XXX
 *
 * Digestion of food is handled in "dungeon.c", in which, normally,
 * the player digests about 20 food units per 100 game turns, more
 * when "fast", more when "regenerating", less with "slow digestion",
 * but when the player is "gorged", he digests 100 food units per 10
 * game turns, or a full 1000 food units per 100 game turns.
 *
 * Note that the player's speed is reduced by 10 units while gorged,
 * so if the player eats a single food ration (5000 food units) when
 * full (15000 food units), he will be gorged for (5000/100)*10 = 500
 * game turns, or 500/(100/5) = 25 player turns (if nothing else is
 * affecting the player speed).
 */
bool set_food(int v)
{
	int old_aux, new_aux;

	bool notice = FALSE;

	/* Hack -- Force good values */
	v = MIN(v, PY_FOOD_UPPER);
	v = MAX(v, 0);

	/* Fainting / Starving */
	if (p_ptr->food < PY_FOOD_FAINT)
	{
		old_aux = 0;
	}

	/* Weak */
	else if (p_ptr->food < PY_FOOD_WEAK)
	{
		old_aux = 1;
	}

	/* Hungry */
	else if (p_ptr->food < PY_FOOD_ALERT)
	{
		old_aux = 2;
	}

	/* Normal */
	else if (p_ptr->food < PY_FOOD_FULL)
	{
		old_aux = 3;
	}

	/* Full */
	else if (p_ptr->food < PY_FOOD_MAX)
	{
		old_aux = 4;
	}

	/* Gorged */
	else
	{
		old_aux = 5;
	}

	/* Fainting / Starving */
	if (v < PY_FOOD_FAINT)
	{
		new_aux = 0;
	}

	/* Weak */
	else if (v < PY_FOOD_WEAK)
	{
		new_aux = 1;
	}

	/* Hungry */
	else if (v < PY_FOOD_ALERT)
	{
		new_aux = 2;
	}

	/* Normal */
	else if (v < PY_FOOD_FULL)
	{
		new_aux = 3;
	}

	/* Full */
	else if (v < PY_FOOD_MAX)
	{
		new_aux = 4;
	}

	/* Gorged */
	else
//.........这里部分代码省略.........

/**
 * g2d_blend - blend image data in source and destination buffers.
 *
 * @ctx: a pointer to g2d_context structure.
 * @src: a pointer to g2d_image structure including image and buffer
 *	information to source.
 * @dst: a pointer to g2d_image structure including image and buffer
 *	information to destination.
 * @src_x: x start position to source buffer.
 * @src_y: y start position to source buffer.
 * @dst_x: x start position to destination buffer.
 * @dst_y: y start position to destination buffer.
 * @w: width value to source and destination buffer.
 * @h: height value to source and destination buffer.
 * @op: blend operation type.
 */
int
g2d_blend(struct g2d_context *ctx, struct g2d_image *src,
		struct g2d_image *dst, unsigned int src_x,
		unsigned int src_y, unsigned int dst_x, unsigned int dst_y,
		unsigned int w, unsigned int h, enum e_g2d_op op)
{
	union g2d_point_val pt;
	union g2d_bitblt_cmd_val bitblt;
	union g2d_blend_func_val blend;
	unsigned int src_w = 0, src_h = 0, dst_w = 0, dst_h = 0;

	bitblt.val = 0;
	blend.val = 0;

	if (op == G2D_OP_SRC || op == G2D_OP_CLEAR)
		g2d_add_cmd(ctx, DST_SELECT_REG, G2D_SELECT_MODE_BGCOLOR);
	else
		g2d_add_cmd(ctx, DST_SELECT_REG, G2D_SELECT_MODE_NORMAL);

	g2d_add_cmd(ctx, DST_COLOR_MODE_REG, dst->color_mode);
	g2d_add_base_addr(ctx, dst, g2d_dst);
	g2d_add_cmd(ctx, DST_STRIDE_REG, dst->stride);

	g2d_add_cmd(ctx, SRC_SELECT_REG, src->select_mode);
	g2d_add_cmd(ctx, SRC_COLOR_MODE_REG, src->color_mode);

	switch (src->select_mode) {
	case G2D_SELECT_MODE_NORMAL:
		g2d_add_base_addr(ctx, src, g2d_src);
		g2d_add_cmd(ctx, SRC_STRIDE_REG, src->stride);
		break;
	case G2D_SELECT_MODE_FGCOLOR:
		g2d_add_cmd(ctx, FG_COLOR_REG, src->color);
		break;
	case G2D_SELECT_MODE_BGCOLOR:
		g2d_add_cmd(ctx, BG_COLOR_REG, src->color);
		break;
	default:
		fprintf(stderr , "failed to set src.\n");
		return -EINVAL;
	}

	src_w = w;
	src_h = h;
	if (src_x + w > src->width)
		src_w = src->width - src_x;
	if (src_y + h > src->height)
		src_h = src->height - src_y;

	dst_w = w;
	dst_h = h;
	if (dst_x + w > dst->width)
		dst_w = dst->width - dst_x;
	if (dst_y + h > dst->height)
		dst_h = dst->height - dst_y;

	w = MIN(src_w, dst_w);
	h = MIN(src_h, dst_h);

	if (w <= 0 || h <= 0) {
		fprintf(stderr, "invalid width or height.\n");
		g2d_reset(ctx);
		return -EINVAL;
	}

	bitblt.data.alpha_blend_mode = G2D_ALPHA_BLEND_MODE_ENABLE;
	blend.val = g2d_get_blend_op(op);
	g2d_add_cmd(ctx, BITBLT_COMMAND_REG, bitblt.val);
	g2d_add_cmd(ctx, BLEND_FUNCTION_REG, blend.val);

	pt.val = 0;
	pt.data.x = src_x;
	pt.data.y = src_y;
	g2d_add_cmd(ctx, SRC_LEFT_TOP_REG, pt.val);
	pt.val = 0;
	pt.data.x = src_x + w;
	pt.data.y = src_y + h;
	g2d_add_cmd(ctx, SRC_RIGHT_BOTTOM_REG, pt.val);

	pt.val = 0;
	pt.data.x = dst_x;
	pt.data.y = dst_y;
	g2d_add_cmd(ctx, DST_LEFT_TOP_REG, pt.val);
	pt.val = 0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
/* ==================================================================== */
    int		iScanline;

    for( iScanline = 0; iScanline < nYSize; iScanline++ )
    {
        int	nLastRedError, nLastGreenError, nLastBlueError, i;

/* -------------------------------------------------------------------- */
/*      Report progress                                                 */
/* -------------------------------------------------------------------- */
        if( !pfnProgress( iScanline / (double) nYSize, NULL, pProgressArg ) )
        {
            CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" );
            err = CE_Failure;
            goto end_and_cleanup;
        }

/* -------------------------------------------------------------------- */
/*      Read source data.                                               */
/* -------------------------------------------------------------------- */
        GDALRasterIO( hRed, GF_Read, 0, iScanline, nXSize, 1, 
                      pabyRed, nXSize, 1, GDT_Byte, 0, 0 );
        GDALRasterIO( hGreen, GF_Read, 0, iScanline, nXSize, 1, 
                      pabyGreen, nXSize, 1, GDT_Byte, 0, 0 );
        GDALRasterIO( hBlue, GF_Read, 0, iScanline, nXSize, 1, 
                      pabyBlue, nXSize, 1, GDT_Byte, 0, 0 );

/* -------------------------------------------------------------------- */
/*	Apply the error from the previous line to this one.		*/
/* -------------------------------------------------------------------- */
        for( i = 0; i < nXSize; i++ )
        {
            pabyRed[i] = (GByte)
                MAX(0,MIN(255,(pabyRed[i]   + panError[i*3+0+3])));
            pabyGreen[i] = (GByte)
                MAX(0,MIN(255,(pabyGreen[i] + panError[i*3+1+3])));
            pabyBlue[i] =  (GByte)
                MAX(0,MIN(255,(pabyBlue[i]  + panError[i*3+2+3])));
        }

        memset( panError, 0, sizeof(int) * (nXSize+2) * 3 );

/* -------------------------------------------------------------------- */
/*	Figure out the nearest color to the RGB value.			*/
/* -------------------------------------------------------------------- */
        nLastRedError = 0;
        nLastGreenError = 0;
        nLastBlueError = 0;

        for( i = 0; i < nXSize; i++ )
        {
            int		iIndex, nError, nSixth, iRed, iGreen, iBlue;
            int		nRedValue, nGreenValue, nBlueValue;

            nRedValue =   MAX(0,MIN(255, pabyRed[i]   + nLastRedError));
            nGreenValue = MAX(0,MIN(255, pabyGreen[i] + nLastGreenError));
            nBlueValue =  MAX(0,MIN(255, pabyBlue[i]  + nLastBlueError));

            iRed   = nRedValue *   C_LEVELS   / 256;
            iGreen = nGreenValue * C_LEVELS / 256;
            iBlue  = nBlueValue *  C_LEVELS  / 256;
            
            iIndex = pabyColorMap[iRed + iGreen * C_LEVELS 
                                 + iBlue * C_LEVELS * C_LEVELS];
	
            pabyIndex[i] = (GByte) iIndex;

/// calc_desired_velocity - updates desired velocity (i.e. feed forward) with pilot requested acceleration and fake wind resistance
///		updated velocity sent directly to position controller
void AC_Loiter::calc_desired_velocity(float nav_dt)
{
    float ekfGndSpdLimit, ekfNavVelGainScaler;
    AP::ahrs_navekf().getEkfControlLimits(ekfGndSpdLimit, ekfNavVelGainScaler);

    // calculate a loiter speed limit which is the minimum of the value set by the LOITER_SPEED
    // parameter and the value set by the EKF to observe optical flow limits
    float gnd_speed_limit_cms = MIN(_speed_cms, ekfGndSpdLimit*100.0f);
    gnd_speed_limit_cms = MAX(gnd_speed_limit_cms, LOITER_SPEED_MIN);

    float pilot_acceleration_max = GRAVITY_MSS*100.0f * tanf(radians(get_angle_max_cd()*0.01f));

    // range check nav_dt
    if (nav_dt < 0) {
        return;
    }

    _pos_control.set_max_speed_xy(gnd_speed_limit_cms);
    _pos_control.set_max_accel_xy(_accel_cmss);
    _pos_control.set_leash_length_xy(LOITER_POS_CORRECTION_MAX);

    // get loiters desired velocity from the position controller where it is being stored.
    const Vector3f &desired_vel_3d = _pos_control.get_desired_velocity();
    Vector2f desired_vel(desired_vel_3d.x,desired_vel_3d.y);

    // update the desired velocity using our predicted acceleration
    desired_vel.x += _predicted_accel.x * nav_dt;
    desired_vel.y += _predicted_accel.y * nav_dt;

    Vector2f loiter_accel_brake;
    float desired_speed = desired_vel.length();
    if (!is_zero(desired_speed)) {
        Vector2f desired_vel_norm = desired_vel/desired_speed;

        // TODO: consider using a velocity squared relationship like
        // pilot_acceleration_max*(desired_speed/gnd_speed_limit_cms)^2;
        // the drag characteristic of a multirotor should be examined to generate a curve
        // we could add a expo function here to fine tune it

        // calculate a drag acceleration based on the desired speed.
        float drag_decel = pilot_acceleration_max*desired_speed/gnd_speed_limit_cms;

        // calculate a braking acceleration if sticks are at zero
        float loiter_brake_accel = 0.0f;
        if (_desired_accel.is_zero()) {
            if ((AP_HAL::millis()-_brake_timer) > _brake_delay * 1000.0f) {
                float brake_gain = _pos_control.get_vel_xy_pid().kP() * 0.5f;
                loiter_brake_accel = constrain_float(AC_AttitudeControl::sqrt_controller(desired_speed, brake_gain, _brake_jerk_max_cmsss, nav_dt), 0.0f, _brake_accel_cmss);
            }
        } else {
            loiter_brake_accel = 0.0f;
            _brake_timer = AP_HAL::millis();
        }
        _brake_accel += constrain_float(loiter_brake_accel-_brake_accel, -_brake_jerk_max_cmsss*nav_dt, _brake_jerk_max_cmsss*nav_dt);
        loiter_accel_brake = desired_vel_norm*_brake_accel;

        // update the desired velocity using the drag and braking accelerations
        desired_speed = MAX(desired_speed-(drag_decel+_brake_accel)*nav_dt,0.0f);
        desired_vel = desired_vel_norm*desired_speed;
    }

    // add braking to the desired acceleration
    _desired_accel -= loiter_accel_brake;

    // Apply EKF limit to desired velocity -  this limit is calculated by the EKF and adjusted as required to ensure certain sensor limits are respected (eg optical flow sensing)
    float horizSpdDem = desired_vel.length();
    if (horizSpdDem > gnd_speed_limit_cms) {
        desired_vel.x = desired_vel.x * gnd_speed_limit_cms / horizSpdDem;
        desired_vel.y = desired_vel.y * gnd_speed_limit_cms / horizSpdDem;
    }

    // Limit the velocity to prevent fence violations
    // TODO: We need to also limit the _desired_accel
    if (_avoid != nullptr) {
        _avoid->adjust_velocity(_pos_control.get_pos_xy_p().kP(), _accel_cmss, desired_vel, nav_dt);
    }

    // send adjusted feed forward acceleration and velocity back to the Position Controller
    _pos_control.set_desired_accel_xy(_desired_accel.x, _desired_accel.y);
    _pos_control.set_desired_velocity_xy(desired_vel.x, desired_vel.y);
}

//.........这里部分代码省略.........
	for (j=0; j<n; j++) { D[j] = z[j]/x[j]; }
	for (i=0; i<m; i++) { E[i] = w[i]/y[i]; }

	ldltfac(n, m, kAt, iAt, At, E, D, kA, iA, A, v);

	for (j=0; j<n; j++) { fx[j] = -sigma[j]; }
	for (i=0; i<m; i++) { fy[i] =  rho[i]; }

	forwardbackward(E, D, fy, fx);

	for (j=0; j<n; j++) { gx[j] = -c[j]; }
	for (i=0; i<m; i++) { gy[i] = -b[i]; }

	forwardbackward(E, D, gy, gx);

	dphi = (dotprod(c,fx,n)-dotprod(b,fy,m)+gamma)/
	       (dotprod(c,gx,n)-dotprod(b,gy,m)-psi/phi);

	for (j=0; j<n; j++) { dx[j] = fx[j] - gx[j]*dphi; }
	for (i=0; i<m; i++) { dy[i] = fy[i] - gy[i]*dphi; }

	for (j=0; j<n; j++) { dz[j] = delta*mu/x[j] - z[j] - D[j]*dx[j]; }
	for (i=0; i<m; i++) { dw[i] = delta*mu/y[i] - w[i] - E[i]*dy[i]; }
	dpsi = delta*mu/phi - psi - (psi/phi)*dphi;

        /*************************************************************
	* STEP 5: Compute step length.
        *************************************************************/

	theta = 1.0;
	for (j=0; j<n; j++) { 
	    theta 
	    = 
	    MIN(theta, linesearch(x[j],z[j],dx[j],dz[j],beta,delta,mu));
	}
	for (i=0; i<m; i++) { 
	    theta 
	    = 
	    MIN(theta,linesearch(y[i],w[i],dy[i],dw[i],beta,delta,mu));
	}
	theta = MIN(theta,linesearch(phi,psi,dphi,dpsi,beta,delta,mu));
	/*
	if (theta < 4*beta/(n+m+1)) {
		printf("ratio = %10.3e \n", theta*(n+m+1)/(4*beta));
		status = 7;  
		break;
	}
	*/
	if (theta < 1.0) theta *= 0.9999;

        /*************************************************************
	* STEP 6: Step to new point
        *************************************************************/

	for (j=0; j<n; j++) { 
	    x[j] = x[j] + theta*dx[j];
	    z[j] = z[j] + theta*dz[j];
	}
	for (i=0; i<m; i++) { 
	    y[i] = y[i] + theta*dy[i];
	    w[i] = w[i] + theta*dw[i];
	}
	phi = phi + theta*dphi;
	psi = psi + theta*dpsi;
    }  	

//.........这里部分代码省略.........
    if (thumbstr_utf8 != NULL)
	default_thumbstr = FALSE;

    error = FT_Init_FreeType (&library);
    if (error) {
	g_printerr("Could not initialise freetype: %s\n", get_ft_error (error));
	goto out;
    }

    totem_resources_monitor_start (arguments[0], 30 * G_USEC_PER_SEC);

    file = g_file_new_for_commandline_arg (arguments[0]);
    uri = g_file_get_uri (file);
    g_object_unref (file);

    face = sushi_new_ft_face_from_uri (library, uri, &contents, &gerror);
    if (gerror) {
	g_printerr ("Could not load face '%s': %s\n", uri,
		    gerror->message);
        g_free (uri);
        g_error_free (gerror);
	goto out;
    }

    g_free (uri);

    if (default_thumbstr) {
        if (check_font_contain_text (face, "Aa"))
            str = g_strdup ("Aa");
        else
            str = build_fallback_thumbstr (face);
    } else {
        str = thumbstr_utf8;
    }

    surface = cairo_image_surface_create (CAIRO_FORMAT_ARGB32,
                                          thumb_size, thumb_size);
    cr = cairo_create (surface);

    gdk_cairo_set_source_rgba (cr, &white);
    cairo_paint (cr);

    font = cairo_ft_font_face_create_for_ft_face (face, 0);
    cairo_set_font_face (cr, font);
    cairo_font_face_destroy (font);

    font_size = thumb_size - 2 * PADDING_VERTICAL;
    cairo_set_font_size (cr, font_size);
    cairo_text_extents (cr, str, &text_extents);

    if ((text_extents.width) > (thumb_size - 2 * PADDING_HORIZONTAL)) {
        scale_x = (gdouble) (thumb_size - 2 * PADDING_HORIZONTAL) / (text_extents.width);
    } else {
        scale_x = 1.0;
    }

    if ((text_extents.height) > (thumb_size - 2 * PADDING_VERTICAL)) {
        scale_y = (gdouble) (thumb_size - 2 * PADDING_VERTICAL) / (text_extents.height);
    } else {
        scale_y = 1.0;
    }

    scale = MIN (scale_x, scale_y);
    cairo_scale (cr, scale, scale);
    cairo_translate (cr, 
                     PADDING_HORIZONTAL - text_extents.x_bearing + (thumb_size - scale * text_extents.width) / 2.0,
                     PADDING_VERTICAL - text_extents.y_bearing + (thumb_size - scale * text_extents.height) / 2.0);

    gdk_cairo_set_source_rgba (cr, &black);
    cairo_show_text (cr, str);
    cairo_destroy (cr);

    cairo_surface_write_to_png (surface, arguments[1]);
    cairo_surface_destroy (surface);

    totem_resources_monitor_stop ();

    error = FT_Done_Face (face);
    if (error) {
	g_printerr("Could not unload face: %s\n", get_ft_error (error));
	goto out;
    }

    error = FT_Done_FreeType (library);
    if (error) {
	g_printerr ("Could not finalize freetype library: %s\n",
		   get_ft_error (error));
	goto out;
    }

    rv = 0; /* success */

  out:

    g_strfreev (arguments);
    g_free (str);
    g_free (contents);

    return rv;
}

/*
 * Returns true iff the specified sector is present in the disk image. Drivers
 * not implementing the functionality are assumed to not support backing files,
 * hence all their sectors are reported as allocated.
 *
 * If 'sector_num' is beyond the end of the disk image the return value is 0
 * and 'pnum' is set to 0.
 *
 * 'pnum' is set to the number of sectors (including and immediately following
 * the specified sector) that are known to be in the same
 * allocated/unallocated state.
 *
 * 'nb_sectors' is the max value 'pnum' should be set to.  If nb_sectors goes
 * beyond the end of the disk image it will be clamped.
 */
static int coroutine_fn raw_co_is_allocated(BlockDriverState *bs,
                                            int64_t sector_num,
                                            int nb_sectors, int *pnum)
{
    off_t start, data, hole;
    int ret;

    ret = fd_open(bs);
    if (ret < 0) {
        return ret;
    }

    start = sector_num * BDRV_SECTOR_SIZE;

#ifdef CONFIG_FIEMAP

    BDRVRawState *s = bs->opaque;
    struct {
        struct fiemap fm;
        struct fiemap_extent fe;
    } f;

    f.fm.fm_start = start;
    f.fm.fm_length = (int64_t)nb_sectors * BDRV_SECTOR_SIZE;
    f.fm.fm_flags = 0;
    f.fm.fm_extent_count = 1;
    f.fm.fm_reserved = 0;
    if (ioctl(s->fd, FS_IOC_FIEMAP, &f) == -1) {
        /* Assume everything is allocated.  */
        *pnum = nb_sectors;
        return 1;
    }

    if (f.fm.fm_mapped_extents == 0) {
        /* No extents found, data is beyond f.fm.fm_start + f.fm.fm_length.
         * f.fm.fm_start + f.fm.fm_length must be clamped to the file size!
         */
        off_t length = lseek(s->fd, 0, SEEK_END);
        hole = f.fm.fm_start;
        data = MIN(f.fm.fm_start + f.fm.fm_length, length);
    } else {
        data = f.fe.fe_logical;
        hole = f.fe.fe_logical + f.fe.fe_length;
    }

#elif defined SEEK_HOLE && defined SEEK_DATA

    BDRVRawState *s = bs->opaque;

    hole = lseek(s->fd, start, SEEK_HOLE);
    if (hole == -1) {
        /* -ENXIO indicates that sector_num was past the end of the file.
         * There is a virtual hole there.  */
        assert(errno != -ENXIO);

        /* Most likely EINVAL.  Assume everything is allocated.  */
        *pnum = nb_sectors;
        return 1;
    }

    if (hole > start) {
        data = start;
    } else {
        /* On a hole.  We need another syscall to find its end.  */
        data = lseek(s->fd, start, SEEK_DATA);
        if (data == -1) {
            data = lseek(s->fd, 0, SEEK_END);
        }
    }
#else
    *pnum = nb_sectors;
    return 1;
#endif

    if (data <= start) {
        /* On a data extent, compute sectors to the end of the extent.  */
        *pnum = MIN(nb_sectors, (hole - start) / BDRV_SECTOR_SIZE);
        return 1;
    } else {
        /* On a hole, compute sectors to the beginning of the next extent.  */
        *pnum = MIN(nb_sectors, (data - start) / BDRV_SECTOR_SIZE);
        return 0;
    }
}

size_t ZSTDMT_compressCCtx(ZSTDMT_CCtx* mtctx,
                           void* dst, size_t dstCapacity,
                     const void* src, size_t srcSize,
                           int compressionLevel)
{
    ZSTD_parameters params = ZSTD_getParams(compressionLevel, srcSize, 0);
    U32 const overlapLog = (compressionLevel >= ZSTD_maxCLevel()) ? 0 : 3;
    size_t const overlapSize = (size_t)1 << (params.cParams.windowLog - overlapLog);
    size_t const chunkTargetSize = (size_t)1 << (params.cParams.windowLog + 2);
    unsigned const nbChunksMax = (unsigned)(srcSize / chunkTargetSize) + 1;
    unsigned nbChunks = MIN(nbChunksMax, mtctx->nbThreads);
    size_t const proposedChunkSize = (srcSize + (nbChunks-1)) / nbChunks;
    size_t const avgChunkSize = ((proposedChunkSize & 0x1FFFF) < 0xFFFF) ? proposedChunkSize + 0xFFFF : proposedChunkSize;   /* avoid too small last block */
    size_t remainingSrcSize = srcSize;
    const char* const srcStart = (const char*)src;
    unsigned const compressWithinDst = (dstCapacity >= ZSTD_compressBound(srcSize)) ? nbChunks : (unsigned)(dstCapacity / ZSTD_compressBound(avgChunkSize));  /* presumes avgChunkSize >= 256 KB, which should be the case */
    size_t frameStartPos = 0, dstBufferPos = 0;

    DEBUGLOG(3, "windowLog : %2u => chunkTargetSize : %u bytes  ", params.cParams.windowLog, (U32)chunkTargetSize);
    DEBUGLOG(2, "nbChunks  : %2u   (chunkSize : %u bytes)   ", nbChunks, (U32)avgChunkSize);
    params.fParams.contentSizeFlag = 1;

    if (nbChunks==1) {   /* fallback to single-thread mode */
        ZSTD_CCtx* const cctx = mtctx->cctxPool->cctx[0];
        return ZSTD_compressCCtx(cctx, dst, dstCapacity, src, srcSize, compressionLevel);
    }

    {   unsigned u;
        for (u=0; u<nbChunks; u++) {
            size_t const chunkSize = MIN(remainingSrcSize, avgChunkSize);
            size_t const dstBufferCapacity = ZSTD_compressBound(chunkSize);
            buffer_t const dstAsBuffer = { (char*)dst + dstBufferPos, dstBufferCapacity };
            buffer_t const dstBuffer = u < compressWithinDst ? dstAsBuffer : ZSTDMT_getBuffer(mtctx->buffPool, dstBufferCapacity);
            ZSTD_CCtx* const cctx = ZSTDMT_getCCtx(mtctx->cctxPool);
            size_t dictSize = u ? overlapSize : 0;

            if ((cctx==NULL) || (dstBuffer.start==NULL)) {
                mtctx->jobs[u].cSize = ERROR(memory_allocation);   /* job result */
                mtctx->jobs[u].jobCompleted = 1;
                nbChunks = u+1;
                break;   /* let's wait for previous jobs to complete, but don't start new ones */
            }

            mtctx->jobs[u].srcStart = srcStart + frameStartPos - dictSize;
            mtctx->jobs[u].dictSize = dictSize;
            mtctx->jobs[u].srcSize = chunkSize;
            mtctx->jobs[u].fullFrameSize = srcSize;
            mtctx->jobs[u].params = params;
            mtctx->jobs[u].dstBuff = dstBuffer;
            mtctx->jobs[u].cctx = cctx;
            mtctx->jobs[u].firstChunk = (u==0);
            mtctx->jobs[u].lastChunk = (u==nbChunks-1);
            mtctx->jobs[u].jobCompleted = 0;
            mtctx->jobs[u].jobCompleted_mutex = &mtctx->jobCompleted_mutex;
            mtctx->jobs[u].jobCompleted_cond = &mtctx->jobCompleted_cond;

            DEBUGLOG(3, "posting job %u   (%u bytes)", u, (U32)chunkSize);
            DEBUG_PRINTHEX(3, mtctx->jobs[u].srcStart, 12);
            POOL_add(mtctx->factory, ZSTDMT_compressChunk, &mtctx->jobs[u]);

            frameStartPos += chunkSize;
            dstBufferPos += dstBufferCapacity;
            remainingSrcSize -= chunkSize;
    }   }
    /* note : since nbChunks <= nbThreads, all jobs should be running immediately in parallel */

    {   unsigned chunkID;
        size_t error = 0, dstPos = 0;
        for (chunkID=0; chunkID<nbChunks; chunkID++) {
            DEBUGLOG(3, "waiting for chunk %u ", chunkID);
            PTHREAD_MUTEX_LOCK(&mtctx->jobCompleted_mutex);
            while (mtctx->jobs[chunkID].jobCompleted==0) {
                DEBUGLOG(4, "waiting for jobCompleted signal from chunk %u", chunkID);
                pthread_cond_wait(&mtctx->jobCompleted_cond, &mtctx->jobCompleted_mutex);
            }
            pthread_mutex_unlock(&mtctx->jobCompleted_mutex);
            DEBUGLOG(3, "ready to write chunk %u ", chunkID);

            ZSTDMT_releaseCCtx(mtctx->cctxPool, mtctx->jobs[chunkID].cctx);
            mtctx->jobs[chunkID].cctx = NULL;
            mtctx->jobs[chunkID].srcStart = NULL;
            {   size_t const cSize = mtctx->jobs[chunkID].cSize;
                if (ZSTD_isError(cSize)) error = cSize;
                if ((!error) && (dstPos + cSize > dstCapacity)) error = ERROR(dstSize_tooSmall);
                if (chunkID) {   /* note : chunk 0 is already written directly into dst */
                    if (!error)
                        memmove((char*)dst + dstPos, mtctx->jobs[chunkID].dstBuff.start, cSize);  /* may overlap if chunk decompressed within dst */
                    if (chunkID >= compressWithinDst)   /* otherwise, it decompresses within dst */
                        ZSTDMT_releaseBuffer(mtctx->buffPool, mtctx->jobs[chunkID].dstBuff);
                    mtctx->jobs[chunkID].dstBuff = g_nullBuffer;
                }
                dstPos += cSize ;
            }
        }
        if (!error) DEBUGLOG(3, "compressed size : %u  ", (U32)dstPos);
        return error ? error : dstPos;
    }

}

static size_t ZSTDMT_createCompressionJob(ZSTDMT_CCtx* zcs, size_t srcSize, unsigned endFrame)
{
    size_t const dstBufferCapacity = ZSTD_compressBound(srcSize);
    buffer_t const dstBuffer = ZSTDMT_getBuffer(zcs->buffPool, dstBufferCapacity);
    ZSTD_CCtx* const cctx = ZSTDMT_getCCtx(zcs->cctxPool);
    unsigned const jobID = zcs->nextJobID & zcs->jobIDMask;

    if ((cctx==NULL) || (dstBuffer.start==NULL)) {
        zcs->jobs[jobID].jobCompleted = 1;
        zcs->nextJobID++;
        ZSTDMT_waitForAllJobsCompleted(zcs);
        ZSTDMT_releaseAllJobResources(zcs);
        return ERROR(memory_allocation);
    }

    DEBUGLOG(4, "preparing job %u to compress %u bytes with %u preload ", zcs->nextJobID, (U32)srcSize, (U32)zcs->dictSize);
    zcs->jobs[jobID].src = zcs->inBuff.buffer;
    zcs->jobs[jobID].srcStart = zcs->inBuff.buffer.start;
    zcs->jobs[jobID].srcSize = srcSize;
    zcs->jobs[jobID].dictSize = zcs->dictSize;   /* note : zcs->inBuff.filled is presumed >= srcSize + dictSize */
    zcs->jobs[jobID].params = zcs->params;
    if (zcs->nextJobID) zcs->jobs[jobID].params.fParams.checksumFlag = 0;  /* do not calculate checksum within sections, just keep it in header for first section */
    zcs->jobs[jobID].cdict = zcs->nextJobID==0 ? zcs->cdict : NULL;
    zcs->jobs[jobID].fullFrameSize = zcs->frameContentSize;
    zcs->jobs[jobID].dstBuff = dstBuffer;
    zcs->jobs[jobID].cctx = cctx;
    zcs->jobs[jobID].firstChunk = (zcs->nextJobID==0);
    zcs->jobs[jobID].lastChunk = endFrame;
    zcs->jobs[jobID].jobCompleted = 0;
    zcs->jobs[jobID].dstFlushed = 0;
    zcs->jobs[jobID].jobCompleted_mutex = &zcs->jobCompleted_mutex;
    zcs->jobs[jobID].jobCompleted_cond = &zcs->jobCompleted_cond;

    /* get a new buffer for next input */
    if (!endFrame) {
        size_t const newDictSize = MIN(srcSize + zcs->dictSize, zcs->targetDictSize);
        zcs->inBuff.buffer = ZSTDMT_getBuffer(zcs->buffPool, zcs->inBuffSize);
        if (zcs->inBuff.buffer.start == NULL) {   /* not enough memory to allocate next input buffer */
            zcs->jobs[jobID].jobCompleted = 1;
            zcs->nextJobID++;
            ZSTDMT_waitForAllJobsCompleted(zcs);
            ZSTDMT_releaseAllJobResources(zcs);
            return ERROR(memory_allocation);
        }
        DEBUGLOG(5, "inBuff filled to %u", (U32)zcs->inBuff.filled);
        zcs->inBuff.filled -= srcSize + zcs->dictSize - newDictSize;
        DEBUGLOG(5, "new job : filled to %u, with %u dict and %u src", (U32)zcs->inBuff.filled, (U32)newDictSize, (U32)(zcs->inBuff.filled - newDictSize));
        memmove(zcs->inBuff.buffer.start, (const char*)zcs->jobs[jobID].srcStart + zcs->dictSize + srcSize - newDictSize, zcs->inBuff.filled);
        DEBUGLOG(5, "new inBuff pre-filled");
        zcs->dictSize = newDictSize;
    } else {
        zcs->inBuff.buffer = g_nullBuffer;
        zcs->inBuff.filled = 0;
        zcs->dictSize = 0;
        zcs->frameEnded = 1;
        if (zcs->nextJobID == 0)
            zcs->params.fParams.checksumFlag = 0;   /* single chunk : checksum is calculated directly within worker thread */
    }

    DEBUGLOG(3, "posting job %u : %u bytes  (end:%u) (note : doneJob = %u=>%u)", zcs->nextJobID, (U32)zcs->jobs[jobID].srcSize, zcs->jobs[jobID].lastChunk, zcs->doneJobID, zcs->doneJobID & zcs->jobIDMask);
    POOL_add(zcs->factory, ZSTDMT_compressChunk, &zcs->jobs[jobID]);   /* this call is blocking when thread worker pool is exhausted */
    zcs->nextJobID++;
    return 0;
}

int cfasta_gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
    {
/*TREATMENT OF THE TERMINAL GAP PENALTIES*/
/*TG_MODE=0---> gop and gep*/
/*TG_MODE=1---> ---     gep*/
/*TG_MODE=2---> ---     ---*/


	int maximise;

/*VARIABLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
	int **tot_diag;

	int *diag;
	int ktup;
	static int n_groups;
	static char **group_list;
	int score, new_score;
        int n_chosen_diag=20;
        int step;
	int max_n_chosen_diag;
	int l1, l2;
        /********Prepare Penalties******/


	maximise=CL->maximise;
	ktup=CL->ktup;

	/********************************/




	if ( !group_list)
	   {

	       group_list=make_group_aa (&n_groups, CL->matrix_for_aa_group);
	   }

	l1=strlen (A->seq_al[l_s[0][0]]);
	l2=strlen (A->seq_al[l_s[1][0]]);

	if ( !CL->fasta_step)
	    {
	    step=MIN(l1,l2);
	    step=(int) log ((double)MAX(step, 1));
	    step=MAX(step, 20);
	    }
	else
	    {
		step=CL->fasta_step;
	    }


	tot_diag=evaluate_diagonals ( A, ns, l_s, CL, maximise,n_groups,group_list, ktup);


	max_n_chosen_diag=strlen (A->seq_al[l_s[0][0]])+strlen (A->seq_al[l_s[1][0]])-1;
	/*max_n_chosen_diag=(int)log10((double)(l1+l2))*10;*/

	n_chosen_diag+=step;
	n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);


	diag=extract_N_diag (strlen (A->seq_al[l_s[0][0]]),strlen (A->seq_al[l_s[1][0]]), tot_diag, n_chosen_diag, 0);


	score    =make_fasta_gotoh_pair_wise ( A, ns, l_s, CL, diag);

	new_score=0;
	vfree ( diag);


	while (new_score!=score && n_chosen_diag< max_n_chosen_diag    )
	  {


	    score=new_score;

	    ungap_sub_aln ( A, ns[0], l_s[0]);
	    ungap_sub_aln ( A, ns[1], l_s[1]);


	    n_chosen_diag+=step;
	    n_chosen_diag=MIN(n_chosen_diag, max_n_chosen_diag);


	    diag     =extract_N_diag (strlen (A->seq_al[l_s[0][0]]),strlen (A->seq_al[l_s[1][0]]), tot_diag, n_chosen_diag, 0);
	    new_score=make_fasta_gotoh_pair_wise (  A, ns, l_s, CL, diag);

	    vfree ( diag);

	  }

	score=new_score;
	free_int (tot_diag, -1);

	return score;
    }

/*
 * The audit system call. Trust what the user has sent down and save it
 * away in the audit file. User passes a complete audit record and its
 * length.  We will fill in the time stamp, check the header and the length
 * Put a trailer and a sequence token if policy requires.
 * In the future length might become size_t instead of an int.
 *
 * The call is valid whether or not AUDIT_PERZONE is set (think of
 * login to a zone).  When the local audit state (auk_auditstate) is
 * AUC_INIT_AUDIT, records are accepted even though auditd isn't
 * running.
 */
int
audit(caddr_t record, int length)
{
	char	c;
	int	count, l;
	token_t	*m, *n, *s, *ad;
	int	hdrlen, delta;
	adr_t	hadr;
	adr_t	sadr;
	int	size;	/* 0: 32 bit utility  1: 64 bit utility */
	int	host_len;
	size_t	zlen;
	au_kcontext_t	*kctx = GET_KCTX_PZ;

	/* if auditing not enabled, then don't generate an audit record */
	if (kctx->auk_auditstate != AUC_AUDITING &&
	    kctx->auk_auditstate != AUC_INIT_AUDIT)
		return (0);

	/* Only privileged processes can audit */
	if (secpolicy_audit_modify(CRED()) != 0)
		return (EPERM);

	/* Max user record size is 32K */
	if (length > AUDIT_REC_SIZE)
		return (E2BIG);

	/*
	 * The specified length must be at least as big as the smallest
	 * possible header token. Later after beginning to scan the
	 * header we'll determine the true minimum length according to
	 * the header type and attributes.
	 */
#define	AU_MIN_HEADER_LEN	(sizeof (char) + sizeof (int32_t) + \
	sizeof (char) + sizeof (short) + sizeof (short) + \
	(sizeof (int32_t) * 2))

	if (length < AU_MIN_HEADER_LEN)
		return (EINVAL);

	/* Read in user's audit record */
	count = length;
	m = n = s = ad = NULL;
	while (count) {
		m = au_getclr();
		if (!s)
			s = n = m;
		else {
			n->next_buf = m;
			n = m;
		}
		l = MIN(count, AU_BUFSIZE);
		if (copyin(record, memtod(m, caddr_t), (size_t)l)) {
			/* copyin failed release au_membuf */
			au_free_rec(s);
			return (EFAULT);
		}
		record += l;
		count -= l;
		m->len = (uchar_t)l;
	}

	/* Now attach the entire thing to ad */
	au_write((caddr_t *)&(ad), s);

	/* validate header token type. trust everything following it */
	adr_start(&hadr, memtod(s, char *));
	(void) adr_getchar(&hadr, &c);
	switch (c) {
	case AUT_HEADER32:
		/* size vers+event_ID+event_modifier fields */
		delta = 1 + 2 + 2;
		hdrlen = 1 + 4 + delta + (sizeof (int32_t) * 2);
		size = HEADER_SIZE32;
		break;

#ifdef _LP64
	case AUT_HEADER64:
		/* size vers+event_ID+event_modifier fields */
		delta = 1 + 2 + 2;
		hdrlen = 1 + 4 + delta + (sizeof (int64_t) * 2);
		size = HEADER_SIZE64;
		break;
#endif

	case AUT_HEADER32_EX:
		/*
		 * Skip over the length/version/type/mod fields and
//.........这里部分代码省略.........

static int
privcmd_ioctl(struct cdev *dev, unsigned long cmd, caddr_t arg,
	      int mode, struct thread *td)
{
	int error, i;

	switch (cmd) {
	case IOCTL_PRIVCMD_HYPERCALL: {
		struct ioctl_privcmd_hypercall *hcall;

		hcall = (struct ioctl_privcmd_hypercall *)arg;
#ifdef __amd64__
		/*
		 * The hypervisor page table walker will refuse to access
		 * user-space pages if SMAP is enabled, so temporary disable it
		 * while performing the hypercall.
		 */
		if (cpu_stdext_feature & CPUID_STDEXT_SMAP)
			stac();
#endif
		error = privcmd_hypercall(hcall->op, hcall->arg[0],
		    hcall->arg[1], hcall->arg[2], hcall->arg[3], hcall->arg[4]);
#ifdef __amd64__
		if (cpu_stdext_feature & CPUID_STDEXT_SMAP)
			clac();
#endif
		if (error >= 0) {
			hcall->retval = error;
			error = 0;
		} else {
			error = xen_translate_error(error);
			hcall->retval = 0;
		}
		break;
	}
	case IOCTL_PRIVCMD_MMAPBATCH: {
		struct ioctl_privcmd_mmapbatch *mmap;
		vm_map_t map;
		vm_map_entry_t entry;
		vm_object_t mem;
		vm_pindex_t pindex;
		vm_prot_t prot;
		boolean_t wired;
		struct xen_add_to_physmap_range add;
		xen_ulong_t *idxs;
		xen_pfn_t *gpfns;
		int *errs, index;
		struct privcmd_map *umap;
		uint16_t num;

		mmap = (struct ioctl_privcmd_mmapbatch *)arg;

		if ((mmap->num == 0) ||
		    ((mmap->addr & PAGE_MASK) != 0)) {
			error = EINVAL;
			break;
		}

		map = &td->td_proc->p_vmspace->vm_map;
		error = vm_map_lookup(&map, mmap->addr, VM_PROT_NONE, &entry,
		    &mem, &pindex, &prot, &wired);
		if (error != KERN_SUCCESS) {
			error = EINVAL;
			break;
		}
		if ((entry->start != mmap->addr) ||
		    (entry->end != mmap->addr + (mmap->num * PAGE_SIZE))) {
			vm_map_lookup_done(map, entry);
			error = EINVAL;
			break;
		}
		vm_map_lookup_done(map, entry);
		if ((mem->type != OBJT_MGTDEVICE) ||
		    (mem->un_pager.devp.ops != &privcmd_pg_ops)) {
			error = EINVAL;
			break;
		}
		umap = mem->handle;

		add.domid = DOMID_SELF;
		add.space = XENMAPSPACE_gmfn_foreign;
		add.foreign_domid = mmap->dom;

		/*
		 * The 'size' field in the xen_add_to_physmap_range only
		 * allows for UINT16_MAX mappings in a single hypercall.
		 */
		num = MIN(mmap->num, UINT16_MAX);

		idxs = malloc(sizeof(*idxs) * num, M_PRIVCMD, M_WAITOK);
		gpfns = malloc(sizeof(*gpfns) * num, M_PRIVCMD, M_WAITOK);
		errs = malloc(sizeof(*errs) * num, M_PRIVCMD, M_WAITOK);

		set_xen_guest_handle(add.idxs, idxs);
		set_xen_guest_handle(add.gpfns, gpfns);
		set_xen_guest_handle(add.errs, errs);

		/* Allocate a bitset to store broken page mappings. */
		umap->err = BITSET_ALLOC(mmap->num, M_PRIVCMD,
		    M_WAITOK | M_ZERO);
//.........这里部分代码省略.........

void DrawRegion( Region key, float scale ) {
	
	if ( key == NULL ) return;
	
	int stable = key->stable;
	
	char name[256];
	sprintf(name,"/tmp/T%03d.ppm",stable);
	Image out = ReadPPMFile(name);
	
	static int count = 0;
	
	if ( !ImageIsGood(out) ) {
		out = ConvertImage1(CopyImage(key->image));
		sprintf(name,"/tmp/R%05d.ppm",count++);
	} else sprintf(name,"/tmp/T%03d.ppm",stable);
	
	fprintf(stderr,".");
	
	int rv = RandomNumber(0,255);
	int gv = RandomNumber(0,rv);
	int bv = RandomNumber(0,gv);
	
	int color = PIX3(rv,gv,bv);
	
	DrawPolygon(key->border,out,color);
	
	Ellipse e1 = NewEllipse(key->row,key->col,key->maj*scale,key->min*scale,key->phi);
	DrawEllipse(e1,out,color); free(e1);
	Image patch = CreateImage(41*sqrt(2),41*sqrt(2));
	RegionToPatch(key,key->image,patch,scale);

	FVec hist = GenerateOrientationHistogram(patch);
	GaussianBlur1D(hist->values,hist->l,hist->r,2);
	DrawFVec(hist,10,10,200,400,PIX3(0,0,250),out);
	FVecFree(hist);
	
	if ( PolygonIsGood(key->sizes) ) {
		
		struct PointSt p1 = key->sizes->vertices[0];
		struct PointSt p2 = key->sizes->vertices[key->sizes->numberOfVertices-1];

		int i;
		hist = FVecNew(0,255);
		Point p;
		while ( ( p = NextPolygonVertex(key->sizes) ) != NULL ) FVecSetAt(hist,p->y,p->x);
		if ( p1.y < p2.y ) {
			for(i=p1.y;i<=p2.y;i++) if ( hist->values[i] == 0.0 ) FVecAddAt(hist,i,1);
		} else {
			for(i=p2.y;i>=p1.y;i--) if ( hist->values[i] == 0.0 ) FVecAddAt(hist,i,1);
		}
		
		hist->l = MIN(p1.y,p2.y);
		hist->r = MAX(p2.y-1,p1.y-1);
		
		DrawSizeFVec(hist,497,0,1021,1023,color,stable,out);
		DrawSizeFVec(hist,498,0,1022,1023,color,stable,out);
		DrawSizeFVec(hist,499,0,1023,1023,color,stable,out);
		
	}
	
	WritePPM(name,out);
	FreeImage(out);

}