/******************************************************************************
 ** Filename:    intfx.c
 ** Purpose:     Integer character normalization & feature extraction
 ** Author:      Robert Moss, rays@google.com (Ray Smith)
 **
 ** (c) Copyright Hewlett-Packard Company, 1988.
 ** Licensed under the Apache License, Version 2.0 (the "License");
 ** you may not use this file except in compliance with the License.
 ** You may obtain a copy of the License at
 ** http://www.apache.org/licenses/LICENSE-2.0
 ** Unless required by applicable law or agreed to in writing, software
 ** distributed under the License is distributed on an "AS IS" BASIS,
 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 ** See the License for the specific language governing permissions and
 ** limitations under the License.
 *****************************************************************************/
/**----------------------------------------------------------------------------
          Include Files and Type Defines
----------------------------------------------------------------------------**/

#define _USE_MATH_DEFINES // for M_PI

#include "intfx.h"

#include "classify.h"
#include "intmatcher.h"
#include "linlsq.h"
#include "normalis.h"
#include "statistc.h"
#include "trainingsample.h"

#include "helpers.h"

#include <allheaders.h>

#include <cmath> // for M_PI
#include <mutex> // for std::mutex

namespace tesseract {

/**----------------------------------------------------------------------------
        Global Data Definitions and Declarations
----------------------------------------------------------------------------**/
// Look up table for cos and sin to turn the intfx feature angle to a vector.
// Protected by atan_table_mutex.
// The entries are in binary degrees where a full circle is 256 binary degrees.
static float cos_table[INT_CHAR_NORM_RANGE];
static float sin_table[INT_CHAR_NORM_RANGE];

/**----------------------------------------------------------------------------
            Public Code
----------------------------------------------------------------------------**/

void InitIntegerFX() {
  // Guards write access to AtanTable so we don't create it more than once.
  static std::mutex atan_table_mutex;
  static bool atan_table_init = false;
  std::lock_guard<std::mutex> guard(atan_table_mutex);
  if (!atan_table_init) {
    for (int i = 0; i < INT_CHAR_NORM_RANGE; ++i) {
      cos_table[i] = cos(i * 2 * M_PI / INT_CHAR_NORM_RANGE + M_PI);
      sin_table[i] = sin(i * 2 * M_PI / INT_CHAR_NORM_RANGE + M_PI);
    }
    atan_table_init = true;
  }
}

// Returns a vector representing the direction of a feature with the given
// theta direction in an INT_FEATURE_STRUCT.
FCOORD FeatureDirection(uint8_t theta) {
  return FCOORD(cos_table[theta], sin_table[theta]);
}

// Generates a TrainingSample from a TBLOB. Extracts features and sets
// the bounding box, so classifiers that operate on the image can work.
// TODO(rays) Make BlobToTrainingSample a member of Classify now that
// the FlexFx and FeatureDescription code have been removed and LearnBlob
// is now a member of Classify.
TrainingSample *BlobToTrainingSample(const TBLOB &blob, bool nonlinear_norm,
                                     INT_FX_RESULT_STRUCT *fx_info,
                                     std::vector<INT_FEATURE_STRUCT> *bl_features) {
  std::vector<INT_FEATURE_STRUCT> cn_features;
  Classify::ExtractFeatures(blob, nonlinear_norm, bl_features, &cn_features, fx_info, nullptr);
  // TODO(rays) Use blob->PreciseBoundingBox() instead.
  TBOX box = blob.bounding_box();
  TrainingSample *sample = nullptr;
  int num_features = fx_info->NumCN;
  if (num_features > 0) {
    sample = TrainingSample::CopyFromFeatures(*fx_info, box, &cn_features[0], num_features);
  }
  if (sample != nullptr) {
    // Set the bounding box (in original image coordinates) in the sample.
    TPOINT topleft, botright;
    topleft.x = box.left();
    topleft.y = box.top();
    botright.x = box.right();
    botright.y = box.bottom();
    TPOINT original_topleft, original_botright;
    blob.denorm().DenormTransform(nullptr, topleft, &original_topleft);
    blob.denorm().DenormTransform(nullptr, botright, &original_botright);
    sample->set_bounding_box(
        TBOX(original_topleft.x, original_botright.y, original_botright.x, original_topleft.y));
  }
  return sample;
}

// Computes the DENORMS for bl(baseline) and cn(character) normalization
// during feature extraction. The input denorm describes the current state
// of the blob, which is usually a baseline-normalized word.
// The Transforms setup are as follows:
// Baseline Normalized (bl) Output:
//   We center the grapheme by aligning the x-coordinate of its centroid with
//   x=128 and leaving the already-baseline-normalized y as-is.
//
// Character Normalized (cn) Output:
//   We align the grapheme's centroid at the origin and scale it
//   asymmetrically in x and y so that the 2nd moments are a standard value
//   (51.2) ie the result is vaguely square.
// If classify_nonlinear_norm is true:
//   A non-linear normalization is setup that attempts to evenly distribute
//   edges across x and y.
//
// Some of the fields of fx_info are also setup:
// Length: Total length of outline.
// Rx:     Rounded y second moment. (Reversed by convention.)
// Ry:     rounded x second moment.
// Xmean:  Rounded x center of mass of the blob.
// Ymean:  Rounded y center of mass of the blob.
void Classify::SetupBLCNDenorms(const TBLOB &blob, bool nonlinear_norm, DENORM *bl_denorm,
                                DENORM *cn_denorm, INT_FX_RESULT_STRUCT *fx_info) {
  // Compute 1st and 2nd moments of the original outline.
  FCOORD center, second_moments;
  int length = blob.ComputeMoments(&center, &second_moments);
  if (fx_info != nullptr) {
    fx_info->Length = length;
    fx_info->Rx = IntCastRounded(second_moments.y());
    fx_info->Ry = IntCastRounded(second_moments.x());

    fx_info->Xmean = IntCastRounded(center.x());
    fx_info->Ymean = IntCastRounded(center.y());
  }
  // Setup the denorm for Baseline normalization.
  bl_denorm->SetupNormalization(nullptr, nullptr, &blob.denorm(), center.x(), 128.0f, 1.0f, 1.0f,
                                128.0f, 128.0f);
  // Setup the denorm for character normalization.
  if (nonlinear_norm) {
    std::vector<std::vector<int>> x_coords;
    std::vector<std::vector<int>> y_coords;
    TBOX box;
    blob.GetPreciseBoundingBox(&box);
    box.pad(1, 1);
    blob.GetEdgeCoords(box, x_coords, y_coords);
    cn_denorm->SetupNonLinear(&blob.denorm(), box, UINT8_MAX, UINT8_MAX, 0.0f, 0.0f, x_coords,
                              y_coords);
  } else {
    cn_denorm->SetupNormalization(nullptr, nullptr, &blob.denorm(), center.x(), center.y(),
                                  51.2f / second_moments.x(), 51.2f / second_moments.y(), 128.0f,
                                  128.0f);
  }
}

// Helper normalizes the direction, assuming that it is at the given
// unnormed_pos, using the given denorm, starting at the root_denorm.
static uint8_t NormalizeDirection(uint8_t dir, const FCOORD &unnormed_pos, const DENORM &denorm,
                                  const DENORM *root_denorm) {
  // Convert direction to a vector.
  FCOORD unnormed_end;
  unnormed_end.from_direction(dir);
  unnormed_end += unnormed_pos;
  FCOORD normed_pos, normed_end;
  denorm.NormTransform(root_denorm, unnormed_pos, &normed_pos);
  denorm.NormTransform(root_denorm, unnormed_end, &normed_end);
  normed_end -= normed_pos;
  return normed_end.to_direction();
}

// Helper returns the mean direction vector from the given stats. Use the
// mean direction from dirs if there is information available, otherwise, use
// the fit_vector from point_diffs.
static FCOORD MeanDirectionVector(const LLSQ &point_diffs, const LLSQ &dirs, const FCOORD &start_pt,
                                  const FCOORD &end_pt) {
  FCOORD fit_vector;
  if (dirs.count() > 0) {
    // There were directions, so use them. To avoid wrap-around problems, we
    // have 2 accumulators in dirs: x for normal directions and y for
    // directions offset by 128. We will use the one with the least variance.
    FCOORD mean_pt = dirs.mean_point();
    double mean_dir = 0.0;
    if (dirs.x_variance() <= dirs.y_variance()) {
      mean_dir = mean_pt.x();
    } else {
      mean_dir = mean_pt.y() + 128;
    }
    fit_vector.from_direction(Modulo(IntCastRounded(mean_dir), 256));
  } else {
    // There were no directions, so we rely on the vector_fit to the points.
    // Since the vector_fit is 180 degrees ambiguous, we align with the
    // supplied feature_dir by making the scalar product non-negative.
    FCOORD feature_dir(end_pt - start_pt);
    fit_vector = point_diffs.vector_fit();
    if (fit_vector.x() == 0.0f && fit_vector.y() == 0.0f) {
      // There was only a single point. Use feature_dir directly.
      fit_vector = feature_dir;
    } else {
      // Sometimes the least mean squares fit is wrong, due to the small sample
      // of points and scaling. Use a 90 degree rotated vector if that matches
      // feature_dir better.
      FCOORD fit_vector2 = !fit_vector;
      // The fit_vector is 180 degrees ambiguous, so resolve the ambiguity by
      // insisting that the scalar product with the feature_dir should be +ve.
      if (fit_vector % feature_dir < 0.0) {
        fit_vector = -fit_vector;
      }
      if (fit_vector2 % feature_dir < 0.0) {
        fit_vector2 = -fit_vector2;
      }
      // Even though fit_vector2 has a higher mean squared error, it might be
      // a better fit, so use it if the dot product with feature_dir is bigger.
      if (fit_vector2 % feature_dir > fit_vector % feature_dir) {
        fit_vector = fit_vector2;
      }
    }
  }
  return fit_vector;
}

// Helper computes one or more features corresponding to the given points.
// Emitted features are on the line defined by:
// start_pt + lambda * (end_pt - start_pt) for scalar lambda.
// Features are spaced at feature_length intervals.
static int ComputeFeatures(const FCOORD &start_pt, const FCOORD &end_pt, double feature_length,
                           std::vector<INT_FEATURE_STRUCT> *features) {
  FCOORD feature_vector(end_pt - start_pt);
  if (feature_vector.x() == 0.0f && feature_vector.y() == 0.0f) {
    return 0;
  }
  // Compute theta for the feature based on its direction.
  uint8_t theta = feature_vector.to_direction();
  // Compute the number of features and lambda_step.
  double target_length = feature_vector.length();
  int num_features = IntCastRounded(target_length / feature_length);
  if (num_features == 0) {
    return 0;
  }
  // Divide the length evenly into num_features pieces.
  double lambda_step = 1.0 / num_features;
  double lambda = lambda_step / 2.0;
  for (int f = 0; f < num_features; ++f, lambda += lambda_step) {
    FCOORD feature_pt(start_pt);
    feature_pt += feature_vector * lambda;
    INT_FEATURE_STRUCT feature(feature_pt, theta);
    features->push_back(feature);
  }
  return num_features;
}

// Gathers outline points and their directions from start_index into dirs by
// stepping along the outline and normalizing the coordinates until the
// required feature_length has been collected or end_index is reached.
// On input pos must point to the position corresponding to start_index and on
// return pos is updated to the current raw position, and pos_normed is set to
// the normed version of pos.
// Since directions wrap-around, they need special treatment to get the mean.
// Provided the cluster of directions doesn't straddle the wrap-around point,
// the simple mean works. If they do, then, unless the directions are wildly
// varying, the cluster rotated by 180 degrees will not straddle the wrap-
// around point, so mean(dir + 180 degrees) - 180 degrees will work. Since
// LLSQ conveniently stores the mean of 2 variables, we use it to store
// dir and dir+128 (128 is 180 degrees) and then use the resulting mean
// with the least variance.
static int GatherPoints(const C_OUTLINE *outline, double feature_length, const DENORM &denorm,
                        const DENORM *root_denorm, int start_index, int end_index, ICOORD *pos,
                        FCOORD *pos_normed, LLSQ *points, LLSQ *dirs) {
  int step_length = outline->pathlength();
  ICOORD step = outline->step(start_index % step_length);
  // Prev_normed is the start point of this collection and will be set on the
  // first iteration, and on later iterations used to determine the length
  // that has been collected.
  FCOORD prev_normed;
  points->clear();
  dirs->clear();
  int num_points = 0;
  int index;
  for (index = start_index; index <= end_index; ++index, *pos += step) {
    step = outline->step(index % step_length);
    int edge_weight = outline->edge_strength_at_index(index % step_length);
    if (edge_weight == 0) {
      // This point has conflicting gradient and step direction, so ignore it.
      continue;
    }
    // Get the sub-pixel precise location and normalize.
    FCOORD f_pos = outline->sub_pixel_pos_at_index(*pos, index % step_length);
    denorm.NormTransform(root_denorm, f_pos, pos_normed);
    if (num_points == 0) {
      // The start of this segment.
      prev_normed = *pos_normed;
    } else {
      FCOORD offset = *pos_normed - prev_normed;
      float length = offset.length();
      if (length > feature_length) {
        // We have gone far enough from the start. We will use this point in
        // the next set so return what we have so far.
        return index;
      }
    }
    points->add(pos_normed->x(), pos_normed->y(), edge_weight);
    int direction = outline->direction_at_index(index % step_length);
    if (direction >= 0) {
      direction = NormalizeDirection(direction, f_pos, denorm, root_denorm);
      // Use both the direction and direction +128 so we are not trying to
      // take the mean of something straddling the wrap-around point.
      dirs->add(direction, Modulo(direction + 128, 256));
    }
    ++num_points;
  }
  return index;
}

// Extracts Tesseract features and appends them to the features vector.
// Startpt to lastpt, inclusive, MUST have the same src_outline member,
// which may be nullptr. The vector from lastpt to its next is included in
// the feature extraction. Hidden edges should be excluded by the caller.
// If force_poly is true, the features will be extracted from the polygonal
// approximation even if more accurate data is available.
static void ExtractFeaturesFromRun(const EDGEPT *startpt, const EDGEPT *lastpt,
                                   const DENORM &denorm, double feature_length, bool force_poly,
                                   std::vector<INT_FEATURE_STRUCT> *features) {
  const EDGEPT *endpt = lastpt->next;
  const C_OUTLINE *outline = startpt->src_outline;
  if (outline != nullptr && !force_poly) {
    // Detailed information is available. We have to normalize only from
    // the root_denorm to denorm.
    const DENORM *root_denorm = denorm.RootDenorm();
    int total_features = 0;
    // Get the features from the outline.
    int step_length = outline->pathlength();
    int start_index = startpt->start_step;
    // pos is the integer coordinates of the binary image steps.
    ICOORD pos = outline->position_at_index(start_index);
    // We use an end_index that allows us to use a positive increment, but that
    // may be beyond the bounds of the outline steps/ due to wrap-around, to
    // so we use % step_length everywhere, except for start_index.
    int end_index = lastpt->start_step + lastpt->step_count;
    if (end_index <= start_index) {
      end_index += step_length;
    }
    LLSQ prev_points;
    LLSQ prev_dirs;
    FCOORD prev_normed_pos = outline->sub_pixel_pos_at_index(pos, start_index);
    denorm.NormTransform(root_denorm, prev_normed_pos, &prev_normed_pos);
    LLSQ points;
    LLSQ dirs;
    FCOORD normed_pos(0.0f, 0.0f);
    int index = GatherPoints(outline, feature_length, denorm, root_denorm, start_index, end_index,
                             &pos, &normed_pos, &points, &dirs);
    while (index <= end_index) {
      // At each iteration we nominally have 3 accumulated sets of points and
      // dirs: prev_points/dirs, points/dirs, next_points/dirs and sum them
      // into sum_points/dirs, but we don't necessarily get any features out,
      // so if that is the case, we keep accumulating instead of rotating the
      // accumulators.
      LLSQ next_points;
      LLSQ next_dirs;
      FCOORD next_normed_pos(0.0f, 0.0f);
      index = GatherPoints(outline, feature_length, denorm, root_denorm, index, end_index, &pos,
                           &next_normed_pos, &next_points, &next_dirs);
      LLSQ sum_points(prev_points);
      // TODO(rays) find out why it is better to use just dirs and next_dirs
      // in sum_dirs, instead of using prev_dirs as well.
      LLSQ sum_dirs(dirs);
      sum_points.add(points);
      sum_points.add(next_points);
      sum_dirs.add(next_dirs);
      bool made_features = false;
      // If we have some points, we can try making some features.
      if (sum_points.count() > 0) {
        // We have gone far enough from the start. Make a feature and restart.
        FCOORD fit_pt = sum_points.mean_point();
        FCOORD fit_vector = MeanDirectionVector(sum_points, sum_dirs, prev_normed_pos, normed_pos);
        // The segment to which we fit features is the line passing through
        // fit_pt in direction of fit_vector that starts nearest to
        // prev_normed_pos and ends nearest to normed_pos.
        FCOORD start_pos = prev_normed_pos.nearest_pt_on_line(fit_pt, fit_vector);
        FCOORD end_pos = normed_pos.nearest_pt_on_line(fit_pt, fit_vector);
        // Possible correction to match the adjacent polygon segment.
        if (total_features == 0 && startpt != endpt) {
          FCOORD poly_pos(startpt->pos.x, startpt->pos.y);
          denorm.LocalNormTransform(poly_pos, &start_pos);
        }
        if (index > end_index && startpt != endpt) {
          FCOORD poly_pos(endpt->pos.x, endpt->pos.y);
          denorm.LocalNormTransform(poly_pos, &end_pos);
        }
        int num_features = ComputeFeatures(start_pos, end_pos, feature_length, features);
        if (num_features > 0) {
          // We made some features so shuffle the accumulators.
          prev_points = points;
          prev_dirs = dirs;
          prev_normed_pos = normed_pos;
          points = next_points;
          dirs = next_dirs;
          made_features = true;
          total_features += num_features;
        }
        // The end of the next set becomes the end next time around.
        normed_pos = next_normed_pos;
      }
      if (!made_features) {
        // We didn't make any features, so keep the prev accumulators and
        // add the next ones into the current.
        points.add(next_points);
        dirs.add(next_dirs);
      }
    }
  } else {
    // There is no outline, so we are forced to use the polygonal approximation.
    const EDGEPT *pt = startpt;
    do {
      FCOORD start_pos(pt->pos.x, pt->pos.y);
      FCOORD end_pos(pt->next->pos.x, pt->next->pos.y);
      denorm.LocalNormTransform(start_pos, &start_pos);
      denorm.LocalNormTransform(end_pos, &end_pos);
      ComputeFeatures(start_pos, end_pos, feature_length, features);
    } while ((pt = pt->next) != endpt);
  }
}

// Extracts sets of 3-D features of length kStandardFeatureLength (=12.8), as
// (x,y) position and angle as measured counterclockwise from the vector
// <-1, 0>, from blob using two normalizations defined by bl_denorm and
// cn_denorm. See SetpuBLCNDenorms for definitions.
// If outline_cn_counts is not nullptr, on return it contains the cumulative
// number of cn features generated for each outline in the blob (in order).
// Thus after the first outline, there were (*outline_cn_counts)[0] features,
// after the second outline, there were (*outline_cn_counts)[1] features etc.
void Classify::ExtractFeatures(const TBLOB &blob, bool nonlinear_norm,
                               std::vector<INT_FEATURE_STRUCT> *bl_features,
                               std::vector<INT_FEATURE_STRUCT> *cn_features,
                               INT_FX_RESULT_STRUCT *results,
                               std::vector<int> *outline_cn_counts) {
  DENORM bl_denorm, cn_denorm;
  tesseract::Classify::SetupBLCNDenorms(blob, nonlinear_norm, &bl_denorm, &cn_denorm, results);
  if (outline_cn_counts != nullptr) {
    outline_cn_counts->clear();
  }
  // Iterate the outlines.
  for (TESSLINE *ol = blob.outlines; ol != nullptr; ol = ol->next) {
    // Iterate the polygon.
    EDGEPT *loop_pt = ol->FindBestStartPt();
    EDGEPT *pt = loop_pt;
    if (pt == nullptr) {
      continue;
    }
    do {
      if (pt->IsHidden()) {
        continue;
      }
      // Find a run of equal src_outline.
      EDGEPT *last_pt = pt;
      do {
        last_pt = last_pt->next;
      } while (last_pt != loop_pt && !last_pt->IsHidden() &&
               last_pt->src_outline == pt->src_outline);
      last_pt = last_pt->prev;
      // Until the adaptive classifier can be weaned off polygon segments,
      // we have to force extraction from the polygon for the bl_features.
      ExtractFeaturesFromRun(pt, last_pt, bl_denorm, kStandardFeatureLength, true, bl_features);
      ExtractFeaturesFromRun(pt, last_pt, cn_denorm, kStandardFeatureLength, false, cn_features);
      pt = last_pt;
    } while ((pt = pt->next) != loop_pt);
    if (outline_cn_counts != nullptr) {
      outline_cn_counts->push_back(cn_features->size());
    }
  }
  results->NumBL = bl_features->size();
  results->NumCN = cn_features->size();
  results->YBottom = blob.bounding_box().bottom();
  results->YTop = blob.bounding_box().top();
  results->Width = blob.bounding_box().width();
}

} // namespace tesseract