/* SPDX-License-Identifier: GPL-2.0-or-later
 * Copyright Blender Foundation. All rights reserved. */

/** \file
 * \ingroup bke
 */

#include <cfloat>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <ctime>
#include <memory.h>

#include "DNA_gpencil_modifier_types.h"
#include "DNA_mesh_types.h"
#include "DNA_meshdata_types.h"
#include "DNA_modifier_types.h"
#include "DNA_object_types.h"

#include "BLI_math.h"
#include "BLI_math_solvers.h"
#include "BLI_task.h"
#include "BLI_utildefines.h"

#include "BKE_DerivedMesh.h"
#include "BKE_cdderivedmesh.h"
#include "BKE_context.h"
#include "BKE_lattice.h"
#include "BKE_lib_id.h"
#include "BKE_modifier.h"
#include "BKE_shrinkwrap.h"

#include "BKE_deform.h"
#include "BKE_editmesh.h"
#include "BKE_mesh.h" /* for OMP limits. */
#include "BKE_mesh_runtime.h"
#include "BKE_mesh_wrapper.h"
#include "BKE_subsurf.h"

#include "DEG_depsgraph_query.h"

#include "MEM_guardedalloc.h"

#include "BLI_strict_flags.h"

/* for timing... */
#if 0
#  include "PIL_time_utildefines.h"
#else
#  define TIMEIT_BENCH(expr, id) (expr)
#endif

/* Util macros */
#define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n"))

struct ShrinkwrapCalcData {
  ShrinkwrapModifierData *smd; /* shrinkwrap modifier data */

  Object *ob; /* object we are applying shrinkwrap to */

  MVert *vert; /* Array of verts being projected. */
  const float (*vert_normals)[3];
  /* Vertices being shrink-wrapped. */
  float (*vertexCos)[3];
  int numVerts;

  const MDeformVert *dvert; /* Pointer to mdeform array */
  int vgroup;               /* Vertex group num */
  bool invert_vgroup;       /* invert vertex group influence */

  Mesh *target;                /* mesh we are shrinking to */
  SpaceTransform local2target; /* transform to move between local and target space */
  ShrinkwrapTreeData *tree;    /* mesh BVH tree data */

  Object *aux_target;

  float keepDist; /* Distance to keep above target surface (units are in local space) */
};

struct ShrinkwrapCalcCBData {
  ShrinkwrapCalcData *calc;

  ShrinkwrapTreeData *tree;
  ShrinkwrapTreeData *aux_tree;

  float *proj_axis;
  SpaceTransform *local2aux;
};

bool BKE_shrinkwrap_needs_normals(int shrinkType, int shrinkMode)
{
  return (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) ||
         (shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX &&
          shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE);
}

bool BKE_shrinkwrap_init_tree(
    ShrinkwrapTreeData *data, Mesh *mesh, int shrinkType, int shrinkMode, bool force_normals)
{
  memset(data, 0, sizeof(*data));

  if (mesh == nullptr) {
    return false;
  }

  /* We could create a BVH tree from the edit mesh,
   * however accessing normals from the face/loop normals gets more involved.
   * Convert mesh data since this isn't typically used in edit-mode. */
  BKE_mesh_wrapper_ensure_mdata(mesh);

  if (mesh->totvert <= 0) {
    return false;
  }

  data->mesh = mesh;
  data->polys = BKE_mesh_polys(mesh);

  if (shrinkType == MOD_SHRINKWRAP_NEAREST_VERTEX) {
    data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_VERTS, 2);

    return data->bvh != nullptr;
  }

  if (mesh->totpoly <= 0) {
    return false;
  }

  data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_LOOPTRI, 4);

  if (data->bvh == nullptr) {
    return false;
  }

  if (force_normals || BKE_shrinkwrap_needs_normals(shrinkType, shrinkMode)) {
    data->pnors = BKE_mesh_poly_normals_ensure(mesh);
    if ((mesh->flag & ME_AUTOSMOOTH) != 0) {
      data->clnors = static_cast<const float(*)[3]>(CustomData_get_layer(&mesh->ldata, CD_NORMAL));
    }
  }

  if (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
    data->boundary = mesh->runtime.shrinkwrap_data;
  }

  return true;
}

void BKE_shrinkwrap_free_tree(ShrinkwrapTreeData *data)
{
  free_bvhtree_from_mesh(&data->treeData);
}

void BKE_shrinkwrap_discard_boundary_data(Mesh *mesh)
{
  ShrinkwrapBoundaryData *data = mesh->runtime.shrinkwrap_data;

  if (data != nullptr) {
    MEM_freeN((void *)data->edge_is_boundary);
    MEM_freeN((void *)data->looptri_has_boundary);
    MEM_freeN((void *)data->vert_boundary_id);
    MEM_freeN((void *)data->boundary_verts);

    MEM_freeN(data);
  }

  mesh->runtime.shrinkwrap_data = nullptr;
}

/* Accumulate edge for average boundary edge direction. */
static void merge_vert_dir(ShrinkwrapBoundaryVertData *vdata,
                           signed char *status,
                           int index,
                           const float edge_dir[3],
                           signed char side)
{
  BLI_assert(index >= 0);
  float *direction = vdata[index].direction;

  /* Invert the direction vector if either:
   * - This is the second edge and both edges have the vertex as start or end.
   * - For third and above, if it points in the wrong direction.
   */
  if (status[index] >= 0 ? status[index] == side : dot_v3v3(direction, edge_dir) < 0) {
    sub_v3_v3(direction, edge_dir);
  }
  else {
    add_v3_v3(direction, edge_dir);
  }

  status[index] = (status[index] == 0) ? side : -1;
}

