From 3e26a5208c0158f2652bda44e26fe9ea9838e585 Mon Sep 17 00:00:00 2001
From: Andrew Jewett <jewett.aij@gmail.com>
Date: Wed, 9 Jun 2021 11:59:07 -0700
Subject: [PATCH] tensor voting now works with masked images (bugfix).

---
 bin/filter_mrc/handlers.cpp | 42 ++++++++++++++++++++-----------------
 lib/visfd/feature.hpp       | 21 +++++++++++--------
 2 files changed, 35 insertions(+), 28 deletions(-)

diff --git a/bin/filter_mrc/handlers.cpp b/bin/filter_mrc/handlers.cpp
index 7a2a124..82c1c7f 100644
--- a/bin/filter_mrc/handlers.cpp
+++ b/bin/filter_mrc/handlers.cpp
@@ -1341,7 +1341,7 @@ HandleRidgeDetector(Settings settings,
 
         // DEBUG: REMOVE THE NEXT IF STATMENT AFTER DEBUGGING IS FINISHED
         #ifndef NDEBUG
-        if ((ix==image_size[0]/2) && //if ((ix==78) &&
+        if ((ix == image_size[0] / 2) &&
             (iy == image_size[1] / 2) &&
             (iz == image_size[2] / 2))
         {
@@ -1474,17 +1474,19 @@ HandleRidgeDetector(Settings settings,
     for(int iz=0; iz < image_size[2]; iz++) {
       for(int iy=0; iy < image_size[1]; iy++) {
         for(int ix=0; ix < image_size[0]; ix++) {
-          float diagonalized_hessian[6];
-          DiagonalizeFlatSym3(tmp_tensor.aaaafI[iz][iy][ix],
-                              diagonalized_hessian,
-                              eival_order);
-          float score = ScoreTensorPlanar(diagonalized_hessian);
-          float peak_height = 1.0;
-          if (tomo_background.aaafI)
-            peak_height = (tomo_in.aaafI[iz][iy][ix] -
+          if (tmp_tensor.aaaafI[iz][iy][ix]) {
+            float diagonalized_hessian[6];
+            DiagonalizeFlatSym3(tmp_tensor.aaaafI[iz][iy][ix],
+                                diagonalized_hessian,
+                                eival_order);
+            float score = ScoreTensorPlanar(diagonalized_hessian);
+            float peak_height = 1.0;
+            if (tomo_background.aaafI)
+              peak_height = (tomo_in.aaafI[iz][iy][ix] -
                            tomo_background.aaafI[iz][iy][ix]);
-          score *= peak_height;
-          tomo_out.aaafI[iz][iy][ix] = score;
+            score *= peak_height;
+            tomo_out.aaafI[iz][iy][ix] = score;
+          }
         }
       }
     }
@@ -1502,14 +1504,16 @@ HandleRidgeDetector(Settings settings,
     for(int iz=0; iz < image_size[2]; iz++) {
       for(int iy=0; iy < image_size[1]; iy++) {
         for(int ix=0; ix < image_size[0]; ix++) {
-          float eivals[3];
-          float eivects[3][3];
-          ConvertFlatSym2Evects3(tmp_tensor.aaaafI[iz][iy][ix],
-                                 eivals,
-                                 eivects,
-                                 eival_order);
-          for (int d=0; d<3; d++)
-            aaaafDirection[iz][iy][ix][d] = eivects[0][d];
+          if (tmp_tensor.aaaafI[iz][iy][ix]) {
+            float eivals[3];
+            float eivects[3][3];
+            ConvertFlatSym2Evects3(tmp_tensor.aaaafI[iz][iy][ix],
+                                   eivals,
+                                   eivects,
+                                   eival_order);
+            for (int d=0; d<3; d++)
+              aaaafDirection[iz][iy][ix][d] = eivects[0][d];
+          }
         }
       }
     }
diff --git a/lib/visfd/feature.hpp b/lib/visfd/feature.hpp
index bdf6604..36c9ea4 100644
--- a/lib/visfd/feature.hpp
+++ b/lib/visfd/feature.hpp
@@ -2281,11 +2281,13 @@ class TV3D {
             for (Integer ix=0; ix<image_size[0]; ix++) {
               if ((! aaafMaskDest) || (aaafMaskDest[iz][iy][ix] == 0))
                 continue;
-              assert(aaafDenominator[iz][iy][ix] > 0.0);
-              for (int di = 0; di < 3; di++) {
-                for (int dj = di; dj < 3; dj++) {
-                  aaaafDest[iz][iy][ix][ MapIndices_3x3_to_linear[di][dj] ]
-                    /= aaafDenominator[iz][iy][ix];
+              if (aaafDenominator[iz][iy][ix] > 0.0) {
+                assert(aaaafDest[iz][iy][ix]);
+                for (int di = 0; di < 3; di++) {
+                  for (int dj = di; dj < 3; dj++) {
+                    aaaafDest[iz][iy][ix][ MapIndices_3x3_to_linear[di][dj] ]
+                      /= aaafDenominator[iz][iy][ix];
+                  }
                 }
               }
             }
@@ -2414,7 +2416,6 @@ class TV3D {
     assert(aaaafV);
     assert(aaaafDest);
 
-
     //optional: count the number of voxels which can
     //          cast votes (useful for benchmarking)
     if (pReportProgress) {
@@ -2444,9 +2445,10 @@ class TV3D {
     for (Integer iz=0; iz<image_size[2]; iz++)
       for (Integer iy=0; iy<image_size[1]; iy++)
         for (Integer ix=0; ix<image_size[0]; ix++)
-          for (int di=0; di<3; di++)
-            for (int dj=0; dj<3; dj++)
-              aaaafDest[iz][iy][ix][ MapIndices_3x3_to_linear[di][dj] ] = 0.0;
+          if (aaaafDest[iz][iy][ix])
+            for (int di=0; di<3; di++)
+              for (int dj=0; dj<3; dj++)
+                aaaafDest[iz][iy][ix][ MapIndices_3x3_to_linear[di][dj] ] = 0.0;
 
     if (aaafDenominator) {
       // aaafDenominator[][][] keeps track of how much of the sum of
@@ -2583,6 +2585,7 @@ class TV3D {
 
           if ((aaafMaskDest) && (aaafMaskDest[iz_jz][iy_jy][ix_jx] == 0.0))
             continue;
+          assert(aaaafDest[iz_jz][iy_jy][ix_jx]);
 
           // The function describing how the vote-strength falls off with
           // distance has been precomputed and is stored in