diff --git a/interface/resources/qml/Stats.qml b/interface/resources/qml/Stats.qml
index f8ff461f6c..d5cfcff1e2 100644
--- a/interface/resources/qml/Stats.qml
+++ b/interface/resources/qml/Stats.qml
@@ -119,6 +119,22 @@ Item {
                         visible: root.expanded
                         text: "Avatars NOT Updated: " + root.notUpdatedAvatarCount
                     }
+                    StatText {
+                        visible: root.expanded
+                        text: "Total picks:\n    " +
+                                    root.stylusPicksCount + " styluses\n    " +
+                                    root.rayPicksCount + " rays\n    " +
+                                    root.parabolaPicksCount + " parabolas\n    " +
+                                    root.collisionPicksCount + " colliders"
+                    }
+                    StatText {
+                        visible: root.expanded
+                        text: "Intersection calls: Entities/Overlays/Avatars/HUD\n    " +
+                                    "Styluses:\t" + root.stylusPicksUpdated.x + "/" + root.stylusPicksUpdated.y + "/" + root.stylusPicksUpdated.z + "/" + root.stylusPicksUpdated.w + "\n    " +
+                                    "Rays:\t" + root.rayPicksUpdated.x + "/" + root.rayPicksUpdated.y + "/" + root.rayPicksUpdated.z + "/" + root.rayPicksUpdated.w + "\n    " +
+                                    "Parabolas:\t" + root.parabolaPicksUpdated.x + "/" + root.parabolaPicksUpdated.y + "/" + root.parabolaPicksUpdated.z + "/" + root.parabolaPicksUpdated.w + "\n    " +
+                                    "Colliders:\t" + root.collisionPicksUpdated.x + "/" + root.collisionPicksUpdated.y + "/" + root.collisionPicksUpdated.z + "/" + root.collisionPicksUpdated.w
+                    }
                 }
             }
 