static ShrinkwrapBoundaryData *shrinkwrap_build_boundary_data(Mesh *mesh)
{
  const MVert *mvert = BKE_mesh_verts(mesh);
  const MEdge *medge = BKE_mesh_edges(mesh);
  const MLoop *mloop = BKE_mesh_loops(mesh);

  /* Count faces per edge (up to 2). */
  char *edge_mode = static_cast<char *>(
      MEM_calloc_arrayN(size_t(mesh->totedge), sizeof(char), __func__));

  for (int i = 0; i < mesh->totloop; i++) {
    uint eidx = mloop[i].e;

    if (edge_mode[eidx] < 2) {
      edge_mode[eidx]++;
    }
  }

  /* Build the boundary edge bitmask. */
  BLI_bitmap *edge_is_boundary = BLI_BITMAP_NEW(mesh->totedge,
                                                "ShrinkwrapBoundaryData::edge_is_boundary");
  uint num_boundary_edges = 0;

  for (int i = 0; i < mesh->totedge; i++) {
    edge_mode[i] = (edge_mode[i] == 1);

    if (edge_mode[i]) {
      BLI_BITMAP_ENABLE(edge_is_boundary, i);
      num_boundary_edges++;
    }
  }

  /* If no boundary, return nullptr. */
  if (num_boundary_edges == 0) {
    MEM_freeN(edge_is_boundary);
    MEM_freeN(edge_mode);
    return nullptr;
  }

  /* Allocate the data object. */
  ShrinkwrapBoundaryData *data = MEM_cnew<ShrinkwrapBoundaryData>(__func__);

  data->edge_is_boundary = edge_is_boundary;

  /* Build the boundary looptri bitmask. */
  const MLoopTri *mlooptri = BKE_mesh_runtime_looptri_ensure(mesh);
  int totlooptri = BKE_mesh_runtime_looptri_len(mesh);

  BLI_bitmap *looptri_has_boundary = BLI_BITMAP_NEW(totlooptri,
                                                    "ShrinkwrapBoundaryData::looptri_is_boundary");

  for (int i = 0; i < totlooptri; i++) {
    int edges[3];
    BKE_mesh_looptri_get_real_edges(mesh, &mlooptri[i], edges);

    for (int j = 0; j < 3; j++) {
      if (edges[j] >= 0 && edge_mode[edges[j]]) {
        BLI_BITMAP_ENABLE(looptri_has_boundary, i);
        break;
      }
    }
  }

  data->looptri_has_boundary = looptri_has_boundary;

  /* Find boundary vertices and build a mapping table for compact storage of data. */
  int *vert_boundary_id = static_cast<int *>(
      MEM_calloc_arrayN(size_t(mesh->totvert), sizeof(int), __func__));

  for (int i = 0; i < mesh->totedge; i++) {
    if (edge_mode[i]) {
      const MEdge *edge = &medge[i];

      vert_boundary_id[edge->v1] = 1;
      vert_boundary_id[edge->v2] = 1;
    }
  }

  uint num_boundary_verts = 0;

  for (int i = 0; i < mesh->totvert; i++) {
    vert_boundary_id[i] = (vert_boundary_id[i] != 0) ? int(num_boundary_verts++) : -1;
  }

  data->vert_boundary_id = vert_boundary_id;
  data->num_boundary_verts = num_boundary_verts;

  /* Compute average directions. */
  ShrinkwrapBoundaryVertData *boundary_verts = static_cast<ShrinkwrapBoundaryVertData *>(
      MEM_calloc_arrayN(num_boundary_verts, sizeof(*boundary_verts), __func__));

  signed char *vert_status = static_cast<signed char *>(
      MEM_calloc_arrayN(num_boundary_verts, sizeof(char), __func__));

  for (int i = 0; i < mesh->totedge; i++) {
    if (edge_mode[i]) {
      const MEdge *edge = &medge[i];

      float dir[3];
      sub_v3_v3v3(dir, mvert[edge->v2].co, mvert[edge->v1].co);
      normalize_v3(dir);

      merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v1], dir, 1);
      merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v2], dir, 2);
    }
  }

  MEM_freeN(vert_status);

  /* Finalize average direction and compute normal. */
  const float(*vert_normals)[3] = BKE_mesh_vertex_normals_ensure(mesh);
  for (int i = 0; i < mesh->totvert; i++) {
    int bidx = vert_boundary_id[i];

    if (bidx >= 0) {
      ShrinkwrapBoundaryVertData *vdata = &boundary_verts[bidx];
      float tmp[3];

      normalize_v3(vdata->direction);

      cross_v3_v3v3(tmp, vert_normals[i], vdata->direction);
      cross_v3_v3v3(vdata->normal_plane, tmp, vert_normals[i]);
      normalize_v3(vdata->normal_plane);
    }
  }

  data->boundary_verts = boundary_verts;

  MEM_freeN(edge_mode);
  return data;
}

void BKE_shrinkwrap_compute_boundary_data(Mesh *mesh)
{
  BKE_shrinkwrap_discard_boundary_data(mesh);

  mesh->runtime.shrinkwrap_data = shrinkwrap_build_boundary_data(mesh);
}

/**
 * Shrink-wrap to the nearest vertex
 *
 * it builds a BVH-tree of vertices we can attach to and then
 * for each vertex performs a nearest vertex search on the tree.
 */
static void shrinkwrap_calc_nearest_vertex_cb_ex(void *__restrict userdata,
                                                 const int i,
                                                 const TaskParallelTLS *__restrict tls)
{
  ShrinkwrapCalcCBData *data = static_cast<ShrinkwrapCalcCBData *>(userdata);

  ShrinkwrapCalcData *calc = data->calc;
  BVHTreeFromMesh *treeData = &data->tree->treeData;
  BVHTreeNearest *nearest = static_cast<BVHTreeNearest *>(tls->userdata_chunk);

  float *co = calc->vertexCos[i];
  float tmp_co[3];
  float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);

  if (calc->invert_vgroup) {
    weight = 1.0f - weight;
  }

  if (weight == 0.0f) {
    return;
  }

  /* Convert the vertex to tree coordinates */
  if (calc->vert) {
    copy_v3_v3(tmp_co, calc->vert[i].co);
  }
  else {
    copy_v3_v3(tmp_co, co);
  }
  BLI_space_transform_apply(&calc->local2target, tmp_co);

  /* Use local proximity heuristics (to reduce the nearest search)
   *
   * If we already had an hit before.. we assume this vertex is going to have a close hit to that
   * other vertex so we can initiate the "nearest.dist" with the expected value to that last hit.
   * This will lead in pruning of the search tree. */
  if (nearest->index != -1) {
    nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
  }
  else {
    nearest->dist_sq = FLT_MAX;
  }

  BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData);

  /* Found the nearest vertex */
  if (nearest->index != -1) {
    /* Adjusting the vertex weight,
     * so that after interpolating it keeps a certain distance from the nearest position */
    if (nearest->dist_sq > FLT_EPSILON) {
      const float dist = sqrtf(nearest->dist_sq);
      weight *= (dist - calc->keepDist) / dist;
    }

    /* Convert the coordinates back to mesh coordinates */
    copy_v3_v3(tmp_co, nearest->co);
    BLI_space_transform_invert(&calc->local2target, tmp_co);

    interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */
  }
}

static void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
{
  BVHTreeNearest nearest = NULL_BVHTreeNearest;

  /* Setup nearest */
  nearest.index = -1;
  nearest.dist_sq = FLT_MAX;

  ShrinkwrapCalcCBData data{};
  data.calc = calc;
  data.tree = calc->tree;
  TaskParallelSettings settings;
  BLI_parallel_range_settings_defaults(&settings);
  settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
  settings.userdata_chunk = &nearest;
  settings.userdata_chunk_size = sizeof(nearest);
  BLI_task_parallel_range(
      0, calc->numVerts, &data, shrinkwrap_calc_nearest_vertex_cb_ex, &settings);
}

