/* Last edited on 2007-05-09 08:14:43 by stolfi */
typedef struct irange_t { int end[2]; } irange_t;
  /* An interval with integer endpoints.
    When used to describe image domains, it represents all pixel
    indices in {end[0]..end[1]-1}. */

typedef struct plan_t 
  { i2_t SSize;     /* Dimensions of scaled-down image {S}. */
    i2_t ESize;     /* Dimensions of extracted images {E[r,s]}. */
    i2_t EPos;      /* Shift of {E[0,0]} relative to {S}. */
    irange_t tr[2]; /* Tile index along {ax} ranges in {tr[ax].end[0]..tr[ax].end[1]}. */
    double score;   /* Score of plan (negative if unacceptable). */
  } plan_t;
  /* An image-extraction plan. */
  
typedef enum
  { edge_L = 0, /* Left (min X) edge. */
    edge_B = 1, /* Bottom (min Y) edge. */
    edge_R = 2, /* Right (max X) edge. */
    edge_T = 3  /* Top (max Y) edge. */
  } edge_t;
  /* Index of an edge of a rectangle. Note that bit 0 is the 
    perpendicular axis, and bit 1 specifes min or max. */


plan_t fse_find_best_plan(i2_t *OSize, bool_t paddable[], options_t *o);
  /* Finds the best image extraction plan for an image {O} whose size
    (after marging stripping) is {OSize}. If {paddable[e]} is TRUE,
    for any {e} in {0..3}, assumes that the image {O} can be extended
    in the direction of edge {e} by padding with some appropriate
    color. Plans are enumerated and compared according to the
    parameters in {o}. In particular, considers for the extracted
    images {E[r,s]} all sizes in the series that begins with
    {o.minSize} or {o.midESize} and proceeds by doubling both
    dimensions, without exceeding {o.maxSize}. */

void fse_eval_plans_for_esize(i2_t *OSize, i2_t *ESize, double maxScale, double maxStretch, double maxPad, double maxCut, bool_t paddable[], plan_t *pl_best);
  /* Considers the plans {pl} that entail scaling down an image {O}
    whose size is {OSize} to produce some image {S}, and then
    extracting from {S} a set of images {E[r,s]} whose size is
    {ESize}. Whenever it finds a plan {pl} with better score 
    than {pl_best}, updates {pl_best} with it. The client must
    initialize {pl_best} suitable before the call.
    
    For the image {S}, considers all interesting dimensions {SX × SY}
    such that the scale factor from {O} to {S} is at most {maxScale},
    and the longest side of {S} is stretched by at a factor not exceeding
    {maxStretch} relative to the size that would preserve the aspect of {O}. 

    Assumes that edge {e} in {0..3]} of {O} (hence {S}) can be padded
    iff {paddable[e]} is TRUE. The amont by which an image {E[r,s]}
    extends outside {S} along either axis should not exceed {maxPad}
    times the size of {E} along that axis. The amount by which the
    union of all {E[r,s]} falls short of {S} along either axis must
    not exceed {maxCut} times the extent of {S} along that axis. */

void fse_eval_plans_for_ssize_esize(i2_t *OSize, i2_t *SSize, i2_t *ESize, double maxScale, double maxStretch, double maxPad, double maxCut, bool_t paddable[], plan_t *pl_best);
  /* Considers the plans {pl} that entail scaling down an image {O}
    whose size is {OSize} to size {SSize}, and then extracting from
    {S} a set of images {E[r,s]} whose size is {ESize}. Whenever it
    finds a plan {pl} with better score than than {pl_best}, updates
    {pl_best} with {pl}. The client must initialize {pl_best} suitable
    before the call.
    
    The plans considered differ only on the number of tiles,
    namely the ranges for the indices {r} and {s}.
    
    The parameters {maxScale}, {maxStretch}, {maxPad}, {maxCut}, and
    {paddable} are interpreted as in {fse_eval_plans_for_esize}. */
    
void fse_write_description(i2_t *ISize, irange_t OBox[], i2_t *SSize, int r, int s, i2_t *ESize, i2_t *EPos, char *prefix);
  /* Writes the description of extracted image {E[r,s]}
    to a text file called "{prefix}-{R}{S}.txt" where {R} is the digit
    {r+5} and {S} is the digit {s+5}.
    
    Assumes that the original image had dimensions {ISize}, was
    margin-stripped down to the rectangle {OBox}, and then scaled down
    to dimensions {SSize}. Also assumes that the sub-image {E[r,s]}
    starts at position {EPos + (r*DX,s*DY)} where
    {DX=floor(ESize.c[0]/2)} and {DY=floor(ESiz.c[1]/2)}. 
    
    When calling this procedure, the displacements refer to the LOWER
    LEFT corners with Y axis pointing UP. In the output file, the
    displacements refer to the UPPER LEFT corners, with the Y axis
    pointing DOWN.
    
    Fails if {|r|} or {|s|} exceed 4. */

