Finding paths between triangles efficiently in 3D geometry #2

Question

This post is an update of the one from here. I've updated the code and a couple pieces of the post itself.

I've been writing some functions used to find paths between two types of triangles - alphas and betas. Alphas are triangles that have been in a zone we consider important, have an "interesting" value above a given threshold, and are "active". Betas are essentially anything that isn't an Alpha.

The position of the zone and the geometry of the model can change between invocations. Thus, both alphas and betas change almost every invocation to some extent. This requires a complete re-computation of the paths between them.

This is written in C++03, compiled into a MEX file (.mexa64) to be executed by MATLAB R2016B on a Linux machine. Those are all hard limits I do not have control over. I cannot make it C++11 or later.

This code uses a good deal of functions and data from an external libraries and objects. However, most of the methods used are very simple array lookups, nothing performance-hindering.

Everything works correctly so far in testing, but performance has become a significant problem.

Note on globals: They're needed because some information must persist between calls to the MEX file. This is the only way to do it, other than perhaps by writing them out to a file, which would be much slower. I know globals aren't ideal - it's just what I have to work with.

What the MATLAB Script is providing:

I should be clear on this as well. The MATLAB script or the Raytracer (both of which I cannot modify) provide nodeIValues, elemFace, and anything from m_nrt or the CRTWrapper that I use. I can't touch those.

The code:

// Doxygen block goes here

// Various includes

// Only needed because ultimately the MATLAB script needs an error code, not a
// C++ exception
#define SUCCESS 0
#define PTHREAD_ERR 1

typedef std::pair<unsigned int, unsigned int> ABPair;

// Useful for multithreading
struct ThreadData
{
  CRayTracer* rt;
  pthread_t threadId;                          // ID returned by pthread_create
  unsigned uThreadID;                          // Index
  std::vector<ABPair> validPathsThread;        // valid pairs that this thread
                                               // found
  unsigned int numTris;                        // Input argument, the number of
                                               // triangles in the mesh
  double distThreshold;                        // Input argument, the maximum
                                               // distance between triangles
};

// Exception for experimentation
class PThreadException: public std::exception
{
  virtual const char* what() const throw()
  {
    return "Exception occured in a pthread_attr_init or pthread_create\n";
  }
};

// Data about each individual tri, could be brought intro a vector of structs
// Needed to check if geometry has changed since last invokation
std::vector<bool> triActive;
// Needed to check if alphas have changed since last invokation
std::vector<bool> validAlphaIndex;
// Needed to keep history of what tris have ever been in the beam, for alphas
std::vector<bool> hasBeenInBeam;

// A "map" from a given face to the element it resides in. Many faces may share
// the same element.
std::vector<unsigned int> faceToElementMap;

// Not neccesary as a global if it's just getting re-generated each time.
// However, if we do decide to append and cull instead of regenerating, this
// needs to stay.
std::vector<unsigned int> validAlphas;

// All valid paths. Must be maintained, because we don't know if
// findPaths() will be called. It may not be if geometry hasnt changed.
std::vector<ThermalPair> validPaths;
unsigned int prevPathNum = 0;

// Useful everywhere
CRTWrapper* rayTracer = NULL;
NanoRTWrapper* m_nrt = NULL;

// Function declarations
// Not required, but prevents warnings depending on how functions are ordered
// and call each other
// (Including the mexFunction here would be redundant, as it exists in mex.h)
void exitFcn();
bool isTriInZoneRadius(const unsigned int iTri);
bool checkForModelChanges(const unsigned int numTris,
                          const float* nodeIValues,
                          const double iValueThreshold
                          );
void initialize(const float* elemFace,
                const unsigned int numElems,
                const unsigned int facePerElMax,
                unsigned int* numTri,
                unsigned int* numFace
                );
void* findPathsThread(void *data);
void findPathsThreadSpooler(const unsigned int numTris,
                            const double distThreshold
                            );
void mapFacesToElements(const float* elemFace,
                        const unsigned int numElems,
                        const unsigned int facePerElMax
                        );
bool checkPairValid(const unsigned int i,
                    const unsigned int j,
                    const double distThreshold
                    );
bool isTriAlpha(const unsigned int iTri,
                const float* nodeIValues,
                const double iValueThreshold
                );
int mainFunc(some args gohere);

