Program Splt_CD5(Input, Output) ;
{ Split CD5 files into something Star can process }
{ Input has already been processed with ProcCD5 }
{ These files have 2032 words per N1 chunk }
{ The output files need N1 768 + fill = 800 to match MKIII }
{ The other way goes all the way to 2032 }
{ So it needs 3 images to cover }
{ No! Of course nothing is a multiple of FITS record 1440 words }
Uses
   WinAPI, { No - it's not Windows. It's BP7.0 32-bit stuff }
   X_U, { My routines for 32-bit addressing }
   FITSX_U, { Random assortment of FITS read/write stuff }
   SORT_U ; { Quicksort }
Const
{ Update the first 4 and header items in main program }
  InDir   = 'C:\BP\TASS\CD5\IMAGES\' ;
  LogName = 'C:\BP\TASS\CD5\SPLT_CD5.TXT' ;
  OutDir  = 'C:\BP\TASS\CD5\' ;
{ Bright source purging tests pixel for value > CritVal            }
{ It then finds local median MedVal and replaces all neighbouring  }
{ pixels within the purge box (see below) above MedVal + MinOff    }
{ with MedVal }
  MinOff  = 200 ; { << saturation: > noise }
  CritVal = 5000 ; { -> 37768; was 8000 -> 40768; This is crude!  }
{ Images with clouds can have median values near, or even above,  }
{ "saturation" which muck up the purging process. Throw 'em away! }
  MaxVal  = 0 ; { -> 32768 Largest allowed median }
{ Mapping }
  NOut1 = 800 ; { Old MkIII used for output }
  NOut2 = 2032 ;
  N1_1 = 14 ; { Matches MkIII }
  N1_2 = 782 ;
  N1_Off1 =  -14 ; { giving    0.. 767 }
  N1_Off2 =  618 ; { giving  632..1399 with an overlap of 136 }
  N1_Off3 = 1250 ; { giving 1264..2031 }
{ Medians }
  Box = 40 ; { half side of median box }
  MStep = 50 ; { Computation interval for medians }
  NMed = (NOut2+MStep-1) Div MStep ;
  DepthLim = 5000 ; { Maximum recursion depth for ProcessPoint }
  Extra = 5000 ; { Overflow space for ProcessPoint }
Type
  ExtraT = Array[0..Extra-1] of Integer ;
  ExtraP = ^ExtraT ;
Var
  Data, Sel : FITSInfoT ; { Input, output control structures     }
  Log : Text ; { File for output list of "saturated" sources etc }
  N1_M, N2_M : Longint ; { Actual sizes of input images }
  MedVal, MinVal : Integer ; { Values for purging (see above) }
  MedDat : Array[0..NMed-1,0..NMed-1] of Integer ;
  PurgeCt : Longint ;
  Depth : Integer ; { Recursion depth for ProcessPoint }
  MaxDepth : Integer ;
  XP1, XP2 : ExtraP ;
  ExtraCt : Integer ;
  MaxExtra : Integer ;
  A : Longint ; { Accumulated pixel sum for bright sources }
  Sum1, Sum2 : Extended ; { Accumulated sums for positions }
  StarNo : Integer ; { Running count of purged bright sources }

Procedure Purge(X1, X2 : Longint) ;
{ Accumulate statistics of purged pixel; replace pixel with MedVal }
Var
  T : Longint ;
  ET : Extended ;
Begin
{ The Longint typecasts force the compiler to do what it ought }
{ to have been doing anyway ... }
  T := Longint(IPtr(Data.Image, X2*N1_M+X1)^) - Longint(MedVal) ;
{ This mess is subverting BP7.0 into doing:                    }
{ Inc(A, Image[X2, X1])   properly, with 32-bit addressing     }
  A := A + T ; { Split over two lines to help debugging! }
  ET := T ; { Gets round compiler integer problem }
{ No! I NEVER said Borland Pascal 7.0 was the best way to do this }
  Sum1 := Sum1 + ET*X1 ; { Sum for Axis1 position }
  Sum2 := Sum2 + ET*X2 ; { Sum for Axis2 position }
  IPtr(Data.Image, X2*N1_M+X1)^ := MedVal ; { Purge }
  Inc(PurgeCt) ;
End ; { Purge }

Function Med(I1, I2 : Integer) : Integer ;
{ Get median value in box }
Const
  Factor = 0.5 ; { 0.5 for median but may be biased }
  Step1 = 4 ; { Arbitrary decimation factors }
  Step2 = 4 ;
  NMax = ((2*Box+Step1) Div Step1)*((2*Box+Step2) Div Step2) ;
Var
  X1, X2 : Longint ;
  Count : Integer ;
  L, Index : Array[0..NMax-1] of Integer ;
  Min1, Max1, Min2, Max2 : Integer ; { Box boundaries }
Begin
{ Set up box for Med }
  Min1 := I1 - Box ; If (Min1 < 0) Then Min1 := 0 ;
  Max1 := I1 + Box ; If (Max1 >= N1_M) Then Max1 := N1_M - 1 ;
  Min2 := I2 - Box ; If (Min2 < 0) then Min2 := 0 ;
  Max2 := I2 + Box ; If (Max2 >= N2_M) Then Max2 := N2_M - 1 ;
  X2 := Min2 ; Count := 0 ;
  Repeat X1 := Min1 ;
    Repeat
      L[Count] := IPtr(Data.Image, X2*N1_M+X1)^ ;
      Inc(Count) ; Inc(X1, Step1) ;
    Until (X1 >= Max1) ;
    Inc(X2, Step2) ;
  Until (X2 >= Max2) ;
  ISort(L, 0, Count-1, Index) ;
  Med := L[Index[Trunc(Factor*Count)]] ;
End ; {Med }

Procedure ProcessPoint(X1, X2 : Integer) ; Forward ;
Procedure TryPoint(X1, X2 : Integer) ;
{ = ProcessPoint but for saved Extra points }
Begin
  MedVal := MedDat[X2 Div MStep, X1 Div MStep] ;
  MinVal := MedVal + MinOff ;
  Inc(Depth) ; { Count recursions. No need to check overflow here }
{ Check all 4 directions but not diagonals }
  If (X1+1 < N1_M) And
  (IPtr(Data.Image, Longint(X2)*N1_M+Longint(X1)+1)^ > MinVal) Then
    ProcessPoint(X1+1, X2) ;
  If (X2+1 < N2_M) And
  (IPtr(Data.Image, (Longint(X2)+1)*N1_M+Longint(X1))^ > MinVal) Then
    ProcessPoint(X1, X2+1) ;
  If (X1-1 >= 0) And
  (IPtr(Data.Image, Longint(X2)*N1_M+Longint(X1)-1)^ > MinVal) Then
    ProcessPoint(X1-1, X2) ;
  If (X2-1 >= 0) And
  (IPtr(Data.Image, (Longint(X2)-1)*N1_M+Longint(X1))^ > MinVal) Then
    ProcessPoint(X1, X2-1) ;
  Dec(Depth) ;
End ; { TryPoint }

Procedure ProcessPoint(X1, X2 : Integer) ;
{ Iterative routine to seek and destroy pixels above MinVal }
{ around a bright source }
Begin
  MedVal := MedDat[X2 Div MStep, X1 Div MStep] ;
  MinVal := MedVal + MinOff ;
  Purge(X1, X2) ;
{ Check depth of recursion }
  If (Depth >= DepthLim) Then Begin
{ Off the end; save the parameters }
    If (ExtraCt >= Extra) Then Begin
{ Not reached ... so far! }
      Writeln('You blew it  ') ; Readln ; Halt ;
    End Else Begin
      XP1^[ExtraCt] := X1 ; XP2^[ExtraCt] := X2 ;
      Inc(ExtraCt) ;
      If (ExtraCt > MaxExtra) Then MaxExtra := ExtraCt ;
    End ;
  End Else Begin
    Inc(Depth) ; { Count recursions }
    If (Depth > MaxDepth) Then MaxDepth := Depth ;
{ Check all 4 directions but not diagonals }
    If (X1+1 < N1_M) And
    (IPtr(Data.Image, Longint(X2)*N1_M+Longint(X1)+1)^ > MinVal) Then
      ProcessPoint(X1+1, X2) ;
    If (X2+1 < N2_M) And
    (IPtr(Data.Image, (Longint(X2)+1)*N1_M+Longint(X1))^ > MinVal) Then
      ProcessPoint(X1, X2+1) ;
    If (X1-1 >= 0) And
    (IPtr(Data.Image, Longint(X2)*N1_M+Longint(X1)-1)^ > MinVal) Then
      ProcessPoint(X1-1, X2) ;
    If (X2-1 >= 0) And
    (IPtr(Data.Image, (Longint(X2)-1)*N1_M+Longint(X1))^ > MinVal) Then
      ProcessPoint(X1, X2-1) ;
    Dec(Depth) ;
    While ((Depth <= 0) and (ExtraCt > 0)) Do Begin
      Dec(ExtraCt) ;
      TryPoint(XP1^[ExtraCt], XP2^[ExtraCt]) ;
    End ;
  End ;
End ; { ProcessPoint }

Procedure Process(InName, OutName : String) ;
{ Process one image, removing bright sources and splitting }
Var
  I2, I1 : Longint ;
  X1, X2 : Integer ;
  BigFlag : Boolean ; { Flag for median too big }
Begin
  StarNo := 0 ;
  Append(Log) ;
  Writeln(Log, '"SPLT_CD5","', InName, '"') ;
  Close(Log) ; { Another fix for BP7.0; file gets written each time      }
{ that it is closed so if something blows up one at least knows how      }
{ far one had got. The file is reopened and reclosed every line. This    }
{ costs time and makes the disk sound pretty demented but it's worth it! }
  Write('Reading data ', InName, '  ') ;
  ReadSome(Data, InDir + InName, 0, N1_M, 0, N2_M) ;
  SetHead(Sel, NOut1, N2_M) ; { Used to allocate space; also sets up header }
  Move(Data.Head^, Sel.Head^, LRecB) ; { Overwrite the new header with    }
{ the one read in with the image. Most of it can be kept ... more or less }
  SetHeadCardI(Sel.Head^.Hdr[3], 'NAXIS1', NOut1) ;
  Write('Medians  ') ;
{ Set up array of medians }
  BigFlag := False ;
  I2 := MStep Div 2 ; X2 := 0 ;
  Repeat
    I1 := MStep Div 2 ; X1 := 0 ;
    Repeat
      MedDat[X2, X1] := Med(I1, I2) ; ;
      If (MedDat[X2, X1] + MinOff > MaxVal) Then BigFlag := True ;
      Inc(I1, MStep) ; Inc(X1) ;
    Until (I1 >= N1_M) ;
    Inc(I2, MStep) ; Inc(X2) ;
  Until (I2 >= N2_M) ;
{ Actual purge loop }
  Writeln('Purging') ;
  If BigFlag Then Begin
    Append(Log) ;
    Writeln(Log, '"Background too high"') ;
    Close(Log) ;
  End Else Begin
    For I2 := 0 To N2_M-1 Do Begin
      For I1 := 0 To N1_M-1  Do
        If (IPtr(Data.Image, I2*N1_M+I1)^ > Critval) Then Begin
        Inc(StarNo) ;
        PurgeCt := 0 ;
        Depth := 0 ; ExtraCt := 0 ;
        MaxDepth := 0 ; MaxExtra := 0 ;
{ Initialise bright source statistics }
        A := 0 ; Sum1 := 0; Sum2 := 0 ;
        ProcessPoint(I1, I2) ;
        Append(Log) ;
        Writeln(Log, StarNo, ',', Sum2/A:9:2, ',', Sum1/A:9:2, ',',
          A, ',', PurgeCt, ',', MaxDepth, ',', MaxExtra) ;
        Close(Log) ;
      End ; { If }
    End ; { For }
{ Copy to Sel }
    Write('Writing strip 1  ') ;
    For I2 := 0 To N2_M-1 Do Begin
      For I1 :=    0 To N1_1-1  Do IPtr(Sel.Image, I2*NOut1+I1)^ := 0 ;
      For I1 := N1_1 To N1_2-1  Do
        IPtr(Sel.Image, I2*NOut1+I1)^ := IPtr(Data.Image, I2*N1_M+I1+N1_Off1)^ ;
      For I1 := N1_2 To NOut1-1 Do IPtr(Sel.Image, I2*NOut1+I1)^ := 0 ;
    End ;
    OutName[7] := '1' ;
    WriteFITS(Sel, OutDir + OutName) ;
    Write('2  ') ;
    For I2 := 0 To N2_M-1 Do Begin
      For I1 :=    0 To N1_1-1  Do IPtr(Sel.Image, I2*NOut1+I1)^ := 0 ;
      For I1 := N1_1 To N1_2-1  Do
        IPtr(Sel.Image, I2*NOut1+I1)^ := IPtr(Data.Image, I2*N1_M+I1+N1_Off2)^ ;
      For I1 := N1_2 To NOut1-1 Do IPtr(Sel.Image, I2*NOut1+I1)^ := 0 ;
    End ;
    OutName[7] := '2' ;
    WriteFITS(Sel, OutDir + OutName) ;
    Writeln('3') ;
    For I2 := 0 To N2_M-1 Do Begin
      For I1 :=    0 To N1_1-1  Do IPtr(Sel.Image, I2*NOut1+I1)^ := 0 ;
      For I1 := N1_1 To N1_2-1  Do
        IPtr(Sel.Image, I2*NOut1+I1)^ := IPtr(Data.Image, I2*N1_M+I1+N1_Off3)^ ;
      For I1 := N1_2 To NOut1-1 Do IPtr(Sel.Image, I2*NOut1+I1)^ := 0 ;
    End ;
    OutName[7] := '3' ;
    WriteFITS(Sel, OutDir + OutName) ;
  End ; { If BigFlag }
End ; { Process }

Begin
  Assign(Log, LogName) ;
  Rewrite(Log) ;
  Writeln(Log, '"Splt_CD5"') ;
  Close(Log) ;
  New(XP1) ; New(XP2) ; { Arrays to hold purging overflow }
{ Compute sizes }
  N2_M := 2032 ;
  N1_M := 2032 ;
  Process('H3R1438.878', 'P1__38_.878') ;
  Process('H3R1438.882', 'P1__38_.882') ;
  Process('H3R1438.886', 'P1__38_.886') ;
  Process('H3R1438.890', 'P1__38_.890') ;
  Process('H3R1438.895', 'P1__38_.895') ;
  Process('H3R1438.899', 'P1__38_.899') ;
  Process('H3R1438.903', 'P1__38_.903') ;
  Process('H3R1439.870', 'P1__39_.870') ;
  Process('H3R1439.875', 'P1__39_.875') ;
  Process('H3R1439.880', 'P1__39_.880') ;
  Process('H3R1439.885', 'P1__39_.885') ;
  Process('H3R1439.891', 'P1__39_.891') ;
  Process('H3R1439.896', 'P1__39_.896') ;
  Process('H3R1439.901', 'P1__39_.901') ;
  Process('H3R1440.871', 'P1__40_.871') ;
  Process('H3R1440.875', 'P1__40_.875') ;
  Process('H3R1440.879', 'P1__40_.879') ;
  Process('H3R1440.884', 'P1__40_.884') ;
  Process('H3R1440.888', 'P1__40_.888') ;
  Process('H3R1440.892', 'P1__40_.892') ;
  Process('H3R1440.897', 'P1__40_.897') ;
  Process('H3R1442.863', 'P1__42_.863') ;
  Process('H3R1442.867', 'P1__42_.867') ;
  Process('H3R1442.871', 'P1__42_.871') ;
  Process('H3R1442.876', 'P1__42_.876') ;
  Process('H3R1442.880', 'P1__42_.880') ;
  Process('H3R1442.885', 'P1__42_.885') ;
  Process('H3R1442.889', 'P1__42_.889') ;
  Process('H3R1446.853', 'P1__46_.853') ;
  Process('H3R1446.858', 'P1__46_.858') ;
  Process('H3R1446.862', 'P1__46_.862') ;
  Process('H3R1446.866', 'P1__46_.866') ;
  Process('H3R1446.871', 'P1__46_.871') ;
  Process('H3R1446.875', 'P1__46_.875') ;
  Process('H3R1446.880', 'P1__46_.880') ;
  Process('H4R1438.878', 'P2__38_.878') ;
  Process('H4R1438.882', 'P2__38_.882') ;
  Process('H4R1438.886', 'P2__38_.886') ;
  Process('H4R1438.890', 'P2__38_.890') ;
  Process('H4R1438.895', 'P2__38_.895') ;
  Process('H4R1438.899', 'P2__38_.899') ;
  Process('H4R1438.903', 'P2__38_.903') ;
  Process('H4R1439.870', 'P2__39_.870') ;
  Process('H4R1439.875', 'P2__39_.875') ;
  Process('H4R1439.880', 'P2__39_.880') ;
  Process('H4R1439.885', 'P2__39_.885') ;
  Process('H4R1439.891', 'P2__39_.891') ;
  Process('H4R1439.896', 'P2__39_.896') ;
  Process('H4R1439.901', 'P2__39_.901') ;
  Process('H4R1440.871', 'P2__40_.871') ;
  Process('H4R1440.875', 'P2__40_.875') ;
  Process('H4R1440.879', 'P2__40_.879') ;
  Process('H4R1440.884', 'P2__40_.884') ;
  Process('H4R1440.888', 'P2__40_.888') ;
  Process('H4R1440.892', 'P2__40_.892') ;
  Process('H4R1440.897', 'P2__40_.897') ;
  Process('H4R1442.863', 'P2__42_.863') ;
  Process('H4R1442.867', 'P2__42_.867') ;
  Process('H4R1442.871', 'P2__42_.871') ;
  Process('H4R1442.876', 'P2__42_.876') ;
  Process('H4R1442.880', 'P2__42_.880') ;
  Process('H4R1442.885', 'P2__42_.885') ;
  Process('H4R1442.889', 'P2__42_.889') ;
  Process('H4R1446.853', 'P2__46_.853') ;
  Process('H4R1446.858', 'P2__46_.858') ;
  Process('H4R1446.862', 'P2__46_.862') ;
  Process('H4R1446.866', 'P2__46_.866') ;
  Process('H4R1446.871', 'P2__46_.871') ;
  Process('H4R1446.875', 'P2__46_.875') ;
  Process('H4R1446.880', 'P2__46_.880') ;
End. { Splt_CD5 }