plan_t fse_find_best_plan(i2_t *OSize, bool_t paddable[], options_t *o)
  {
    plan_t pl_best;
    
    pl_best.score = -INF;
    bool_t mid;  /* FALSE for main size series, TRUE for intermediate series. */
    bool_t turn; /* FALSE for unturned sizes, TRUE for turned sizes. */
    for (mid = FALSE; mid <= TRUE; mid++)
      for (turn = FALSE; turn <= o->rotate; turn++)
        { i2_t minESize = (mid ? o->midSize : o->minSize);
          i2_t maxESize = o->maxSize;
          if (turn) 
            { /* Swap axes in {minESize}, {fin{Size}: */
              int t;
              t = minESize.c[0]; minESize.c[0] = minESize.c[1]; minESize.c[1] = t;
              t = maxESize.c[0]; maxESize.c[0] = maxESize.c[1]; maxESize.c[1] = t;
            }
          /* Enumerate valid image sizes {EX,EY} and consider each in turn: */
          i2_t ESize = minESize;
          while ((ESize.c[0] <= maxESize.c[0]) && (ESize.c[1] <= maxESize.c[1]))
            { /* Consider plans that extract images of size {ESize}: */
              fse_eval_plans_for_esize(OSize, &ESize, o->maxScale, o->maxPad, o->maxCut, paddable, &pl_best);
              /* Get next size pair in same series: */
              ESize.c[0] *= 2; ESize.c[1] *= 2;
            }
        }
    return pl_best;
  }

void fse_eval_plans_for_esize(i2_t *OSize, i2_t *ESize, double maxScale, double maxStretch, double maxPad, double maxCut, bool_t paddable[], plan_t *pl_best)
  {
    double supScale = maxScale; /* Maximum linear factor for {O->S} scaling. */
    double infScale = 0.0;      /* Minimum linear factor for {O->S} scaling. */
    bool_t axis_pad[2];         /* Tells whether {O} can be padded along each axis. */ 
    /* Tighten {infScale} by considering {maxPad}: */
    int ax;
    for (ax = 0; ax < 2; ax++)
      { /* Get the two edges perpendicular to {ax}: */
        edge_t eLo = (ax == 0 ? edge_L : edge_B);
        edge_t eHi = (ax == 0 ? edge_R : edge_T);
        /* Compute the scale factor {fitScale} for exact fit along this axis: */
        double fitScale = ((double)ESize.c[ax])/((double)OSize.c[ax]);
        if ((! paddable[eLo]) && (! paddable[eHi]))
          { /* Image {O} is not paddable in the direction {ax}: */ 
            axis_pad[ax] = TRUE;
            /* Then {SSize} must be at least {ESize} along axis {ax}: */
            double r = 0.999999*fitScale;
            if (r > infScale) { infScale = r; }
          }
        else
          { /* Image {O} is paddable in the direction {ax}: */
            axis_pad[ax] = TRUE;
            /* Then the padded fraction of {E[0,0]} along {ax} must not exceed {maxPad}: */
            double r = 0.999999*(1 - maxPax)*fitScale;
            if (r > infScale) { infScale = r; }
          }
      }
    /* Find the longest and shortest axes of {OSize}: */
    int axMax = (OSize.c[0] >= OSize.c[1] ? 0 : 1);
    int axMin = 1 - axMax;
    /* Range of variation for the longest dimension of {SSize}: */
    int tad = 1; /* Vary the longest side by this much below and above the ideal value. */
    /* Compute the maximum and minimum sizes for {S}: */
    i2_t minSSize, maxSSize;
    for (ax = 0; ax < 2; ax++)
      { /* The minimum size is defined by {OSize} and {infScale}: */
        minSSize.c[ax] = (int)floor(infScale*OSize.c[ax]);
        if (minSSize.c[ax] < 1) { minSSize.c[ax] = 1; }
        /* The maximum size is defined by {OSize} and {supScale}: */
        maxSSize.c[ax] = (int)ceil(supScale*OSize.c[ax]);
        if (maxSSize.c[ax] > OSize.c[ax]) { maxSSize.c[ax] = OSize.c[ax]; }
        /* For the longest axis, reduce {minSize} by a tad: */
        if (ax == axMax) { maxSSize.c[ax] -= tad; }
      }

    /* Consider all alternatives for the scaled image size {SSize}: */
    i2_t SSize = minSSize;
    while ((SSize.c[0] <= maxSSize.c[0]) && (SSize.c[1] <= maxSSize.c[1]))
      { /* Evaluate the plan: */
        fse_eval_plans_for_ssize_esize(OSize, &SSize, ESize, maxScale, maxPad, maxCut, paddable, pl_best);
        /* Compute upper limits for {SSize.c[axMax]} given {SSize.c[axMin]}: */
        int max1 = maxSSize.c[axMax];
        int max2 = (int)ceil(((double)OSize.c[axMax]*SSize.c[axMin])/((double)OSize.c[axMin])) + tad;
        /* Increment {SSize} along longest axis, if possible: */
        SSize.c[axMax]++;
        if ((SSize.c[axMax] > max1) || (SSize.c[axMax] > max2))
          { /* Cannot increment {SSize.c[axMax]} any more. */
            /* Increment {SSize.c[axMin]} instead: */
            SSize.c[axMin]++;
            /* Compute lower limits for {SSize.c[axMax]} given {SSize.c[axMin]}: */
            int min1 = minSSize.c[axMax]
            int min2 = (int)floor(((double)OSize.c[axMax]*SSize.c[axMin])/((double)OSize.c[axMin])) - tad;
            SSize.c[axMax] = (min1 > min2 ? min1 : min2);
          }
      }
  }

