Real LL[LLVsize+4][LLHsize+4]; // LL subband symmetrically extended
int row, col;

for (i=0; i<coeff[0]->nData; i++)
	LL[floor(i/LLHsize)+2][i%LLHsize+2] = coeff[0]->data[i];

// symmetric extension for corners
LL[0][0] = LL[3][3]; LL[0][1] = LL[2][3];
LL[1][0] = LL[3][2]; LL[1][1] = LL[2][2];

LL[0][LLHsize+2] = LL[2][LLHsize];
LL[0][LLHsize+3] = LL[3][LLHsize];
LL[1][LLHsize+2] = LL[2][LLHsize+1];
LL[1][LLHsize+3] = LL[3][LLHsize+1];

LL[LLVsize+2][0] = LL[LLVsize][2];
LL[LLVsize+3][0] = LL[LLVsize][3];
LL[LLVsize+2][1] = LL[LLVsize+1][2];
LL[LLVsize+3][1] = LL[LLVsize+1][3];

LL[LLVsize+2][LLHsize+2] = LL[LLVsize+1][LLHsize+1];
LL[LLVsize+2][LLHsize+3] = LL[LLVsize][LLHsize+1];
LL[LLVsize+3][LLHsize+2] = LL[LLVsize+1][LLHsize];
LL[LLVsize+3][LLHsize+3] = LL[LLVsize][LLHsize];

// symmetric extension for boundaries
for (i=0; i<LLVsize+4; i++)
{
  for (j=0; j<2; j++)
   // ...
}

// calculate quantities for bicubic interpolation
for (i=0; i<coeff[0]->nData; i++)
	if (missing[coeff[0]->packet(0,i)-1])
 	{
		// missing coefficient - process its four neighbors
		row = floor(i/LLHsize);
		col = i%LLHsize;
		// upper left
		if (row>1 && col>1)
		{
			y[3] = LL[row-1][col-1];
			// vertical derivative at [row-1][col-1]
			y1[3] = (LL[row][col-1] - LL[row-2][col-1])/2;
			// horizontal derivative at [row-1][col-1]
			y2[3] = (LL[row-1][col] - LL[row-1][col-2])/2;
			// cross derivative at [row-1][col-1]
			y12[3] = (LL[row][col] - LL[row][col-1] - LL[row-2][col] +
					      LL[row-2][col-2])/4;
		}
		// ...
	}
			