/**
 * @brief exitFcn - Cleans up malloc'd or calloc'd memory if someone in the
 * MATLAB script calls "clear mexFileName" or "clear all". Does nothing ATM.
 */
void exitFcn()
{
  // mexPrintf("exitFcn() called\n");
  // Does nothing right now, since I don't malloc/calloc anything
}

/**
 * @brief Checks if a given tri is currently in the zone's external radius.
 * @param iTri - The index of the triangle to check
 * @return True if in the radius, false if not
 */
bool isTriInZoneRadius(const unsigned int iTri)
{
  // Omitted, relies on some funky external stuff that'd be hard to explain
  // hasBeenInZone[] gets updated here
}

/**
 * @brief Checks if the model has changed (either in terms of alphas or
 * geometry) and re-generates the vector of alphas
 * @param numTris -       The number of triangles in the mesh
 * @param nodeIValues -   The iValue at each node
 * @param iValueThreshold - The iValue threshold beyond which an alpha
 * is interesting enough to be valid
 * @return True if the list of alphas or the geometry has changed, false if
 * neither have
 */
bool checkForModelChanges(const unsigned int numTris,
                          const float* nodeIValues,
                          const double iValueThreshold
                          )
{
  bool modelChanged = false;
  bool isAlpha;
  bool currentlyActive;

  // Two checks need to happen - geometry changes and for the list of valid
  // alphas to change
  // Also regenerates the vector of valid alphas from scratch as it goes

  for(unsigned int i = 0; i < numTris; i++)
  {
    // Active means it has 1 exposed face, not 2 (internal) or 0 (eroded)
    currentlyActive = m_nrt->getTriActive(i);

    // Has the geometry changed?
    if(currentlyActive != triActive[i])
    {
      modelChanged = true;
      triActive[i] = currentlyActive;
    }

    // Get whether this triangle is an alpha:
    isAlpha = isTriAlpha(i, nodeIValues, iValueThreshold);

    // Triangle is a valid alpha now, but wasn't before
    if((isAlpha == true) && (validAlphaIndex[i] == false))
    {
      validAlphaIndex[i] = true;
      modelChanged = true;
    }
    // Was valid before, is no longer valid now
    else if((isAlpha == false) && (validAlphaIndex[i] == true))
    {
      validAlphaIndex[i] = false;
      modelChanged = true;
      //cullalphasFlag = true;
    }

    // Generating the set of all valid alphas
    if(isAlpha)
    {
      validAlphas.push_back(i);
    }
  }

  return modelChanged;
}

/**
 * @brief Initializes this MEX file for its first run
 * @param rt -            A pointer to the raytracer object
 * @param numTris -       The total number of triangles in the mesh
 * @param numFaces -      The total number of faces in the mesh
 * @param elemFace -      The map of elements to the faces that they have
 * @param numElems -      The number of elements in the mesh
 * @param facePerElMax -  The maximum number of faces per element
 */
void initialize(const float* elemFace,
                const unsigned int numElems,
                const unsigned int facePerElMax,
                unsigned int* numTri,
                unsigned int* numFace
                )
{
  // Fetch number of tris and faces
  // Must be done every time, since we're storing locally and not globally
  // However:
  // They're never modified
  // They never change between calls to rtThermalCalculate()
  // They're used frequently in many functions
  // I think that's a strong candidate for being a global

  unsigned int numTris = m_nrt->getTriCount();
  *numTri = numTris;

  unsigned int numFaces = m_nrt->getFaceCount();
  *numFace = numFaces;

  /*
   * Allocate some space for things we need to be persistent between runs of
   * this MEX file.
   */
  if(triActive.empty())
  {
    triActive.resize(numTris, false);
  }
  if(hasBeenInZone.empty())
  {
    hasBeenInZone.resize(numTris, false);
  }
  if(validAlphaIndex.empty())
  {
    validAlphaIndex.resize(numTris, false);
  }
  if(faceToElementMap.empty())
  {
    faceToElementMap.resize(numFaces);
    mapFacesToElements(elemFace, numElems, facePerElMax);
  }

  return;
}

/**
 * @brief Is something that can be used by pthread_create(). Threads will skip
 * over some of the work, and do isValidPair on others. Thus...multithreading.
 * @param data - The data structure that will hold the results and arguments
 */