bool BKE_shrinkwrap_project_normal(char options,
                                   const float vert[3],
                                   const float dir[3],
                                   const float ray_radius,
                                   const SpaceTransform *transf,
                                   ShrinkwrapTreeData *tree,
                                   BVHTreeRayHit *hit)
{
  /* don't use this because this dist value could be incompatible
   * this value used by the callback for comparing previous/new dist values.
   * also, at the moment there is no need to have a corrected 'dist' value */
  // #define USE_DIST_CORRECT

  float tmp_co[3], tmp_no[3];
  const float *co, *no;
  BVHTreeRayHit hit_tmp;

  /* Copy from hit (we need to convert hit rays from one space coordinates to the other */
  memcpy(&hit_tmp, hit, sizeof(hit_tmp));

  /* Apply space transform (TODO readjust dist) */
  if (transf) {
    copy_v3_v3(tmp_co, vert);
    BLI_space_transform_apply(transf, tmp_co);
    co = tmp_co;

    copy_v3_v3(tmp_no, dir);
    BLI_space_transform_apply_normal(transf, tmp_no);
    no = tmp_no;

#ifdef USE_DIST_CORRECT
    hit_tmp.dist *= mat4_to_scale(((SpaceTransform *)transf)->local2target);
#endif
  }
  else {
    co = vert;
    no = dir;
  }

  hit_tmp.index = -1;

  BLI_bvhtree_ray_cast(
      tree->bvh, co, no, ray_radius, &hit_tmp, tree->treeData.raycast_callback, &tree->treeData);

  if (hit_tmp.index != -1) {
    /* invert the normal first so face culling works on rotated objects */
    if (transf) {
      BLI_space_transform_invert_normal(transf, hit_tmp.no);
    }

    if (options & MOD_SHRINKWRAP_CULL_TARGET_MASK) {
      /* Apply back-face. */
      const float dot = dot_v3v3(dir, hit_tmp.no);
      if (((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f) ||
          ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f)) {
        return false; /* Ignore hit */
      }
    }

    if (transf) {
      /* Inverting space transform (TODO: make coherent with the initial dist readjust). */
      BLI_space_transform_invert(transf, hit_tmp.co);
#ifdef USE_DIST_CORRECT
      hit_tmp.dist = len_v3v3(vert, hit_tmp.co);
#endif
    }

    BLI_assert(hit_tmp.dist <= hit->dist);

    memcpy(hit, &hit_tmp, sizeof(hit_tmp));
    return true;
  }
  return false;
}

static void shrinkwrap_calc_normal_projection_cb_ex(void *__restrict userdata,
                                                    const int i,
                                                    const TaskParallelTLS *__restrict tls)
{
  ShrinkwrapCalcCBData *data = static_cast<ShrinkwrapCalcCBData *>(userdata);

  ShrinkwrapCalcData *calc = data->calc;
  ShrinkwrapTreeData *tree = data->tree;
  ShrinkwrapTreeData *aux_tree = data->aux_tree;

  float *proj_axis = data->proj_axis;
  SpaceTransform *local2aux = data->local2aux;

  BVHTreeRayHit *hit = static_cast<BVHTreeRayHit *>(tls->userdata_chunk);

  const float proj_limit_squared = calc->smd->projLimit * calc->smd->projLimit;
  float *co = calc->vertexCos[i];
  float tmp_co[3], tmp_no[3];
  float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);

  if (calc->invert_vgroup) {
    weight = 1.0f - weight;
  }

  if (weight == 0.0f) {
    return;
  }

  if (calc->vert != nullptr && calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
    /* calc->vert contains verts from evaluated mesh. */
    /* These coordinates are deformed by vertexCos only for normal projection
     * (to get correct normals) for other cases calc->verts contains undeformed coordinates and
     * vertexCos should be used */
    copy_v3_v3(tmp_co, calc->vert[i].co);
    copy_v3_v3(tmp_no, calc->vert_normals[i]);
  }
  else {
    copy_v3_v3(tmp_co, co);
    copy_v3_v3(tmp_no, proj_axis);
  }

  hit->index = -1;

  /* TODO: we should use FLT_MAX here, but sweep-sphere code isn't prepared for that. */
  hit->dist = BVH_RAYCAST_DIST_MAX;

  bool is_aux = false;

  /* Project over positive direction of axis. */
  if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR) {
    if (aux_tree) {
      if (BKE_shrinkwrap_project_normal(0, tmp_co, tmp_no, 0.0, local2aux, aux_tree, hit)) {
        is_aux = true;
      }
    }

    if (BKE_shrinkwrap_project_normal(
            calc->smd->shrinkOpts, tmp_co, tmp_no, 0.0, &calc->local2target, tree, hit)) {
      is_aux = false;
    }
  }

  /* Project over negative direction of axis */
  if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR) {
    float inv_no[3];
    negate_v3_v3(inv_no, tmp_no);

    char options = calc->smd->shrinkOpts;

    if ((options & MOD_SHRINKWRAP_INVERT_CULL_TARGET) &&
        (options & MOD_SHRINKWRAP_CULL_TARGET_MASK)) {
      options ^= MOD_SHRINKWRAP_CULL_TARGET_MASK;
    }

    if (aux_tree) {
      if (BKE_shrinkwrap_project_normal(0, tmp_co, inv_no, 0.0, local2aux, aux_tree, hit)) {
        is_aux = true;
      }
    }

    if (BKE_shrinkwrap_project_normal(
            options, tmp_co, inv_no, 0.0, &calc->local2target, tree, hit)) {
      is_aux = false;
    }
  }

  /* don't set the initial dist (which is more efficient),
   * because its calculated in the targets space, we want the dist in our own space */
  if (proj_limit_squared != 0.0f) {
    if (hit->index != -1 && len_squared_v3v3(hit->co, co) > proj_limit_squared) {
      hit->index = -1;
    }
  }

  if (hit->index != -1) {
    if (is_aux) {
      BKE_shrinkwrap_snap_point_to_surface(aux_tree,
                                           local2aux,
                                           calc->smd->shrinkMode,
                                           hit->index,
                                           hit->co,
                                           hit->no,
                                           calc->keepDist,
                                           tmp_co,
                                           hit->co);
    }
    else {
      BKE_shrinkwrap_snap_point_to_surface(tree,
                                           &calc->local2target,
                                           calc->smd->shrinkMode,
                                           hit->index,
                                           hit->co,
                                           hit->no,
                                           calc->keepDist,
                                           tmp_co,
                                           hit->co);
    }

    interp_v3_v3v3(co, co, hit->co, weight);
  }
}