diff --git a/interface/src/Application.cpp b/interface/src/Application.cpp
index 324ce6d833..c9912a979c 100644
--- a/interface/src/Application.cpp
+++ b/interface/src/Application.cpp
@@ -1760,11 +1760,6 @@ Application::Application(int& argc, char** argv, QElapsedTimer& startupTimer, bo
     QTimer* settingsTimer = new QTimer();
     moveToNewNamedThread(settingsTimer, "Settings Thread", [this, settingsTimer]{
         connect(qApp, &Application::beforeAboutToQuit, [this, settingsTimer]{
-            bool autoLogout = Setting::Handle<bool>(AUTO_LOGOUT_SETTING_NAME, false).get();
-            if (autoLogout) {
-                auto accountManager = DependencyManager::get<AccountManager>();
-                accountManager->logout();
-            }
             // Disconnect the signal from the save settings
             QObject::disconnect(settingsTimer, &QTimer::timeout, this, &Application::saveSettings);
             // Stop the settings timer
@@ -2516,6 +2511,11 @@ void Application::cleanupBeforeQuit() {
     }
     DependencyManager::destroy<ScriptEngines>();
 
+    bool autoLogout = Setting::Handle<bool>(AUTO_LOGOUT_SETTING_NAME, false).get();
+    if (autoLogout) {
+        DependencyManager::get<AccountManager>()->removeAccountFromFile();
+    }
+
     _displayPlugin.reset();
     PluginManager::getInstance()->shutdown();
 
diff --git a/interface/src/AvatarBookmarks.cpp b/interface/src/AvatarBookmarks.cpp
index 4119b843e1..5c79bedc9a 100644
--- a/interface/src/AvatarBookmarks.cpp
+++ b/interface/src/AvatarBookmarks.cpp
@@ -145,20 +145,9 @@ void AvatarBookmarks::removeBookmark(const QString& bookmarkName) {
     emit bookmarkDeleted(bookmarkName);
 }
 
-bool isWearableEntity(const EntityItemPointer& entity) {
-    return entity->isVisible() && (entity->getParentJointIndex() != INVALID_JOINT_INDEX || (entity->getType() == EntityTypes::Model && (std::static_pointer_cast<ModelEntityItem>(entity))->getRelayParentJoints()))
-        && (entity->getParentID() == DependencyManager::get<NodeList>()->getSessionUUID() || entity->getParentID() == DependencyManager::get<AvatarManager>()->getMyAvatar()->getSelfID());
-}
-
 void AvatarBookmarks::updateAvatarEntities(const QVariantList &avatarEntities) {
     auto myAvatar = DependencyManager::get<AvatarManager>()->getMyAvatar();
-    auto treeRenderer = DependencyManager::get<EntityTreeRenderer>();
-    EntityTreePointer entityTree = treeRenderer ? treeRenderer->getTree() : nullptr;
-    myAvatar->removeAvatarEntities([&](const QUuid& entityID) {
-        auto entity = entityTree->findEntityByID(entityID);
-        return entity && isWearableEntity(entity);
-    });
-
+    myAvatar->removeWearableAvatarEntities();
     addAvatarEntities(avatarEntities);
 }
 
@@ -183,10 +172,7 @@ void AvatarBookmarks::loadBookmark(const QString& bookmarkName) {
             auto myAvatar = DependencyManager::get<AvatarManager>()->getMyAvatar();
             auto treeRenderer = DependencyManager::get<EntityTreeRenderer>();
             EntityTreePointer entityTree = treeRenderer ? treeRenderer->getTree() : nullptr;
-            myAvatar->removeAvatarEntities([&](const QUuid& entityID) {
-                auto entity = entityTree->findEntityByID(entityID);
-                return entity && isWearableEntity(entity);
-            });
+            myAvatar->removeWearableAvatarEntities();
             const QString& avatarUrl = bookmark.value(ENTRY_AVATAR_URL, "").toString();
             myAvatar->useFullAvatarURL(avatarUrl);
             qCDebug(interfaceapp) << "Avatar On " << avatarUrl;
diff --git a/interface/src/Menu.cpp b/interface/src/Menu.cpp
index a7d7dcf8b3..a6ba983ab5 100644
--- a/interface/src/Menu.cpp
+++ b/interface/src/Menu.cpp
@@ -46,6 +46,7 @@
 #include "InterfaceLogging.h"
 #include "LocationBookmarks.h"
 #include "DeferredLightingEffect.h"
+#include "PickManager.h"
 
 #include "AmbientOcclusionEffect.h"
 #include "RenderShadowTask.h"
@@ -691,6 +692,11 @@ Menu::Menu() {
     addCheckableActionToQMenuAndActionHash(physicsOptionsMenu, MenuOption::PhysicsShowBulletConstraints, 0, false, qApp, SLOT(setShowBulletConstraints(bool)));
     addCheckableActionToQMenuAndActionHash(physicsOptionsMenu, MenuOption::PhysicsShowBulletConstraintLimits, 0, false, qApp, SLOT(setShowBulletConstraintLimits(bool)));
 
+    // Developer > Picking >>>
+    MenuWrapper* pickingOptionsMenu = developerMenu->addMenu("Picking");
+    addCheckableActionToQMenuAndActionHash(pickingOptionsMenu, MenuOption::ForceCoarsePicking, 0, false,
+        DependencyManager::get<PickManager>().data(), SLOT(setForceCoarsePicking(bool)));
+
     // Developer > Display Crash Options
     addCheckableActionToQMenuAndActionHash(developerMenu, MenuOption::DisplayCrashOptions, 0, true);
     // Developer > Crash >>>
diff --git a/interface/src/Menu.h b/interface/src/Menu.h
index c55bffdcb6..c4ea3734f5 100644
--- a/interface/src/Menu.h
+++ b/interface/src/Menu.h
@@ -221,6 +221,7 @@ namespace MenuOption {
     const QString NotificationSounds = "play_notification_sounds";
     const QString NotificationSoundsSnapshot = "play_notification_sounds_snapshot";
     const QString NotificationSoundsTablet = "play_notification_sounds_tablet";
+    const QString ForceCoarsePicking = "Force Coarse Picking";
     const QString ComputeBlendshapes = "Compute Blendshapes";
 }
 
diff --git a/interface/src/avatar/AvatarManager.cpp b/interface/src/avatar/AvatarManager.cpp
index 9f45482e84..9a7d8ef0c8 100644
--- a/interface/src/avatar/AvatarManager.cpp
+++ b/interface/src/avatar/AvatarManager.cpp
@@ -583,8 +583,14 @@ RayToAvatarIntersectionResult AvatarManager::findRayIntersectionVector(const Pic
         return result;
     }
 
-    glm::vec3 normDirection = glm::normalize(ray.direction);
+    // It's better to intersect the ray against the avatar's actual mesh, but this is currently difficult to
+    // do, because the transformed mesh data only exists over in GPU-land.  As a compromise, this code
+    // intersects against the avatars capsule and then against the (T-pose) mesh.  The end effect is that picking
+    // against the avatar is sort-of right, but you likely wont be able to pick against the arms.
 
+    // TODO -- find a way to extract transformed avatar mesh data from the rendering engine.
+
+    std::vector<SortedAvatar> sortedAvatars;
     auto avatarHashCopy = getHashCopy();
     for (auto avatarData : avatarHashCopy) {
         auto avatar = std::static_pointer_cast<Avatar>(avatarData);
@@ -593,52 +599,65 @@ RayToAvatarIntersectionResult AvatarManager::findRayIntersectionVector(const Pic
             continue;
         }
 
-        float distance;
-        BoxFace face;
-        glm::vec3 surfaceNormal;
-
-        SkeletonModelPointer avatarModel = avatar->getSkeletonModel();
-
-        // It's better to intersect the ray against the avatar's actual mesh, but this is currently difficult to
-        // do, because the transformed mesh data only exists over in GPU-land.  As a compromise, this code
-        // intersects against the avatars capsule and then against the (T-pose) mesh.  The end effect is that picking
-        // against the avatar is sort-of right, but you likely wont be able to pick against the arms.
-
-        // TODO -- find a way to extract transformed avatar mesh data from the rendering engine.
-
+        float distance = FLT_MAX;
+#if 0
         // if we weren't picking against the capsule, we would want to pick against the avatarBounds...
-        // AABox avatarBounds = avatarModel->getRenderableMeshBound();
-        // if (!avatarBounds.findRayIntersection(ray.origin, normDirection, distance, face, surfaceNormal)) {
-        //     // ray doesn't intersect avatar's bounding-box
-        //     continue;
-        // }
-
+        SkeletonModelPointer avatarModel = avatar->getSkeletonModel();
+        AABox avatarBounds = avatarModel->getRenderableMeshBound();
+        if (avatarBounds.contains(ray.origin)) {
+            distance = 0.0f;
+        } else {
+            float boundDistance = FLT_MAX;
+            BoxFace face;
+            glm::vec3 surfaceNormal;
+            if (avatarBounds.findRayIntersection(ray.origin, ray.direction, boundDistance, face, surfaceNormal)) {
+                distance = boundDistance;
+            }
+        }
+#else
         glm::vec3 start;
         glm::vec3 end;
         float radius;
         avatar->getCapsule(start, end, radius);
-        bool intersects = findRayCapsuleIntersection(ray.origin, normDirection, start, end, radius, distance);
-        if (!intersects) {
-            // ray doesn't intersect avatar's capsule
-            continue;
+        findRayCapsuleIntersection(ray.origin, ray.direction, start, end, radius, distance);
+#endif
+
+        if (distance < FLT_MAX) {
+            sortedAvatars.emplace_back(distance, avatar);
+        }
+    }
+
+    if (sortedAvatars.size() > 1) {
+        static auto comparator = [](const SortedAvatar& left, const SortedAvatar& right) { return left.first < right.first; };
+        std::sort(sortedAvatars.begin(), sortedAvatars.end(), comparator);
+    }
+
+    for (auto it = sortedAvatars.begin(); it != sortedAvatars.end(); ++it) {
+        const SortedAvatar& sortedAvatar = *it;
+        // We can exit once avatarCapsuleDistance > bestDistance
+        if (sortedAvatar.first > result.distance) {
+            break;
         }
 
+        float distance = FLT_MAX;
+        BoxFace face;
+        glm::vec3 surfaceNormal;
         QVariantMap extraInfo;
-        intersects = avatarModel->findRayIntersectionAgainstSubMeshes(ray.origin, normDirection,
-                                                                      distance, face, surfaceNormal, extraInfo, true);
-
-        if (intersects && (!result.intersects || distance < result.distance)) {
-            result.intersects = true;
-            result.avatarID = avatar->getID();
-            result.distance = distance;
-            result.face = face;
-            result.surfaceNormal = surfaceNormal;
-            result.extraInfo = extraInfo;
+        SkeletonModelPointer avatarModel = sortedAvatar.second->getSkeletonModel();
+        if (avatarModel->findRayIntersectionAgainstSubMeshes(ray.origin, ray.direction, distance, face, surfaceNormal, extraInfo, true)) {
+            if (distance < result.distance) {
+                result.intersects = true;
+                result.avatarID = sortedAvatar.second->getID();
+                result.distance = distance;
+                result.face = face;
+                result.surfaceNormal = surfaceNormal;
+                result.extraInfo = extraInfo;
+            }
         }
     }
 
     if (result.intersects) {
-        result.intersection = ray.origin + normDirection * result.distance;
+        result.intersection = ray.origin + ray.direction * result.distance;
     }
 
     return result;
@@ -657,6 +676,14 @@ ParabolaToAvatarIntersectionResult AvatarManager::findParabolaIntersectionVector
         return result;
     }
 
+    // It's better to intersect the ray against the avatar's actual mesh, but this is currently difficult to
+    // do, because the transformed mesh data only exists over in GPU-land.  As a compromise, this code
+    // intersects against the avatars capsule and then against the (T-pose) mesh.  The end effect is that picking
+    // against the avatar is sort-of right, but you likely wont be able to pick against the arms.
+
+    // TODO -- find a way to extract transformed avatar mesh data from the rendering engine.
+
+    std::vector<SortedAvatar> sortedAvatars;
     auto avatarHashCopy = getHashCopy();
     for (auto avatarData : avatarHashCopy) {
         auto avatar = std::static_pointer_cast<Avatar>(avatarData);
@@ -665,47 +692,60 @@ ParabolaToAvatarIntersectionResult AvatarManager::findParabolaIntersectionVector
             continue;
         }
 
-        float parabolicDistance;
-        BoxFace face;
-        glm::vec3 surfaceNormal;
-
-        SkeletonModelPointer avatarModel = avatar->getSkeletonModel();
-
-        // It's better to intersect the parabola against the avatar's actual mesh, but this is currently difficult to
-        // do, because the transformed mesh data only exists over in GPU-land.  As a compromise, this code
-        // intersects against the avatars capsule and then against the (T-pose) mesh.  The end effect is that picking
-        // against the avatar is sort-of right, but you likely wont be able to pick against the arms.
-
-        // TODO -- find a way to extract transformed avatar mesh data from the rendering engine.
-
+        float distance = FLT_MAX;
+#if 0
         // if we weren't picking against the capsule, we would want to pick against the avatarBounds...
-        // AABox avatarBounds = avatarModel->getRenderableMeshBound();
-        // if (!avatarBounds.findParabolaIntersection(pick.origin, pick.velocity, pick.acceleration, parabolicDistance, face, surfaceNormal)) {
-        //     // parabola doesn't intersect avatar's bounding-box
-        //     continue;
-        // }
-
+        SkeletonModelPointer avatarModel = avatar->getSkeletonModel();
+        AABox avatarBounds = avatarModel->getRenderableMeshBound();
+        if (avatarBounds.contains(pick.origin)) {
+            distance = 0.0f;
+        } else {
+            float boundDistance = FLT_MAX;
+            BoxFace face;
+            glm::vec3 surfaceNormal;
+            if (avatarBounds.findParabolaIntersection(pick.origin, pick.velocity, pick.acceleration, boundDistance, face, surfaceNormal)) {
+                distance = boundDistance;
+            }
+        }
+#else
         glm::vec3 start;
         glm::vec3 end;
         float radius;
         avatar->getCapsule(start, end, radius);
-        bool intersects = findParabolaCapsuleIntersection(pick.origin, pick.velocity, pick.acceleration, start, end, radius, avatar->getWorldOrientation(), parabolicDistance);
-        if (!intersects) {
-            // ray doesn't intersect avatar's capsule
-            continue;
+        findParabolaCapsuleIntersection(pick.origin, pick.velocity, pick.acceleration, start, end, radius, avatar->getWorldOrientation(), distance);
+#endif
+
+        if (distance < FLT_MAX) {
+            sortedAvatars.emplace_back(distance, avatar);
+        }
+    }
+
+    if (sortedAvatars.size() > 1) {
+        static auto comparator = [](const SortedAvatar& left, const SortedAvatar& right) { return left.first < right.first; };
+        std::sort(sortedAvatars.begin(), sortedAvatars.end(), comparator);
+    }
+
+    for (auto it = sortedAvatars.begin(); it != sortedAvatars.end(); ++it) {
+        const SortedAvatar& sortedAvatar = *it;
+        // We can exit once avatarCapsuleDistance > bestDistance
+        if (sortedAvatar.first > result.parabolicDistance) {
+            break;
         }
 
+        float parabolicDistance = FLT_MAX;
+        BoxFace face;
+        glm::vec3 surfaceNormal;
         QVariantMap extraInfo;
-        intersects = avatarModel->findParabolaIntersectionAgainstSubMeshes(pick.origin, pick.velocity, pick.acceleration,
-                                                                           parabolicDistance, face, surfaceNormal, extraInfo, true);
-
-        if (intersects && (!result.intersects || parabolicDistance < result.parabolicDistance)) {
-            result.intersects = true;
-            result.avatarID = avatar->getID();
-            result.parabolicDistance = parabolicDistance;
-            result.face = face;
-            result.surfaceNormal = surfaceNormal;
-            result.extraInfo = extraInfo;
+        SkeletonModelPointer avatarModel = sortedAvatar.second->getSkeletonModel();
+        if (avatarModel->findParabolaIntersectionAgainstSubMeshes(pick.origin, pick.velocity, pick.acceleration, parabolicDistance, face, surfaceNormal, extraInfo, true)) {
+            if (parabolicDistance < result.parabolicDistance) {
+                result.intersects = true;
+                result.avatarID = sortedAvatar.second->getID();
+                result.parabolicDistance = parabolicDistance;
+                result.face = face;
+                result.surfaceNormal = surfaceNormal;
+                result.extraInfo = extraInfo;
+            }
         }
     }
 
diff --git a/interface/src/avatar/AvatarManager.h b/interface/src/avatar/AvatarManager.h
index 36df2f0aaf..bcdfc064bd 100644
--- a/interface/src/avatar/AvatarManager.h
+++ b/interface/src/avatar/AvatarManager.h
@@ -31,6 +31,8 @@
 #include "MyAvatar.h"
 #include "OtherAvatar.h"
 
+using SortedAvatar = std::pair<float, std::shared_ptr<Avatar>>;
+
 /**jsdoc 
  * The AvatarManager API has properties and methods which manage Avatars within the same domain.
  *
diff --git a/interface/src/avatar/MyAvatar.cpp b/interface/src/avatar/MyAvatar.cpp
index 923d071cc0..7bf18f468c 100755
--- a/interface/src/avatar/MyAvatar.cpp
+++ b/interface/src/avatar/MyAvatar.cpp
@@ -1703,18 +1703,50 @@ void MyAvatar::setSkeletonModelURL(const QUrl& skeletonModelURL) {
     emit skeletonChanged();
 }
 
-void MyAvatar::removeAvatarEntities(const std::function<bool(const QUuid& entityID)>& condition) {
+bool isWearableEntity(const EntityItemPointer& entity) {
+    return entity->isVisible()
+        && (entity->getParentJointIndex() != INVALID_JOINT_INDEX
+            || (entity->getType() == EntityTypes::Model && (std::static_pointer_cast<ModelEntityItem>(entity))->getRelayParentJoints()))
+        && (entity->getParentID() == DependencyManager::get<NodeList>()->getSessionUUID()
+            || entity->getParentID() == AVATAR_SELF_ID);
+}
+
+void MyAvatar::clearAvatarEntities() {
     auto treeRenderer = DependencyManager::get<EntityTreeRenderer>();
     EntityTreePointer entityTree = treeRenderer ? treeRenderer->getTree() : nullptr;
-    if (entityTree) {
-        entityTree->withWriteLock([&] {
-            AvatarEntityMap avatarEntities = getAvatarEntityData();
-            for (auto entityID : avatarEntities.keys()) {
-                if (!condition || condition(entityID)) {
-                    entityTree->deleteEntity(entityID, true, true);
-                }
-            }
+
+    AvatarEntityMap avatarEntities = getAvatarEntityData();
+    for (auto entityID : avatarEntities.keys()) {
+        entityTree->withWriteLock([&entityID, &entityTree] {
+            // remove this entity first from the entity tree
+            entityTree->deleteEntity(entityID, true, true);
         });
+
+        // remove the avatar entity from our internal list
+        // (but indicate it doesn't need to be pulled from the tree)
+        clearAvatarEntity(entityID, false);
+    }
+}
+
+void MyAvatar::removeWearableAvatarEntities() {
+    auto treeRenderer = DependencyManager::get<EntityTreeRenderer>();
+    EntityTreePointer entityTree = treeRenderer ? treeRenderer->getTree() : nullptr;
+    
+    if (entityTree) {
+        AvatarEntityMap avatarEntities = getAvatarEntityData();
+        for (auto entityID : avatarEntities.keys()) {
+            auto entity = entityTree->findEntityByID(entityID);
+            if (entity && isWearableEntity(entity)) {
+                entityTree->withWriteLock([&entityID, &entityTree] {
+                    // remove this entity first from the entity tree
+                    entityTree->deleteEntity(entityID, true, true);
+                });
+
+                // remove the avatar entity from our internal list
+                // (but indicate it doesn't need to be pulled from the tree)
+                clearAvatarEntity(entityID, false);
+            }
+        }
     }
 }
 
@@ -2116,7 +2148,7 @@ void MyAvatar::setAttachmentData(const QVector<AttachmentData>& attachmentData)
     }
 
     // clear any existing avatar entities
-    setAvatarEntityData(AvatarEntityMap());
+    clearAvatarEntities();
 
     for (auto& properties : newEntitiesProperties) {
         DependencyManager::get<EntityScriptingInterface>()->addEntity(properties, true);
diff --git a/interface/src/avatar/MyAvatar.h b/interface/src/avatar/MyAvatar.h
index ba6348cc22..bb0df1af62 100644
--- a/interface/src/avatar/MyAvatar.h
+++ b/interface/src/avatar/MyAvatar.h
@@ -931,7 +931,8 @@ public:
     * @returns {object[]}
     */
     Q_INVOKABLE QVariantList getAvatarEntitiesVariant();
-    void removeAvatarEntities(const std::function<bool(const QUuid& entityID)>& condition = {});
+    void clearAvatarEntities();
+    void removeWearableAvatarEntities();
 
     /**jsdoc
      * @function MyAvatar.isFlying
@@ -1782,4 +1783,6 @@ void audioListenModeFromScriptValue(const QScriptValue& object, AudioListenerMod
 QScriptValue driveKeysToScriptValue(QScriptEngine* engine, const MyAvatar::DriveKeys& driveKeys);
 void driveKeysFromScriptValue(const QScriptValue& object, MyAvatar::DriveKeys& driveKeys);
 
+bool isWearableEntity(const EntityItemPointer& entity);
+
 #endif // hifi_MyAvatar_h
diff --git a/interface/src/raypick/LaserPointer.cpp b/interface/src/raypick/LaserPointer.cpp
index 2382a95105..40a5e5cb69 100644
--- a/interface/src/raypick/LaserPointer.cpp
+++ b/interface/src/raypick/LaserPointer.cpp
@@ -42,6 +42,9 @@ glm::vec3 LaserPointer::getPickOrigin(const PickResultPointer& pickResult) const
 
 glm::vec3 LaserPointer::getPickEnd(const PickResultPointer& pickResult, float distance) const {
     auto rayPickResult = std::static_pointer_cast<RayPickResult>(pickResult);
+    if (!rayPickResult) {
+        return glm::vec3(0.0f);
+    }
     if (distance > 0.0f) {
         PickRay pick = PickRay(rayPickResult->pickVariant);
         return pick.origin + distance * pick.direction;
diff --git a/interface/src/raypick/ParabolaPick.cpp b/interface/src/raypick/ParabolaPick.cpp
index b3e3f16345..f4db1df25a 100644
--- a/interface/src/raypick/ParabolaPick.cpp
+++ b/interface/src/raypick/ParabolaPick.cpp
@@ -13,11 +13,13 @@
 #include "avatar/AvatarManager.h"
 #include "scripting/HMDScriptingInterface.h"
 #include "DependencyManager.h"
+#include "PickManager.h"
 
 PickResultPointer ParabolaPick::getEntityIntersection(const PickParabola& pick) {
     if (glm::length2(pick.acceleration) > EPSILON && glm::length2(pick.velocity) > EPSILON) {
+        bool precisionPicking = !(getFilter().doesPickCoarse() || DependencyManager::get<PickManager>()->getForceCoarsePicking());
         ParabolaToEntityIntersectionResult entityRes =
-            DependencyManager::get<EntityScriptingInterface>()->findParabolaIntersectionVector(pick, !getFilter().doesPickCoarse(),
+            DependencyManager::get<EntityScriptingInterface>()->findParabolaIntersectionVector(pick, precisionPicking,
                 getIncludeItemsAs<EntityItemID>(), getIgnoreItemsAs<EntityItemID>(), !getFilter().doesPickInvisible(), !getFilter().doesPickNonCollidable());
         if (entityRes.intersects) {
             return std::make_shared<ParabolaPickResult>(IntersectionType::ENTITY, entityRes.entityID, entityRes.distance, entityRes.parabolicDistance, entityRes.intersection, pick, entityRes.surfaceNormal, entityRes.extraInfo);
@@ -28,8 +30,9 @@ PickResultPointer ParabolaPick::getEntityIntersection(const PickParabola& pick)
 
 PickResultPointer ParabolaPick::getOverlayIntersection(const PickParabola& pick) {
     if (glm::length2(pick.acceleration) > EPSILON && glm::length2(pick.velocity) > EPSILON) {
+        bool precisionPicking = !(getFilter().doesPickCoarse() || DependencyManager::get<PickManager>()->getForceCoarsePicking());
         ParabolaToOverlayIntersectionResult overlayRes =
-            qApp->getOverlays().findParabolaIntersectionVector(pick, !getFilter().doesPickCoarse(),
+            qApp->getOverlays().findParabolaIntersectionVector(pick, precisionPicking,
                 getIncludeItemsAs<OverlayID>(), getIgnoreItemsAs<OverlayID>(), !getFilter().doesPickInvisible(), !getFilter().doesPickNonCollidable());
         if (overlayRes.intersects) {
             return std::make_shared<ParabolaPickResult>(IntersectionType::OVERLAY, overlayRes.overlayID, overlayRes.distance, overlayRes.parabolicDistance, overlayRes.intersection, pick, overlayRes.surfaceNormal, overlayRes.extraInfo);
diff --git a/interface/src/raypick/ParabolaPointer.cpp b/interface/src/raypick/ParabolaPointer.cpp
index 097c98340c..57d57e11c4 100644
--- a/interface/src/raypick/ParabolaPointer.cpp
+++ b/interface/src/raypick/ParabolaPointer.cpp
@@ -67,6 +67,9 @@ glm::vec3 ParabolaPointer::getPickOrigin(const PickResultPointer& pickResult) co
 
 glm::vec3 ParabolaPointer::getPickEnd(const PickResultPointer& pickResult, float distance) const {
     auto parabolaPickResult = std::static_pointer_cast<ParabolaPickResult>(pickResult);
+    if (!parabolaPickResult) {
+        return glm::vec3(0.0f);
+    }
     if (distance > 0.0f) {
         PickParabola pick = PickParabola(parabolaPickResult->pickVariant);
         return pick.origin + pick.velocity * distance + 0.5f * pick.acceleration * distance * distance;
diff --git a/interface/src/raypick/RayPick.cpp b/interface/src/raypick/RayPick.cpp
index 96b41dcc72..fde1da3f87 100644
--- a/interface/src/raypick/RayPick.cpp
+++ b/interface/src/raypick/RayPick.cpp
@@ -13,10 +13,12 @@
 #include "avatar/AvatarManager.h"
 #include "scripting/HMDScriptingInterface.h"
 #include "DependencyManager.h"
+#include "PickManager.h"
 
 PickResultPointer RayPick::getEntityIntersection(const PickRay& pick) {
+    bool precisionPicking = !(getFilter().doesPickCoarse() || DependencyManager::get<PickManager>()->getForceCoarsePicking());
     RayToEntityIntersectionResult entityRes =
-        DependencyManager::get<EntityScriptingInterface>()->findRayIntersectionVector(pick, !getFilter().doesPickCoarse(),
+        DependencyManager::get<EntityScriptingInterface>()->findRayIntersectionVector(pick, precisionPicking,
             getIncludeItemsAs<EntityItemID>(), getIgnoreItemsAs<EntityItemID>(), !getFilter().doesPickInvisible(), !getFilter().doesPickNonCollidable());
     if (entityRes.intersects) {
         return std::make_shared<RayPickResult>(IntersectionType::ENTITY, entityRes.entityID, entityRes.distance, entityRes.intersection, pick, entityRes.surfaceNormal, entityRes.extraInfo);
@@ -26,8 +28,9 @@ PickResultPointer RayPick::getEntityIntersection(const PickRay& pick) {
 }
 
 PickResultPointer RayPick::getOverlayIntersection(const PickRay& pick) {
+    bool precisionPicking = !(getFilter().doesPickCoarse() || DependencyManager::get<PickManager>()->getForceCoarsePicking());
     RayToOverlayIntersectionResult overlayRes =
-        qApp->getOverlays().findRayIntersectionVector(pick, !getFilter().doesPickCoarse(),
+        qApp->getOverlays().findRayIntersectionVector(pick, precisionPicking,
             getIncludeItemsAs<OverlayID>(), getIgnoreItemsAs<OverlayID>(), !getFilter().doesPickInvisible(), !getFilter().doesPickNonCollidable());
     if (overlayRes.intersects) {
         return std::make_shared<RayPickResult>(IntersectionType::OVERLAY, overlayRes.overlayID, overlayRes.distance, overlayRes.intersection, pick, overlayRes.surfaceNormal, overlayRes.extraInfo);
diff --git a/interface/src/ui/Stats.cpp b/interface/src/ui/Stats.cpp
index f9a96077f5..3fa712d267 100644
--- a/interface/src/ui/Stats.cpp
+++ b/interface/src/ui/Stats.cpp
@@ -26,6 +26,7 @@
 #include <OffscreenUi.h>
 #include <PerfStat.h>
 #include <plugins/DisplayPlugin.h>
+#include <PickManager.h>
 
 #include <gl/Context.h>
 
@@ -147,6 +148,20 @@ void Stats::updateStats(bool force) {
     }
     STAT_UPDATE(gameLoopRate, (int)qApp->getGameLoopRate());
 
+    auto pickManager = DependencyManager::get<PickManager>();
+    if (pickManager && (_expanded || force)) {
+        std::vector<int> totalPicks = pickManager->getTotalPickCounts();
+        STAT_UPDATE(stylusPicksCount, totalPicks[PickQuery::Stylus]);
+        STAT_UPDATE(rayPicksCount, totalPicks[PickQuery::Ray]);
+        STAT_UPDATE(parabolaPicksCount, totalPicks[PickQuery::Parabola]);
+        STAT_UPDATE(collisionPicksCount, totalPicks[PickQuery::Collision]);
+        std::vector<QVector4D> updatedPicks = pickManager->getUpdatedPickCounts();
+        STAT_UPDATE(stylusPicksUpdated, updatedPicks[PickQuery::Stylus]);
+        STAT_UPDATE(rayPicksUpdated, updatedPicks[PickQuery::Ray]);
+        STAT_UPDATE(parabolaPicksUpdated, updatedPicks[PickQuery::Parabola]);
+        STAT_UPDATE(collisionPicksUpdated, updatedPicks[PickQuery::Collision]);
+    }
+
     auto bandwidthRecorder = DependencyManager::get<BandwidthRecorder>();
     STAT_UPDATE(packetInCount, (int)bandwidthRecorder->getCachedTotalAverageInputPacketsPerSecond());
     STAT_UPDATE(packetOutCount, (int)bandwidthRecorder->getCachedTotalAverageOutputPacketsPerSecond());
@@ -286,7 +301,7 @@ void Stats::updateStats(bool force) {
         //    downloads << (int)(resource->getProgress() * 100.0f) << "% ";
         //}
         //downloads << "(" <<  << " pending)";
-    } // expanded avatar column
+    }
 
     // Fourth column, octree stats
     int serverCount = 0;
diff --git a/interface/src/ui/Stats.h b/interface/src/ui/Stats.h
index 998690f302..ed88f4fc6b 100644
--- a/interface/src/ui/Stats.h
+++ b/interface/src/ui/Stats.h
@@ -22,7 +22,6 @@ public: \
 private: \
     type _##name{ initialValue };
 
-
 /**jsdoc
  * @namespace Stats
  *
@@ -169,6 +168,15 @@ private: \
  * @property {number} implicitHeight
  *
  * @property {object} layer - <em>Read-only.</em>
+
+ * @property {number} stylusPicksCount - <em>Read-only.</em>
+ * @property {number} rayPicksCount - <em>Read-only.</em>
+ * @property {number} parabolaPicksCount - <em>Read-only.</em>
+ * @property {number} collisionPicksCount - <em>Read-only.</em>
+ * @property {Vec4} stylusPicksUpdated - <em>Read-only.</em>
+ * @property {Vec4} rayPicksUpdated - <em>Read-only.</em>
+ * @property {Vec4} parabolaPicksUpdated - <em>Read-only.</em>
+ * @property {Vec4} collisionPicksUpdated - <em>Read-only.</em>
  */
 // Properties from x onwards are QQuickItem properties.
 
@@ -288,6 +296,15 @@ class Stats : public QQuickItem {
     STATS_PROPERTY(float, avatarSimulationTime, 0)
     Q_PROPERTY(QStringList animStackNames READ animStackNames NOTIFY animStackNamesChanged)
 
+    STATS_PROPERTY(int, stylusPicksCount, 0)
+    STATS_PROPERTY(int, rayPicksCount, 0)
+    STATS_PROPERTY(int, parabolaPicksCount, 0)
+    STATS_PROPERTY(int, collisionPicksCount, 0)
+    STATS_PROPERTY(QVector4D, stylusPicksUpdated, QVector4D(0, 0, 0, 0))
+    STATS_PROPERTY(QVector4D, rayPicksUpdated, QVector4D(0, 0, 0, 0))
+    STATS_PROPERTY(QVector4D, parabolaPicksUpdated, QVector4D(0, 0, 0, 0))
+    STATS_PROPERTY(QVector4D, collisionPicksUpdated, QVector4D(0, 0, 0, 0))
+
 public:
     static Stats* getInstance();
 
@@ -1263,6 +1280,62 @@ signals:
      * @function Stats.update
      */
 
+    /**jsdoc
+     * Triggered when the value of the <code>stylusPicksCount</code> property changes.
+     * @function Stats.stylusPicksCountChanged
+     * @returns {Signal}
+     */
+    void stylusPicksCountChanged();
+
+    /**jsdoc
+     * Triggered when the value of the <code>rayPicksCount</code> property changes.
+     * @function Stats.rayPicksCountChanged
+     * @returns {Signal}
+     */
+    void rayPicksCountChanged();
+
+    /**jsdoc
+     * Triggered when the value of the <code>parabolaPicksCount</code> property changes.
+     * @function Stats.parabolaPicksCountChanged
+     * @returns {Signal}
+     */
+    void parabolaPicksCountChanged();
+
+    /**jsdoc
+     * Triggered when the value of the <code>collisionPicksCount</code> property changes.
+     * @function Stats.collisionPicksCountChanged
+     * @returns {Signal}
+     */
+    void collisionPicksCountChanged();
+
+    /**jsdoc
+     * Triggered when the value of the <code>stylusPicksUpdated</code> property changes.
+     * @function Stats.stylusPicksUpdatedChanged
+     * @returns {Signal}
+     */
+    void stylusPicksUpdatedChanged();
+
+    /**jsdoc
+     * Triggered when the value of the <code>rayPicksUpdated</code> property changes.
+     * @function Stats.rayPicksUpdatedChanged
+     * @returns {Signal}
+     */
+    void rayPicksUpdatedChanged();
+
+    /**jsdoc
+     * Triggered when the value of the <code>parabolaPicksUpdated</code> property changes.
+     * @function Stats.parabolaPicksUpdatedChanged
+     * @returns {Signal}
+     */
+    void parabolaPicksUpdatedChanged();
+
+    /**jsdoc
+     * Triggered when the value of the <code>collisionPicksUpdated</code> property changes.
+     * @function Stats.collisionPicksUpdatedChanged
+     * @returns {Signal}
+     */
+    void collisionPicksUpdatedChanged();
+
 private:
     int _recentMaxPackets{ 0 } ; // recent max incoming voxel packets to process
     bool _resetRecentMaxPacketsSoon{ true };
diff --git a/interface/src/ui/overlays/ContextOverlayInterface.cpp b/interface/src/ui/overlays/ContextOverlayInterface.cpp
index 789b1f9969..ba9a1f9fc9 100644
--- a/interface/src/ui/overlays/ContextOverlayInterface.cpp
+++ b/interface/src/ui/overlays/ContextOverlayInterface.cpp
@@ -180,7 +180,7 @@ bool ContextOverlayInterface::createOrDestroyContextOverlay(const EntityItemID&
                 float distance;
                 BoxFace face;
                 glm::vec3 normal;
-                boundingBox.findRayIntersection(cameraPosition, direction, distance, face, normal);
+                boundingBox.findRayIntersection(cameraPosition, direction, 1.0f / direction, distance, face, normal);
                 float offsetAngle = -CONTEXT_OVERLAY_OFFSET_ANGLE;
                 if (event.getID() == 1) { // "1" is left hand
                     offsetAngle *= -1.0f;
diff --git a/interface/src/ui/overlays/Volume3DOverlay.cpp b/interface/src/ui/overlays/Volume3DOverlay.cpp
index c87650a77b..a307d445c0 100644
--- a/interface/src/ui/overlays/Volume3DOverlay.cpp
+++ b/interface/src/ui/overlays/Volume3DOverlay.cpp
@@ -88,7 +88,7 @@ bool Volume3DOverlay::findRayIntersection(const glm::vec3& origin, const glm::ve
 
     // we can use the AABox's ray intersection by mapping our origin and direction into the overlays frame
     // and testing intersection there.
-    bool hit = _localBoundingBox.findRayIntersection(overlayFrameOrigin, overlayFrameDirection, distance, face, surfaceNormal);
+    bool hit = _localBoundingBox.findRayIntersection(overlayFrameOrigin, overlayFrameDirection, 1.0f / overlayFrameDirection, distance, face, surfaceNormal);
 
     if (hit) {
         surfaceNormal = transform.getRotation() * surfaceNormal;
diff --git a/libraries/avatars/src/AvatarData.h b/libraries/avatars/src/AvatarData.h
index bb3b1fe2d3..0474d07acd 100644
--- a/libraries/avatars/src/AvatarData.h
+++ b/libraries/avatars/src/AvatarData.h
@@ -1565,7 +1565,7 @@ class RayToAvatarIntersectionResult {
 public:
     bool intersects { false };
     QUuid avatarID;
-    float distance { 0.0f };
+    float distance { FLT_MAX };
     BoxFace face;
     glm::vec3 intersection;
     glm::vec3 surfaceNormal;
diff --git a/libraries/entities-renderer/src/RenderablePolyVoxEntityItem.cpp b/libraries/entities-renderer/src/RenderablePolyVoxEntityItem.cpp
index c11ccb70a0..2da5c30dc0 100644
--- a/libraries/entities-renderer/src/RenderablePolyVoxEntityItem.cpp
+++ b/libraries/entities-renderer/src/RenderablePolyVoxEntityItem.cpp
@@ -571,7 +571,6 @@ bool RenderablePolyVoxEntityItem::findDetailedRayIntersection(const glm::vec3& o
     }
 
     glm::mat4 wtvMatrix = worldToVoxelMatrix();
-    glm::mat4 vtwMatrix = voxelToWorldMatrix();
     glm::vec3 normDirection = glm::normalize(direction);
 
     // the PolyVox ray intersection code requires a near and far point.
@@ -584,8 +583,6 @@ bool RenderablePolyVoxEntityItem::findDetailedRayIntersection(const glm::vec3& o
     glm::vec4 originInVoxel = wtvMatrix * glm::vec4(origin, 1.0f);
     glm::vec4 farInVoxel = wtvMatrix * glm::vec4(farPoint, 1.0f);
 
-    glm::vec4 directionInVoxel = glm::normalize(farInVoxel - originInVoxel);
-
     glm::vec4 result = glm::vec4(0.0f, 0.0f, 0.0f, 0.0f);
     PolyVox::RaycastResult raycastResult = doRayCast(originInVoxel, farInVoxel, result);
     if (raycastResult == PolyVox::RaycastResults::Completed) {
@@ -599,14 +596,9 @@ bool RenderablePolyVoxEntityItem::findDetailedRayIntersection(const glm::vec3& o
     voxelBox += result3 - Vectors::HALF;
     voxelBox += result3 + Vectors::HALF;
 
-    float voxelDistance;
-    bool hit = voxelBox.findRayIntersection(glm::vec3(originInVoxel), glm::vec3(directionInVoxel),
-                                            voxelDistance, face, surfaceNormal);
-
-    glm::vec4 voxelIntersectionPoint = glm::vec4(glm::vec3(originInVoxel) + glm::vec3(directionInVoxel) * voxelDistance, 1.0);
-    glm::vec4 intersectionPoint = vtwMatrix * voxelIntersectionPoint;
-    distance = glm::distance(origin, glm::vec3(intersectionPoint));
-    return hit;
+    glm::vec3 directionInVoxel = vec3(wtvMatrix * glm::vec4(direction, 0.0f));
+    return voxelBox.findRayIntersection(glm::vec3(originInVoxel), directionInVoxel, 1.0f / directionInVoxel,
+                                        distance, face, surfaceNormal);
 }
 
 bool RenderablePolyVoxEntityItem::findDetailedParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity,
diff --git a/libraries/entities/src/EntityTree.cpp b/libraries/entities/src/EntityTree.cpp
index 377e192bb1..a7c88ddc7d 100644
--- a/libraries/entities/src/EntityTree.cpp
+++ b/libraries/entities/src/EntityTree.cpp
@@ -48,6 +48,7 @@ public:
     // Inputs
     glm::vec3 origin;
     glm::vec3 direction;
+    glm::vec3 invDirection;
     const QVector<EntityItemID>& entityIdsToInclude;
     const QVector<EntityItemID>& entityIdsToDiscard;
     bool visibleOnly;
@@ -825,28 +826,51 @@ bool findRayIntersectionOp(const OctreeElementPointer& element, void* extraData)
     RayArgs* args = static_cast<RayArgs*>(extraData);
     bool keepSearching = true;
     EntityTreeElementPointer entityTreeElementPointer = std::static_pointer_cast<EntityTreeElement>(element);
-    EntityItemID entityID = entityTreeElementPointer->findRayIntersection(args->origin, args->direction, keepSearching,
+    EntityItemID entityID = entityTreeElementPointer->findRayIntersection(args->origin, args->direction,
         args->element, args->distance, args->face, args->surfaceNormal, args->entityIdsToInclude,
         args->entityIdsToDiscard, args->visibleOnly, args->collidableOnly, args->extraInfo, args->precisionPicking);
     if (!entityID.isNull()) {
         args->entityID = entityID;
+        // We recurse OctreeElements in order, so if we hit something, we can stop immediately
+        keepSearching = false;
     }
     return keepSearching;
 }
 
+float findRayIntersectionSortingOp(const OctreeElementPointer& element, void* extraData) {
+    RayArgs* args = static_cast<RayArgs*>(extraData);
+    EntityTreeElementPointer entityTreeElementPointer = std::static_pointer_cast<EntityTreeElement>(element);
+    float distance = FLT_MAX;
+    // If origin is inside the cube, always check this element first
+    if (entityTreeElementPointer->getAACube().contains(args->origin)) {
+        distance = 0.0f;
+    } else {
+        float boundDistance = FLT_MAX;
+        BoxFace face;
+        glm::vec3 surfaceNormal;
+        if (entityTreeElementPointer->getAACube().findRayIntersection(args->origin, args->direction, args->invDirection, boundDistance, face, surfaceNormal)) {
+            // Don't add this cell if it's already farther than our best distance so far
+            if (boundDistance < args->distance) {
+                distance = boundDistance;
+            }
+        }
+    }
+    return distance;
+}
+
 EntityItemID EntityTree::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction,
                                     QVector<EntityItemID> entityIdsToInclude, QVector<EntityItemID> entityIdsToDiscard,
                                     bool visibleOnly, bool collidableOnly, bool precisionPicking, 
                                     OctreeElementPointer& element, float& distance,
                                     BoxFace& face, glm::vec3& surfaceNormal, QVariantMap& extraInfo,
                                     Octree::lockType lockType, bool* accurateResult) {
-    RayArgs args = { origin, direction, entityIdsToInclude, entityIdsToDiscard,
+    RayArgs args = { origin, direction, 1.0f / direction, entityIdsToInclude, entityIdsToDiscard,
             visibleOnly, collidableOnly, precisionPicking, element, distance, face, surfaceNormal, extraInfo, EntityItemID() };
     distance = FLT_MAX;
 
     bool requireLock = lockType == Octree::Lock;
     bool lockResult = withReadLock([&]{
-        recurseTreeWithOperation(findRayIntersectionOp, &args);
+        recurseTreeWithOperationSorted(findRayIntersectionOp, findRayIntersectionSortingOp, &args);
     }, requireLock);
 
     if (accurateResult) {
@@ -860,15 +884,38 @@ bool findParabolaIntersectionOp(const OctreeElementPointer& element, void* extra
     ParabolaArgs* args = static_cast<ParabolaArgs*>(extraData);
     bool keepSearching = true;
     EntityTreeElementPointer entityTreeElementPointer = std::static_pointer_cast<EntityTreeElement>(element);
-    EntityItemID entityID = entityTreeElementPointer->findParabolaIntersection(args->origin, args->velocity, args->acceleration, keepSearching,
+    EntityItemID entityID = entityTreeElementPointer->findParabolaIntersection(args->origin, args->velocity, args->acceleration,
         args->element, args->parabolicDistance, args->face, args->surfaceNormal, args->entityIdsToInclude,
         args->entityIdsToDiscard, args->visibleOnly, args->collidableOnly, args->extraInfo, args->precisionPicking);
     if (!entityID.isNull()) {
         args->entityID = entityID;
+        // We recurse OctreeElements in order, so if we hit something, we can stop immediately
+        keepSearching = false;
     }
     return keepSearching;
 }
 
+float findParabolaIntersectionSortingOp(const OctreeElementPointer& element, void* extraData) {
+    ParabolaArgs* args = static_cast<ParabolaArgs*>(extraData);
+    EntityTreeElementPointer entityTreeElementPointer = std::static_pointer_cast<EntityTreeElement>(element);
+    float distance = FLT_MAX;
+    // If origin is inside the cube, always check this element first
+    if (entityTreeElementPointer->getAACube().contains(args->origin)) {
+        distance = 0.0f;
+    } else {
+        float boundDistance = FLT_MAX;
+        BoxFace face;
+        glm::vec3 surfaceNormal;
+        if (entityTreeElementPointer->getAACube().findParabolaIntersection(args->origin, args->velocity, args->acceleration, boundDistance, face, surfaceNormal)) {
+            // Don't add this cell if it's already farther than our best distance so far
+            if (boundDistance < args->parabolicDistance) {
+                distance = boundDistance;
+            }
+        }
+    }
+    return distance;
+}
+
 EntityItemID EntityTree::findParabolaIntersection(const PickParabola& parabola,
                                     QVector<EntityItemID> entityIdsToInclude, QVector<EntityItemID> entityIdsToDiscard,
                                     bool visibleOnly, bool collidableOnly, bool precisionPicking,
@@ -882,7 +929,7 @@ EntityItemID EntityTree::findParabolaIntersection(const PickParabola& parabola,
 
     bool requireLock = lockType == Octree::Lock;
     bool lockResult = withReadLock([&] {
-        recurseTreeWithOperation(findParabolaIntersectionOp, &args);
+        recurseTreeWithOperationSorted(findParabolaIntersectionOp, findParabolaIntersectionSortingOp, &args);
     }, requireLock);
 
     if (accurateResult) {
diff --git a/libraries/entities/src/EntityTreeElement.cpp b/libraries/entities/src/EntityTreeElement.cpp
index 5974fce6c5..efe5dafccf 100644
--- a/libraries/entities/src/EntityTreeElement.cpp
+++ b/libraries/entities/src/EntityTreeElement.cpp
@@ -140,31 +140,18 @@ bool EntityTreeElement::bestFitBounds(const glm::vec3& minPoint, const glm::vec3
 }
 
 EntityItemID EntityTreeElement::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction,
-    bool& keepSearching, OctreeElementPointer& element, float& distance,
-    BoxFace& face, glm::vec3& surfaceNormal, const QVector<EntityItemID>& entityIdsToInclude,
-    const QVector<EntityItemID>& entityIdsToDiscard, bool visibleOnly, bool collidableOnly,
-    QVariantMap& extraInfo, bool precisionPicking) {
+    OctreeElementPointer& element, float& distance, BoxFace& face, glm::vec3& surfaceNormal,
+    const QVector<EntityItemID>& entityIdsToInclude, const QVector<EntityItemID>& entityIdsToDiscard,
+    bool visibleOnly, bool collidableOnly, QVariantMap& extraInfo, bool precisionPicking) {
 
     EntityItemID result;
-    float distanceToElementCube = std::numeric_limits<float>::max();
     BoxFace localFace;
     glm::vec3 localSurfaceNormal;
 
-    // if the ray doesn't intersect with our cube OR the distance to element is less than current best distance
-    // we can stop searching!
-    bool hit = _cube.findRayIntersection(origin, direction, distanceToElementCube, localFace, localSurfaceNormal);
-    if (!hit || (!_cube.contains(origin) && distanceToElementCube > distance)) {
-        keepSearching = false; // no point in continuing to search
-        return result; // we did not intersect
-    }
-
-    // by default, we only allow intersections with leaves with content
     if (!canPickIntersect()) {
-        return result; // we don't intersect with non-leaves, and we keep searching
+        return result;
     }
 
-    // if the distance to the element cube is not less than the current best distance, then it's not possible
-    // for any details inside the cube to be closer so we don't need to consider them.
     QVariantMap localExtraInfo;
     float distanceToElementDetails = distance;
     EntityItemID entityID = findDetailedRayIntersection(origin, direction, element, distanceToElementDetails,
@@ -228,7 +215,7 @@ EntityItemID EntityTreeElement::findDetailedRayIntersection(const glm::vec3& ori
         float localDistance;
         BoxFace localFace;
         glm::vec3 localSurfaceNormal;
-        if (entityFrameBox.findRayIntersection(entityFrameOrigin, entityFrameDirection, localDistance,
+        if (entityFrameBox.findRayIntersection(entityFrameOrigin, entityFrameDirection, 1.0f / entityFrameDirection, localDistance,
                                                 localFace, localSurfaceNormal)) {
             if (entityFrameBox.contains(entityFrameOrigin) || localDistance < distance) {
                 // now ask the entity if we actually intersect
@@ -289,31 +276,19 @@ bool EntityTreeElement::findSpherePenetration(const glm::vec3& center, float rad
 }
 
 EntityItemID EntityTreeElement::findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity,
-    const glm::vec3& acceleration, bool& keepSearching, OctreeElementPointer& element, float& parabolicDistance,
+    const glm::vec3& acceleration, OctreeElementPointer& element, float& parabolicDistance,
     BoxFace& face, glm::vec3& surfaceNormal, const QVector<EntityItemID>& entityIdsToInclude,
     const QVector<EntityItemID>& entityIdsToDiscard, bool visibleOnly, bool collidableOnly,
     QVariantMap& extraInfo, bool precisionPicking) {
 
     EntityItemID result;
-    float distanceToElementCube = std::numeric_limits<float>::max();
     BoxFace localFace;
     glm::vec3 localSurfaceNormal;
 
-    // if the parabola doesn't intersect with our cube OR the distance to element is less than current best distance
-    // we can stop searching!
-    bool hit = _cube.findParabolaIntersection(origin, velocity, acceleration, distanceToElementCube, localFace, localSurfaceNormal);
-    if (!hit || (!_cube.contains(origin) && distanceToElementCube > parabolicDistance)) {
-        keepSearching = false; // no point in continuing to search
-        return result; // we did not intersect
-    }
-
-    // by default, we only allow intersections with leaves with content
     if (!canPickIntersect()) {
-        return result; // we don't intersect with non-leaves, and we keep searching
+        return result;
     }
 
-    // if the distance to the element cube is not less than the current best distance, then it's not possible
-    // for any details inside the cube to be closer so we don't need to consider them.
     QVariantMap localExtraInfo;
     float distanceToElementDetails = parabolicDistance;
     // We can precompute the world-space parabola normal and reuse it for the parabola plane intersects AABox sphere check
diff --git a/libraries/entities/src/EntityTreeElement.h b/libraries/entities/src/EntityTreeElement.h
index d6f9db08d6..793340c9a4 100644
--- a/libraries/entities/src/EntityTreeElement.h
+++ b/libraries/entities/src/EntityTreeElement.h
@@ -136,10 +136,9 @@ public:
 
     virtual bool canPickIntersect() const override { return hasEntities(); }
     virtual EntityItemID findRayIntersection(const glm::vec3& origin, const glm::vec3& direction,
-        bool& keepSearching, OctreeElementPointer& element, float& distance,
-        BoxFace& face, glm::vec3& surfaceNormal, const QVector<EntityItemID>& entityIdsToInclude,
-        const QVector<EntityItemID>& entityIdsToDiscard, bool visibleOnly, bool collidableOnly,
-        QVariantMap& extraInfo, bool precisionPicking = false);
+        OctreeElementPointer& element, float& distance, BoxFace& face, glm::vec3& surfaceNormal,
+        const QVector<EntityItemID>& entityIdsToInclude, const QVector<EntityItemID>& entityIdsToDiscard,
+        bool visibleOnly, bool collidableOnly, QVariantMap& extraInfo, bool precisionPicking = false);
     virtual EntityItemID findDetailedRayIntersection(const glm::vec3& origin, const glm::vec3& direction,
                          OctreeElementPointer& element, float& distance,
                          BoxFace& face, glm::vec3& surfaceNormal, const QVector<EntityItemID>& entityIdsToInclude,
@@ -149,7 +148,7 @@ public:
                         glm::vec3& penetration, void** penetratedObject) const override;
 
     virtual EntityItemID findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity,
-        const glm::vec3& acceleration, bool& keepSearching, OctreeElementPointer& element, float& parabolicDistance,
+        const glm::vec3& acceleration, OctreeElementPointer& element, float& parabolicDistance,
         BoxFace& face, glm::vec3& surfaceNormal, const QVector<EntityItemID>& entityIdsToInclude,
         const QVector<EntityItemID>& entityIdsToDiscard, bool visibleOnly, bool collidableOnly,
         QVariantMap& extraInfo, bool precisionPicking = false);
diff --git a/libraries/entities/src/ShapeEntityItem.cpp b/libraries/entities/src/ShapeEntityItem.cpp
index e4ea1470c1..773a7059dc 100644
--- a/libraries/entities/src/ShapeEntityItem.cpp
+++ b/libraries/entities/src/ShapeEntityItem.cpp
@@ -262,20 +262,18 @@ bool ShapeEntityItem::findDetailedRayIntersection(const glm::vec3& origin, const
     glm::mat4 entityToWorldMatrix = getEntityToWorldMatrix();
     glm::mat4 worldToEntityMatrix = glm::inverse(entityToWorldMatrix);
     glm::vec3 entityFrameOrigin = glm::vec3(worldToEntityMatrix * glm::vec4(origin, 1.0f));
-    glm::vec3 entityFrameDirection = glm::normalize(glm::vec3(worldToEntityMatrix * glm::vec4(direction, 0.0f)));
+    glm::vec3 entityFrameDirection = glm::vec3(worldToEntityMatrix * glm::vec4(direction, 0.0f));
 
-    float localDistance;
     // NOTE: unit sphere has center of 0,0,0 and radius of 0.5
-    if (findRaySphereIntersection(entityFrameOrigin, entityFrameDirection, glm::vec3(0.0f), 0.5f, localDistance)) {
-        // determine where on the unit sphere the hit point occured
-        glm::vec3 entityFrameHitAt = entityFrameOrigin + (entityFrameDirection * localDistance);
-        // then translate back to work coordinates
-        glm::vec3 hitAt = glm::vec3(entityToWorldMatrix * glm::vec4(entityFrameHitAt, 1.0f));
-        distance = glm::distance(origin, hitAt);
+    if (findRaySphereIntersection(entityFrameOrigin, entityFrameDirection, glm::vec3(0.0f), 0.5f, distance)) {
         bool success;
-        // FIXME: this is only correct for uniformly scaled spheres
-        surfaceNormal = glm::normalize(hitAt - getCenterPosition(success));
-        if (!success) {
+        glm::vec3 center = getCenterPosition(success);
+        if (success) {
+            // FIXME: this is only correct for uniformly scaled spheres
+            // determine where on the unit sphere the hit point occured
+            glm::vec3 hitAt = origin + (direction * distance);
+            surfaceNormal = glm::normalize(hitAt - center);
+        } else {
             return false;
         }
         return true;
@@ -297,9 +295,11 @@ bool ShapeEntityItem::findDetailedParabolaIntersection(const glm::vec3& origin,
     // NOTE: unit sphere has center of 0,0,0 and radius of 0.5
     if (findParabolaSphereIntersection(entityFrameOrigin, entityFrameVelocity, entityFrameAcceleration, glm::vec3(0.0f), 0.5f, parabolicDistance)) {
         bool success;
-        // FIXME: this is only correct for uniformly scaled spheres
-        surfaceNormal = glm::normalize((origin + velocity * parabolicDistance + 0.5f * acceleration * parabolicDistance * parabolicDistance) - getCenterPosition(success));
-        if (!success) {
+        glm::vec3 center = getCenterPosition(success);
+        if (success) {
+            // FIXME: this is only correct for uniformly scaled spheres
+            surfaceNormal = glm::normalize((origin + velocity * parabolicDistance + 0.5f * acceleration * parabolicDistance * parabolicDistance) - center);
+        } else {
             return false;
         }
         return true;
diff --git a/libraries/networking/src/AccountManager.h b/libraries/networking/src/AccountManager.h
index a79b69fe2b..b122115dd0 100644
--- a/libraries/networking/src/AccountManager.h
+++ b/libraries/networking/src/AccountManager.h
@@ -96,6 +96,8 @@ public:
 
     QUrl getMetaverseServerURL() { return NetworkingConstants::METAVERSE_SERVER_URL(); }
 
+    void removeAccountFromFile();
+
 public slots:
     void requestAccessToken(const QString& login, const QString& password);
     void requestAccessTokenWithSteam(QByteArray authSessionTicket);
@@ -133,7 +135,6 @@ private:
     void operator=(AccountManager const& other) = delete;
 
     void persistAccountToFile();
-    void removeAccountFromFile();
 
     void passSuccessToCallback(QNetworkReply* reply);
     void passErrorToCallback(QNetworkReply* reply);
diff --git a/libraries/octree/src/Octree.cpp b/libraries/octree/src/Octree.cpp
index 9bb0e25982..3d387e0956 100644
--- a/libraries/octree/src/Octree.cpp
+++ b/libraries/octree/src/Octree.cpp
@@ -68,55 +68,12 @@ Octree::~Octree() {
     eraseAllOctreeElements(false);
 }
 
-
-// Inserts the value and key into three arrays sorted by the key array, the first array is the value,
-// the second array is a sorted key for the value, the third array is the index for the value in it original
-// non-sorted array
-// returns -1 if size exceeded
-// originalIndexArray is optional
-int insertOctreeElementIntoSortedArrays(const OctreeElementPointer& value, float key, int originalIndex,
-                                        OctreeElementPointer* valueArray, float* keyArray, int* originalIndexArray,
-                                        int currentCount, int maxCount) {
-
-    if (currentCount < maxCount) {
-        int i = 0;
-        if (currentCount > 0) {
-            while (i < currentCount && key > keyArray[i]) {
-                i++;
-            }
-            // i is our desired location
-            // shift array elements to the right
-            if (i < currentCount && i+1 < maxCount) {
-                for (int j = currentCount - 1; j > i; j--) {
-                    valueArray[j] = valueArray[j - 1];
-                    keyArray[j] = keyArray[j - 1];
-                }
-            }
-        }
-        // place new element at i
-        valueArray[i] = value;
-        keyArray[i] = key;
-        if (originalIndexArray) {
-            originalIndexArray[i] = originalIndex;
-        }
-        return currentCount + 1;
-    }
-    return -1; // error case
-}
-
-
-
 // Recurses voxel tree calling the RecurseOctreeOperation function for each element.
 // stops recursion if operation function returns false.
 void Octree::recurseTreeWithOperation(const RecurseOctreeOperation& operation, void* extraData) {
     recurseElementWithOperation(_rootElement, operation, extraData);
 }
 
-// Recurses voxel tree calling the RecurseOctreePostFixOperation function for each element in post-fix order.
-void Octree::recurseTreeWithPostOperation(const RecurseOctreeOperation& operation, void* extraData) {
-    recurseElementWithPostOperation(_rootElement, operation, extraData);
-}
-
 // Recurses voxel element with an operation function
 void Octree::recurseElementWithOperation(const OctreeElementPointer& element, const RecurseOctreeOperation& operation, void* extraData,
                         int recursionCount) {
@@ -129,71 +86,53 @@ void Octree::recurseElementWithOperation(const OctreeElementPointer& element, co
         for (int i = 0; i < NUMBER_OF_CHILDREN; i++) {
             OctreeElementPointer child = element->getChildAtIndex(i);
             if (child) {
-                recurseElementWithOperation(child, operation, extraData, recursionCount+1);
+                recurseElementWithOperation(child, operation, extraData, recursionCount + 1);
             }
         }
     }
 }
 
-// Recurses voxel element with an operation function
-void Octree::recurseElementWithPostOperation(const OctreeElementPointer& element, const RecurseOctreeOperation& operation,
-                                             void* extraData, int recursionCount) {
+void Octree::recurseTreeWithOperationSorted(const RecurseOctreeOperation& operation, const RecurseOctreeSortingOperation& sortingOperation, void* extraData) {
+    recurseElementWithOperationSorted(_rootElement, operation, sortingOperation, extraData);
+}
+
+// Recurses voxel element with an operation function, calling operation on its children in a specific order
+bool Octree::recurseElementWithOperationSorted(const OctreeElementPointer& element, const RecurseOctreeOperation& operation,
+                                               const RecurseOctreeSortingOperation& sortingOperation, void* extraData, int recursionCount) {
     if (recursionCount > DANGEROUSLY_DEEP_RECURSION) {
-        HIFI_FCDEBUG(octree(), "Octree::recurseElementWithPostOperation() reached DANGEROUSLY_DEEP_RECURSION, bailing!");
-        return;
+        HIFI_FCDEBUG(octree(), "Octree::recurseElementWithOperationSorted() reached DANGEROUSLY_DEEP_RECURSION, bailing!");
+        // If we go too deep, we want to keep searching other paths
+        return true;
     }
 
+    bool keepSearching = operation(element, extraData);
+
+    std::vector<SortedChild> sortedChildren;
     for (int i = 0; i < NUMBER_OF_CHILDREN; i++) {
         OctreeElementPointer child = element->getChildAtIndex(i);
         if (child) {
-            recurseElementWithPostOperation(child, operation, extraData, recursionCount+1);
-        }
-    }
-    operation(element, extraData);
-}
-
-// Recurses voxel tree calling the RecurseOctreeOperation function for each element.
-// stops recursion if operation function returns false.
-void Octree::recurseTreeWithOperationDistanceSorted(const RecurseOctreeOperation& operation,
-                                                       const glm::vec3& point, void* extraData) {
-
-    recurseElementWithOperationDistanceSorted(_rootElement, operation, point, extraData);
-}
-
-// Recurses voxel element with an operation function
-void Octree::recurseElementWithOperationDistanceSorted(const OctreeElementPointer& element, const RecurseOctreeOperation& operation,
-                                                       const glm::vec3& point, void* extraData, int recursionCount) {
-
-    if (recursionCount > DANGEROUSLY_DEEP_RECURSION) {
-        HIFI_FCDEBUG(octree(), "Octree::recurseElementWithOperationDistanceSorted() reached DANGEROUSLY_DEEP_RECURSION, bailing!");
-        return;
-    }
-
-    if (operation(element, extraData)) {
-        // determine the distance sorted order of our children
-        OctreeElementPointer sortedChildren[NUMBER_OF_CHILDREN] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };
-        float distancesToChildren[NUMBER_OF_CHILDREN] = { 0, 0, 0, 0, 0, 0, 0, 0 };
-        int indexOfChildren[NUMBER_OF_CHILDREN] = { 0, 0, 0, 0, 0, 0, 0, 0 };
-        int currentCount = 0;
-
-        for (int i = 0; i < NUMBER_OF_CHILDREN; i++) {
-            OctreeElementPointer childElement = element->getChildAtIndex(i);
-            if (childElement) {
-                // chance to optimize, doesn't need to be actual distance!! Could be distance squared
-                float distanceSquared = childElement->distanceSquareToPoint(point);
-                currentCount = insertOctreeElementIntoSortedArrays(childElement, distanceSquared, i,
-                                                                   sortedChildren, (float*)&distancesToChildren,
-                                                                   (int*)&indexOfChildren, currentCount, NUMBER_OF_CHILDREN);
-            }
-        }
-
-        for (int i = 0; i < currentCount; i++) {
-            OctreeElementPointer childElement = sortedChildren[i];
-            if (childElement) {
-                recurseElementWithOperationDistanceSorted(childElement, operation, point, extraData);
+            float priority = sortingOperation(child, extraData);
+            if (priority < FLT_MAX) {
+                sortedChildren.emplace_back(priority, child);
             }
         }
     }
+
+    if (sortedChildren.size() > 1) {
+        static auto comparator = [](const SortedChild& left, const SortedChild& right) { return left.first < right.first; };
+        std::sort(sortedChildren.begin(), sortedChildren.end(), comparator);
+    }
+
+    for (auto it = sortedChildren.begin(); it != sortedChildren.end(); ++it) {
+        const SortedChild& sortedChild = *it;
+        // Our children were sorted, so if one hits something, we don't need to check the others
+        if (!recurseElementWithOperationSorted(sortedChild.second, operation, sortingOperation, extraData, recursionCount + 1)) {
+            return false;
+        }
+    }
+    // We checked all our children and didn't find anything.
+    // Stop if we hit something in this element.  Continue if we didn't.
+    return keepSearching;
 }
 
 void Octree::recurseTreeWithOperator(RecurseOctreeOperator* operatorObject) {
diff --git a/libraries/octree/src/Octree.h b/libraries/octree/src/Octree.h
index b827ed7cd0..a2b2f227cb 100644
--- a/libraries/octree/src/Octree.h
+++ b/libraries/octree/src/Octree.h
@@ -49,6 +49,9 @@ public:
 
 // Callback function, for recuseTreeWithOperation
 using RecurseOctreeOperation = std::function<bool(const OctreeElementPointer&, void*)>;
+// Function for sorting octree children during recursion.  If return value == FLT_MAX, child is discarded
+using RecurseOctreeSortingOperation = std::function<float(const OctreeElementPointer&, void*)>;
+using SortedChild = std::pair<float, OctreeElementPointer>;
 typedef QHash<uint, AACube> CubeList;
 
 const bool NO_EXISTS_BITS         = false;
@@ -163,17 +166,10 @@ public:
     OctreeElementPointer getOrCreateChildElementContaining(const AACube& box);
 
     void recurseTreeWithOperation(const RecurseOctreeOperation& operation, void* extraData = NULL);
-    void recurseTreeWithPostOperation(const RecurseOctreeOperation& operation, void* extraData = NULL);
-
-    /// \param operation type of operation
-    /// \param point point in world-frame (meters)
-    /// \param extraData hook for user data to be interpreted by special context
-    void recurseTreeWithOperationDistanceSorted(const RecurseOctreeOperation& operation,
-                                                const glm::vec3& point, void* extraData = NULL);
+    void recurseTreeWithOperationSorted(const RecurseOctreeOperation& operation, const RecurseOctreeSortingOperation& sortingOperation, void* extraData = NULL);
 
     void recurseTreeWithOperator(RecurseOctreeOperator* operatorObject);
 
-
     bool isDirty() const { return _isDirty; }
     void clearDirtyBit() { _isDirty = false; }
     void setDirtyBit() { _isDirty = true; }
@@ -227,14 +223,8 @@ public:
 
     void recurseElementWithOperation(const OctreeElementPointer& element, const RecurseOctreeOperation& operation,
                 void* extraData, int recursionCount = 0);
-
-    /// Traverse child nodes of node applying operation in post-fix order
-    ///
-    void recurseElementWithPostOperation(const OctreeElementPointer& element, const RecurseOctreeOperation& operation,
-                void* extraData, int recursionCount = 0);
-
-    void recurseElementWithOperationDistanceSorted(const OctreeElementPointer& element, const RecurseOctreeOperation& operation,
-                const glm::vec3& point, void* extraData, int recursionCount = 0);
+    bool recurseElementWithOperationSorted(const OctreeElementPointer& element, const RecurseOctreeOperation& operation,
+        const RecurseOctreeSortingOperation& sortingOperation, void* extraData, int recursionCount = 0);
 
     bool recurseElementWithOperator(const OctreeElementPointer& element, RecurseOctreeOperator* operatorObject, int recursionCount = 0);
 
diff --git a/libraries/pointers/src/Pick.cpp b/libraries/pointers/src/Pick.cpp
index 4315409bdb..7ac53a2643 100644
--- a/libraries/pointers/src/Pick.cpp
+++ b/libraries/pointers/src/Pick.cpp
@@ -72,11 +72,15 @@ PickResultPointer PickQuery::getPrevPickResult() const {
 void PickQuery::setIgnoreItems(const QVector<QUuid>& ignoreItems) {
     withWriteLock([&] {
         _ignoreItems = ignoreItems;
+        // We sort these items here so the PickCacheOptimizer can catch cases where two picks have the same ignoreItems in a different order
+        std::sort(_ignoreItems.begin(), _ignoreItems.end(), std::less<QUuid>());
     });
 }
 
 void PickQuery::setIncludeItems(const QVector<QUuid>& includeItems) {
     withWriteLock([&] {
         _includeItems = includeItems;
+        // We sort these items here so the PickCacheOptimizer can catch cases where two picks have the same includeItems in a different order
+        std::sort(_includeItems.begin(), _includeItems.end(), std::less<QUuid>());
     });
 }
\ No newline at end of file
diff --git a/libraries/pointers/src/PickCacheOptimizer.h b/libraries/pointers/src/PickCacheOptimizer.h
index 49a039935c..7a52cfc410 100644
--- a/libraries/pointers/src/PickCacheOptimizer.h
+++ b/libraries/pointers/src/PickCacheOptimizer.h
@@ -37,7 +37,7 @@ template<typename T>
 class PickCacheOptimizer {
 
 public:
-    void update(std::unordered_map<uint32_t, std::shared_ptr<PickQuery>>& picks, uint32_t& nextToUpdate, uint64_t expiry, bool shouldPickHUD);
+    QVector4D update(std::unordered_map<uint32_t, std::shared_ptr<PickQuery>>& picks, uint32_t& nextToUpdate, uint64_t expiry, bool shouldPickHUD);
 
 protected:
     typedef std::unordered_map<T, std::unordered_map<PickCacheKey, PickResultPointer>> PickCache;
@@ -67,8 +67,9 @@ void PickCacheOptimizer<T>::cacheResult(const bool intersects, const PickResultP
 }
 
 template<typename T>
-void PickCacheOptimizer<T>::update(std::unordered_map<uint32_t, std::shared_ptr<PickQuery>>& picks,
+QVector4D PickCacheOptimizer<T>::update(std::unordered_map<uint32_t, std::shared_ptr<PickQuery>>& picks,
         uint32_t& nextToUpdate, uint64_t expiry, bool shouldPickHUD) {
+    QVector4D numIntersectionsComputed;
     PickCache results;
     const uint32_t INVALID_PICK_ID = 0;
     auto itr = picks.begin();
@@ -91,6 +92,7 @@ void PickCacheOptimizer<T>::update(std::unordered_map<uint32_t, std::shared_ptr<
                 PickCacheKey entityKey = { pick->getFilter().getEntityFlags(), pick->getIncludeItems(), pick->getIgnoreItems() };
                 if (!checkAndCompareCachedResults(mathematicalPick, results, res, entityKey)) {
                     PickResultPointer entityRes = pick->getEntityIntersection(mathematicalPick);
+                    numIntersectionsComputed[0]++;
                     if (entityRes) {
                         cacheResult(entityRes->doesIntersect(), entityRes, entityKey, res, mathematicalPick, results, pick);
                     }
@@ -101,6 +103,7 @@ void PickCacheOptimizer<T>::update(std::unordered_map<uint32_t, std::shared_ptr<
                 PickCacheKey overlayKey = { pick->getFilter().getOverlayFlags(), pick->getIncludeItems(), pick->getIgnoreItems() };
                 if (!checkAndCompareCachedResults(mathematicalPick, results, res, overlayKey)) {
                     PickResultPointer overlayRes = pick->getOverlayIntersection(mathematicalPick);
+                    numIntersectionsComputed[1]++;
                     if (overlayRes) {
                         cacheResult(overlayRes->doesIntersect(), overlayRes, overlayKey, res, mathematicalPick, results, pick);
                     }
@@ -111,6 +114,7 @@ void PickCacheOptimizer<T>::update(std::unordered_map<uint32_t, std::shared_ptr<
                 PickCacheKey avatarKey = { pick->getFilter().getAvatarFlags(), pick->getIncludeItems(), pick->getIgnoreItems() };
                 if (!checkAndCompareCachedResults(mathematicalPick, results, res, avatarKey)) {
                     PickResultPointer avatarRes = pick->getAvatarIntersection(mathematicalPick);
+                    numIntersectionsComputed[2]++;
                     if (avatarRes) {
                         cacheResult(avatarRes->doesIntersect(), avatarRes, avatarKey, res, mathematicalPick, results, pick);
                     }
@@ -122,6 +126,7 @@ void PickCacheOptimizer<T>::update(std::unordered_map<uint32_t, std::shared_ptr<
                 PickCacheKey hudKey = { pick->getFilter().getHUDFlags(), QVector<QUuid>(), QVector<QUuid>() };
                 if (!checkAndCompareCachedResults(mathematicalPick, results, res, hudKey)) {
                     PickResultPointer hudRes = pick->getHUDIntersection(mathematicalPick);
+                    numIntersectionsComputed[3]++;
                     if (hudRes) {
                         cacheResult(true, hudRes, hudKey, res, mathematicalPick, results, pick);
                     }
@@ -145,6 +150,7 @@ void PickCacheOptimizer<T>::update(std::unordered_map<uint32_t, std::shared_ptr<
             break;
         }
     }
+    return numIntersectionsComputed;
 }
 
 #endif // hifi_PickCacheOptimizer_h
diff --git a/libraries/pointers/src/PickManager.cpp b/libraries/pointers/src/PickManager.cpp
index 470d88a46c..a384eb661f 100644
--- a/libraries/pointers/src/PickManager.cpp
+++ b/libraries/pointers/src/PickManager.cpp
@@ -20,6 +20,7 @@ unsigned int PickManager::addPick(PickQuery::PickType type, const std::shared_pt
             id = _nextPickID++;
             _picks[type][id] = pick;
             _typeMap[id] = type;
+            _totalPickCounts[type]++;
         }
     });
     return id;
@@ -41,6 +42,7 @@ void PickManager::removePick(unsigned int uid) {
         if (type != _typeMap.end()) {
             _picks[type->second].erase(uid);
             _typeMap.erase(uid);
+            _totalPickCounts[type->second]--;
         }
     });
 }
@@ -96,12 +98,12 @@ void PickManager::update() {
     });
 
     bool shouldPickHUD = _shouldPickHUDOperator();
-    // we pass the same expiry to both updates, but the stylus updates are relatively cheap
-    // and the rayPicks updae will ALWAYS update at least one ray even when there is no budget
-    _stylusPickCacheOptimizer.update(cachedPicks[PickQuery::Stylus], _nextPickToUpdate[PickQuery::Stylus], expiry, false);
-    _rayPickCacheOptimizer.update(cachedPicks[PickQuery::Ray], _nextPickToUpdate[PickQuery::Ray], expiry, shouldPickHUD);
-    _parabolaPickCacheOptimizer.update(cachedPicks[PickQuery::Parabola], _nextPickToUpdate[PickQuery::Parabola], expiry, shouldPickHUD);
-    _collisionPickCacheOptimizer.update(cachedPicks[PickQuery::Collision], _nextPickToUpdate[PickQuery::Collision], expiry, false);
+    // FIXME: give each type its own expiry
+    // Each type will update at least one pick, regardless of the expiry
+    _updatedPickCounts[PickQuery::Stylus] = _stylusPickCacheOptimizer.update(cachedPicks[PickQuery::Stylus], _nextPickToUpdate[PickQuery::Stylus], expiry, false);
+    _updatedPickCounts[PickQuery::Ray] = _rayPickCacheOptimizer.update(cachedPicks[PickQuery::Ray], _nextPickToUpdate[PickQuery::Ray], expiry, shouldPickHUD);
+    _updatedPickCounts[PickQuery::Parabola] = _parabolaPickCacheOptimizer.update(cachedPicks[PickQuery::Parabola], _nextPickToUpdate[PickQuery::Parabola], expiry, shouldPickHUD);
+    _updatedPickCounts[PickQuery::Collision] = _collisionPickCacheOptimizer.update(cachedPicks[PickQuery::Collision], _nextPickToUpdate[PickQuery::Collision], expiry, false);
 }
 
 bool PickManager::isLeftHand(unsigned int uid) {
diff --git a/libraries/pointers/src/PickManager.h b/libraries/pointers/src/PickManager.h
index 242550d837..e8d32771e0 100644
--- a/libraries/pointers/src/PickManager.h
+++ b/libraries/pointers/src/PickManager.h
@@ -16,7 +16,10 @@
 
 #include <NumericalConstants.h>
 
-class PickManager : public Dependency, protected ReadWriteLockable {
+#include <QObject>
+
+class PickManager : public QObject, public Dependency, protected ReadWriteLockable {
+    Q_OBJECT
     SINGLETON_DEPENDENCY
 
 public:
@@ -53,7 +56,19 @@ public:
     unsigned int getPerFrameTimeBudget() const { return _perFrameTimeBudget; }
     void setPerFrameTimeBudget(unsigned int numUsecs) { _perFrameTimeBudget = numUsecs; }
 
+    bool getForceCoarsePicking() { return _forceCoarsePicking; }
+
+    const std::vector<QVector4D>& getUpdatedPickCounts() { return _updatedPickCounts; }
+    const std::vector<int>& getTotalPickCounts() { return _totalPickCounts; }
+
+public slots:
+    void setForceCoarsePicking(bool forceCoarsePicking) { _forceCoarsePicking = forceCoarsePicking; }
+
 protected:
+    std::vector<QVector4D> _updatedPickCounts { PickQuery::NUM_PICK_TYPES };
+    std::vector<int> _totalPickCounts { 0, 0, 0, 0 };
+
+    bool _forceCoarsePicking { false };
     std::function<bool()> _shouldPickHUDOperator;
     std::function<glm::vec2(const glm::vec3&)> _calculatePos2DFromHUDOperator;
 
diff --git a/libraries/render-utils/src/Model.cpp b/libraries/render-utils/src/Model.cpp
index 71250c6483..d19328aa8a 100644
--- a/libraries/render-utils/src/Model.cpp
+++ b/libraries/render-utils/src/Model.cpp
@@ -388,17 +388,20 @@ bool Model::findRayIntersectionAgainstSubMeshes(const glm::vec3& origin, const g
 
     // we can use the AABox's intersection by mapping our origin and direction into the model frame
     // and testing intersection there.
-    if (modelFrameBox.findRayIntersection(modelFrameOrigin, modelFrameDirection, distance, face, surfaceNormal)) {
+    if (modelFrameBox.findRayIntersection(modelFrameOrigin, modelFrameDirection, 1.0f / modelFrameDirection, distance, face, surfaceNormal)) {
         QMutexLocker locker(&_mutex);
 
-        float bestDistance = std::numeric_limits<float>::max();
+        float bestDistance = FLT_MAX;
+        BoxFace bestFace;
         Triangle bestModelTriangle;
         Triangle bestWorldTriangle;
+        glm::vec3 bestWorldIntersectionPoint;
+        glm::vec3 bestMeshIntersectionPoint;
+        int bestPartIndex = 0;
+        int bestShapeID = 0;
         int bestSubMeshIndex = 0;
 
-        int subMeshIndex = 0;
         const FBXGeometry& geometry = getFBXGeometry();
-
         if (!_triangleSetsValid) {
             calculateTriangleSets(geometry);
         }
@@ -409,41 +412,78 @@ bool Model::findRayIntersectionAgainstSubMeshes(const glm::vec3& origin, const g
 
         glm::vec3 meshFrameOrigin = glm::vec3(worldToMeshMatrix * glm::vec4(origin, 1.0f));
         glm::vec3 meshFrameDirection = glm::vec3(worldToMeshMatrix * glm::vec4(direction, 0.0f));
+        glm::vec3 meshFrameInvDirection = 1.0f / meshFrameDirection;
 
         int shapeID = 0;
+        int subMeshIndex = 0;
+
+        std::vector<SortedTriangleSet> sortedTriangleSets;
         for (auto& meshTriangleSets : _modelSpaceMeshTriangleSets) {
             int partIndex = 0;
-            for (auto &partTriangleSet : meshTriangleSets) {
-                float triangleSetDistance;
-                BoxFace triangleSetFace;
-                Triangle triangleSetTriangle;
-                if (partTriangleSet.findRayIntersection(meshFrameOrigin, meshFrameDirection, triangleSetDistance, triangleSetFace, triangleSetTriangle, pickAgainstTriangles, allowBackface)) {
-                    glm::vec3 meshIntersectionPoint = meshFrameOrigin + (meshFrameDirection * triangleSetDistance);
-                    glm::vec3 worldIntersectionPoint = glm::vec3(meshToWorldMatrix * glm::vec4(meshIntersectionPoint, 1.0f));
-                    float worldDistance = glm::distance(origin, worldIntersectionPoint);
-
-                    if (worldDistance < bestDistance) {
-                        bestDistance = worldDistance;
-                        intersectedSomething = true;
-                        face = triangleSetFace;
-                        bestModelTriangle = triangleSetTriangle;
-                        bestWorldTriangle = triangleSetTriangle * meshToWorldMatrix;
-                        extraInfo["worldIntersectionPoint"] = vec3toVariant(worldIntersectionPoint);
-                        extraInfo["meshIntersectionPoint"] = vec3toVariant(meshIntersectionPoint);
-                        extraInfo["partIndex"] = partIndex;
-                        extraInfo["shapeID"] = shapeID;
-                        bestSubMeshIndex = subMeshIndex;
+            for (auto& partTriangleSet : meshTriangleSets) {
+                float priority = FLT_MAX;
+                if (partTriangleSet.getBounds().contains(meshFrameOrigin)) {
+                    priority = 0.0f;
+                } else {
+                    float partBoundDistance = FLT_MAX;
+                    BoxFace partBoundFace;
+                    glm::vec3 partBoundNormal;
+                    if (partTriangleSet.getBounds().findRayIntersection(meshFrameOrigin, meshFrameDirection, meshFrameInvDirection,
+                                                                        partBoundDistance, partBoundFace, partBoundNormal)) {
+                        priority = partBoundDistance;
                     }
                 }
+
+                if (priority < FLT_MAX) {
+                    sortedTriangleSets.emplace_back(priority, &partTriangleSet, partIndex, shapeID, subMeshIndex);
+                }
                 partIndex++;
                 shapeID++;
             }
             subMeshIndex++;
         }
 
+        if (sortedTriangleSets.size() > 1) {
+            static auto comparator = [](const SortedTriangleSet& left, const SortedTriangleSet& right) { return left.distance < right.distance; };
+            std::sort(sortedTriangleSets.begin(), sortedTriangleSets.end(), comparator);
+        }
+
+        for (auto it = sortedTriangleSets.begin(); it != sortedTriangleSets.end(); ++it) {
+            const SortedTriangleSet& sortedTriangleSet = *it;
+            // We can exit once triangleSetDistance > bestDistance
+            if (sortedTriangleSet.distance > bestDistance) {
+                break;
+            }
+            float triangleSetDistance = FLT_MAX;
+            BoxFace triangleSetFace;
+            Triangle triangleSetTriangle;
+            if (sortedTriangleSet.triangleSet->findRayIntersection(meshFrameOrigin, meshFrameDirection, meshFrameInvDirection, triangleSetDistance, triangleSetFace,
+                                                                   triangleSetTriangle, pickAgainstTriangles, allowBackface)) {
+                if (triangleSetDistance < bestDistance) {
+                    bestDistance = triangleSetDistance;
+                    intersectedSomething = true;
+                    bestFace = triangleSetFace;
+                    bestModelTriangle = triangleSetTriangle;
+                    bestWorldTriangle = triangleSetTriangle * meshToWorldMatrix;
+                    glm::vec3 meshIntersectionPoint = meshFrameOrigin + (meshFrameDirection * triangleSetDistance);
+                    glm::vec3 worldIntersectionPoint = glm::vec3(meshToWorldMatrix * glm::vec4(meshIntersectionPoint, 1.0f));
+                    bestWorldIntersectionPoint = worldIntersectionPoint;
+                    bestMeshIntersectionPoint = meshIntersectionPoint;
+                    bestPartIndex = sortedTriangleSet.partIndex;
+                    bestShapeID = sortedTriangleSet.shapeID;
+                    bestSubMeshIndex = sortedTriangleSet.subMeshIndex;
+                }
+            }
+        }
+
         if (intersectedSomething) {
             distance = bestDistance;
+            face = bestFace;
             surfaceNormal = bestWorldTriangle.getNormal();
+            extraInfo["worldIntersectionPoint"] = vec3toVariant(bestWorldIntersectionPoint);
+            extraInfo["meshIntersectionPoint"] = vec3toVariant(bestMeshIntersectionPoint);
+            extraInfo["partIndex"] = bestPartIndex;
+            extraInfo["shapeID"] = bestShapeID;
             if (pickAgainstTriangles) {
                 extraInfo["subMeshIndex"] = bestSubMeshIndex;
                 extraInfo["subMeshName"] = geometry.getModelNameOfMesh(bestSubMeshIndex);
@@ -495,13 +535,16 @@ bool Model::findParabolaIntersectionAgainstSubMeshes(const glm::vec3& origin, co
         QMutexLocker locker(&_mutex);
 
         float bestDistance = FLT_MAX;
+        BoxFace bestFace;
         Triangle bestModelTriangle;
         Triangle bestWorldTriangle;
+        glm::vec3 bestWorldIntersectionPoint;
+        glm::vec3 bestMeshIntersectionPoint;
+        int bestPartIndex = 0;
+        int bestShapeID = 0;
         int bestSubMeshIndex = 0;
 
-        int subMeshIndex = 0;
         const FBXGeometry& geometry = getFBXGeometry();
-
         if (!_triangleSetsValid) {
             calculateTriangleSets(geometry);
         }
@@ -515,40 +558,79 @@ bool Model::findParabolaIntersectionAgainstSubMeshes(const glm::vec3& origin, co
         glm::vec3 meshFrameAcceleration = glm::vec3(worldToMeshMatrix * glm::vec4(acceleration, 0.0f));
 
         int shapeID = 0;
+        int subMeshIndex = 0;
+
+        std::vector<SortedTriangleSet> sortedTriangleSets;
         for (auto& meshTriangleSets : _modelSpaceMeshTriangleSets) {
             int partIndex = 0;
-            for (auto &partTriangleSet : meshTriangleSets) {
-                float triangleSetDistance;
-                BoxFace triangleSetFace;
-                Triangle triangleSetTriangle;
-                if (partTriangleSet.findParabolaIntersection(meshFrameOrigin, meshFrameVelocity, meshFrameAcceleration,
-                    triangleSetDistance, triangleSetFace, triangleSetTriangle, pickAgainstTriangles, allowBackface)) {
-                    if (triangleSetDistance < bestDistance) {
-                        bestDistance = triangleSetDistance;
-                        intersectedSomething = true;
-                        face = triangleSetFace;
-                        bestModelTriangle = triangleSetTriangle;
-                        bestWorldTriangle = triangleSetTriangle * meshToWorldMatrix;
-                        glm::vec3 meshIntersectionPoint = meshFrameOrigin + meshFrameVelocity * triangleSetDistance +
-                            0.5f * meshFrameAcceleration * triangleSetDistance * triangleSetDistance;
-                        glm::vec3 worldIntersectionPoint = origin + velocity * triangleSetDistance +
-                            0.5f * acceleration * triangleSetDistance * triangleSetDistance;
-                        extraInfo["worldIntersectionPoint"] = vec3toVariant(worldIntersectionPoint);
-                        extraInfo["meshIntersectionPoint"] = vec3toVariant(meshIntersectionPoint);
-                        extraInfo["partIndex"] = partIndex;
-                        extraInfo["shapeID"] = shapeID;
-                        bestSubMeshIndex = subMeshIndex;
+            for (auto& partTriangleSet : meshTriangleSets) {
+                float priority = FLT_MAX;
+                if (partTriangleSet.getBounds().contains(meshFrameOrigin)) {
+                    priority = 0.0f;
+                } else {
+                    float partBoundDistance = FLT_MAX;
+                    BoxFace partBoundFace;
+                    glm::vec3 partBoundNormal;
+                    if (partTriangleSet.getBounds().findParabolaIntersection(meshFrameOrigin, meshFrameVelocity, meshFrameAcceleration,
+                                                                             partBoundDistance, partBoundFace, partBoundNormal)) {
+                        priority = partBoundDistance;
                     }
                 }
+
+                if (priority < FLT_MAX) {
+                    sortedTriangleSets.emplace_back(priority, &partTriangleSet, partIndex, shapeID, subMeshIndex);
+                }
                 partIndex++;
                 shapeID++;
             }
             subMeshIndex++;
         }
 
+        if (sortedTriangleSets.size() > 1) {
+            static auto comparator = [](const SortedTriangleSet& left, const SortedTriangleSet& right) { return left.distance < right.distance; };
+            std::sort(sortedTriangleSets.begin(), sortedTriangleSets.end(), comparator);
+        }
+
+        for (auto it = sortedTriangleSets.begin(); it != sortedTriangleSets.end(); ++it) {
+            const SortedTriangleSet& sortedTriangleSet = *it;
+            // We can exit once triangleSetDistance > bestDistance
+            if (sortedTriangleSet.distance > bestDistance) {
+                break;
+            }
+            float triangleSetDistance = FLT_MAX;
+            BoxFace triangleSetFace;
+            Triangle triangleSetTriangle;
+            if (sortedTriangleSet.triangleSet->findParabolaIntersection(meshFrameOrigin, meshFrameVelocity, meshFrameAcceleration,
+                                                                        triangleSetDistance, triangleSetFace, triangleSetTriangle,
+                                                                        pickAgainstTriangles, allowBackface)) {
+                if (triangleSetDistance < bestDistance) {
+                    bestDistance = triangleSetDistance;
+                    intersectedSomething = true;
+                    bestFace = triangleSetFace;
+                    bestModelTriangle = triangleSetTriangle;
+                    bestWorldTriangle = triangleSetTriangle * meshToWorldMatrix;
+                    glm::vec3 meshIntersectionPoint = meshFrameOrigin + meshFrameVelocity * triangleSetDistance +
+                        0.5f * meshFrameAcceleration * triangleSetDistance * triangleSetDistance;
+                    glm::vec3 worldIntersectionPoint = origin + velocity * triangleSetDistance +
+                        0.5f * acceleration * triangleSetDistance * triangleSetDistance;
+                    bestWorldIntersectionPoint = worldIntersectionPoint;
+                    bestMeshIntersectionPoint = meshIntersectionPoint;
+                    bestPartIndex = sortedTriangleSet.partIndex;
+                    bestShapeID = sortedTriangleSet.shapeID;
+                    bestSubMeshIndex = sortedTriangleSet.subMeshIndex;
+                    // These sets can overlap, so we can't exit early if we find something
+                }
+            }
+        }
+
         if (intersectedSomething) {
             parabolicDistance = bestDistance;
+            face = bestFace;
             surfaceNormal = bestWorldTriangle.getNormal();
+            extraInfo["worldIntersectionPoint"] = vec3toVariant(bestWorldIntersectionPoint);
+            extraInfo["meshIntersectionPoint"] = vec3toVariant(bestMeshIntersectionPoint);
+            extraInfo["partIndex"] = bestPartIndex;
+            extraInfo["shapeID"] = bestShapeID;
             if (pickAgainstTriangles) {
                 extraInfo["subMeshIndex"] = bestSubMeshIndex;
                 extraInfo["subMeshName"] = geometry.getModelNameOfMesh(bestSubMeshIndex);
diff --git a/libraries/render-utils/src/Model.h b/libraries/render-utils/src/Model.h
index 677be61261..9a7f9f3848 100644
--- a/libraries/render-utils/src/Model.h
+++ b/libraries/render-utils/src/Model.h
@@ -64,6 +64,16 @@ class Model;
 using ModelPointer = std::shared_ptr<Model>;
 using ModelWeakPointer = std::weak_ptr<Model>;
 
+struct SortedTriangleSet {
+    SortedTriangleSet(float distance, TriangleSet* triangleSet, int partIndex, int shapeID, int subMeshIndex) :
+        distance(distance), triangleSet(triangleSet), partIndex(partIndex), shapeID(shapeID), subMeshIndex(subMeshIndex) {}
+
+    float distance;
+    TriangleSet* triangleSet;
+    int partIndex;
+    int shapeID;
+    int subMeshIndex;
+};
 
 /// A generic 3D model displaying geometry loaded from a URL.
 class Model : public QObject, public std::enable_shared_from_this<Model>, public scriptable::ModelProvider {
diff --git a/libraries/render-utils/src/PickItemsJob.cpp b/libraries/render-utils/src/PickItemsJob.cpp
index 860262a969..f41222d358 100644
--- a/libraries/render-utils/src/PickItemsJob.cpp
+++ b/libraries/render-utils/src/PickItemsJob.cpp
@@ -33,6 +33,7 @@ void PickItemsJob::run(const render::RenderContextPointer& renderContext, const
 render::ItemBound PickItemsJob::findNearestItem(const render::RenderContextPointer& renderContext, const render::ItemBounds& inputs, float& minIsectDistance) const {
     const glm::vec3 rayOrigin = renderContext->args->getViewFrustum().getPosition();
     const glm::vec3 rayDirection = renderContext->args->getViewFrustum().getDirection();
+    const glm::vec3 rayInvDirection = 1.0f / rayDirection;
     BoxFace face;
     glm::vec3 normal;
     float isectDistance;
@@ -42,7 +43,7 @@ render::ItemBound PickItemsJob::findNearestItem(const render::RenderContextPoint
     render::ItemKey itemKey;
 
     for (const auto& itemBound : inputs) {
-        if (!itemBound.bound.contains(rayOrigin) && itemBound.bound.findRayIntersection(rayOrigin, rayDirection, isectDistance, face, normal)) {
+        if (!itemBound.bound.contains(rayOrigin) && itemBound.bound.findRayIntersection(rayOrigin, rayDirection, rayInvDirection, isectDistance, face, normal)) {
             auto& item = renderContext->_scene->getItem(itemBound.id);
             itemKey = item.getKey();
             if (itemKey.isWorldSpace() && isectDistance>minDistance && isectDistance < minIsectDistance && isectDistance<maxDistance
diff --git a/libraries/shared/src/AABox.cpp b/libraries/shared/src/AABox.cpp
index 994df551fe..b4384c494f 100644
--- a/libraries/shared/src/AABox.cpp
+++ b/libraries/shared/src/AABox.cpp
@@ -192,9 +192,9 @@ bool AABox::expandedIntersectsSegment(const glm::vec3& start, const glm::vec3& e
                 isWithin(start.x + axisDistance*direction.x, expandedCorner.x, expandedSize.x));
 }
 
-bool AABox::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, float& distance,
-                                BoxFace& face, glm::vec3& surfaceNormal) const {
-    return findRayAABoxIntersection(origin, direction, _corner, _scale, distance, face, surfaceNormal);
+bool AABox::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
+                                float& distance, BoxFace& face, glm::vec3& surfaceNormal) const {
+    return findRayAABoxIntersection(origin, direction, invDirection, _corner, _scale, distance, face, surfaceNormal);
 }
 
 bool AABox::findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
diff --git a/libraries/shared/src/AABox.h b/libraries/shared/src/AABox.h
index fbc90cff47..daad01d7c7 100644
--- a/libraries/shared/src/AABox.h
+++ b/libraries/shared/src/AABox.h
@@ -69,7 +69,7 @@ public:
 
     bool expandedContains(const glm::vec3& point, float expansion) const;
     bool expandedIntersectsSegment(const glm::vec3& start, const glm::vec3& end, float expansion) const;
-    bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, float& distance,
+    bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection, float& distance,
                              BoxFace& face, glm::vec3& surfaceNormal) const;
     bool findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
                                   float& parabolicDistance, BoxFace& face, glm::vec3& surfaceNormal) const;
diff --git a/libraries/shared/src/AACube.cpp b/libraries/shared/src/AACube.cpp
index dc1003215d..af2e531917 100644
--- a/libraries/shared/src/AACube.cpp
+++ b/libraries/shared/src/AACube.cpp
@@ -187,9 +187,9 @@ bool AACube::expandedIntersectsSegment(const glm::vec3& start, const glm::vec3&
                 isWithin(start.x + axisDistance*direction.x, expandedCorner.x, expandedSize.x));
 }
 
-bool AACube::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, float& distance, 
-                                    BoxFace& face, glm::vec3& surfaceNormal) const {
-    return findRayAABoxIntersection(origin, direction, _corner, glm::vec3(_scale), distance, face, surfaceNormal);
+bool AACube::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
+                                 float& distance, BoxFace& face, glm::vec3& surfaceNormal) const {
+    return findRayAABoxIntersection(origin, direction, invDirection, _corner, glm::vec3(_scale), distance, face, surfaceNormal);
 }
 
 bool AACube::findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
diff --git a/libraries/shared/src/AACube.h b/libraries/shared/src/AACube.h
index 72aed31999..66b29e3185 100644
--- a/libraries/shared/src/AACube.h
+++ b/libraries/shared/src/AACube.h
@@ -56,10 +56,10 @@ public:
     bool touches(const AABox& otherBox) const;
     bool expandedContains(const glm::vec3& point, float expansion) const;
     bool expandedIntersectsSegment(const glm::vec3& start, const glm::vec3& end, float expansion) const;
-    bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, float& distance,
-                                BoxFace& face, glm::vec3& surfaceNormal) const;
+    bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
+                             float& distance, BoxFace& face, glm::vec3& surfaceNormal) const;
     bool findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
-                                float& parabolicDistance, BoxFace& face, glm::vec3& surfaceNormal) const;
+                                  float& parabolicDistance, BoxFace& face, glm::vec3& surfaceNormal) const;
     bool touchesSphere(const glm::vec3& center, float radius) const;
     bool findSpherePenetration(const glm::vec3& center, float radius, glm::vec3& penetration) const;
     bool findCapsulePenetration(const glm::vec3& start, const glm::vec3& end, float radius, glm::vec3& penetration) const;
diff --git a/libraries/shared/src/GeometryUtil.cpp b/libraries/shared/src/GeometryUtil.cpp
index 6fb06eb624..3f2f7cd7fb 100644
--- a/libraries/shared/src/GeometryUtil.cpp
+++ b/libraries/shared/src/GeometryUtil.cpp
@@ -214,65 +214,39 @@ bool findInsideOutIntersection(float origin, float direction, float corner, floa
     return false;
 }
 
-bool findRayAABoxIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& corner, const glm::vec3& scale, float& distance,
-                              BoxFace& face, glm::vec3& surfaceNormal) {
-    // handle the trivial case where the box contains the origin
-    if (aaBoxContains(origin, corner, scale)) {
-        // We still want to calculate the distance from the origin to the inside out plane
-        float axisDistance;
-        if ((findInsideOutIntersection(origin.x, direction.x, corner.x, scale.x, axisDistance) && axisDistance >= 0 &&
-            isWithin(origin.y + axisDistance * direction.y, corner.y, scale.y) &&
-            isWithin(origin.z + axisDistance * direction.z, corner.z, scale.z))) {
-            distance = axisDistance;
-            face = direction.x > 0 ? MAX_X_FACE : MIN_X_FACE;
-            surfaceNormal = glm::vec3(direction.x > 0 ? 1.0f : -1.0f, 0.0f, 0.0f);
-            return true;
-        }
-        if ((findInsideOutIntersection(origin.y, direction.y, corner.y, scale.y, axisDistance) && axisDistance >= 0 &&
-            isWithin(origin.x + axisDistance * direction.x, corner.x, scale.x) &&
-            isWithin(origin.z + axisDistance * direction.z, corner.z, scale.z))) {
-            distance = axisDistance;
-            face = direction.y > 0 ? MAX_Y_FACE : MIN_Y_FACE;
-            surfaceNormal = glm::vec3(0.0f, direction.y > 0 ? 1.0f : -1.0f, 0.0f);
-            return true;
-        }
-        if ((findInsideOutIntersection(origin.z, direction.z, corner.z, scale.z, axisDistance) && axisDistance >= 0 &&
-            isWithin(origin.y + axisDistance * direction.y, corner.y, scale.y) &&
-            isWithin(origin.x + axisDistance * direction.x, corner.x, scale.x))) {
-            distance = axisDistance;
-            face = direction.z > 0 ? MAX_Z_FACE : MIN_Z_FACE;
-            surfaceNormal = glm::vec3(0.0f, 0.0f, direction.z > 0 ? 1.0f : -1.0f);
-            return true;
-        }
-        // This case is unexpected, but mimics the previous behavior for inside out intersections
-        distance = 0;
-        return true;
+// https://tavianator.com/fast-branchless-raybounding-box-intersections/
+bool findRayAABoxIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
+                              const glm::vec3& corner, const glm::vec3& scale, float& distance, BoxFace& face, glm::vec3& surfaceNormal) {
+    float t1, t2, newTmin, newTmax, tmin = -INFINITY, tmax = INFINITY;
+    int minAxis = -1, maxAxis = -1;
+
+    for (int i = 0; i < 3; ++i) {
+        t1 = (corner[i] - origin[i]) * invDirection[i];
+        t2 = (corner[i] + scale[i] - origin[i]) * invDirection[i];
+
+        newTmin = glm::min(t1, t2);
+        newTmax = glm::max(t1, t2);
+
+        minAxis = newTmin > tmin ? i : minAxis;
+        tmin = glm::max(tmin, newTmin);
+        maxAxis = newTmax < tmax ? i : maxAxis;
+        tmax = glm::min(tmax, newTmax);
     }
 
-    // check each axis
-    float axisDistance;
-    if ((findIntersection(origin.x, direction.x, corner.x, scale.x, axisDistance) && axisDistance >= 0 &&
-        isWithin(origin.y + axisDistance * direction.y, corner.y, scale.y) &&
-        isWithin(origin.z + axisDistance * direction.z, corner.z, scale.z))) {
-        distance = axisDistance;
-        face = direction.x > 0 ? MIN_X_FACE : MAX_X_FACE;
-        surfaceNormal = glm::vec3(direction.x > 0 ? -1.0f : 1.0f, 0.0f, 0.0f);
-        return true;
-    }
-    if ((findIntersection(origin.y, direction.y, corner.y, scale.y, axisDistance) && axisDistance >= 0 &&
-        isWithin(origin.x + axisDistance * direction.x, corner.x, scale.x) &&
-        isWithin(origin.z + axisDistance * direction.z, corner.z, scale.z))) {
-        distance = axisDistance;
-        face = direction.y > 0 ? MIN_Y_FACE : MAX_Y_FACE;
-        surfaceNormal = glm::vec3(0.0f, direction.y > 0 ? -1.0f : 1.0f, 0.0f);
-        return true;
-    }
-    if ((findIntersection(origin.z, direction.z, corner.z, scale.z, axisDistance) && axisDistance >= 0 &&
-        isWithin(origin.y + axisDistance * direction.y, corner.y, scale.y) &&
-        isWithin(origin.x + axisDistance * direction.x, corner.x, scale.x))) {
-        distance = axisDistance;
-        face = direction.z > 0 ? MIN_Z_FACE : MAX_Z_FACE;
-        surfaceNormal = glm::vec3(0.0f, 0.0f, direction.z > 0 ? -1.0f : 1.0f);
+    if (tmax >= glm::max(tmin, 0.0f)) {
+        if (tmin < 0.0f) {
+            distance = tmax;
+            bool positiveDirection = direction[maxAxis] > 0.0f;
+            surfaceNormal = glm::vec3(0.0f);
+            surfaceNormal[maxAxis] = positiveDirection ? -1.0f : 1.0f;
+            face = positiveDirection ? BoxFace(2 * maxAxis + 1) : BoxFace(2 * maxAxis);
+        } else {
+            distance = tmin;
+            bool positiveDirection = direction[minAxis] > 0.0f;
+            surfaceNormal = glm::vec3(0.0f);
+            surfaceNormal[minAxis] = positiveDirection ? -1.0f : 1.0f;
+            face = positiveDirection ? BoxFace(2 * minAxis) : BoxFace(2 * minAxis + 1);
+        }
         return true;
     }
     return false;
@@ -286,12 +260,13 @@ bool findRaySphereIntersection(const glm::vec3& origin, const glm::vec3& directi
         distance = 0.0f;
         return true; // starts inside the sphere
     }
-    float b = glm::dot(direction, relativeOrigin);
-    float radicand = b * b - c;
+    float b = 2.0f * glm::dot(direction, relativeOrigin);
+    float a = glm::dot(direction, direction);
+    float radicand = b * b - 4.0f * a * c;
     if (radicand < 0.0f) {
         return false; // doesn't hit the sphere
     }
-    float t = -b - sqrtf(radicand);
+    float t = 0.5f * (-b - sqrtf(radicand)) / a;
     if (t < 0.0f) {
         return false; // doesn't hit the sphere
     }
@@ -391,24 +366,34 @@ Triangle Triangle::operator*(const glm::mat4& transform) const {
     };
 }
 
+// https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
 bool findRayTriangleIntersection(const glm::vec3& origin, const glm::vec3& direction,
         const glm::vec3& v0, const glm::vec3& v1, const glm::vec3& v2, float& distance, bool allowBackface) {
-    glm::vec3 firstSide = v0 - v1;
-    glm::vec3 secondSide = v2 - v1;
-    glm::vec3 normal = glm::cross(secondSide, firstSide);
-    float dividend = glm::dot(normal, v1) - glm::dot(origin, normal);
-    if (!allowBackface && dividend > 0.0f) {
-        return false; // origin below plane
-    }
-    float divisor = glm::dot(normal, direction);
-    if (divisor >= 0.0f) {
+    glm::vec3 firstSide = v1 - v0;
+    glm::vec3 secondSide = v2 - v0;
+    glm::vec3 P = glm::cross(direction, secondSide);
+    float det = glm::dot(firstSide, P);
+    if (!allowBackface && det < EPSILON) {
+        return false;
+    } else if (fabsf(det) < EPSILON) {
         return false;
     }
-    float t = dividend / divisor;
-    glm::vec3 point = origin + direction * t;
-    if (glm::dot(normal, glm::cross(point - v1, firstSide)) > 0.0f &&
-            glm::dot(normal, glm::cross(secondSide, point - v1)) > 0.0f &&
-            glm::dot(normal, glm::cross(point - v0, v2 - v0)) > 0.0f) {
+
+    float invDet = 1.0f / det;
+    glm::vec3 T = origin - v0;
+    float u = glm::dot(T, P) * invDet;
+    if (u < 0.0f || u > 1.0f) {
+        return false;
+    }
+
+    glm::vec3 Q = glm::cross(T, firstSide);
+    float v = glm::dot(direction, Q) * invDet;
+    if (v < 0.0f || u + v > 1.0f) {
+        return false;
+    }
+
+    float t = glm::dot(secondSide, Q) * invDet;
+    if (t > EPSILON) {
         distance = t;
         return true;
     }
diff --git a/libraries/shared/src/GeometryUtil.h b/libraries/shared/src/GeometryUtil.h
index 54f9062469..8ec75f71bd 100644
--- a/libraries/shared/src/GeometryUtil.h
+++ b/libraries/shared/src/GeometryUtil.h
@@ -76,8 +76,8 @@ glm::vec3 addPenetrations(const glm::vec3& currentPenetration, const glm::vec3&
 
 bool findIntersection(float origin, float direction, float corner, float size, float& distance);
 bool findInsideOutIntersection(float origin, float direction, float corner, float size, float& distance);
-bool findRayAABoxIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& corner, const glm::vec3& scale, float& distance,
-                              BoxFace& face, glm::vec3& surfaceNormal);
+bool findRayAABoxIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
+                              const glm::vec3& corner, const glm::vec3& scale, float& distance, BoxFace& face, glm::vec3& surfaceNormal);
 
 bool findRaySphereIntersection(const glm::vec3& origin, const glm::vec3& direction,
     const glm::vec3& center, float radius, float& distance);
diff --git a/libraries/shared/src/TriangleSet.cpp b/libraries/shared/src/TriangleSet.cpp
index cde9c20cab..02a8d458a9 100644
--- a/libraries/shared/src/TriangleSet.cpp
+++ b/libraries/shared/src/TriangleSet.cpp
@@ -13,6 +13,8 @@
 
 #include "GLMHelpers.h"
 
+#include <list>
+
 void TriangleSet::insert(const Triangle& t) {
     _isBalanced = false;
 
@@ -27,48 +29,7 @@ void TriangleSet::clear() {
     _bounds.clear();
     _isBalanced = false;
 
-    _triangleOctree.clear();
-}
-
-bool TriangleSet::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction,
-         float& distance, BoxFace& face, Triangle& triangle, bool precision, bool allowBackface) {
-
-    // reset our distance to be the max possible, lower level tests will store best distance here
-    distance = std::numeric_limits<float>::max();
-
-    if (!_isBalanced) {
-        balanceOctree();
-    }
-
-    int trianglesTouched = 0;
-    auto result = _triangleOctree.findRayIntersection(origin, direction, distance, face, triangle, precision, trianglesTouched, allowBackface);
-
-    #if WANT_DEBUGGING
-    if (precision) {
-        qDebug() << "trianglesTouched :" << trianglesTouched << "out of:" << _triangleOctree._population << "_triangles.size:" << _triangles.size();
-    }
-    #endif
-    return result;
-}
-
-bool TriangleSet::findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
-    float& parabolicDistance, BoxFace& face, Triangle& triangle, bool precision, bool allowBackface) {
-    // reset our distance to be the max possible, lower level tests will store best distance here
-    parabolicDistance = FLT_MAX;
-
-    if (!_isBalanced) {
-        balanceOctree();
-    }
-
-    int trianglesTouched = 0;
-    auto result = _triangleOctree.findParabolaIntersection(origin, velocity, acceleration, parabolicDistance, face, triangle, precision, trianglesTouched, allowBackface);
-
-#if WANT_DEBUGGING
-    if (precision) {
-        qDebug() << "trianglesTouched :" << trianglesTouched << "out of:" << _triangleOctree._population << "_triangles.size:" << _triangles.size();
-    }
-#endif
-    return result;
+    _triangleTree.clear();
 }
 
 bool TriangleSet::convexHullContains(const glm::vec3& point) const {
@@ -92,33 +53,158 @@ void TriangleSet::debugDump() {
     qDebug() << __FUNCTION__;
     qDebug() << "bounds:" << getBounds();
     qDebug() << "triangles:" << size() << "at top level....";
-    qDebug() << "----- _triangleOctree -----";
-    _triangleOctree.debugDump();
+    qDebug() << "----- _triangleTree -----";
+    _triangleTree.debugDump();
 }
 
-void TriangleSet::balanceOctree() {
-    _triangleOctree.reset(_bounds, 0);
+void TriangleSet::balanceTree() {
+    _triangleTree.reset(_bounds);
 
     // insert all the triangles
     for (size_t i = 0; i < _triangles.size(); i++) {
-        _triangleOctree.insert(i);
+        _triangleTree.insert(i);
     }
 
     _isBalanced = true;
 
-    #if WANT_DEBUGGING
+#if WANT_DEBUGGING
     debugDump();
-    #endif
+#endif
 }
 
+// With an octree: 8 ^ MAX_DEPTH = 4096 leaves
+//static const int MAX_DEPTH = 4;
+// With a k-d tree: 2 ^ MAX_DEPTH = 4096 leaves
+static const int MAX_DEPTH = 12;
+
+TriangleSet::TriangleTreeCell::TriangleTreeCell(std::vector<Triangle>& allTriangles, const AABox& bounds, int depth) :
+    _allTriangles(allTriangles)
+{
+    reset(bounds, depth);
+}
+
+void TriangleSet::TriangleTreeCell::clear() {
+    _population = 0;
+    _triangleIndices.clear();
+    _bounds.clear();
+    _children.first.reset();
+    _children.second.reset();
+}
+
+void TriangleSet::TriangleTreeCell::reset(const AABox& bounds, int depth) {
+    clear();
+    _bounds = bounds;
+    _depth = depth;
+}
+
+void TriangleSet::TriangleTreeCell::debugDump() {
+    qDebug() << __FUNCTION__;
+    qDebug() << "               bounds:" << getBounds();
+    qDebug() << "                depth:" << _depth;
+    qDebug() << "           population:" << _population << "this level or below"
+             << " ---- triangleIndices:" << _triangleIndices.size() << "in this cell";
+
+    int numChildren = 0;
+    if (_children.first) {
+        numChildren++;
+    } else if (_children.second) {
+        numChildren++;
+    }
+    qDebug() << "child cells:" << numChildren;
+    if (_depth < MAX_DEPTH) {
+        if (_children.first) {
+            qDebug() << "child: 0";
+            _children.first->debugDump();
+        }
+        if (_children.second) {
+            qDebug() << "child: 1";
+            _children.second->debugDump();
+        }
+    }
+}
+
+std::pair<AABox, AABox> TriangleSet::TriangleTreeCell::getTriangleTreeCellChildBounds() {
+    std::pair<AABox, AABox> toReturn;
+    int axis = 0;
+    // find biggest axis
+    glm::vec3 dimensions = _bounds.getDimensions();
+    for (int i = 0; i < 3; i++) {
+        if (dimensions[i] >= dimensions[(i + 1) % 3] && dimensions[i] >= dimensions[(i + 2) % 3]) {
+            axis = i;
+            break;
+        }
+    }
+
+    // The new boxes are half the size in the largest dimension
+    glm::vec3 newDimensions = dimensions;
+    newDimensions[axis] *= 0.5f;
+    toReturn.first.setBox(_bounds.getCorner(), newDimensions);
+    glm::vec3 offset = glm::vec3(0.0f);
+    offset[axis] = newDimensions[axis];
+    toReturn.second.setBox(_bounds.getCorner() + offset, newDimensions);
+    return toReturn;
+}
+
+void TriangleSet::TriangleTreeCell::insert(size_t triangleIndex) {
+    _population++;
+
+    // if we're not yet at the max depth, then check which child the triangle fits in
+    if (_depth < MAX_DEPTH) {
+        const Triangle& triangle = _allTriangles[triangleIndex];
+        auto childBounds = getTriangleTreeCellChildBounds();
+
+        auto insertOperator = [&](const AABox& childBound, std::shared_ptr<TriangleTreeCell>& child) {
+            // if the child AABox would contain the triangle...
+            if (childBound.contains(triangle)) {
+                // if the child cell doesn't yet exist, create it...
+                if (!child) {
+                    child = std::make_shared<TriangleTreeCell>(_allTriangles, childBound, _depth + 1);
+                }
+
+                // insert the triangleIndex in the child cell
+                child->insert(triangleIndex);
+                return true;
+            }
+            return false;
+        };
+        if (insertOperator(childBounds.first, _children.first) || insertOperator(childBounds.second, _children.second)) {
+            return;
+        }
+    }
+    // either we're at max depth, or the triangle doesn't fit in one of our
+    // children and so we want to just record it here
+    _triangleIndices.push_back(triangleIndex);
+}
+
+bool TriangleSet::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection, float& distance,
+                                      BoxFace& face, Triangle& triangle, bool precision, bool allowBackface) {
+    if (!_isBalanced) {
+        balanceTree();
+    }
+
+    float localDistance = distance;
+    int trianglesTouched = 0;
+    bool hit = _triangleTree.findRayIntersection(origin, direction, invDirection, localDistance, face, triangle, precision, trianglesTouched, allowBackface);
+    if (hit) {
+        distance = localDistance;
+    }
+
+#if WANT_DEBUGGING
+    if (precision) {
+        qDebug() << "trianglesTouched :" << trianglesTouched << "out of:" << _triangleTree._population << "_triangles.size:" << _triangles.size();
+    }
+#endif
+    return hit;
+}
 
 // Determine of the given ray (origin/direction) in model space intersects with any triangles
 // in the set. If an intersection occurs, the distance and surface normal will be provided.
-bool TriangleSet::TriangleOctreeCell::findRayIntersectionInternal(const glm::vec3& origin, const glm::vec3& direction,
+bool TriangleSet::TriangleTreeCell::findRayIntersectionInternal(const glm::vec3& origin, const glm::vec3& direction,
                                                                   float& distance, BoxFace& face, Triangle& triangle, bool precision,
                                                                   int& trianglesTouched, bool allowBackface) {
     bool intersectedSomething = false;
     float bestDistance = FLT_MAX;
+    Triangle bestTriangle;
 
     if (precision) {
         for (const auto& triangleIndex : _triangleIndices) {
@@ -128,8 +214,8 @@ bool TriangleSet::TriangleOctreeCell::findRayIntersectionInternal(const glm::vec
             if (findRayTriangleIntersection(origin, direction, thisTriangle, thisTriangleDistance, allowBackface)) {
                 if (thisTriangleDistance < bestDistance) {
                     bestDistance = thisTriangleDistance;
+                    bestTriangle = thisTriangle;
                     intersectedSomething = true;
-                    triangle = thisTriangle;
                 }
             }
         }
@@ -140,17 +226,133 @@ bool TriangleSet::TriangleOctreeCell::findRayIntersectionInternal(const glm::vec
 
     if (intersectedSomething) {
         distance = bestDistance;
+        triangle = bestTriangle;
     }
 
     return intersectedSomething;
 }
 
-bool TriangleSet::TriangleOctreeCell::findParabolaIntersectionInternal(const glm::vec3& origin, const glm::vec3& velocity,
+bool TriangleSet::TriangleTreeCell::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
+                                                        float& distance, BoxFace& face, Triangle& triangle, bool precision, int& trianglesTouched,
+                                                        bool allowBackface) {
+    if (_population < 1) {
+        return false; // no triangles below here, so we can't intersect
+    }
+
+    float bestLocalDistance = FLT_MAX;
+    BoxFace bestLocalFace;
+    Triangle bestLocalTriangle;
+    bool intersects = false;
+
+    // Check our local triangle set first
+    // The distance passed in here is the distance to our bounding box.  If !precision, that distance is used
+    {
+        float internalDistance = distance;
+        BoxFace internalFace;
+        Triangle internalTriangle;
+        if (findRayIntersectionInternal(origin, direction, internalDistance, internalFace, internalTriangle, precision, trianglesTouched, allowBackface)) {
+            bestLocalDistance = internalDistance;
+            bestLocalFace = internalFace;
+            bestLocalTriangle = internalTriangle;
+            intersects = true;
+        }
+    }
+
+    // if we're not yet at the max depth, then check our children
+    if (_depth < MAX_DEPTH) {
+        std::list<SortedTriangleCell> sortedTriangleCells;
+        auto sortingOperator = [&](std::shared_ptr<TriangleTreeCell>& child) {
+            if (child) {
+                float priority = FLT_MAX;
+                if (child->getBounds().contains(origin)) {
+                    priority = 0.0f;
+                } else {
+                    float childBoundDistance = FLT_MAX;
+                    BoxFace childBoundFace;
+                    glm::vec3 childBoundNormal;
+                    if (child->getBounds().findRayIntersection(origin, direction, invDirection, childBoundDistance, childBoundFace, childBoundNormal)) {
+                        // We only need to add this cell if it's closer than the local triangle set intersection (if there was one)
+                        if (childBoundDistance < bestLocalDistance) {
+                            priority = childBoundDistance;
+                        }
+                    }
+                }
+
+                if (priority < FLT_MAX) {
+                    if (sortedTriangleCells.size() > 0 && priority < sortedTriangleCells.front().first) {
+                        sortedTriangleCells.emplace_front(priority, child);
+                    } else {
+                        sortedTriangleCells.emplace_back(priority, child);
+                    }
+                }
+            }
+        };
+        sortingOperator(_children.first);
+        sortingOperator(_children.second);
+
+        for (auto it = sortedTriangleCells.begin(); it != sortedTriangleCells.end(); ++it) {
+            const SortedTriangleCell& sortedTriangleCell = *it;
+            float childDistance = sortedTriangleCell.first;
+            // We can exit once childDistance > bestLocalDistance
+            if (childDistance > bestLocalDistance) {
+                break;
+            }
+            // If we're inside the child cell and !precision, we need the actual distance to the cell bounds
+            if (!precision && childDistance < EPSILON) {
+                BoxFace childBoundFace;
+                glm::vec3 childBoundNormal;
+                sortedTriangleCell.second->getBounds().findRayIntersection(origin, direction, invDirection, childDistance, childBoundFace, childBoundNormal);
+            }
+            BoxFace childFace;
+            Triangle childTriangle;
+            if (sortedTriangleCell.second->findRayIntersection(origin, direction, invDirection, childDistance, childFace, childTriangle, precision, trianglesTouched)) {
+                if (childDistance < bestLocalDistance) {
+                    bestLocalDistance = childDistance;
+                    bestLocalFace = childFace;
+                    bestLocalTriangle = childTriangle;
+                    intersects = true;
+                    break;
+                }
+            }
+        }
+    }
+
+    if (intersects) {
+        distance = bestLocalDistance;
+        face = bestLocalFace;
+        triangle = bestLocalTriangle;
+    }
+    return intersects;
+}
+
+bool TriangleSet::findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
+                                           float& parabolicDistance, BoxFace& face, Triangle& triangle, bool precision, bool allowBackface) {
+    if (!_isBalanced) {
+        balanceTree();
+    }
+
+    float localDistance = parabolicDistance;
+    int trianglesTouched = 0;
+    bool hit = _triangleTree.findParabolaIntersection(origin, velocity, acceleration, localDistance, face, triangle, precision, trianglesTouched, allowBackface);
+    if (hit) {
+        parabolicDistance = localDistance;
+    }
+
+#if WANT_DEBUGGING
+    if (precision) {
+        qDebug() << "trianglesTouched :" << trianglesTouched << "out of:" << _triangleTree._population << "_triangles.size:" << _triangles.size();
+    }
+#endif
+    return hit;
+}
+
+bool TriangleSet::TriangleTreeCell::findParabolaIntersectionInternal(const glm::vec3& origin, const glm::vec3& velocity,
                                                                        const glm::vec3& acceleration, float& parabolicDistance,
                                                                        BoxFace& face, Triangle& triangle, bool precision,
                                                                        int& trianglesTouched, bool allowBackface) {
     bool intersectedSomething = false;
     float bestDistance = FLT_MAX;
+    Triangle bestTriangle;
 
     if (precision) {
         for (const auto& triangleIndex : _triangleIndices) {
@@ -160,8 +362,8 @@ bool TriangleSet::TriangleOctreeCell::findParabolaIntersectionInternal(const glm
             if (findParabolaTriangleIntersection(origin, velocity, acceleration, thisTriangle, thisTriangleDistance, allowBackface)) {
                 if (thisTriangleDistance < bestDistance) {
                     bestDistance = thisTriangleDistance;
+                    bestTriangle = thisTriangle;
                     intersectedSomething = true;
-                    triangle = thisTriangle;
                 }
             }
         }
@@ -172,146 +374,13 @@ bool TriangleSet::TriangleOctreeCell::findParabolaIntersectionInternal(const glm
 
     if (intersectedSomething) {
         parabolicDistance = bestDistance;
+        triangle = bestTriangle;
     }
 
     return intersectedSomething;
 }
 
-static const int MAX_DEPTH = 4; // for now
-static const int MAX_CHILDREN = 8;
-
-TriangleSet::TriangleOctreeCell::TriangleOctreeCell(std::vector<Triangle>& allTriangles, const AABox& bounds, int depth) :
-    _allTriangles(allTriangles)
-{
-    reset(bounds, depth);
-}
-
-void TriangleSet::TriangleOctreeCell::clear() {
-    _population = 0;
-    _triangleIndices.clear();
-    _bounds.clear();
-    _children.clear();
-}
-
-void TriangleSet::TriangleOctreeCell::reset(const AABox& bounds, int depth) {
-    clear();
-    _bounds = bounds;
-    _depth = depth;
-}
-
-void TriangleSet::TriangleOctreeCell::debugDump() {
-    qDebug() << __FUNCTION__;
-    qDebug() << "bounds:" << getBounds();
-    qDebug() << "depth:" << _depth;
-    qDebug() << "population:" << _population << "this level or below" 
-             << " ---- triangleIndices:" << _triangleIndices.size() << "in this cell";
-
-    qDebug() << "child cells:" << _children.size();
-    if (_depth < MAX_DEPTH) {
-        int childNum = 0;
-        for (auto& child : _children) {
-            qDebug() << "child:" << childNum;
-            child.second.debugDump();
-            childNum++;
-        }
-    }
-}
-
-void TriangleSet::TriangleOctreeCell::insert(size_t triangleIndex) {
-    const Triangle& triangle = _allTriangles[triangleIndex];
-    _population++;
-    // if we're not yet at the max depth, then check which child the triangle fits in
-    if (_depth < MAX_DEPTH) {
-
-        for (int child = 0; child < MAX_CHILDREN; child++) {
-            AABox childBounds = getBounds().getOctreeChild((AABox::OctreeChild)child);
-
-
-            // if the child AABox would contain the triangle...
-            if (childBounds.contains(triangle)) {
-                // if the child cell doesn't yet exist, create it...
-                if (_children.find((AABox::OctreeChild)child) == _children.end()) {
-                    _children.insert(
-                        std::pair<AABox::OctreeChild, TriangleOctreeCell>
-                        ((AABox::OctreeChild)child, TriangleOctreeCell(_allTriangles, childBounds, _depth + 1)));
-                }
-
-                // insert the triangleIndex in the child cell
-                _children.at((AABox::OctreeChild)child).insert(triangleIndex);
-                return;
-            }
-        }
-    }
-    // either we're at max depth, or the triangle doesn't fit in one of our
-    // children and so we want to just record it here
-    _triangleIndices.push_back(triangleIndex);
-}
-
-bool TriangleSet::TriangleOctreeCell::findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, float& distance,
-                                                          BoxFace& face, Triangle& triangle, bool precision, int& trianglesTouched,
-                                                          bool allowBackface) {
-    if (_population < 1) {
-        return false; // no triangles below here, so we can't intersect
-    }
-
-    float bestLocalDistance = FLT_MAX;
-    BoxFace bestLocalFace;
-    Triangle bestLocalTriangle;
-    glm::vec3 bestLocalNormal;
-    bool intersects = false;
-
-    float boxDistance = FLT_MAX;
-    // if the pick intersects our bounding box, then continue
-    if (getBounds().findRayIntersection(origin, direction, boxDistance, bestLocalFace, bestLocalNormal)) {
-        // if the intersection with our bounding box, is greater than the current best distance (the distance passed in)
-        // then we know that none of our triangles can represent a better intersection and we can return
-        if (boxDistance > distance) {
-            return false;
-        }
-
-        // if we're not yet at the max depth, then check which child the triangle fits in
-        if (_depth < MAX_DEPTH) {
-            float bestChildDistance = FLT_MAX;
-            for (auto& child : _children) {
-                // check each child, if there's an intersection, it will return some distance that we need
-                // to compare against the other results, because there might be multiple intersections and
-                // we will always choose the best (shortest) intersection
-                float childDistance = bestChildDistance;
-                BoxFace childFace;
-                Triangle childTriangle;
-                if (child.second.findRayIntersection(origin, direction, childDistance, childFace, childTriangle, precision, trianglesTouched)) {
-                    if (childDistance < bestLocalDistance) {
-                        bestLocalDistance = childDistance;
-                        bestChildDistance = childDistance;
-                        bestLocalFace = childFace;
-                        bestLocalTriangle = childTriangle;
-                        intersects = true;
-                    }
-                }
-            }
-        }
-        // also check our local triangle set
-        float internalDistance = boxDistance;
-        BoxFace internalFace;
-        Triangle internalTriangle;
-        if (findRayIntersectionInternal(origin, direction, internalDistance, internalFace, internalTriangle, precision, trianglesTouched, allowBackface)) {
-            if (internalDistance < bestLocalDistance) {
-                bestLocalDistance = internalDistance;
-                bestLocalFace = internalFace;
-                bestLocalTriangle = internalTriangle;
-                intersects = true;
-            }
-        }
-    }
-    if (intersects) {
-        distance = bestLocalDistance;
-        face = bestLocalFace;
-        triangle = bestLocalTriangle;
-    }
-    return intersects;
-}
-
-bool TriangleSet::TriangleOctreeCell::findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity,
+bool TriangleSet::TriangleTreeCell::findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity,
                                                                const glm::vec3& acceleration, float& parabolicDistance,
                                                                BoxFace& face, Triangle& triangle, bool precision,
                                                                int& trianglesTouched, bool allowBackface) {
@@ -322,52 +391,81 @@ bool TriangleSet::TriangleOctreeCell::findParabolaIntersection(const glm::vec3&
     float bestLocalDistance = FLT_MAX;
     BoxFace bestLocalFace;
     Triangle bestLocalTriangle;
-    glm::vec3 bestLocalNormal;
     bool intersects = false;
 
-    float boxDistance = FLT_MAX;
-    // if the pick intersects our bounding box, then continue
-    if (getBounds().findParabolaIntersection(origin, velocity, acceleration, boxDistance, bestLocalFace, bestLocalNormal)) {
-        // if the intersection with our bounding box, is greater than the current best distance (the distance passed in)
-        // then we know that none of our triangles can represent a better intersection and we can return
-        if (boxDistance > parabolicDistance) {
-            return false;
-        }
-
-        // if we're not yet at the max depth, then check which child the triangle fits in
-        if (_depth < MAX_DEPTH) {
-            float bestChildDistance = FLT_MAX;
-            for (auto& child : _children) {
-                // check each child, if there's an intersection, it will return some distance that we need
-                // to compare against the other results, because there might be multiple intersections and
-                // we will always choose the best (shortest) intersection
-                float childDistance = bestChildDistance;
-                BoxFace childFace;
-                Triangle childTriangle;
-                if (child.second.findParabolaIntersection(origin, velocity, acceleration, childDistance, childFace, childTriangle, precision, trianglesTouched)) {
-                    if (childDistance < bestLocalDistance) {
-                        bestLocalDistance = childDistance;
-                        bestChildDistance = childDistance;
-                        bestLocalFace = childFace;
-                        bestLocalTriangle = childTriangle;
-                        intersects = true;
-                    }
-                }
-            }
-        }
-        // also check our local triangle set
-        float internalDistance = boxDistance;
+    // Check our local triangle set first
+    // The distance passed in here is the distance to our bounding box.  If !precision, that distance is used
+    {
+        float internalDistance = parabolicDistance;
         BoxFace internalFace;
         Triangle internalTriangle;
         if (findParabolaIntersectionInternal(origin, velocity, acceleration, internalDistance, internalFace, internalTriangle, precision, trianglesTouched, allowBackface)) {
-            if (internalDistance < bestLocalDistance) {
-                bestLocalDistance = internalDistance;
-                bestLocalFace = internalFace;
-                bestLocalTriangle = internalTriangle;
-                intersects = true;
+            bestLocalDistance = internalDistance;
+            bestLocalFace = internalFace;
+            bestLocalTriangle = internalTriangle;
+            intersects = true;
+        }
+    }
+
+    // if we're not yet at the max depth, then check our children
+    if (_depth < MAX_DEPTH) {
+        std::list<SortedTriangleCell> sortedTriangleCells;
+        auto sortingOperator = [&](std::shared_ptr<TriangleTreeCell>& child) {
+            if (child) {
+                float priority = FLT_MAX;
+                if (child->getBounds().contains(origin)) {
+                    priority = 0.0f;
+                } else {
+                    float childBoundDistance = FLT_MAX;
+                    BoxFace childBoundFace;
+                    glm::vec3 childBoundNormal;
+                    if (child->getBounds().findParabolaIntersection(origin, velocity, acceleration, childBoundDistance, childBoundFace, childBoundNormal)) {
+                        // We only need to add this cell if it's closer than the local triangle set intersection (if there was one)
+                        if (childBoundDistance < bestLocalDistance) {
+                            priority = childBoundDistance;
+                        }
+                    }
+                }
+
+                if (priority < FLT_MAX) {
+                    if (sortedTriangleCells.size() > 0 && priority < sortedTriangleCells.front().first) {
+                        sortedTriangleCells.emplace_front(priority, child);
+                    } else {
+                        sortedTriangleCells.emplace_back(priority, child);
+                    }
+                }
+            }
+        };
+        sortingOperator(_children.first);
+        sortingOperator(_children.second);
+
+        for (auto it = sortedTriangleCells.begin(); it != sortedTriangleCells.end(); ++it) {
+            const SortedTriangleCell& sortedTriangleCell = *it;
+            float childDistance = sortedTriangleCell.first;
+            // We can exit once childDistance > bestLocalDistance
+            if (childDistance > bestLocalDistance) {
+                break;
+            }
+            // If we're inside the child cell and !precision, we need the actual distance to the cell bounds
+            if (!precision && childDistance < EPSILON) {
+                BoxFace childBoundFace;
+                glm::vec3 childBoundNormal;
+                sortedTriangleCell.second->getBounds().findParabolaIntersection(origin, velocity, acceleration, childDistance, childBoundFace, childBoundNormal);
+            }
+            BoxFace childFace;
+            Triangle childTriangle;
+            if (sortedTriangleCell.second->findParabolaIntersection(origin, velocity, acceleration, childDistance, childFace, childTriangle, precision, trianglesTouched)) {
+                if (childDistance < bestLocalDistance) {
+                    bestLocalDistance = childDistance;
+                    bestLocalFace = childFace;
+                    bestLocalTriangle = childTriangle;
+                    intersects = true;
+                    break;
+                }
             }
         }
     }
+
     if (intersects) {
         parabolicDistance = bestLocalDistance;
         face = bestLocalFace;
diff --git a/libraries/shared/src/TriangleSet.h b/libraries/shared/src/TriangleSet.h
index 0b0d0a9ac5..3417b36b4a 100644
--- a/libraries/shared/src/TriangleSet.h
+++ b/libraries/shared/src/TriangleSet.h
@@ -12,23 +12,23 @@
 #pragma once
 
 #include <vector>
+#include <memory>
 
 #include "AABox.h"
 #include "GeometryUtil.h"
 
 class TriangleSet {
 
-    class TriangleOctreeCell {
+    class TriangleTreeCell {
     public:
-        TriangleOctreeCell(std::vector<Triangle>& allTriangles) :
-            _allTriangles(allTriangles)
-        { }
+        TriangleTreeCell(std::vector<Triangle>& allTriangles) : _allTriangles(allTriangles) {}
+        TriangleTreeCell(std::vector<Triangle>& allTriangles, const AABox& bounds, int depth);
 
         void insert(size_t triangleIndex);
         void reset(const AABox& bounds, int depth = 0);
         void clear();
 
-        bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction,
+        bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
             float& distance, BoxFace& face, Triangle& triangle, bool precision, int& trianglesTouched,
             bool allowBackface = false);
         bool findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
@@ -40,8 +40,6 @@ class TriangleSet {
         void debugDump();
 
     protected:
-        TriangleOctreeCell(std::vector<Triangle>& allTriangles, const AABox& bounds, int depth);
-
         // checks our internal list of triangles
         bool findRayIntersectionInternal(const glm::vec3& origin, const glm::vec3& direction,
             float& distance, BoxFace& face, Triangle& triangle, bool precision, int& trianglesTouched,
@@ -50,31 +48,33 @@ class TriangleSet {
             float& parabolicDistance, BoxFace& face, Triangle& triangle, bool precision, int& trianglesTouched,
             bool allowBackface = false);
 
+        std::pair<AABox, AABox> getTriangleTreeCellChildBounds();
+
         std::vector<Triangle>& _allTriangles;
-        std::map<AABox::OctreeChild, TriangleOctreeCell> _children;
-        int _depth{ 0 };
-        int _population{ 0 };
+        std::pair<std::shared_ptr<TriangleTreeCell>, std::shared_ptr<TriangleTreeCell>> _children;
+        int _depth { 0 };
+        int _population { 0 };
         AABox _bounds;
         std::vector<size_t> _triangleIndices;
 
         friend class TriangleSet;
     };
 
+    using SortedTriangleCell = std::pair<float, std::shared_ptr<TriangleTreeCell>>;
+
 public:
-    TriangleSet() :
-        _triangleOctree(_triangles)
-    {}
+    TriangleSet() : _triangleTree(_triangles) {}
 
     void debugDump();
 
     void insert(const Triangle& t);
 
-    bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction,
+    bool findRayIntersection(const glm::vec3& origin, const glm::vec3& direction, const glm::vec3& invDirection,
         float& distance, BoxFace& face, Triangle& triangle, bool precision, bool allowBackface = false);
     bool findParabolaIntersection(const glm::vec3& origin, const glm::vec3& velocity, const glm::vec3& acceleration,
         float& parabolicDistance, BoxFace& face, Triangle& triangle, bool precision, bool allowBackface = false);
 
-    void balanceOctree();
+    void balanceTree();
 
     void reserve(size_t size) { _triangles.reserve(size); } // reserve space in the datastructure for size number of triangles
     size_t size() const { return _triangles.size(); }
@@ -87,9 +87,8 @@ public:
     const AABox& getBounds() const { return _bounds; }
 
 protected:
-
-    bool _isBalanced{ false };
+    bool _isBalanced { false };
     std::vector<Triangle> _triangles;
-    TriangleOctreeCell _triangleOctree;
+    TriangleTreeCell _triangleTree;
     AABox _bounds;
 };
diff --git a/scripts/system/commerce/wallet.js b/scripts/system/commerce/wallet.js
index 8730273e7c..5939b36438 100644
--- a/scripts/system/commerce/wallet.js
+++ b/scripts/system/commerce/wallet.js
@@ -407,8 +407,10 @@
     }
 
     function onUsernameChanged() {
-        Settings.setValue("wallet/autoLogout", false);
-        Settings.setValue("wallet/savedUsername", "");
+        if (Account.username !== Settings.getValue("wallet/savedUsername")) {
+            Settings.setValue("wallet/autoLogout", false);
+            Settings.setValue("wallet/savedUsername", "");
+        }
     }
 
     // Function Name: fromQml()