void* findPathsThread(void *data)
{
  struct ThreadData* thisThreadsData = static_cast<struct ThreadData*>(data);
  const unsigned uThreadID = thisThreadsData->uThreadID;
  const unsigned uNumThreads = rayTracer->m_uNumThreads;
  const double distThreshold = thisThreadsData->distThreshold;
  const unsigned int numTris = thisThreadsData->numTris;
  unsigned int validI;

  std::vector<ABPair>& validPathsThread = thisThreadsData->validPathsThread;

  // Loop over all valid alphas
  for(unsigned int i = uThreadID; i < validAlphas.size(); i += uNumThreads)
  {
    // Get this to avoid needing to index into the array 4 times total
    // Each time
    validI = validAlphas[i];

    // Loop over all triangles (potential betas)
    for(unsigned int j = 0; j < numTris; j++)
    {
      // Do the easy checks first to avoid function call overhead
      if(!validAlphaIndex[j] && triActive[j])
      {
        if(checkPairValid(validI, j, distThreshold))
        {
          validPathsThread.push_back(std::make_pair(validI, j));
        }
      }
    }
  }
  return NULL;
}

/**
 * @brief Uses the raytracer object's current state as well as arguments to
 * generate pairs of unobstructed paths between alphas and betas. Creates
 * as many threads as the system has available, and then uses pthread_create()
 * to dish out the work of findPaths()
 * @param numTris - The number of triangles in the mesh
 * @param distThreshold - The maximum distance an alpha and beta can be
 * apart
 */
void findPathsThreadSpooler(const unsigned int numTris,
                            const double distThreshold
                            )
{
  std::vector<ThreadData> threadData(rayTracer->m_nProc);
  pthread_attr_t attr;
  int rc;

  // I think this is checking to make sure something doesn't already exist,
  // not sure what though
  if((rc = pthread_attr_init(&attr)))
  {
    throw PThreadException();
  }

  // We know how many threads the system supports
  // So all this does is walk through an array of them and start them up
  for(unsigned uThread = 0; uThread < rayTracer->m_uNumThreads; uThread++)
  {
    ThreadData& data = threadData[uThread];

    data.rt = rayTracer;
    data.uThreadID = uThread;
    data.numTris = numTris;
    data.distThreshold = distThreshold;

    if(rayTracer->m_uNumThreads > 1)
    {
      if((rc = pthread_create(&data.threadId, &attr, &findPathsThread, &data)))
      {
        throw PThreadException();
      }
    }
    else
    {
      findPathsThread(&data);
    }
  }

  // Join all threads
  for(unsigned uThread = 0; uThread < rayTracer->m_uNumThreads; uThread++)
  {
    std::vector<ABPair>& validPathsThread =
        threadData[uThread].validPathsThread;

    if(rayTracer->m_uNumThreads > 1)
    {
      void* res;

      if((rc = pthread_join(threadData[uThread].threadId, &res)))
      {
        throw PThreadException();
      }
    }

    // validPathsThread is the set of ABPairs that this thread found
    // while validPaths is the globally maintained set of valid paths
    // Take each thread's results and merge it into the overall results
    validPaths.insert(validPaths.end(),
                      validPathsThread.begin(),
                      validPathsThread.end());
  }

  // Useful for preallocation next time
  prevPathNum = validPaths.size();

  return;
}


/*
void cullalphas()
{
  for(unsigned int i = 0; i < validAlphas.size(); i++)
  {
    if(!isValidalpha(validAlphas[i]))
    {
      validAlphas.erase(i);
    }
  }
}
*/

/**
 * @brief Determines the elements that each face belongs to
 * @details the MATLAB script maintains a map of all faces per element.
 * This is the opposite of what we want. Accessing it linearly
 * walks by column, not by row. Also, MATLAB stores everything 1-indexed.
 * Finally, the MATLAB script left them stored as the default, which are
 * singles.
 * @param elemFace - A MATLAB facePerElMax by numElems array, storing which
 * faces belong to each element (elements being the row number)
 * @param numElems - The total number of elements (rows) in the array
 * @param facePerElMax - The max number of faces per element (the number of
 * columns)
 */