static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
{
  /* Options about projection direction */
  float proj_axis[3] = {0.0f, 0.0f, 0.0f};

  /* Ray-cast and tree stuff. */

  /** \note 'hit.dist' is kept in the targets space, this is only used
   * for finding the best hit, to get the real dist,
   * measure the len_v3v3() from the input coord to hit.co */
  BVHTreeRayHit hit;

  /* auxiliary target */
  Mesh *auxMesh = nullptr;
  ShrinkwrapTreeData *aux_tree = nullptr;
  ShrinkwrapTreeData aux_tree_stack;
  SpaceTransform local2aux;

  /* If the user doesn't allows to project in any direction of projection axis
   * then there's nothing todo. */
  if ((calc->smd->shrinkOpts &
       (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0) {
    return;
  }

  /* Prepare data to retrieve the direction in which we should project each vertex */
  if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
    if (calc->vert == nullptr) {
      return;
    }
  }
  else {
    /* The code supports any axis that is a combination of X,Y,Z
     * although currently UI only allows to set the 3 different axis */
    if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) {
      proj_axis[0] = 1.0f;
    }
    if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) {
      proj_axis[1] = 1.0f;
    }
    if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) {
      proj_axis[2] = 1.0f;
    }

    normalize_v3(proj_axis);

    /* Invalid projection direction */
    if (len_squared_v3(proj_axis) < FLT_EPSILON) {
      return;
    }
  }

  if (calc->aux_target) {
    auxMesh = BKE_modifier_get_evaluated_mesh_from_evaluated_object(calc->aux_target);
    if (!auxMesh) {
      return;
    }
    BLI_SPACE_TRANSFORM_SETUP(&local2aux, calc->ob, calc->aux_target);
  }

  if (BKE_shrinkwrap_init_tree(
          &aux_tree_stack, auxMesh, calc->smd->shrinkType, calc->smd->shrinkMode, false)) {
    aux_tree = &aux_tree_stack;
  }

  /* After successfully build the trees, start projection vertices. */
  ShrinkwrapCalcCBData data{};
  data.calc = calc;
  data.tree = calc->tree;
  data.aux_tree = aux_tree;
  data.proj_axis = proj_axis;
  data.local2aux = &local2aux;
  TaskParallelSettings settings;
  BLI_parallel_range_settings_defaults(&settings);
  settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
  settings.userdata_chunk = &hit;
  settings.userdata_chunk_size = sizeof(hit);
  BLI_task_parallel_range(
      0, calc->numVerts, &data, shrinkwrap_calc_normal_projection_cb_ex, &settings);

  /* free data structures */
  if (aux_tree) {
    BKE_shrinkwrap_free_tree(aux_tree);
  }
}

/*
 * Shrinkwrap Target Surface Project mode
 *
 * It uses Newton's method to find a surface location with its
 * smooth normal pointing at the original point.
 *
 * The equation system on barycentric weights and normal multiplier:
 *
 *   (w0*V0 + w1*V1 + w2*V2) + l * (w0*N0 + w1*N1 + w2*N2) - CO = 0
 *   w0 + w1 + w2 = 1
 *
 * The actual solution vector is [ w0, w1, l ], with w2 eliminated.
 */

//#define TRACE_TARGET_PROJECT

struct TargetProjectTriData {
  const float **vtri_co;
  const float (*vtri_no)[3];
  const float *point_co;

  float n0_minus_n2[3], n1_minus_n2[3];
  float c0_minus_c2[3], c1_minus_c2[3];

  /* Current interpolated position and normal. */
  float co_interp[3], no_interp[3];
};

/* Computes the deviation of the equation system from goal. */
static void target_project_tri_deviation(void *userdata, const float x[3], float r_delta[3])
{
  TargetProjectTriData *data = static_cast<TargetProjectTriData *>(userdata);

  const float w[3] = {x[0], x[1], 1.0f - x[0] - x[1]};
  interp_v3_v3v3v3(data->co_interp, data->vtri_co[0], data->vtri_co[1], data->vtri_co[2], w);
  interp_v3_v3v3v3(data->no_interp, data->vtri_no[0], data->vtri_no[1], data->vtri_no[2], w);

  madd_v3_v3v3fl(r_delta, data->co_interp, data->no_interp, x[2]);
  sub_v3_v3(r_delta, data->point_co);
}

/* Computes the Jacobian matrix of the equation system. */
static void target_project_tri_jacobian(void *userdata, const float x[3], float r_jacobian[3][3])
{
  TargetProjectTriData *data = static_cast<TargetProjectTriData *>(userdata);

  madd_v3_v3v3fl(r_jacobian[0], data->c0_minus_c2, data->n0_minus_n2, x[2]);
  madd_v3_v3v3fl(r_jacobian[1], data->c1_minus_c2, data->n1_minus_n2, x[2]);

  copy_v3_v3(r_jacobian[2], data->vtri_no[2]);
  madd_v3_v3fl(r_jacobian[2], data->n0_minus_n2, x[0]);
  madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]);
}

/* Clamp barycentric weights to the triangle. */
static void target_project_tri_clamp(float x[3])
{
  if (x[0] < 0.0f) {
    x[0] = 0.0f;
  }
  if (x[1] < 0.0f) {
    x[1] = 0.0f;
  }
  if (x[0] + x[1] > 1.0f) {
    x[0] = x[0] / (x[0] + x[1]);
    x[1] = 1.0f - x[0];
  }
}

/* Correct the Newton's method step to keep the coordinates within the triangle. */
static bool target_project_tri_correct(void * /*userdata*/,
                                       const float x[3],
                                       float step[3],
                                       float x_next[3])
{
  /* Insignificant correction threshold */
  const float epsilon = 1e-5f;
  /* Dot product threshold for checking if step is 'clearly' pointing outside. */
  const float dir_epsilon = 0.5f;
  bool fixed = false, locked = false;

  /* The barycentric coordinate domain is a triangle bounded by
   * the X and Y axes, plus the x+y=1 diagonal. First, clamp the
   * movement against the diagonal. Note that step is subtracted. */
  float sum = x[0] + x[1];
  float sstep = -(step[0] + step[1]);

  if (sum + sstep > 1.0f) {
    float ldist = 1.0f - sum;

    /* If already at the boundary, slide along it. */
    if (ldist < epsilon * float(M_SQRT2)) {
      float step_len = len_v2(step);

      /* Abort if the solution is clearly outside the domain. */
      if (step_len > epsilon && sstep > step_len * dir_epsilon * float(M_SQRT2)) {
        return false;
      }

      /* Project the new position onto the diagonal. */
      add_v2_fl(step, (sum + sstep - 1.0f) * 0.5f);
      fixed = locked = true;
    }
    else {
      /* Scale a significant step down to arrive at the boundary. */
      mul_v3_fl(step, ldist / sstep);
      fixed = true;
    }
  }

  /* Weight 0 and 1 boundary checks - along axis. */
  for (int i = 0; i < 2; i++) {
    if (step[i] > x[i]) {
      /* If already at the boundary, slide along it. */
      if (x[i] < epsilon) {
        float step_len = len_v2(step);

        /* Abort if the solution is clearly outside the domain. */
        if (step_len > epsilon && (locked || step[i] > step_len * dir_epsilon)) {
          return false;
        }

        /* Reset precision errors to stay at the boundary. */
        step[i] = x[i];
        fixed = true;
      }
      else {
        /* Scale a significant step down to arrive at the boundary. */
        mul_v3_fl(step, x[i] / step[i]);
        fixed = true;
      }
    }
  }

  /* Recompute and clamp the new coordinates after step correction. */
  if (fixed) {
    sub_v3_v3v3(x_next, x, step);
    target_project_tri_clamp(x_next);
  }

  return true;
}

