diff --git a/libs/jpegrecoverymap/include/jpegrecoverymap/recoverymapmath.h b/libs/jpegrecoverymap/include/jpegrecoverymap/recoverymapmath.h
index 0fb64d3..2f9799d 100644
--- a/libs/jpegrecoverymap/include/jpegrecoverymap/recoverymapmath.h
+++ b/libs/jpegrecoverymap/include/jpegrecoverymap/recoverymapmath.h
@@ -112,6 +112,52 @@
   return temp /= rhs;
 }
 
+struct ShepardsIDW {
+  ShepardsIDW(int mapScaleFactor) : mMapScaleFactor{mapScaleFactor} {
+    const int size = mMapScaleFactor * mMapScaleFactor * 4;
+    mWeights = new float[size];
+    mWeightsNR = new float[size];
+    mWeightsNB = new float[size];
+    mWeightsC = new float[size];
+    fillShepardsIDW(mWeights, 1, 1);
+    fillShepardsIDW(mWeightsNR, 0, 1);
+    fillShepardsIDW(mWeightsNB, 1, 0);
+    fillShepardsIDW(mWeightsC, 0, 0);
+  }
+  ~ShepardsIDW() {
+    delete[] mWeights;
+    delete[] mWeightsNR;
+    delete[] mWeightsNB;
+    delete[] mWeightsC;
+  }
+
+  int mMapScaleFactor;
+  // Image :-
+  // p00 p01 p02 p03 p04 p05 p06 p07
+  // p10 p11 p12 p13 p14 p15 p16 p17
+  // p20 p21 p22 p23 p24 p25 p26 p27
+  // p30 p31 p32 p33 p34 p35 p36 p37
+  // p40 p41 p42 p43 p44 p45 p46 p47
+  // p50 p51 p52 p53 p54 p55 p56 p57
+  // p60 p61 p62 p63 p64 p65 p66 p67
+  // p70 p71 p72 p73 p74 p75 p76 p77
+
+  // Recovery Map (for 4 scale factor) :-
+  // m00 p01
+  // m10 m11
+
+  // Recovery sample of curr 4x4, right 4x4, bottom 4x4, bottom right 4x4 are used during
+  // reconstruction. hence table weight size is 4.
+  float* mWeights;
+  // TODO: check if its ok to mWeights at places
+  float* mWeightsNR;  // no right
+  float* mWeightsNB;  // no bottom
+  float* mWeightsC;  // no right & bottom
+
+  float euclideanDistance(float x1, float x2, float y1, float y2);
+  void fillShepardsIDW(float *weights, int incR, int incB);
+};
+
 
 ////////////////////////////////////////////////////////////////////////////////
 // sRGB transformations
@@ -282,7 +328,9 @@
  * Sample the recovery value for the map from a given x,y coordinate on a scale
  * that is map scale factor larger than the map size.
  */