void mapFacesToElements(const float* elemFace,
                        const unsigned int numElems,
                        const unsigned int facePerElMax
                        )
{
  unsigned int i;

  // elemFace[0] = 1. We don't know how elemFace will be structured precisely,
  // so we need to keep going until we find a face in it that equals our number
  // of faces, since it's 1-indexed.
  for(i = 0; i < (numElems * facePerElMax); i++)
  {
    faceToElementMap[static_cast<unsigned int>(elemFace[i]) - 1] =
        (i % numElems);

    // Is the next face for that element a NaN? If so, we can skip it. Keep
    // skipping until the next element WON'T be NaN.
    // Don't cast here, as NaN only exists for floating point numbers,
    // not integers.
    while(((i + 1) < (numElems * facePerElMax)) && isnan(elemFace[i + 1]))
    {
      i++;
    }
  }
}

/**
 * @brief checkPairValid - Checks if a pair of an alpha index
 * (of validAlphas), beta index form a valid path
 * @param i -             Index into validAlphas
 * @param j -             Index into all tris (potential beta)
 * @param distThreshold - The max distance the tri's centers can be apart
 * @return Whether the pair forms a valid path
 */
bool checkPairValid(const unsigned int i,
                    const unsigned int j,
                    const double distThreshold
                    )
{
  double pathDist;
  double alphaCoords[3];
  double betaCoords[3];
  nanort::Ray<double> ray;

  alphaCoords[0] = rayTracer->m_vecTriFixedInfo[i].center.x();
  alphaCoords[1] = rayTracer->m_vecTriFixedInfo[i].center.y();
  alphaCoords[2] = rayTracer->m_vecTriFixedInfo[i].center.z();

  betaCoords[0] = rayTracer->m_vecTriFixedInfo[j].center.x();
  betaCoords[1] = rayTracer->m_vecTriFixedInfo[j].center.y();
  betaCoords[2] = rayTracer->m_vecTriFixedInfo[j].center.z();

  // Determine distance squared between alpha and beta
  // sqrt((x2-x1)^2 + (y2-y1)^2 +(z2-z1)^2)
  pathDist = sqrt(pow((betaCoords[0] - alphaCoords[0]), 2)
                + pow((betaCoords[1] - alphaCoords[1]), 2)
                + pow((betaCoords[2] - alphaCoords[2]), 2));

  // Doing this instead of doing the sqrt to save doing the sqrt when not
  // needed for performance
  if(pathDist < distThreshold)
  {
    // Set up a nanort::Ray's origin, direction, and max distance
    ray.org[0] = alphaCoords[0]; // x
    ray.org[1] = alphaCoords[1]; // y
    ray.org[2] = alphaCoords[2]; // z

    ray.dir[0] = (betaCoords[0] - alphaCoords[0]) / pathDist;
    ray.dir[1] = (betaCoords[1] - alphaCoords[1]) / pathDist;
    ray.dir[2] = (betaCoords[2] - alphaCoords[2]) / pathDist;

    // TODO: Subtract some EPSILON here so it doesn't report a hit because it
    // hit the beta itself (assuming that's how it works)
    ray.max_t = pathDist;

    // Call CNmg::ShootRay()'s third form to check if there is a path
    if(!(m_nrt->shootRay(ray)))
    {
      return true;
    }
    else
    {
      // There's no path
      return false;
    }
  }
  else
  {
    // The distance is too far between alpha and beta
    return false;
  }
}

/**
 * @brief Determines if a given triangle is a valid alpha.
 * @param iTri - The triangle index to check
 * @return True if it is an alpha, false if it is not
 */
bool isTriAlpha(const unsigned int iTri,
                const float* nodeIValues,
                const double iValueThreshold
                )
{
  double triAvgIValue;
  const unsigned int* triNodes;

  // Do the simple checks first, as it's more performant to do so
  // alternate consideration for accuracy
  //if(triActive[iTri] && (hasBeenAlpha[iTri] || isTriInZoneRadius(iTri)))
  if(triActive[iTri] && (hasBeenInZone[iTri] || isTriInZoneRadius(iTri)))
  {
    // Retrieve the average iValue of this triangle
    triNodes = m_nrt->getTriNodes(iTri);

    triAvgIValue = (nodeIValues[triNodes[0]]
                  + nodeIValues[triNodes[1]]
                  + nodeIValues[triNodes[2]]) / 3;

    if(triAvgIValue > iValueThreshold)
    {
      return true;
    }
  }

  return false;
}