static bool target_project_solve_point_tri(const float *vtri_co[3],
                                           const float vtri_no[3][3],
                                           const float point_co[3],
                                           const float hit_co[3],
                                           float hit_dist_sq,
                                           float r_hit_co[3],
                                           float r_hit_no[3])
{
  float x[3], tmp[3];
  float dist = sqrtf(hit_dist_sq);
  float magnitude_estimate = dist + len_manhattan_v3(vtri_co[0]) + len_manhattan_v3(vtri_co[1]) +
                             len_manhattan_v3(vtri_co[2]);
  float epsilon = magnitude_estimate * 1.0e-6f;

  /* Initial solution vector: barycentric weights plus distance along normal. */
  interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co);

  interp_v3_v3v3v3(r_hit_no, UNPACK3(vtri_no), x);
  sub_v3_v3v3(tmp, point_co, hit_co);

  x[2] = (dot_v3v3(tmp, r_hit_no) < 0) ? -dist : dist;

  /* Solve the equations iteratively. */
  TargetProjectTriData tri_data{};
  tri_data.vtri_co = vtri_co;
  tri_data.vtri_no = vtri_no;
  tri_data.point_co = point_co;

  sub_v3_v3v3(tri_data.n0_minus_n2, vtri_no[0], vtri_no[2]);
  sub_v3_v3v3(tri_data.n1_minus_n2, vtri_no[1], vtri_no[2]);
  sub_v3_v3v3(tri_data.c0_minus_c2, vtri_co[0], vtri_co[2]);
  sub_v3_v3v3(tri_data.c1_minus_c2, vtri_co[1], vtri_co[2]);

  target_project_tri_clamp(x);

#ifdef TRACE_TARGET_PROJECT
  const bool trace = true;
#else
  const bool trace = false;
#endif

  bool ok = BLI_newton3d_solve(target_project_tri_deviation,
                               target_project_tri_jacobian,
                               target_project_tri_correct,
                               &tri_data,
                               epsilon,
                               20,
                               trace,
                               x,
                               x);

  if (ok) {
    copy_v3_v3(r_hit_co, tri_data.co_interp);
    copy_v3_v3(r_hit_no, tri_data.no_interp);

    return true;
  }

  return false;
}

static bool update_hit(BVHTreeNearest *nearest,
                       int index,
                       const float co[3],
                       const float hit_co[3],
                       const float hit_no[3])
{
  float dist_sq = len_squared_v3v3(hit_co, co);

  if (dist_sq < nearest->dist_sq) {
#ifdef TRACE_TARGET_PROJECT
    printf(
        "#=#=#> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq);
#endif
    nearest->index = index;
    nearest->dist_sq = dist_sq;
    copy_v3_v3(nearest->co, hit_co);
    normalize_v3_v3(nearest->no, hit_no);
    return true;
  }

  return false;
}

/* Target projection on a non-manifold boundary edge -
 * treats it like an infinitely thin cylinder. */
static void target_project_edge(const ShrinkwrapTreeData *tree,
                                int index,
                                const float co[3],
                                BVHTreeNearest *nearest,
                                int eidx)
{
  const BVHTreeFromMesh *data = &tree->treeData;
  const MEdge *edge = &data->edge[eidx];
  const float *vedge_co[2] = {data->vert[edge->v1].co, data->vert[edge->v2].co};

#ifdef TRACE_TARGET_PROJECT
  printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n",
         eidx,
         UNPACK3(vedge_co[0]),
         UNPACK3(vedge_co[1]));
#endif

  /* Retrieve boundary vertex IDs */
  const int *vert_boundary_id = tree->boundary->vert_boundary_id;
  int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2];

  if (bid1 < 0 || bid2 < 0) {
    return;
  }

  /* Retrieve boundary vertex normals and align them to direction. */
  const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts;
  float vedge_dir[2][3], dir[3];

  copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane);
  copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane);

  sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]);

  if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) {
    negate_v3(vedge_dir[0]);
  }
  if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) {
    negate_v3(vedge_dir[1]);
  }

  /* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */
  float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]);
  float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]);
  float d0co = dot_v3v3(vedge_dir[0], co);

  float a = d0v1 - d0v0 + d1v0 - d1v1;
  float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co);
  float c = d0co - d0v0;
  float det = b * b - 4 * a * c;

  if (det >= 0) {
    const float epsilon = 1e-6f;
    float sdet = sqrtf(det);
    float hit_co[3], hit_no[3];

    for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) {
      float x = (-b + (float(i) - 1) * sdet) / (2 * a);

      if (x >= -epsilon && x <= 1.0f + epsilon) {
        CLAMP(x, 0, 1);

        float vedge_no[2][3];
        copy_v3_v3(vedge_no[0], data->vert_normals[edge->v1]);
        copy_v3_v3(vedge_no[1], data->vert_normals[edge->v2]);

        interp_v3_v3v3(hit_co, vedge_co[0], vedge_co[1], x);
        interp_v3_v3v3(hit_no, vedge_no[0], vedge_no[1], x);

        update_hit(nearest, index, co, hit_co, hit_no);
      }
    }
  }
}

