gt:       for (int i2=0; i2<n2; ++i2) {
context: /****************************************************************************
  Copyright (c) 2014, Colorado School of Mines and others. All rights reserved.
  This program and accompanying materials are made available under the terms of
  the Common Public License - v1.0, which accompanies this distribution, and is 
  available at http://www.eclipse.org/legal/cpl-v10.html
  ****************************************************************************/
  package ipf;
  
  import java.util.*;
  
  import edu.mines.jtk.awt.*;
  import edu.mines.jtk.dsp.*;
  import edu.mines.jtk.interp.*;
  import edu.mines.jtk.util.*;
  
  import static edu.mines.jtk.util.ArrayMath.*;
  import static ipf.FaultGeometry.*;
  
  /**
   * Uses image samples alongside fault skins to estimate fault dip-slips.
   * <p>
   * Before estimating slips, the seismic images should be smoothed along
   * reflectors, but not across faults. (Thinned fault likelihoods can be used
   * to stop smoothing at potential fault locations.) This smoothing does what
   * seismic interpreters do visually when estimating fault throws. In effect,
   * such smoothing enables us to use samples of a seismic image located away
   * from faults to estimate slips along those faults.
   * <p>
   * However, even after this smoothing, we may not want to use image samples
   * located immediately adjacent to faults, as these samples may exhibit poor
   * estimates of reflector slopes and smoothing filter artifacts. Therefore,
   * image samples at some small horizontal distance (e.g., two samples away)
   * from fault skins are used to estimate slips. Specified reflector slopes are
   * then used to account for this small offset when computing fault slips.
   *
   * @author Dave Hale, Colorado School of Mines
   * @version 2014.07.05
   */
  public class FaultSlipper {
  
    /**
     * Constructs a fault slipper for the specified seismic image and slopes.
     * @param gs seismic image, smoothed up to (but not across) faults.
     * @param p2 slopes in the 2nd dimension.
     * @param p3 slopes in the 3rd dimension.
     */
    public FaultSlipper(float[][][] gs, float[][][] p2, float[][][] p3) {
      _gs = gs;
      _p2 = p2;
      _p3 = p3;
      _n1 = gs[0][0].length;
      _n2 = gs[0].length;
      _n3 = gs.length;
    }
  
    /**
     * Sets the offset distance used when estimating slips.
     * The default offset is 2 samples.
     * @param offset the offset, in samples.
     */
    public void setOffset(double offset) {
      _offset = (float)abs(offset);
    }
  
    /**
     * Enables or disables a zero-slope reflector assumption. For educational
     * use, only. Ignoring reflector slopes can yield significant errors in
     * estimated dip slips. The default is false.
     * @param zeroSlope true, to assume zero slopes; false, otherwise.
     */
    public void setZeroSlope(boolean zeroSlope) {
      _zeroSlope = zeroSlope;
    }
  
    /**
     * Computes fault dip slips for the specified skins.
     * Specified bounds on throw are used to reduce unnecessary computation.
     * These bounds need not be precise, and are not enforced. If the fault
     * skins represent reverse faults (with negative throws), the specified
     * estimate for minimum fault throw should be negative.
     * @param skins array of skins for which to compute dip slips.
     * @param smin an estimate for minimum fault throw, in samples.
     * @param smax an estimate for maximum fault throw, in samples.
     */
    public void computeDipSlips(FaultSkin[] skins, double smin, double smax) {
      for (FaultSkin skin:skins)
        computeDipSlips(skin,smin,smax);
    }
  
    /**
     * Computes fault dip slips for all cells in the specified skin.
     * Specified bounds on throw are used to reduce unnecessary computation.
     * These bounds need not be precise, and are not enforced. If the fault
     * skins represent reverse faults (with negative throws), the specified
     * estimate for minimum fault throw should be negative.
     * @param skin the skin for which to compute shifts.
     * @param smin an estimate for minimum fault throw, in samples.
     * @param smax an estimate for maximum fault throw, in samples.
     */
    public void computeDipSlips(FaultSkin skin, double smin, double smax) {
      Check.argument(smax>=0.0f,"smax not less than zero");
      FaultCell[][] cab = skin.getCellsAB();
      FaultCell[][] clr = skin.getCellsLR();
      int lmin = round((float)smin-2*_offset); // because shifts != throws
      int lmax = round((float)smax+2*_offset); // TODO: use pmax*_offset?
      DynamicWarping dw = new DynamicWarping(lmin,lmax);
      dw.setStrainMax(0.25,0.25); // TODO: always 0.25? goes with 4 below?
      computeAlignmentErrors(skin,lmin,lmax,_offset,_gs);
      extrapolateAlignmentErrors(lmin,lmax,cab);
      computeShifts(dw,cab,clr);
      clearErrors(skin);
      for (int nsmooth=0; nsmooth<2; ++nsmooth) // TODO: 2?
        smoothShifts(skin);
      computeDipSlips(skin);
    }
  
    /**
     * Gets images of fault dip-slip vectors for an array of fault skins.
     * Assumes that dip-slip vectors have already been computed and stored
     * in all skin cells.
     * <p>
     * Dip-slip vectors are non-zero for only the hanging wall (plus side)
     * of faults. Slips on the footwall (minus side) are zero. Image samples
     * not adjacent to faults are marked using a specified value.
     * @param skins array of fault skins.
     * @param smark the mark for slips not adjacent to a fault.
     * @return array {s1,s2,s3} of components of dip slips.
     */
    public float[][][][] getDipSlips(FaultSkin[] skins, float smark) {
      int n1 = _n1, n2 = _n2, n3 = _n3;
      float[][][] s1 = new float[n3][n2][n1];
      float[][][] s2 = new float[n3][n2][n1];
      float[][][] s3 = new float[n3][n2][n1];
      float[][][] ss = new float[n3][n2][n1];
  
      // Initially set all slip vectors to the specified mark.
      for (int i3=0; i3<n3; ++i3) {
        for (int i2=0; i2<n2; ++i2) {
          for (int i1=0; i1<n1; ++i1) {
            s1[i3][i2][i1] = smark;
            s2[i3][i2][i1] = smark;
            s3[i3][i2][i1] = smark;
          }
        }
      }
  
      // For all cells in all skins, ...
      for (FaultSkin skin:skins) {
        for (FaultCell cell:skin) {
  
          // Get sample indices for the minus and plus sides of the cell.
          int i1 = cell.i1;
          int i2m = cell.i2m;
          int i3m = cell.i3m;
          int i2p = cell.i2p;
          int i3p = cell.i3p;
  
          // If slip on the minus side has not been set, zero it.
          if (s1[i3m][i2m][i1]==smark) {
            s1[i3m][i2m][i1] = 0.0f;
            s2[i3m][i2m][i1] = 0.0f;
            s3[i3m][i2m][i1] = 0.0f;
          }
  
          // Set or accumulate slip on the plus side.
          if (s1[i3p][i2p][i1]==smark) {
            s1[i3p][i2p][i1]  = cell.s1;
            s2[i3p][i2p][i1]  = cell.s2;
            s3[i3p][i2p][i1]  = cell.s3;
            ss[i3p][i2p][i1]  = 1.0f;
          } else {
            s1[i3p][i2p][i1] += cell.s1;
            s2[i3p][i2p][i1] += cell.s2;
            s3[i3p][i2p][i1] += cell.s3;
            ss[i3p][i2p][i1] += 1.0f;
          }
        }
      }
  
      // Where more than one slip was accumulated, compute the average.
      for (int i3=0; i3<n3; ++i3) {
        for (int i2=0; i2<n2; ++i2) {
          for (int i1=0; i1<n1; ++i1) {
            if (ss[i3][i2][i1]>1.0f) {
              float si = 1.0f/ss[i3][i2][i1];
              s1[i3][i2][i1] *= si;
              s2[i3][i2][i1] *= si;
              s3[i3][i2][i1] *= si;
            }
          }
        }
      }
      return new float[][][][]{s1,s2,s3};
    }
  
    /**
     * Interpolates specified dip-slip vectors.
     * @param s array {s1,s2,s3} of dip-slip vectors.
     * @param smark the mark for slips not adjacent to a fault.
     * @return interpolated dip-slip vectors.
     */
    public static float[][][][] interpolateDipSlips(
        float[][][][] s, float smark) {
      int n1 = s[0][0][0].length;
      int n2 = s[0][0].length;
      int n3 = s[0].length;
      float[][][] sp = new float[n3][n2][n1];
      float[][][] st = new float[n3][n2][n1];
      short[][][] k1 = new short[n3][n2][n1];
      short[][][] k2 = new short[n3][n2][n1];
      short[][][] k3 = new short[n3][n2][n1];
      float[][][][] sq = new float[3][n3][n2][n1];
      ClosestPointTransform cpt = new ClosestPointTransform();
      cpt.apply(smark,s[0],st,k1,k2,k3);
      clip(0.0f,100.0f,st,st);
      LocalDiffusionKernel.Stencil stencil = LocalDiffusionKernel.Stencil.D21;
      LocalDiffusionKernel ldk = new LocalDiffusionKernel(stencil);
      BlendedGridder3 bg = new BlendedGridder3();
      bg.setBlendingKernel(ldk);
      bg.setSmoothness(0.5);
      for (int is=0; is<3; ++is) {
        float[][][] si = s[is];
        for (int i3=0; i3<n3; ++i3) {
          for (int i2=0; i2<n2; ++i2) {
            for (int i1=0; i1<n1; ++i1) {
              int j1 = k1[i3][i2][i1];
              int j2 = k2[i3][i2][i1];
              int j3 = k3[i3][i2][i1];
              sp[i3][i2][i1] = si[j3][j2][j1];
            }
          }
        }
        bg.gridBlended(st,sp,sq[is]);
      }
      return sq;
    }
  
    /**
     * Unfaults an image using interpolated dip-slip vectors.
     * @param s array {s1,s2,s3} of interpolated dip-slip vectors.
     * @param g image to be unfaulted.
     * @return unfaulted image.
     */
    public static float[][][] unfault(float[][][][] s, final float[][][] g) {
      final int n1 = g[0][0].length;
      final int n2 = g[0].length;
      final int n3 = g.length;
      final float[][][] s1 = s[0];
      final float[][][] s2 = s[1];
      final float[][][] s3 = s[2];
      final float[][][] gs = new float[n3][n2][n1];
      final SincInterpolator si = new SincInterpolator();
      Parallel.loop(n3,new Parallel.LoopInt() {
      public void compute(int i3) {
  