// Doxygen block, omitted
int mainFunc(args)
{
  // Some local vars, omitted

  // Initialize the program if we're on a first run
  initialize(elemFace, numElems, facePerElMax, &numTris, &numFaces);

  // Need to see if we need to call findPaths
  if(checkForModelChanges(numTris, nodeIValues, iValueThreshold))
  {
    validPaths.clear();
    validPaths.reserve(prevPathNum);

    try
    {
      findPathsThreadSpooler(numTris, distThreshold);
    }
    catch(PThreadException& e)
    {
      return PTHREAD_ERR;
    }
  }

  // Loop over all valid paths, use them to do some more calculations..(omitted)
  // This takes up hundreds of the time findPaths() takes

  // Clear vector of valid alphas, it'll be re-generated from scratch each time
  validAlphas.clear()
}

// Doxygen block goes here, omitted, specific code also omitted as it's
// irrelevant
void mexFunction(int nlhs,
                 mxArray *plhs[],
                 int nrhs,
                 const mxArray *prhs[]
                 )
{
  // Register exit function

  // Prep for writing out results

  // Checking to make sure # of arguments was right from MATLAB

  // Input argument handling to convert from mxArrays to double*, float*, etc

  // *errcode = mainFunc(some args)

  // retrieve execution time in clock cycles, convert to seconds, print

  // Put the outputs in plhs
}

Callgraph(?):

This isn't exactly a callgraph, but it might be useful to get an idea of the flow of the program.

The Problem: Performance

For medium-size models (104k tris, 204k faces, 51k elems) it can take up to a couple seconds for this to complete, even though the worst of it is multi-threaded on a powerful 4C/8T machine. (roughly 100*104k size loop)

For any models where the number of alphas is very large (50K) it can take up to three minutes for a single execution to complete because of how large that double-nested for loop must become. (50k^2 size loop)

Pushing the list of betas onto their own vector can help in cases like that, but seems to significantly hurt performance of more normal cases.

Possible optimizations:

• Creating a sphere around all alphas to use in culling betas that are outside of the range of any alpha could potentially provide benefit, but it's an O(alphas^2) operation, and its benefit is extremely variable on the geometry.

• Creating a vector of Betas and pushing onto it as the alphas are also created seems to only benefit extreme edge cases like the 50k alpha case. In more "normal" cases of small numbers of alphas, it seems to hurt performance significantly.

• Adding to the list of valid alphas and culling it rather than re-building it each time may be an option, however, this will again be dependent on what % are alphas in the geometry.

• As well, it's possible something can be done with nanoRT's BVHs, but I'm not very familiar with BVH's or what they'd let me do in this

Note: How it's being used:

The MATLAB script will likely call this many times. In small models, it may finish its own loop within tenths of a second and then call ours again. In larger ones, there may be half a second between calls. In total, this may be called hundreds of times.

Note: How it's being built:

This isn't built using the MEX command in MATLAB, nor by using Visual Studio. Instead, g++ is used to create an object file (.o) and then g++ is used again to create the .mexw64 file in a method I'm not entirely familiar with. (This is also a hard limit I cannot touch)

I occasionally compiled with very aggressive warnings enabled to catch things like sign conversion, promotions, bad casts, etc.

Profiling:

I would love to be able to profile this code more in depth. However, it seems impossible. MEX files built using MEX command in MATLAB can be done. MEX files compiled in Visual Studio can be profiled. But we're not doing either of those, and so when I try to profile with either MATLAB or Visual Studio, it just doesn't work.

Even if I could, I don't think it would reveal anything surprising. The numbers we're working with are large, so the double-nested loop at the core of it grows very large.

I can (and do) measure per-invocation performance and total runtime after the MATLAB script completes. This is mostly stable, ~1% std dev in runtimes.

Final Note:

While performance is my most major concern, style improvements are always welcome. I'm more familiar with C than C++, and that bleeds into my code sometimes.

Answer 1

This review doesn't cover performance, but writing cleaner code:

---

Global variables are bad because it's hard to reason about them. They can be used anywhere, and (more frustratingly) at any time, making it harder to understand or change the code and its dependencies.