/* Target normal projection BVH callback - based on mesh_looptri_nearest_point. */
static void mesh_looptri_target_project(void *userdata,
                                        int index,
                                        const float co[3],
                                        BVHTreeNearest *nearest)
{
  const ShrinkwrapTreeData *tree = (ShrinkwrapTreeData *)userdata;
  const BVHTreeFromMesh *data = &tree->treeData;
  const MLoopTri *lt = &data->looptri[index];
  const MLoop *loop[3] = {
      &data->loop[lt->tri[0]], &data->loop[lt->tri[1]], &data->loop[lt->tri[2]]};
  const MVert *vtri[3] = {
      &data->vert[loop[0]->v], &data->vert[loop[1]->v], &data->vert[loop[2]->v]};
  const float *vtri_co[3] = {vtri[0]->co, vtri[1]->co, vtri[2]->co};
  float raw_hit_co[3], hit_co[3], hit_no[3], dist_sq, vtri_no[3][3];

  /* First find the closest point and bail out if it's worse than the current solution. */
  closest_on_tri_to_point_v3(raw_hit_co, co, UNPACK3(vtri_co));
  dist_sq = len_squared_v3v3(co, raw_hit_co);

#ifdef TRACE_TARGET_PROJECT
  printf("TRIANGLE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) %g %g\n",
         index,
         UNPACK3(vtri_co[0]),
         UNPACK3(vtri_co[1]),
         UNPACK3(vtri_co[2]),
         dist_sq,
         nearest->dist_sq);
#endif

  if (dist_sq >= nearest->dist_sq) {
    return;
  }

  /* Decode normals */
  copy_v3_v3(vtri_no[0], tree->treeData.vert_normals[loop[0]->v]);
  copy_v3_v3(vtri_no[1], tree->treeData.vert_normals[loop[1]->v]);
  copy_v3_v3(vtri_no[2], tree->treeData.vert_normals[loop[2]->v]);

  /* Solve the equations for the triangle */
  if (target_project_solve_point_tri(vtri_co, vtri_no, co, raw_hit_co, dist_sq, hit_co, hit_no)) {
    update_hit(nearest, index, co, hit_co, hit_no);
  }
  /* Boundary edges */
  else if (tree->boundary && BLI_BITMAP_TEST(tree->boundary->looptri_has_boundary, index)) {
    const BLI_bitmap *is_boundary = tree->boundary->edge_is_boundary;
    int edges[3];

    BKE_mesh_looptri_get_real_edges(tree->mesh, lt, edges);

    for (int i = 0; i < 3; i++) {
      if (edges[i] >= 0 && BLI_BITMAP_TEST(is_boundary, edges[i])) {
        target_project_edge(tree, index, co, nearest, edges[i]);
      }
    }
  }
}

void BKE_shrinkwrap_find_nearest_surface(ShrinkwrapTreeData *tree,
                                         BVHTreeNearest *nearest,
                                         float co[3],
                                         int type)
{
  BVHTreeFromMesh *treeData = &tree->treeData;

  if (type == MOD_SHRINKWRAP_TARGET_PROJECT) {
#ifdef TRACE_TARGET_PROJECT
    printf("\n====== TARGET PROJECT START ======\n");
#endif

    BLI_bvhtree_find_nearest_ex(
        tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER);

#ifdef TRACE_TARGET_PROJECT
    printf("====== TARGET PROJECT END: %d %g ======\n\n", nearest->index, nearest->dist_sq);
#endif

    if (nearest->index < 0) {
      /* fallback to simple nearest */
      BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
    }
  }
  else {
    BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
  }
}

/**
 * Shrink-wrap moving vertices to the nearest surface point on the target.
 *
 * It builds a #BVHTree from the target mesh and then performs a
 * NN matches for each vertex
 */
static void shrinkwrap_calc_nearest_surface_point_cb_ex(void *__restrict userdata,
                                                        const int i,
                                                        const TaskParallelTLS *__restrict tls)
{
  ShrinkwrapCalcCBData *data = static_cast<ShrinkwrapCalcCBData *>(userdata);

  ShrinkwrapCalcData *calc = data->calc;
  BVHTreeNearest *nearest = static_cast<BVHTreeNearest *>(tls->userdata_chunk);

  float *co = calc->vertexCos[i];
  float tmp_co[3];
  float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);

  if (calc->invert_vgroup) {
    weight = 1.0f - weight;
  }

  if (weight == 0.0f) {
    return;
  }

  /* Convert the vertex to tree coordinates */
  if (calc->vert) {
    copy_v3_v3(tmp_co, calc->vert[i].co);
  }
  else {
    copy_v3_v3(tmp_co, co);
  }
  BLI_space_transform_apply(&calc->local2target, tmp_co);

  /* Use local proximity heuristics (to reduce the nearest search)
   *
   * If we already had an hit before.. we assume this vertex is going to have a close hit to that
   * other vertex so we can initiate the "nearest.dist" with the expected value to that last hit.
   * This will lead in pruning of the search tree. */
  if (nearest->index != -1) {
    if (calc->smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
      /* Heuristic doesn't work because of additional restrictions. */
      nearest->index = -1;
      nearest->dist_sq = FLT_MAX;
    }
    else {
      nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
    }
  }
  else {
    nearest->dist_sq = FLT_MAX;
  }

  BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType);

  /* Found the nearest vertex */
  if (nearest->index != -1) {
    BKE_shrinkwrap_snap_point_to_surface(data->tree,
                                         nullptr,
                                         calc->smd->shrinkMode,
                                         nearest->index,
                                         nearest->co,
                                         nearest->no,
                                         calc->keepDist,
                                         tmp_co,
                                         tmp_co);

    /* Convert the coordinates back to mesh coordinates */
    BLI_space_transform_invert(&calc->local2target, tmp_co);
    interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */
  }
}

void BKE_shrinkwrap_compute_smooth_normal(const ShrinkwrapTreeData *tree,
                                          const SpaceTransform *transform,
                                          int looptri_idx,
                                          const float hit_co[3],
                                          const float hit_no[3],
                                          float r_no[3])
{
  const BVHTreeFromMesh *treeData = &tree->treeData;
  const MLoopTri *tri = &treeData->looptri[looptri_idx];
  const float(*vert_normals)[3] = tree->treeData.vert_normals;

  /* Interpolate smooth normals if enabled. */
  if ((tree->polys[tri->poly].flag & ME_SMOOTH) != 0) {
    const uint32_t vert_indices[3] = {treeData->loop[tri->tri[0]].v,
                                      treeData->loop[tri->tri[1]].v,
                                      treeData->loop[tri->tri[2]].v};
    float w[3], no[3][3], tmp_co[3];

    /* Custom and auto smooth split normals. */
    if (tree->clnors) {
      copy_v3_v3(no[0], tree->clnors[tri->tri[0]]);
      copy_v3_v3(no[1], tree->clnors[tri->tri[1]]);
      copy_v3_v3(no[2], tree->clnors[tri->tri[2]]);
    }
    /* Ordinary vertex normals. */
    else {
      copy_v3_v3(no[0], vert_normals[vert_indices[0]]);
      copy_v3_v3(no[1], vert_normals[vert_indices[1]]);
      copy_v3_v3(no[2], vert_normals[vert_indices[2]]);
    }

    /* Barycentric weights from hit point. */
    copy_v3_v3(tmp_co, hit_co);

    if (transform) {
      BLI_space_transform_apply(transform, tmp_co);
    }

    interp_weights_tri_v3(w,
                          treeData->vert[vert_indices[0]].co,
                          treeData->vert[vert_indices[1]].co,
                          treeData->vert[vert_indices[2]].co,
                          tmp_co);

    /* Interpolate using weights. */
    interp_v3_v3v3v3(r_no, no[0], no[1], no[2], w);

    if (transform) {
      BLI_space_transform_invert_normal(transform, r_no);
    }
    else {
      normalize_v3(r_no);
    }
  }
  /* Use the polygon normal if flat. */
  else if (tree->pnors != nullptr) {
    copy_v3_v3(r_no, tree->pnors[tri->poly]);
  }
  /* Finally fallback to the looptri normal. */
  else {
    copy_v3_v3(r_no, hit_no);
  }
}