void fse_eval_plans_for_ssize_esize(i2_t *OSize, i2_t *SSize, i2_t *ESize, double maxScale, double maxStretch, double maxPad, double maxCut, bool_t paddable[], plan_t *pl_best)
  {

    plan_t pl;
    pl.SSize = *SSize;
    pl.ESize = *ESize;

    int ax;
    /* Position of {E[0,0]} is concentric with {S}. */
    /* Round down if {E} is smaller than {S}, up if larger. */
    for (ax = 0; ax < 2; ax++) { pl.EPos.c[ax] = (SSize.c[ax] - ESize.c[ax])/2; }

    /* We consider only one-dimensional tilings. */
    /* Consider each axis in turn: */
    for (ax = 0; ax < 2; ax++) 
      { /* Compute the displacement between tiles along axis {ax}: */
        int D = ESize.c[ax]/2;
        /* Compute a range {v} that is enough to cover {SSize} along {ax}: */
        irange_t v; /* Range that covers the full extent: */
        v.end[0] = - (EPos.c[ax] + D - 1)/D;
        if (v.end[0] > 0) { v.end[0] = 0; }
        v.end[1] = (SSize.c[ax] - EPos.c[ax] - ESize.c[ax] + D - 1)/D;
        if (v.end[1] < 0) { v.end[1] = 0; }
        /* Try that range, and also the trimmed ranges, if possible: */
        int dlo, dhi;
        for (dlo = 0; dlo <= 1; dlo++)
          for (dhi = 0; dhi <= 1; dhi++)
            { /* Compute the adjusted range {w}: */
              irange_t w = (irange_t){{ v.end[0] + dlo, v.end[1]-dhi }};
              if ((w.end[0] <= 0) && (w.end[1] >= 0))
                { /* Avoid considering twice the single pic: */
                  if ((ax == 0) || (w.end[0] < 0) || (w.end[1] > 0))
                    { /* Try the range {w} along {ax}, single along other axis: */
                      pl.tile[ax] = w;
                      pl.tile[1-ax] = (irange_t){{ 0, 0 }};
                      /* Evaluate and update: */
                      pl.score = fse_eval_plan(OSize, &pl, maxScale, maxStretch, maxPad, maxCut, paddable);
                      if (pl.score > pl_best->score) { (*pl_best) = pl; }
                    }
                }
            }
      }
  }