It would be much easier to understand each function if all the dependencies are passed into it (by value or reference as appropriate). e.g. as well as the existing function arguments, isTriAlpha also depends on triActive, hasBeenInZone, and whatever global or state isTriInZoneRadius also depends on.

While it may be necessary to declare variables at file / namespace scope in this I don't think there's a need to actually use them globally. e.g. They can be placed in a "Global" namespace that's only referred to in the top level function, and references to the relevant variables passed down as necessary.

---

This last point might seem like a burden, but grouping data together appropriately will make this much easier. e.g. it looks like all these contain one item per triangle:

// Data about each individual tri, could be brought intro a vector of structs
// Needed to check if geometry has changed since last invokation
std::vector<bool> triActive;
// Needed to check if alphas have changed since last invokation
std::vector<bool> validAlphaIndex;
// Needed to keep history of what tris have ever been in the beam, for alphas
std::vector<bool> hasBeenInBeam;

So perhaps they should all be in a struct TriangleData (or something), and we can pass a reference to it down through the function chain.

---

Prefer references to pointers as function arguments for "passing out" data. e.g. the numTri and numFace arguments to initialize should be references and not pointers. Pointers can be null, whereas references can only be created from a valid object. Since we don't check for a null value before dereferencing the pointers, it looks like references would be more appropriate.

---

It's better to use constant variables than defines. i.e. SUCCESS and PTHREAD_ERR should be:

static const int SUCCESS = 0;
static const int PTHREAD_ERR = 1;

Preprocessor definitions have no scoping, so they can affect your entire project (and any code that may use your project), so are prone to name collisions.

---

Declare variables as close to the point of use as possible and initialize them to the correct value straight away. e.g. in checkForModelChanges, currentlyActive and isAlpha should be declared and initialized inside the loop.

Unless constructing the variables does some very slow resource allocation it's best to let the compiler worry about optimization.

---

Comments should explain why the code does something, not just restate what the code does:

// Get whether this triangle is an alpha:
isAlpha = isTriAlpha(i, nodeIValues, iValueThreshold);

If we have to write a comment that says what the code does because it's not clear from the code itself, we should make the code clearer instead, e.g.:

// Get whether this triangle is an alpha:
isAlpha = isTriAlpha(i, nodeIValues, iValueThreshold);

// Triangle is a valid alpha now, but wasn't before
if((isAlpha == true) && (validAlphaIndex[i] == false))
{
  validAlphaIndex[i] = true;
  modelChanged = true;
}
// Was valid before, is no longer valid now
else if((isAlpha == false) && (validAlphaIndex[i] == true))
{
  validAlphaIndex[i] = false;
  modelChanged = true;
  //cullalphasFlag = true;
}

Could just be:

const bool wasAlpha = validAlphaIndex[i];
const bool isAlpha = isTriAlpha(i, nodeIValues, iValueThreshold);

if (wasAlpha != isAlpha) modelChanged = true;
validAlphaIndex[i] = isAlpha;

---

Don't test booleans by comparing them to true or false, just test the boolean directly:

if (isAlpha) { ... }
if (!isAlpha) { ... }

After all, the == operator returns a bool anyway...

if ((isAlpha == true) == true) { ... } // is it really, definitely true?

---

Similarly, something like this:

if(!(m_nrt->shootRay(ray)))
{
  return true;
}
else
{
  // There's no path
  return false;
}

is 8 lines of code, where we can really use just one:

return !m_nrt->shootRay(ray);

---

Prefer to return early where possible. This allows us to avoid unnecessary indentation and else clauses:

bool isTriAlpha(const unsigned int iTri,
                const float* nodeIValues,
                const double iValueThreshold
                )
{
    if (!triActive[iTri])
        return false;

    if (!hasBeenInZone[iTri] && !isInTriZoneRadius(iTri)))
        return false;

    const unsigned int* triNodes = m_nrt->getTriNodes(iTri);
    double triAvgIValue = (nodeIValues[triNodes[0]] + nodeIValues[triNodes[1]] + nodeIValues[triNodes[2]]) / 3.0;

    return (triAvgValue > iValueThreshold);
}

---

The rc variable here doesn't seem to have any reason to exist. We could just check the result of the function directly.

  int rc;

  if((rc = pthread_attr_init(&attr)))
  {
    throw PThreadException();
  }

