int up = 0;  
      // Define two thresholds, one for edge detection and one for
      // interpolating point manipulation. Edge detection threshold should
      // be lower since it is mainly used for determining whether the edge
      // is horizontal/vertical or diagonal - this info is used in deciding
      // which grid to use. 
      
      Real th = 2.0;   // threshold for interpolating point manipulation
      Real thdet = 2.0; // threshold for edge detection
      int pm1 = -1;
      int p = -1;
      int pp1 = -1;
      int horiz = 0;  // horizontal edge
      int vert = 0;   // vertical edge
      Real HFVhoriz = 0;
      Real HFVvert = 0;
      
      // find if edge is vertical
      // first check if window is in the same row
      if (       
          ((int)(i/LLHsize) == (int)((i-4)/LLHsize)) && 
          ((int)(i/LLHsize) == (int)((i+2)/LLHsize)) &&  
          ((
            min(min(abv(coeff[1]->data[i-2]), abv(coeff[1]->data[i-1])),
                                           abv(coeff[1]->data[i])) >
            thdet*max(
                 max(
                     max(abv(coeff[1]->data[i-4]), abv(coeff[1]->data[i-3])),
                     abv(coeff[1]->data[i+1])
                    ), abv(coeff[1]->data[i+2])
                  )   
           ) ||
           ((
             min(abv(coeff[1]->data[i-2]), abv(coeff[1]->data[i])) >
             thdet*max(
                 max(
                     max(abv(coeff[1]->data[i-4]), abv(coeff[1]->data[i-3])),
                     abv(coeff[1]->data[i+1])
                    ), abv(coeff[1]->data[i+2])
                   )
            ) &&
             (abv(abv(coeff[1]->data[i-2])-abv(coeff[1]->data[i]))/abv(coeff[1]->data[i]) < 0.2)
         )))
           // coeffs found at i-2, i-1, i
          { 
            vert = 1;
            printf("vertical edge found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
          }
      else if 
         (
          ((int)(i/LLHsize) == (int)((i-3)/LLHsize)) && 
          ((int)(i/LLHsize) == (int)((i+3)/LLHsize)) &&
          ((
            min(min(abv(coeff[1]->data[i-1]), abv(coeff[1]->data[i])),
                                              abv(coeff[1]->data[i+1])) >
            thdet*max(
                 max(
                     max(abv(coeff[1]->data[i-3]), abv(coeff[1]->data[i-2])),
                     abv(coeff[1]->data[i+2])
                    ), abv(coeff[1]->data[i+3])
                  ) 
           ) ||
           ((
             min(abv(coeff[1]->data[i-1]), abv(coeff[1]->data[i+1])) >
             thdet*max(
                 max(
                     max(abv(coeff[1]->data[i-3]), abv(coeff[1]->data[i-2])),
                     abv(coeff[1]->data[i+2])
                    ), abv(coeff[1]->data[i+3])
                  )
            ) &&
             (abv(abv(coeff[1]->data[i-1])-abv(coeff[1]->data[i+1]))/abv(coeff[1]->data[i+1]) < 0.2)
         )))
           // coeffs found at i-1, i, i+1
           { 
             vert = 1;
             printf("vertical edge found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
           }
       // now check for HFV in horiz. dir.:  > 3 large coeffs
       else if 
         (
          ((int)(i/LLHsize) == (int)((i-4)/LLHsize)) && 
          ((int)(i/LLHsize) == (int)((i+3)/LLHsize)) &&
          (    //check for 4 large ones
           (
            min(
                min(
                    min(
                        abv(coeff[1]->data[i-1]), 
                        abv(coeff[1]->data[i])
                       ),
                    abv(coeff[1]->data[i+1])
                   ),
                abv(coeff[1]->data[i-2])
               ) >
            thdet*max(
                      max(
                          max(
                              abv(coeff[1]->data[i-4]), 
                              abv(coeff[1]->data[i-3])
                             ),
                          abv(coeff[1]->data[i+2])
                         ), 
                      abv(coeff[1]->data[i+3])
                     ) 
           ) ||   //check for other 4
           (
            min(
                min(
                    min(
                        abv(coeff[1]->data[i-1]), 
                        abv(coeff[1]->data[i])
                       ),
                    abv(coeff[1]->data[i+1])
                   ),
                abv(coeff[1]->data[i+2])
               ) >
            thdet*max(
                      max(
                          max(
                              abv(coeff[1]->data[i-4]), 
                              abv(coeff[1]->data[i-3])
                             ),
                          abv(coeff[1]->data[i-2])
                         ), 
                      abv(coeff[1]->data[i+3])
                     ) 
           ) || 
           (
            min(
                min(
                    min(
                        abv(coeff[1]->data[i-1]), 
                        abv(coeff[1]->data[i])
                       ),
                    abv(coeff[1]->data[i-3])
                   ),
                abv(coeff[1]->data[i-2])
               ) >
            thdet*max(
                      max(
                          max(
                              abv(coeff[1]->data[i-4]), 
                              abv(coeff[1]->data[i+1])
                             ),
                          abv(coeff[1]->data[i+2])
                         ), 
                      abv(coeff[1]->data[i+3])
                     ) 
           ) ||  // check for 5 large ones
           (
            min(
                min(
                    min(
                        min(abv(coeff[1]->data[i-1]), 
                            abv(coeff[1]->data[i-3])
                           ),
                        abv(coeff[1]->data[i])
                       ),
                    abv(coeff[1]->data[i+1])
                   ),
                abv(coeff[1]->data[i-2])
               ) >
            thdet*max(
                      max(
                          abv(coeff[1]->data[i-4]), 
                          abv(coeff[1]->data[i+2])
                         ), 
                      abv(coeff[1]->data[i+3])
                     ) 
           ) ||  // check for 5 large ones
           (
            min(
                min(
                    min(
                        min(abv(coeff[1]->data[i-1]), 
                            abv(coeff[1]->data[i+2])
                           ),
                        abv(coeff[1]->data[i])
                       ),
                    abv(coeff[1]->data[i+1])
                   ),
                abv(coeff[1]->data[i-2])
               ) >
            thdet*max(
                      max(
                          abv(coeff[1]->data[i-4]), 
                          abv(coeff[1]->data[i-3])
                         ), 
                      abv(coeff[1]->data[i+3])
                     ) 
           ) 
          )  
         )
           {
             HFVhoriz = 1;
             printf("horizontal HFV found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
           }
       else
           { 
             vert = 0;
           }

      // find parameters for interpolation point calculation
      // first check if window is in the same row
      if (       
          ((int)(i/LLHsize) == (int)((i-4)/LLHsize)) && 
          ((int)(i/LLHsize) == (int)((i+2)/LLHsize)) &&  
          ((
            min(min(abv(coeff[1]->data[i-2]), abv(coeff[1]->data[i-1])),
                                           abv(coeff[1]->data[i])) >
            th*max(
                 max(
                     max(abv(coeff[1]->data[i-4]), abv(coeff[1]->data[i-3])),
                     abv(coeff[1]->data[i+1])
                    ), abv(coeff[1]->data[i+2])
                  )   
           ) ||
           ((
             min(abv(coeff[1]->data[i-2]), abv(coeff[1]->data[i])) >
             th*max(
                 max(
                     max(abv(coeff[1]->data[i-4]), abv(coeff[1]->data[i-3])),
                     abv(coeff[1]->data[i+1])
                    ), abv(coeff[1]->data[i+2])
                   )
            ) &&
             (abv(abv(coeff[1]->data[i-2])-abv(coeff[1]->data[i]))/abv(coeff[1]->data[i]) < 0.2)
         )))
           // coeffs found at i-2, i-1, i
          { 
            pm1 = i-2; p = i-1; pp1 = i; 
          }
      else if 
         (
          ((int)(i/LLHsize) == (int)((i-3)/LLHsize)) && 
          ((int)(i/LLHsize) == (int)((i+3)/LLHsize)) &&
          ((
            min(min(abv(coeff[1]->data[i-1]), abv(coeff[1]->data[i])),
                                              abv(coeff[1]->data[i+1])) >
            th*max(
                 max(
                     max(abv(coeff[1]->data[i-3]), abv(coeff[1]->data[i-2])),
                     abv(coeff[1]->data[i+2])
                    ), abv(coeff[1]->data[i+3])
                  ) 
           ) ||
           ((
             min(abv(coeff[1]->data[i-1]), abv(coeff[1]->data[i+1])) >
             th*max(
                 max(
                     max(abv(coeff[1]->data[i-3]), abv(coeff[1]->data[i-2])),
                     abv(coeff[1]->data[i+2])
                    ), abv(coeff[1]->data[i+3])
                  )
            ) &&
             (abv(abv(coeff[1]->data[i-1])-abv(coeff[1]->data[i+1]))/abv(coeff[1]->data[i+1]) < 0.2)
         )))
           // coeffs found at i-1, i, i+1
           { 
             pm1 = i-1; p = i; pp1 = i+1; 
           }
       else
           { 
             pm1 = -1; p = -1; pp1 = -1; 
           }
       
       // find value of parameter m, then decide on interpolation point
       if (p == i-1)
       {
         // vertical "up" edge
         up = 1;
         x1 = 0.25; //?
       }
       
       else if (p == i)
       {
         // vertical "down" edge
         x1 = 0.75;  //?  
       }     
         
       up = 0;    
      
      
      // LH subband - HP filtered in vertical direction - horizontal edges

/*    // QMF-like concealment  
      k = i;
      while (k > LLHsize && k > i-2*LLHsize) k -= LLHsize;
      l = i;
      while (l <= LLsize && l < i + 2*LLHsize) l += LLHsize;
      
      m1 = k;
      for (j = k; j < l; j+=LLHsize)
        if (square(coeff[2]->data[j]) > thr*square(coeff[2]->data[m1]))
          m1 = j;
          
      if (m1 == k) m2 = k+LLHsize;
      else m2 = k; 
      for (j = k; j < l; j+=LLHsize)
        if (square(coeff[2]->data[j]) > thr*square(coeff[2]->data[m2])
            && j != m1)
          m2 = j;
      
      if (m1-LLHsize == m2 || m1+LLHsize == m2)
      {
        if (max(m1,m2) == i-LLHsize)  x2 = 0.5;
        else if (max(m1,m2) == i)  x2 = 0.75;
        else if (max(m1,m2) == i+LLHsize)  x2 = 0.5; //??
        //else if (max(m1,m2) == i+LLHsize)  
        //  x2 = 0.5*square((coeff[2]->data[m2])/(coeff[2]->data[m1]));
        else if (max(m1,m2) == i+2*LLHsize)  x2 = 0.25;
        else x2 = 0.5;  
      }      
*/          
      
      // Daubechies 9/7
      
      // find if edge is horizontal
      // first check if the window is entirely in the column
      if (        
          (i-4*LLHsize >= 0) && (i+2*LLHsize <= LLsize) &&
          ((
            min(min(abv(coeff[2]->data[i-2*LLHsize]), abv(coeff[2]->data[i-LLHsize])),
                                           abv(coeff[2]->data[i])) >
            thdet*max(
                 max(
                     max(abv(coeff[2]->data[i-4*LLHsize]), abv(coeff[2]->data[i-3*LLHsize])),
                     abv(coeff[2]->data[i+LLHsize])
                    ), abv(coeff[2]->data[i+2*LLHsize])
                  )
           ) ||
           ((
             min(abv(coeff[2]->data[i-2*LLHsize]), abv(coeff[2]->data[i])) >
             thdet*max(
                 max(
                     max(abv(coeff[2]->data[i-4*LLHsize]), abv(coeff[2]->data[i-3*LLHsize])),
                     abv(coeff[2]->data[i+LLHsize])
                    ), abv(coeff[2]->data[i+2*LLHsize])
                  )
            ) &&
             (abv(abv(coeff[2]->data[i-2*LLHsize])-abv(coeff[2]->data[i]))/abv(coeff[2]->data[i]) < 0.2)
         )))
           // coeffs found at i-2*LLHsize, i-LLHsize, i
           { 
             horiz = 1;
             printf("horizontal edge found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
           }
      else if 
         (
          (i-3*LLHsize >= 0) && (i+3*LLHsize <= LLsize) &&
          ((
            min(min(abv(coeff[2]->data[i-LLHsize]), abv(coeff[2]->data[i])),
                                              abv(coeff[2]->data[i+LLHsize])) >
            thdet*max(
                 max(
                     max(abv(coeff[2]->data[i-3*LLHsize]), abv(coeff[2]->data[i-2*LLHsize])),
                     abv(coeff[2]->data[i+2*LLHsize])
                    ), abv(coeff[2]->data[i+3*LLHsize])
                  )
           ) ||
           ((
             min(abv(coeff[2]->data[i-LLHsize]), abv(coeff[2]->data[i+LLHsize])) >
             thdet*max(
                 max(
                     max(abv(coeff[2]->data[i-3*LLHsize]), abv(coeff[2]->data[i-2*LLHsize])),
                     abv(coeff[2]->data[i+2*LLHsize])
                    ), abv(coeff[2]->data[i+3*LLHsize])
                  )
            ) &&
             (abv(abv(coeff[2]->data[i-LLHsize])-abv(coeff[2]->data[i+LLHsize]))/abv(coeff[2]->data[i+LLHsize]) < 0.2)
         )))
           // coeffs found at i-LLHsize, i, i+LLHsize
          { 
            horiz = 1;
            printf("horizontal edge found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
          }
      
      else if    
        (
          (i-4*LLHsize >= 0) && (i+3*LLHsize <= LLsize) && 
          ((
            min(
                min(
                    min(
                        abv(coeff[2]->data[i-LLHsize]), 
                        abv(coeff[2]->data[i])
                       ),
                    abv(coeff[2]->data[i+LLHsize])
                   ),
                abv(coeff[2]->data[i-2*LLHsize])
               ) >
            thdet*max(
                      max(
                          max(
                              0*abv(coeff[2]->data[i-4*LLHsize]), 
                              abv(coeff[2]->data[i-3*LLHsize])
                             ),
                          abv(coeff[2]->data[i+2*LLHsize])
                         ), 
                      0*abv(coeff[2]->data[i+3*LLHsize])
                     ) 
           ) ||
           (
            min(
                min(
                    min(
                        abv(coeff[2]->data[i-LLHsize]), 
                        abv(coeff[2]->data[i])
                       ),
                    abv(coeff[2]->data[i+LLHsize])
                   ),
                abv(coeff[2]->data[i+2*LLHsize])
               ) >
            thdet*max(
                      max(
                          max(
                              0*abv(coeff[2]->data[i-4*LLHsize]), 
                              abv(coeff[2]->data[i-3*LLHsize])
                             ),
                          abv(coeff[2]->data[i-2*LLHsize])
                         ), 
                      0*abv(coeff[2]->data[i+3*LLHsize])
                     ) 
           ) || 
           (
            min(
                min(
                    min(
                        abv(coeff[2]->data[i-LLHsize]), 
                        abv(coeff[2]->data[i])
                       ),
                    abv(coeff[2]->data[i-3*LLHsize])
                   ),
                abv(coeff[2]->data[i-2*LLHsize])
               ) >
            thdet*max(
                      max(
                          max(
                              0*abv(coeff[2]->data[i-4*LLHsize]), 
                              abv(coeff[2]->data[i+LLHsize])
                             ),
                          abv(coeff[2]->data[i+2*LLHsize])
                         ), 
                      0*abv(coeff[2]->data[i+3*LLHsize])
                     ) 
           ) ||  // check for 5 large ones
           (
            min(
                min(
                    min(
                        min(abv(coeff[2]->data[i-LLHsize]), 
                            abv(coeff[2]->data[i-3*LLHsize])
                           ),
                        abv(coeff[2]->data[i])
                       ),
                    abv(coeff[2]->data[i+LLHsize])
                   ),
                abv(coeff[2]->data[i-2])
               ) >
            thdet*max(
                      max(
                          abv(coeff[2]->data[i-4*LLHsize]), 
                          abv(coeff[2]->data[i+2*LLHsize])
                         ), 
                      abv(coeff[2]->data[i+3*LLHsize])
                     ) 
           ) ||  // check for 5 large ones
           (
            min(
                min(
                    min(
                        min(abv(coeff[2]->data[i-LLHsize]), 
                            abv(coeff[2]->data[i+2*LLHsize])
                           ),
                        abv(coeff[2]->data[i])
                       ),
                    abv(coeff[2]->data[i+LLHsize])
                   ),
                abv(coeff[2]->data[i-2*LLHsize])
               ) >
            thdet*max(
                      max(
                          abv(coeff[2]->data[i-4*LLHsize]), 
                          abv(coeff[2]->data[i-3*LLHsize])
                         ), 
                      abv(coeff[2]->data[i+3*LLHsize])
                     ) 
           )
          )  
         )
           {
             HFVvert = 1; 
             printf("vertical HFV found at (%d,%d)\n",(int)i/LLHsize+1,i%LLHsize+1);
           }
       else
          { 
            horiz = 0;
          }
      
      
      // find parameters for interpolation point calculation
      // first check if the window is entirely in the column
      if (        
          (i-4*LLHsize >= 0) && (i+2*LLHsize <= LLsize) &&
          ((
            min(min(abv(coeff[2]->data[i-2*LLHsize]), abv(coeff[2]->data[i-LLHsize])),
                                           abv(coeff[2]->data[i])) >
            th*max(
                 max(
                     max(abv(coeff[2]->data[i-4*LLHsize]), abv(coeff[2]->data[i-3*LLHsize])),
                     abv(coeff[2]->data[i+LLHsize])
                    ), abv(coeff[2]->data[i+2*LLHsize])
                  )
           ) ||
           ((
             min(abv(coeff[2]->data[i-2*LLHsize]), abv(coeff[2]->data[i])) >
             th*max(
                 max(
                     max(abv(coeff[2]->data[i-4*LLHsize]), abv(coeff[2]->data[i-3*LLHsize])),
                     abv(coeff[2]->data[i+LLHsize])
                    ), abv(coeff[2]->data[i+2*LLHsize])
                  )
            ) &&
             (abv(abv(coeff[2]->data[i-2*LLHsize])-abv(coeff[2]->data[i]))/abv(coeff[2]->data[i]) < 0.2)
         )))
           // coeffs found at i-2*LLHsize, i-LLHsize, i
           { 
             pm1 = i-2*LLHsize; p = i-LLHsize; pp1 = i; 
           }
      else if 
         (
          (i-3*LLHsize >= 0) && (i+3*LLHsize <= LLsize) &&
          ((
            min(min(abv(coeff[2]->data[i-LLHsize]), abv(coeff[2]->data[i])),
                                              abv(coeff[2]->data[i+LLHsize])) >
            th*max(
                 max(
                     max(abv(coeff[2]->data[i-3*LLHsize]), abv(coeff[2]->data[i-2*LLHsize])),
                     abv(coeff[2]->data[i+2*LLHsize])
                    ), abv(coeff[2]->data[i+3*LLHsize])
                  )
           ) ||
           ((
             min(abv(coeff[2]->data[i-LLHsize]), abv(coeff[2]->data[i+LLHsize])) >
             th*max(
                 max(
                     max(abv(coeff[2]->data[i-3*LLHsize]), abv(coeff[2]->data[i-2*LLHsize])),
                     abv(coeff[2]->data[i+2*LLHsize])
                    ), abv(coeff[2]->data[i+3*LLHsize])
                  )
            ) &&
             (abv(abv(coeff[2]->data[i-LLHsize])-abv(coeff[2]->data[i+LLHsize]))/abv(coeff[2]->data[i+LLHsize]) < 0.2)
         )))
           // coeffs found at i-LLHsize, i, i+LLHsize
          { 
            pm1 = i-LLHsize; p = i; pp1 = i+LLHsize; 
          }
       else
          { 
            pm1 = -1; p = -1; pp1 = -1; 
          }
       
       // find value of parameter m, then decide on interpolation point
       if (p == i-LLHsize)
       {
         // horizontal "up" edge
         up = 1;
         x2 = 0.25;  //?
       }
       
       else if (p == i)
       {
         // horizontal "down" edge
         x2 = 0.75;  //?
       }     
      