void fse_eval_plan(i2_t *OSize, plan_t *pl, double maxScale, double maxStretch, double maxPad, double maxCut, bool_t paddable[], plan_t *pl_best)
{
  fprintf(stderr, "plan"); 
  fprintf(stderr, "  OSize = %5d %5d", OSize.c[0], OSize.c[1]);
  fprintf(stderr, "  SSize = %5d %5d", pl.SSize.c[0], pl.SSize.c[1]);
  fprintf(stderr, "  ESize = %5d %5d", pl.ESize.c[0], pl.ESize.c[1]);
  fprintf(stderr, "  EPos = %5d %5d", pl.EPos.c[0], pl.EPos.c[1]);
  fprintf(stderr, "  rrange = %5d %5d", pl.tile[0].end[0], pl.tile[0].end[1]);
  fprintf(stderr, "  srange = %5d %5d", pl.tile[1].end[0], pl.tile[1].end[1]);
  
  /* Compute number of tiles in each diretion: */
  int nx = (pl->tile[0].end[1] - pl->tile[0].end[0] + 1);
  int ny = (pl->tile[1].end[1] - pl->tile[1].end[0] + 1);
  int nt = nx*ny; /* Number of tiles in plan. */
  if ((nx > 1) && (ny > 1))
    { /* Tiling is bidimensional, not allowed: */
      pl->score = -1.0;
    }
  else
    { /* Tiling is singly-dimensional: */
      
      /* Compute displacement {D} between tiles: */
      i2_t D = (i2_t){{ pl->ESize.c[0]/2,  pl->ESize.c[1]/2 }};

      /* Compute partial scores: 1 for best possible, negative if contraint is violated. */

      /* Compute {ScaleScore}: Does {S/O} exceed {maxScale}?  */
      double scaleScore = 1.0; /* Score for {maxScale} compliance. */
      for (ax = 0; ax < 2; ax++) 
        { int MDim = (int) ceil(maxScale*OSize.c[ax]); /* Max size for {S} along {ax}. */
          int SDim = SSize.c[ax]; /* Actual size of {S} along {ax}. */
          double s = (SDim > MDim ? -1.0 : +1.0);
          if (s < scaleScore) { scaleScore = s; }
        }

      /* Compute {shapeScore}: How close is the aspect of {S} to that of {O}? */
      double r = log(((double)OSize.c[0]*SSize.c[1])/((double)OSize.c[1]*SSize.c[0]));
      double shapeScore = 1 - (r*r)/(2*log(maxStretch));

      /* Compute {padScore}: Is {maxPad} honored? How much padding is needed? */
      double padScore = 1.0;
      for (ax = 0; ax < 2; ax++) 
        { /* Compute the real interval {v} covered by all tiles along axis {ax}: */
          irange_t v;
          v.end[0] = pl->EPos.c[ax] + pl->tile[ax].end[0]*D.c[ax]; 
          v.end[1] = pl->EPos.c[ax] + pl->ESize.c[ax] + pl->tile[ax].end[1]*D.c[ax];
          /* Compute the amount of padding required at each end: */
          padLo = (v.end[0] < 0 ? 0 - v.end[0] : 0);
          padHi = (v.end[1] > SSize.c[ax] ? v.end[1] - SSize.c[ax] : 0);
          /* Get the edges {eLo,eHi} perpendicular to {ax}: */
          edge_t eLo = (ax == 0 ? edge_L : edge_B);
          edge_t eHi = (ax == 0 ? edge_R : edge_T);
          /* Compute max padding allwed (in pixels) along this axis: */
          int numPad = (int)ceil(maxPad * pl->ESize.c[ax]); 
          /* Compute pad score {s} for this axis: */
          double s;
          if (((! paddable[eLo]) && (padLo > 0)) || ((! paddable[eHi]) && (padHi > 0)))
            { /* Image does not allow padding as planned: */
              s = -1.0;
            }
          else if (padLo + padHi > numPad)
            { /* Plan requires too much padding: */
              s = -1.0;
            }
          else
            { /* Compute the ratio {r} of planned padding to max padding: */
              double r = ((double)padLo + padHi)/((double)numPad);
              s = 1.0 - r*r;
              if (s < 0) { s = 0.0; }
            }
          /* Update {padScore}: */
          if (s < padScore) { padScore = s; }
        }

      /* Compute {cutScore}: How much of {S} is left uncovered? */
      double cutScore = 1.0;
      for (ax = 0; ax < 2; ax++) 
        { /* Compute the real interval {v} covered by all tiles along axis {ax}: */
          irange_t v;
          v.end[0] = pl->EPos.c[ax] + pl->tile[ax].end[0]*D.c[ax]; 
          v.end[1] = pl->EPos.c[ax] + pl->ESize.c[ax] + pl->tile[ax].end[1]*D.c[ax];
          /* Compute the amount of cropping required at each end: */
          cutLo = (v.end[0] < 0 ? : v.end[0] );
          cutHi = (v.end[1] > SSize.c[ax] ? 0 : SSize.c[ax] - v.end[1] );
          /* Compute the maximum amount of cropping allowed along {ax}: */
          int numCut = (int)ceil(maxCut * pl->SSize.c[ax]); 
          /* Compute the cut score {s} for this axis: */
          double s;
          if (cutLo + cutHi > numCut)
            { /* Plan requires too much cropping: */
              s = -1.0;
            }
          else
            { /* Compute the ratio {r} of planned cropping to max cropping: */
              double r = ((double)cutLo + cutHi)/((double)numCut);
              s = 1.0 - r*r;
              if (s < 0) { s = 0.0; }
            }
          /* Update {cutScore}: */
          if (s < cutScore) { cutScore = s; }
        }

      /* Compute combined score = minimum of partial scores. */
      double score = 1.0;
      if (scaleScore < score) { score = scaleScore; }
      if (shapeScore < score) { score = shapeScore; }
      if (padScore < score) { score = padScore; }
      if (cutScore < score) { score = cutScore; }
      if (tileScore < score) { score = tileScore; }
      tileFactor = 1.0 - 0.5*(nt - 1); 
    }
    
  
  /* If positive, multiply by tileScore: */
  if (score > 0)
    { score * tileScore; }
    
    
      else
        { double r = ((double)SDim)/((double)MDim);
          double s = 1 - r*r;
}