/* Helper for MOD_SHRINKWRAP_INSIDE, MOD_SHRINKWRAP_OUTSIDE and MOD_SHRINKWRAP_OUTSIDE_SURFACE. */
static void shrinkwrap_snap_with_side(float r_point_co[3],
                                      const float point_co[3],
                                      const float hit_co[3],
                                      const float hit_no[3],
                                      float goal_dist,
                                      float forcesign,
                                      bool forcesnap)
{
  float delta[3];
  sub_v3_v3v3(delta, point_co, hit_co);

  float dist = len_v3(delta);

  /* If exactly on the surface, push out along normal */
  if (dist < FLT_EPSILON) {
    if (forcesnap || goal_dist > 0) {
      madd_v3_v3v3fl(r_point_co, hit_co, hit_no, goal_dist * forcesign);
    }
    else {
      copy_v3_v3(r_point_co, hit_co);
    }
  }
  /* Move to the correct side if needed */
  else {
    float dsign = signf(dot_v3v3(delta, hit_no));

    if (forcesign == 0.0f) {
      forcesign = dsign;
    }

    /* If on the wrong side or too close, move to correct */
    if (forcesnap || dsign * dist * forcesign < goal_dist) {
      mul_v3_fl(delta, dsign / dist);

      /* At very small distance, blend in the hit normal to stabilize math. */
      float dist_epsilon = (fabsf(goal_dist) + len_manhattan_v3(hit_co)) * 1e-4f;

      if (dist < dist_epsilon) {
#ifdef TRACE_TARGET_PROJECT
        printf("zero_factor %g = %g / %g\n", dist / dist_epsilon, dist, dist_epsilon);
#endif

        interp_v3_v3v3(delta, hit_no, delta, dist / dist_epsilon);
      }

      madd_v3_v3v3fl(r_point_co, hit_co, delta, goal_dist * forcesign);
    }
    else {
      copy_v3_v3(r_point_co, point_co);
    }
  }
}

void BKE_shrinkwrap_snap_point_to_surface(const ShrinkwrapTreeData *tree,
                                          const SpaceTransform *transform,
                                          int mode,
                                          int hit_idx,
                                          const float hit_co[3],
                                          const float hit_no[3],
                                          float goal_dist,
                                          const float point_co[3],
                                          float r_point_co[3])
{
  float tmp[3];

  switch (mode) {
    /* Offsets along the line between point_co and hit_co. */
    case MOD_SHRINKWRAP_ON_SURFACE:
      if (goal_dist != 0) {
        shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, 0, true);
      }
      else {
        copy_v3_v3(r_point_co, hit_co);
      }
      break;

    case MOD_SHRINKWRAP_INSIDE:
      shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, -1, false);
      break;

    case MOD_SHRINKWRAP_OUTSIDE:
      shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, false);
      break;

    case MOD_SHRINKWRAP_OUTSIDE_SURFACE:
      if (goal_dist != 0) {
        shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, true);
      }
      else {
        copy_v3_v3(r_point_co, hit_co);
      }
      break;

    /* Offsets along the normal */
    case MOD_SHRINKWRAP_ABOVE_SURFACE:
      if (goal_dist != 0) {
        BKE_shrinkwrap_compute_smooth_normal(tree, transform, hit_idx, hit_co, hit_no, tmp);
        madd_v3_v3v3fl(r_point_co, hit_co, tmp, goal_dist);
      }
      else {
        copy_v3_v3(r_point_co, hit_co);
      }
      break;

    default:
      printf("Unknown Shrinkwrap surface snap mode: %d\n", mode);
      copy_v3_v3(r_point_co, hit_co);
  }
}

static void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
{
  BVHTreeNearest nearest = NULL_BVHTreeNearest;

  /* Setup nearest */
  nearest.index = -1;
  nearest.dist_sq = FLT_MAX;

  /* Find the nearest vertex */
  ShrinkwrapCalcCBData data{};
  data.calc = calc;
  data.tree = calc->tree;
  TaskParallelSettings settings;
  BLI_parallel_range_settings_defaults(&settings);
  settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
  settings.userdata_chunk = &nearest;
  settings.userdata_chunk_size = sizeof(nearest);
  BLI_task_parallel_range(
      0, calc->numVerts, &data, shrinkwrap_calc_nearest_surface_point_cb_ex, &settings);
}