-float sampleMap(jr_uncompressed_ptr map, size_t map_scale_factor, size_t x, size_t y);
+float sampleMap(jr_uncompressed_ptr map, float map_scale_factor, size_t x, size_t y);
+float sampleMap(jr_uncompressed_ptr map, size_t map_scale_factor, size_t x, size_t y,
+                ShepardsIDW& weightTables);
 
 /*
  * Convert from Color to RGBA1010102.
diff --git a/libs/jpegrecoverymap/recoverymap.cpp b/libs/jpegrecoverymap/recoverymap.cpp
index ee68043..5224b57 100644
--- a/libs/jpegrecoverymap/recoverymap.cpp
+++ b/libs/jpegrecoverymap/recoverymap.cpp
@@ -773,14 +773,26 @@
       break;
   }
 
+  ShepardsIDW idwTable(kMapDimensionScaleFactor);
+
   for (size_t y = 0; y < height; ++y) {
     for (size_t x = 0; x < width; ++x) {
       Color yuv_gamma_sdr = getYuv420Pixel(uncompressed_yuv_420_image, x, y);
       Color rgb_gamma_sdr = srgbYuvToRgb(yuv_gamma_sdr);
       Color rgb_sdr = srgbInvOetf(rgb_gamma_sdr);
 
+      float recovery;
       // TODO: determine map scaling factor based on actual map dims
-      float recovery = sampleMap(uncompressed_recovery_map, kMapDimensionScaleFactor, x, y);
+      size_t map_scale_factor = kMapDimensionScaleFactor;
+      // TODO: If map_scale_factor is guaranteed to be an integer, then remove the following.
+      // Currently map_scale_factor is of type size_t, but it could be changed to a float
+      // later.
+      if (map_scale_factor != floorf(map_scale_factor)) {
+        recovery = sampleMap(uncompressed_recovery_map, map_scale_factor, x, y);
+      } else {
+        recovery = sampleMap(uncompressed_recovery_map, map_scale_factor, x, y,
+                             idwTable);
+      }
       Color rgb_hdr = applyRecovery(rgb_sdr, recovery, metadata->rangeScalingFactor);
 
       Color rgb_gamma_hdr = hdrOetf(rgb_hdr / metadata->rangeScalingFactor);
diff --git a/libs/jpegrecoverymap/recoverymapmath.cpp b/libs/jpegrecoverymap/recoverymapmath.cpp
index c2937e0..4a6feca 100644
--- a/libs/jpegrecoverymap/recoverymapmath.cpp
+++ b/libs/jpegrecoverymap/recoverymapmath.cpp
@@ -20,6 +20,53 @@
 
 namespace android::recoverymap {
 
+
+// Use Shepard's method for inverse distance weighting. For more information:
+// en.wikipedia.org/wiki/Inverse_distance_weighting#Shepard's_method
+
+float ShepardsIDW::euclideanDistance(float x1, float x2, float y1, float y2) {
+  return sqrt(((y2 - y1) * (y2 - y1)) + (x2 - x1) * (x2 - x1));
+}
+
+void ShepardsIDW::fillShepardsIDW(float *weights, int incR, int incB) {
+  for (int y = 0; y < mMapScaleFactor; y++) {
+    for (int x = 0; x < mMapScaleFactor; x++) {
+      float pos_x = ((float)x) / mMapScaleFactor;
+      float pos_y = ((float)y) / mMapScaleFactor;
+      int curr_x = floor(pos_x);
+      int curr_y = floor(pos_y);
+      int next_x = curr_x + incR;
+      int next_y = curr_y + incB;
+      float e1_distance = euclideanDistance(pos_x, curr_x, pos_y, curr_y);
+      int index = y * mMapScaleFactor * 4 + x * 4;
+      if (e1_distance == 0) {
+        weights[index++] = 1.f;
+        weights[index++] = 0.f;
+        weights[index++] = 0.f;
+        weights[index++] = 0.f;
+      } else {
+        float e1_weight = 1.f / e1_distance;
+
+        float e2_distance = euclideanDistance(pos_x, curr_x, pos_y, next_y);
+        float e2_weight = 1.f / e2_distance;
+
+        float e3_distance = euclideanDistance(pos_x, next_x, pos_y, curr_y);
+        float e3_weight = 1.f / e3_distance;
+
+        float e4_distance = euclideanDistance(pos_x, next_x, pos_y, next_y);
+        float e4_weight = 1.f / e4_distance;
+
+        float total_weight = e1_weight + e2_weight + e3_weight + e4_weight;
+
+        weights[index++] = e1_weight / total_weight;
+        weights[index++] = e2_weight / total_weight;
+        weights[index++] = e3_weight / total_weight;
+        weights[index++] = e4_weight / total_weight;
+      }
+    }
+  }
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 // sRGB transformations
 
@@ -372,6 +419,7 @@
   return sqrt(pow(x_diff, 2.0f) + pow(y_diff, 2.0f));
 }
 
+// TODO: If map_scale_factor is guaranteed to be an integer, then remove the following.
 float sampleMap(jr_uncompressed_ptr map, size_t map_scale_factor, size_t x, size_t y) {
   float x_map = static_cast<float>(x) / static_cast<float>(map_scale_factor);
   float y_map = static_cast<float>(y) / static_cast<float>(map_scale_factor);
@@ -421,6 +469,39 @@
        + e4 * (e4_weight / total_weight);
 }
 
+float sampleMap(jr_uncompressed_ptr map, size_t map_scale_factor, size_t x, size_t y,
+                ShepardsIDW& weightTables) {
+  // TODO: If map_scale_factor is guaranteed to be an integer power of 2, then optimize the
+  // following by computing log2(map_scale_factor) once and then using >> log2(map_scale_factor)
+  int x_lower = x / map_scale_factor;
+  int x_upper = x_lower + 1;
+  int y_lower = y / map_scale_factor;
+  int y_upper = y_lower + 1;
+
+  x_lower = std::min(x_lower, map->width - 1);
+  x_upper = std::min(x_upper, map->width - 1);
+  y_lower = std::min(y_lower, map->height - 1);
+  y_upper = std::min(y_upper, map->height - 1);
+
+  float e1 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_lower + y_lower * map->width]);
+  float e2 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_lower + y_upper * map->width]);
+  float e3 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_upper + y_lower * map->width]);
+  float e4 = mapUintToFloat(reinterpret_cast<uint8_t*>(map->data)[x_upper + y_upper * map->width]);
+
+  // TODO: If map_scale_factor is guaranteed to be an integer power of 2, then optimize the
+  // following by using & (map_scale_factor - 1)
+  int offset_x = x % map_scale_factor;
+  int offset_y = y % map_scale_factor;
+
+  float* weights = weightTables.mWeights;
+  if (x_lower == x_upper && y_lower == y_upper) weights = weightTables.mWeightsC;
+  else if (x_lower == x_upper) weights = weightTables.mWeightsNR;
+  else if (y_lower == y_upper) weights = weightTables.mWeightsNB;
+  weights += offset_y * map_scale_factor * 4 + offset_x * 4;
+
+  return e1 * weights[0] + e2 * weights[1] + e3 * weights[2] + e4 * weights[3];
+}
+
 uint32_t colorToRgba1010102(Color e_gamma) {
   return (0x3ff & static_cast<uint32_t>(e_gamma.r * 1023.0f))
        | ((0x3ff & static_cast<uint32_t>(e_gamma.g * 1023.0f)) << 10)
diff --git a/libs/jpegrecoverymap/tests/recoverymapmath_test.cpp b/libs/jpegrecoverymap/tests/recoverymapmath_test.cpp
index f8dd490..58c139d 100644
--- a/libs/jpegrecoverymap/tests/recoverymapmath_test.cpp
+++ b/libs/jpegrecoverymap/tests/recoverymapmath_test.cpp
@@ -699,6 +699,7 @@
   float (*values)[4] = MapValues();
 
   static const size_t kMapScaleFactor = 2;
+  ShepardsIDW idwTable(kMapScaleFactor);
   for (size_t y = 0; y < 4 * kMapScaleFactor; ++y) {
     for (size_t x = 0; x < 4 * kMapScaleFactor; ++x) {
       size_t x_base = x / kMapScaleFactor;
@@ -725,7 +726,7 @@
       // Instead of reimplementing the sampling algorithm, confirm that the
       // sample output is within the range of the min and max of the nearest
       // points.
-      EXPECT_THAT(sampleMap(&image, kMapScaleFactor, x, y),
+      EXPECT_THAT(sampleMap(&image, kMapScaleFactor, x, y, idwTable),
                   testing::AllOf(testing::Ge(min), testing::Le(max)));
     }
   }