      // Daubechies 9/7
      // Look for three neighboring coeffs with magnitudes "large" wrt
      // neighbors, or two "large" ones with similar magnitudes of opposite
      // sign separated by 1 point. "Similar" here means difference in 
      // magnitude is < 20%
      // looking for a vertical edge
      
      int up = 0;  
      Real th = 1.0; // threshold for edge detection - consider varying it
      int pm1 = 0;
      int p = 0;
      int pp1 = 0;
      k = i;
      while (k%LLHsize > 0 && k > i-4) k--;  // go left 4 points if possible
      l = i;
      while (l%LLHsize == (l+1)%LLHsize && l < i+3) l++; // go right
      if (
          (k < i-2) && (k+6 <= l) &&
          ((
            min(min(abv(coeff[1]->data[k+2]), abv(coeff[1]->data[k+3])),
                                           abv(coeff[1]->data[k+4])) >
            th*max(
                 max(
                     max(abv(coeff[1]->data[k]), abv(coeff[1]->data[k+1])),
                     abv(coeff[1]->data[k+5])
                    ), abv(coeff[1]->data[k+6])
                  )   
           ) ||
           ((
             min(abv(coeff[1]->data[k+2]), abv(coeff[1]->data[k+4])) >
             th*max(
                 max(
                     max(abv(coeff[1]->data[k]), abv(coeff[1]->data[k+1])),
                     abv(coeff[1]->data[k+5])
                    ), abv(coeff[1]->data[k+6])
                   )
            ) &&
             (abv(abv(coeff[1]->data[k+2])-abv(coeff[1]->data[k+4]))/abv(coeff[1]->data[k+2]) < 0.2)
         )))
           // coeffs found at k+2, k+3, k+4
          { pm1 = k+2; p = k+3; pp1 = k+4; }
      else if 
         (
          (k < i-2) && (k+7 <= l) &&
          ((
            min(min(abv(coeff[1]->data[k+3]), abv(coeff[1]->data[k+4])),
                                              abv(coeff[1]->data[k+5])) >
            th*max(
                 max(
                     max(abv(coeff[1]->data[k+1]), abv(coeff[1]->data[k+2])),
                     abv(coeff[1]->data[k+6])
                    ), abv(coeff[1]->data[k+7])
                  ) 
           ) ||
           ((
             min(abv(coeff[1]->data[k+3]), abv(coeff[1]->data[k+5])) >
             th*max(
                 max(
                     max(abv(coeff[1]->data[k+1]), abv(coeff[1]->data[k+2])),
                     abv(coeff[1]->data[k+6])
                    ), abv(coeff[1]->data[k+7])
                  )
            ) &&
             (abv(abv(coeff[1]->data[k+3])-abv(coeff[1]->data[k+5]))/abv(coeff[1]->data[k+3]) < 0.2)
         )))
           // coeffs found at k+3, k+4, k+5
           { pm1 = k+3; p = k+4; pp1 = k+5; }
       else
           { pm1 = 0; p = 0; pp1 = 0; }
       
       // find value of parameter m, then decide on interpolation point
       if (p == i-1)
       {
         // vertical "up" edge
         up = 1;
         x1 = (coeff[1]->data[p])/(abv(coeff[1]->data[pm1])+abv(coeff[1]->data[pp1]));
         x1 = -x1/4.4603;  // -1 <= x1 <= 1
         x1 = (x1+1.0)/2.0; // 0 <= x1 <= 1
         x1 = x1/4.0 + 0.75; // 0.75 <= x1 <= 1
         if (x1 > 1)
           x1 = 1.0;
         else if (x1 < 0.75)
           x1 = 0.75;
       }
       
       if (p == i)
       {
         // vertical "down" edge
         Real x12;
         x12 = (coeff[1]->data[p])/(abv(coeff[1]->data[pm1])+abv(coeff[1]->data[pp1]));
         x12 = x12/4.4603;  // -1 <= x12 <= 1
         x12 = (x12+1.0)/2.0; // 0 <= x12 <= 1
         x12 = x12/4.0; // 0 <= x12 <= 0.25
         if (x12 > 0.25)
           x12 = 0.25;
         else if (x12 < 0)
           x12 = 0;
         if (up == 0)
           x1 = x12;
         else if (abv(x12-0.5) < abv(x1-0.5))  // if both "up" and "down"
           x1 = x12;
       }     
         
       up = 0;    
            
      // LH subband - HP filtered in vertical direction?

      // Daubechies 9/7
      k = i;
      while (k > LLHsize && k > i-4*LLHsize) k -= LLHsize; // go up 4 points
      l = i;
      while (l <= LLsize && l < i + 3*LLHsize) l += LLHsize; // go down 3 pts

      if (
          (k < i-2*LLHsize) && (k+6*LLHsize <= l) &&
          ((
            min(min(abv(coeff[2]->data[k+2*LLHsize]), abv(coeff[2]->data[k+3*LLHsize])),
                                           abv(coeff[2]->data[k+4*LLHsize])) >
            th*max(
                 max(
                     max(abv(coeff[2]->data[k]), abv(coeff[2]->data[k+LLHsize])),
                     abv(coeff[2]->data[k+5*LLHsize])
                    ), abv(coeff[2]->data[k+6*LLHsize])
                  )
           ) ||
           ((
             min(abv(coeff[2]->data[k+2*LLHsize]), abv(coeff[2]->data[k+4*LLHsize])) >
             th*max(
                 max(
                     max(abv(coeff[2]->data[k]), abv(coeff[2]->data[k+LLHsize])),
                     abv(coeff[2]->data[k+5*LLHsize])
                    ), abv(coeff[2]->data[k+6*LLHsize])
                  )
            ) &&
             (abv(abv(coeff[2]->data[k+2*LLHsize])-abv(coeff[2]->data[k+4*LLHsize]))/abv(coeff[2]->data[k+2*LLHsize]) < 0.2)
         )))
           // coeffs found at k+2*LLHsize, k+3*LLHsize, k+4*LLHsize
           { pm1 = k+2*LLHsize; p = k+3*LLHsize; pp1 = k+4*LLHsize; }
      else if 
         (
          (k < i-2*LLHsize) && (k+7*LLHsize <= l) &&
          ((
            min(min(abv(coeff[2]->data[k+3*LLHsize]), abv(coeff[2]->data[k+4*LLHsize])),
                                              abv(coeff[2]->data[k+5*LLHsize])) >
            th*max(
                 max(
                     max(abv(coeff[2]->data[k+LLHsize]), abv(coeff[2]->data[k+2*LLHsize])),
                     abv(coeff[2]->data[k+6*LLHsize])
                    ), abv(coeff[2]->data[k+7*LLHsize])
                  )
           ) ||
           ((
             min(abv(coeff[2]->data[k+3*LLHsize]), abv(coeff[2]->data[k+5*LLHsize])) >
             th*max(
                 max(
                     max(abv(coeff[2]->data[k+LLHsize]), abv(coeff[2]->data[k+2*LLHsize])),
                     abv(coeff[2]->data[k+6*LLHsize])
                    ), abv(coeff[2]->data[k+7*LLHsize])
                  )
            ) &&
             (abv(abv(coeff[2]->data[k+3*LLHsize])-abv(coeff[2]->data[k+5*LLHsize]))/abv(coeff[2]->data[k+3*LLHsize]) < 0.2)
         )))
           // coeffs found at k+3*LLHsize, k+4*LLHsize, k+5*LLHsize
          { pm1 = k+3*LLHsize; p = k+4*LLHsize; pp1 = k+5*LLHsize; }
       else
          { pm1 = 0; p = 0; pp1 = 0; }
       
       // find value of parameter m, then decide on interpolation point
       if (p == i-LLHsize)
       {
         // horizontal "up" edge
         up = 1;
         x2 = (coeff[2]->data[p])/(abv(coeff[2]->data[pm1])+abv(coeff[2]->data[pp1]));
         x2 = x2/4.4603;  // -1 <= x2 <= 1
         x2 = (x2+1.0)/2.0; // 0 <= x2 <= 1
         x2 = x2/4.0;  // 0.75 <= x2 <= 1
         if (x2 > 0.25)
           x2 = 0.25;
         else if (x2 < 0)
           x2 = 0;
       }
       
       if (p == i)
       {
         // horizontal "down" edge
         Real x22;
         x22 = (coeff[2]->data[p])/(abv(coeff[2]->data[pm1])+abv(coeff[2]->data[pp1]));
         x22 = x22/4.4603;  // -1 <= x22 <= 1
         x22 = (x22+1.0)/2.0; // 0 <= x22 <= 1
         x22 = x22/4.0 + 0.75; // 0.75 <= x22 <= 1
         if (x22 > 1)
           x22 = 1;
         else if (x22 < 0.75)
           x22 = 0.75;
         if (up == 0)
           x2 = x22;
         else if (abv(x22-0.5) < abv(x2-0.5))  // if both "up" and "down"
           x2 = x22;
       }     
      
      if ((x1 < 0.5) || (x1 > 0.5)) 
         printf("vertical edge found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
      if ((x2 < 0.5) || (x2 > 0.5)) 
         printf("horizontal edge found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
      
      //if ((fabs(coeff[1]->data[i+1]) > TH || fabs(coeff[2]->data[i+1]) > TH))
          // if (diagonal) edge - bicubic interpolation at (x1, x2)
      
      // if the coefficient is not near the image boundary, then interpolate
      if  (
            ((x1 < 0.5) || (x1 > 0.5) || (x2 < 0.5) || (x2 > 0.5))   
            && (i%LLHsize > 1) && (i%LLHsize < LLHsize-2) 
            && (i > 2*LLHsize) && (i+2*LLHsize < LLsize)
          )
          coeff[0]->data[i] = bicubInt(c, x1, x2);
          