void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd,
                               const ModifierEvalContext *ctx,
                               Scene *scene,
                               Object *ob,
                               Mesh *mesh,
                               const MDeformVert *dvert,
                               const int defgrp_index,
                               float (*vertexCos)[3],
                               int numVerts)
{

  DerivedMesh *ss_mesh = nullptr;
  ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;

  /* remove loop dependencies on derived meshes (TODO should this be done elsewhere?) */
  if (smd->target == ob) {
    smd->target = nullptr;
  }
  if (smd->auxTarget == ob) {
    smd->auxTarget = nullptr;
  }

  /* Configure Shrinkwrap calc data */
  calc.smd = smd;
  calc.ob = ob;
  calc.numVerts = numVerts;
  calc.vertexCos = vertexCos;
  calc.dvert = dvert;
  calc.vgroup = defgrp_index;
  calc.invert_vgroup = (smd->shrinkOpts & MOD_SHRINKWRAP_INVERT_VGROUP) != 0;

  if (smd->target != nullptr) {
    Object *ob_target = DEG_get_evaluated_object(ctx->depsgraph, smd->target);
    calc.target = BKE_modifier_get_evaluated_mesh_from_evaluated_object(ob_target);

    /* TODO: there might be several "bugs" with non-uniform scales matrices
     * because it will no longer be nearest surface, not sphere projection
     * because space has been deformed */
    BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, ob_target);

    /* TODO: smd->keepDist is in global units.. must change to local */
    calc.keepDist = smd->keepDist;
  }
  calc.aux_target = DEG_get_evaluated_object(ctx->depsgraph, smd->auxTarget);

  if (mesh != nullptr && smd->shrinkType == MOD_SHRINKWRAP_PROJECT) {
    /* Setup arrays to get vertices position, normals and deform weights. */
    calc.vert = BKE_mesh_verts_for_write(mesh);
    calc.vert_normals = BKE_mesh_vertex_normals_ensure(mesh);

    /* Using vertices positions/normals as if a subsurface was applied */
    if (smd->subsurfLevels) {
      SubsurfModifierData ssmd = {{nullptr}};
      ssmd.subdivType = ME_CC_SUBSURF;  /* catmull clark */
      ssmd.levels = smd->subsurfLevels; /* levels */

      /* TODO: to be moved to Mesh once we are done with changes in subsurf code. */
      DerivedMesh *dm = CDDM_from_mesh(mesh);

      ss_mesh = subsurf_make_derived_from_derived(
          dm,
          &ssmd,
          scene,
          nullptr,
          (ob->mode & OB_MODE_EDIT) ? SUBSURF_IN_EDIT_MODE : SubsurfFlags(0));

      if (ss_mesh) {
        calc.vert = static_cast<MVert *>(ss_mesh->getVertDataArray(ss_mesh, CD_MVERT));
        if (calc.vert) {
          /* TRICKY: this code assumes subsurface will have the transformed original vertices
           * in their original order at the end of the vert array. */
          calc.vert = calc.vert + ss_mesh->getNumVerts(ss_mesh) - dm->getNumVerts(dm);
        }
      }

      /* Just to make sure we are not leaving any memory behind */
      BLI_assert(ssmd.emCache == nullptr);
      BLI_assert(ssmd.mCache == nullptr);

      dm->release(dm);
    }
  }

  /* Projecting target defined - lets work! */
  ShrinkwrapTreeData tree;

  if (BKE_shrinkwrap_init_tree(&tree, calc.target, smd->shrinkType, smd->shrinkMode, false)) {
    calc.tree = &tree;

    switch (smd->shrinkType) {
      case MOD_SHRINKWRAP_NEAREST_SURFACE:
      case MOD_SHRINKWRAP_TARGET_PROJECT:
        TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface);
        break;

      case MOD_SHRINKWRAP_PROJECT:
        TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project);
        break;

      case MOD_SHRINKWRAP_NEAREST_VERTEX:
        TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), deform_vertex);
        break;
    }

    BKE_shrinkwrap_free_tree(&tree);
  }

  /* free memory */
  if (ss_mesh) {
    ss_mesh->release(ss_mesh);
  }
}

void shrinkwrapGpencilModifier_deform(ShrinkwrapGpencilModifierData *mmd,
                                      Object *ob,
                                      MDeformVert *dvert,
                                      const int defgrp_index,
                                      float (*vertexCos)[3],
                                      int numVerts)
{

  ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
  /* Convert gpencil struct to use the same struct and function used with meshes. */
  ShrinkwrapModifierData smd;
  smd.target = mmd->target;
  smd.auxTarget = mmd->aux_target;
  smd.keepDist = mmd->keep_dist;
  smd.shrinkType = mmd->shrink_type;
  smd.shrinkOpts = mmd->shrink_opts;
  smd.shrinkMode = mmd->shrink_mode;
  smd.projLimit = mmd->proj_limit;
  smd.projAxis = mmd->proj_axis;

  /* Configure Shrinkwrap calc data. */
  calc.smd = &smd;
  calc.ob = ob;
  calc.numVerts = numVerts;
  calc.vertexCos = vertexCos;
  calc.dvert = dvert;
  calc.vgroup = defgrp_index;
  calc.invert_vgroup = (mmd->flag & GP_SHRINKWRAP_INVERT_VGROUP) != 0;

  BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, mmd->target);
  calc.keepDist = mmd->keep_dist;
  calc.tree = mmd->cache_data;

  switch (mmd->shrink_type) {
    case MOD_SHRINKWRAP_NEAREST_SURFACE:
    case MOD_SHRINKWRAP_TARGET_PROJECT:
      TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), gpdeform_surface);
      break;

    case MOD_SHRINKWRAP_PROJECT:
      TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), gpdeform_project);
      break;

    case MOD_SHRINKWRAP_NEAREST_VERTEX:
      TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), gpdeform_vertex);
      break;
  }
}

void BKE_shrinkwrap_mesh_nearest_surface_deform(bContext *C, Object *ob_source, Object *ob_target)
{
  Depsgraph *depsgraph = CTX_data_depsgraph_pointer(C);
  Scene *sce = CTX_data_scene(C);
  ShrinkwrapModifierData ssmd = {{0}};
  ModifierEvalContext ctx = {depsgraph, ob_source, ModifierApplyFlag(0)};
  int totvert;

  ssmd.target = ob_target;
  ssmd.shrinkType = MOD_SHRINKWRAP_NEAREST_SURFACE;
  ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE;
  ssmd.keepDist = 0.0f;

  Mesh *src_me = static_cast<Mesh *>(ob_source->data);
  float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert);

  shrinkwrapModifier_deform(&ssmd, &ctx, sce, ob_source, src_me, nullptr, -1, vertexCos, totvert);

  BKE_mesh_vert_coords_apply(src_me, vertexCos);

  MEM_freeN(vertexCos);
}

void BKE_shrinkwrap_remesh_target_project(Mesh *src_me, Mesh *target_me, Object *ob_target)
{
  ShrinkwrapModifierData ssmd = {{0}};
  int totvert;

  ssmd.target = ob_target;
  ssmd.shrinkType = MOD_SHRINKWRAP_PROJECT;
  ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE;
  ssmd.shrinkOpts = MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR;
  ssmd.keepDist = 0.0f;

  /* Tolerance value to prevent artifacts on sharp edges of a mesh.
   * This constant and based on experimenting with different values. */
  const float projLimitTolerance = 5.0f;
  ssmd.projLimit = target_me->remesh_voxel_size * projLimitTolerance;

  float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert);

  ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;

  calc.smd = &ssmd;
  calc.numVerts = src_me->totvert;
  calc.vertexCos = vertexCos;
  calc.vert_normals = BKE_mesh_vertex_normals_ensure(src_me);
  calc.vgroup = -1;
  calc.target = target_me;
  calc.keepDist = ssmd.keepDist;
  calc.vert = BKE_mesh_verts_for_write(src_me);
  BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob_target, ob_target);

  ShrinkwrapTreeData tree;
  if (BKE_shrinkwrap_init_tree(&tree, calc.target, ssmd.shrinkType, ssmd.shrinkMode, false)) {
    calc.tree = &tree;
    TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project);
    BKE_shrinkwrap_free_tree(&tree);
  }

  BKE_mesh_vert_coords_apply(src_me, vertexCos);

  MEM_freeN(vertexCos);
}