subr subpixelcellc(s1, s2, xoff, yoff, mask)
 ;no support for apodizing mask yet
 nx = dimen(s1,0)	ny = dimen(s1,1)
 ;construct the merit function for a range of sub-pixel shifts
 ;the centered case where we shift both arrays by equal amounts
 s1 = doub(s1)
 s2 = doub(s2)
 x0 = s1(0:(nx-2),0:(ny-2))
 x1 = s1(1:*,0:(ny-2))
 x2 = s1(0:(nx-2),1:*)
 x3 = s1(1:*,1:*)
 y0 = s2(0:(nx-2),0:(ny-2))
 y1 = s2(1:*,0:(ny-2))
 y2 = s2(0:(nx-2),1:*)
 y3 = s2(1:*,1:*)

 a1 = x0 + x1 + x2 + x3 - y0 - y1 - y2 - y3 
 a2 = x0 - x1 + x2 - x3 + y0 - y1 + y2 - y3 
 a3 = x0 + x1 - x2 - x3 + y0 + y1 - y2 - y3 
 a4 = x0 - x1 - x2 + x3 - y0 + y1 + y2 - y3
 
 ;the ten expanded terms
 
 t1 = .0625 * dsumproduct(a1,a1)
 t2 = .25 * dsumproduct(a2,a2)
 t3 = .25 * dsumproduct(a3,a3)
 t4 = dsumproduct(a4,a4)
 t5 = .25 * dsumproduct(a1,a2)
 t6 = .25 * dsumproduct(a1,a3)
 t7 = .50 * dsumproduct(a1,a4)
 t8 = .50 * dsumproduct(a2,a3)
 t9 = dsumproduct(a2,a4)
 t10 = dsumproduct(a3,a4)
 

 lockc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,xmin,ymin
 
 xoff = 2.*xmin
 yoff = 2.*ymin
 
 endsubr 
 ;============================================================================