docs: void constructLinearSystem(Planar<I> input , Se3_F32 g ) {
  		int numBands = imageType.getNumBands();
  
  		// first precompute everything that does not depend on pixel values
  		inboundsPixels = 0;
  		for (int i = 0; i < keypixels.size(); i++) {
  			Pixel p = keypixels.data[i];
  
  			// Apply the known warp
  			SePointOps_F32.transform(g, p.p3, S);
  
  			if( S.z <= 0 ) {
  				p.valid = false;
  				continue;
  			}
  
  			// Compute projected warped pixel coordinate on image I_1
  			p.proj.x = (S.x / S.z) * fx + cx;
  			p.proj.y = (S.y / S.z) * fy + cy;
  
  			// make sure it's in the bounds
  			if (p.proj.x < 0 || p.proj.x > input.width - 1 || p.proj.y < 0 || p.proj.y > input.height - 1) {
  				p.valid = false;
  				continue;
  			} else {
  				p.valid = true;
  			}
  			inboundsPixels++;
  
  			// pi matrix derivative relative to t at S
  			float ZZ = S.z * S.z;
  
  			p.dP11 = fx / S.z;
  			p.dP13 = -S.x * fx / ZZ;
  			p.dP22 = fy / S.z;
  			p.dP23 = -S.y * fy / ZZ;
  		}
  
  		// how compute the components which require
  		errorOptical = 0;
  		int row = 0;
  		for (int band = 0; band < numBands; band++) {
  			interpDX.setImage(derivX.getBand(band));
  			interpDY.setImage(derivY.getBand(band));
  			interpI.setImage(input.getBand(band));
  
  			for (int i = 0; i < keypixels.size(); i++) {
  				Pixel p = keypixels.data[i];
  
  				if( !p.valid )
  					continue;
  
  				// Apply the known warp
  				// TODO precompute?
  				SePointOps_F32.transform(g, p.p3, S);
  
  				// sample pixel values at warped location in I_1
  				// NOTE: This could be highly optimized.  Compute and save interpolation weights once per input
  				//       instead of for each band and image (current,dx,dy)
  				// TODO create a special bilinear class for this?
  				float current = interpI.get( p.proj.x, p.proj.y);
  				float dx      = interpDX.get(p.proj.x, p.proj.y);
  				float dy      = interpDY.get(p.proj.x, p.proj.y);
  
  				// B = grad^T * dPI/dt = shape(1,3)
  				float b1 = dx*p.dP11;
  				float b2 = dy*p.dP22;
  				float b3 = dx*p.dP13 + dy*p.dP23;
  
  				// C * A(S'(x)) = shape(1,6)
  				int indexA = row * 6;
  				A.data[indexA++] = -b2*S.z + b3*S.y;
  				A.data[indexA++] =  b1*S.z - b3*S.x;
  				A.data[indexA++] = -b1*S.y + b2*S.x;
  				A.data[indexA++] = b1;
  				A.data[indexA++] = b2;
  				A.data[indexA]   = b3;
  
  				float error = -(current - p.bands[band]);
  				y.data[row] = error;
  
  				errorOptical += Math.abs(error);
  				row += 1;
  			}
  		}
  		errorOptical /= row;
  
  		A.numRows = row;
  		y.numRows = row;
  	}