---

Maybe split the initialization that needs to be done only once into a separate function from the initialization that is done every time and call these two functions only when appropriate. I'm guessing a lot of these checks are effectively checking the same thing:

 if(triActive.empty())
  {
    triActive.resize(numTris, false);
  }
  if(hasBeenInZone.empty())
  {
    hasBeenInZone.resize(numTris, false);
  }

---

The actual triangle data appears to be floats, but the calculations use a lot of doubles. Are the doubles actually necessary?

---

A decent math library would make stuff like this:

ray.dir[0] = (betaCoords[0] - alphaCoords[0]) / pathDist;
ray.dir[1] = (betaCoords[1] - alphaCoords[1]) / pathDist;
ray.dir[2] = (betaCoords[2] - alphaCoords[2]) / pathDist;

Look more like this:

ray.dir = (betaCoords - alphaCoords) / pathDist;

I wonder if it would be possible to do this with one of the libraries you're already using, instead of manually declaring arrays each time (e.g. double alphaCoords[3]; -> something like vec3<double> alphaCoords).

Answer 2

Here are a number of suggestions and comments that may help you improve your code.

Update your compiler

If you are actually limited to C++03, you're foregoing well over a decade's worth of compiler and language advancement that would have made this code much simpler and probably faster. For instance all of the pthread business could probably be much more elegantly handled using std::async and you'd be able to use references for efficiency and clarity. Without that, your path will be much harder and the code much less elegant and less robust than it should be.

Create and use a 3D point type

There are many instances in which 3 dimensional points are being used, but calculations for each is written out individually. Instead, I'd suggest that the code would be shorter, simpler and easier to read, understand and maintain if it used a Point3D class.

Think carefully about performance

The checkPairValid function is likely to be a performance bottleneck because of its use of floating point operations pow and sqrt. First consider these lines:

// Determine distance squared between alpha and beta
// (x2-x1)^2 + (y2-y1)^2 +(z2-z1)^2
pathDist = sqrt(pow((betaCoords[0] - alphaCoords[0]), 2)
            + pow((betaCoords[1] - alphaCoords[1]), 2)
            + pow((betaCoords[2] - alphaCoords[2]), 2));

The comment and the code don't match. In this case, I'd make them match by omitting sqrt (which should actually be std::sqrt). I'd also suggest that multiplication is likely to be faster than invoking pow (which should be std::pow). I'd use a templated 3D point class (as mentioned above) and define a function like this:

T squaredDist(const Point3D<T>& other) const {
    T dx = loc[0] - other.loc[0];
    T dy = loc[1] - other.loc[1];
    T dz = loc[2] - other.loc[2];
    return dx * dx + dy * dy + dz * dz;
}

Then you can compare with a squared threshold instead of the existing distThreshold for speed.

We also have these three lines:

ray.dir[0] = (betaCoords[0] - alphaCoords[0]) / pathDist;
ray.dir[1] = (betaCoords[1] - alphaCoords[1]) / pathDist;
ray.dir[2] = (betaCoords[2] - alphaCoords[2]) / pathDist;

If this is indeed intended to be a direction vector as the name suggests, it is probably not necessary to divide by the pathDist since it's the same direction either way. That would also save some calculation. In short, here's how I'd rewrite that function:

/**
 * @brief checkPairValid - Checks if a pair of points form a valid path
 * @param alpha - An alpha point
 * @param beta - A beta point
 * @param distThreshold - The square of the max distance apart the 
 * point's centers can be 
 * @param shoot - Function that returns false if there is a path
 * @return Whether the pair forms a valid path
 */
bool checkPairValid(const Point3D<double> &alpha,
                    const Point3D<double> &beta,
                    const double squaredDistThreshold,
                    bool (*shoot)(nanort::Ray<double>)
                    )
{
  double squaredPathDist = alpha.squaredDist(beta);
  if(squaredPathDist < squaredDistThreshold)
  {
    // Set up a nanort::Ray's origin, direction, and max distance
    nanort::Ray<double> ray(alpha, beta-alpha, std::sqrt(squaredPathDist));

    // Call passed shoot function to check for a path
    return !shoot(ray);
  }
  // The distance is too far between alpha and beta
  return false;
}

This is not only easier to read than the original but also no longer has any reliance on global variables.