unit ParticleDetectionUnit;
// -----------------------
// Object detection module
// -----------------------
// 27.05.04
// 27.07.04 No. of objects limit increased to 1000/frame
// 04.09.04 Pre-existing form now closed automatically
//          Object size threshold now in pixel units

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ValidatedEdit, StdCtrls, RangeEdit, ImageFile ;

type
  TParticleDetectionFrm = class(TForm)
    MainGrp: TGroupBox;
    Label2: TLabel;
    cbSourceFile: TComboBox;
    RangeGrp: TGroupBox;
    rbAllFrames: TRadioButton;
    rbRange: TRadioButton;
    edRange: TRangeEdit;
    bOK: TButton;
    bCancel: TButton;
    GroupBox3: TGroupBox;
    edThresholdLo: TValidatedEdit;
    Label1: TLabel;
    Label3: TLabel;
    edThresholdHi: TValidatedEdit;
    SaveFile: TImageFile;
    GroupBox4: TGroupBox;
    edObjectSizeThreshold: TValidatedEdit;
    GroupBox2: TGroupBox;
    rbZProjection: TRadioButton;
    RadioButton1: TRadioButton;
    procedure FormShow(Sender: TObject);
    procedure bOKClick(Sender: TObject);
    procedure cbSourceFileChange(Sender: TObject);
    procedure bCancelClick(Sender: TObject);
  private
    { Private declarations }
    procedure NewFile ;
  public
    { Public declarations }
  end;

var
  ParticleDetectionFrm: TParticleDetectionFrm;

implementation

uses PicViewMain, ViewUnit, maths ;

{$R *.dfm}

procedure TParticleDetectionFrm.FormShow(Sender: TObject);
// -------------------------------------
// Initialisation when form is displayed
// -------------------------------------
var
     i : Integer ;
     Src : Integer ;
begin

     // Get list of available image files
     cbSourceFile.Clear ;
     for i :=  0 to High(MainFrm.ViewFrmsInUse) do begin
         if MainFrm.ViewFrmsInUse[i] then begin
            cbSourceFile.Items.AddObject( ExtractFileName(MainFrm.ViewFrms[i].FileName),
                                          TObject(i)) ;
            end ;
         end ;
     cbSourceFile.ItemIndex := 0 ;

     ClientWidth := MainGrp.Left + MainGrp.Width + 10 ;
     ClientHeight := bOK.Top + bOK.Height + 10 ;

     NewFile ;

     end ;


procedure TParticleDetectionFrm.NewFile ;
// ----------------------
// File selection changed
// ----------------------
var
    Src : Integer ;
begin

     Src := Integer(cbSourceFile.Items.Objects[cbSourceFile.ItemIndex]) ;

     if MainFrm.ViewFrms[Src].NumSectionsPerStack > 1 then begin
        // 4D (stack series) mode
        RangeGrp.Caption := ' Stack Range ' ;
        rbAllFrames.Caption := 'All Stacks' ;
        edRange.HiLimit := MainFrm.ViewFrms[Src].NumStacks ;
        edRange.HiValue := edRange.HiLimit ;
        // Object size threshold units
        edObjectSizeThreshold.Units := MainFrm.ViewFrms[Src].PixelUnits + '^3';
        edObjectSizeThreshold.Scale := MainFrm.ViewFrms[Src].XResolution*
                                       MainFrm.ViewFrms[Src].YResolution*
                                       MainFrm.ViewFrms[Src].ZResolution ;
        end
     else begin
        // 3D (frame series) mode
        RangeGrp.Caption := ' Frame Range ' ;
        rbAllFrames.Caption := 'All Frames' ;
        edRange.HiLimit := MainFrm.ViewFrms[Src].NumFrames ;
        edRange.HiValue := edRange.HiLimit ;

        // Object size threshold units
        edObjectSizeThreshold.Units := MainFrm.ViewFrms[Src].PixelUnits + '^2';
        edObjectSizeThreshold.Scale := MainFrm.ViewFrms[Src].XResolution*
                                       MainFrm.ViewFrms[Src].YResolution ;

        end ;

     edThresholdHi.HiLimit := MainFrm.ViewFrms[cbSourceFile.ItemIndex].GreyMax ;
     edThresholdHi.Value := MainFrm.ViewFrms[cbSourceFile.ItemIndex].GreyMax ;
     edThresholdLo.HiLimit := edThresholdHi.HiLimit ;
     edThresholdLo.Value := 0.5*edThresholdHi.Value ;

     end;


procedure TParticleDetectionFrm.bOKClick(Sender: TObject);
// -------------------------------------------------------
// Create a particle count image from selected source file
// -------------------------------------------------------
const
     MaxRuns = 4096*4096 ;
type
     TRunElement = record
        x : Word ;
        y : Word ;
        z : Word ;
        Count : Word ;
        ObjectID : Word ;
        end ;
     TRuns = Array[0..MaxRuns-1] of TRunElement ;
     PRuns = ^TRuns ;
var
     SrcFileName : String ;  // Source file name
     ParticleFileName : String ; // Projection file name
     Src : Integer ;         // Source image viewfrm index number
     NumComponentsPerPixel : Cardinal ; // No. of colour components/pixel
     NumComponentsPerFrame : Cardinal ; // No. of pixels in source image frame
     FrameWidth : Cardinal ;
     FrameHeight : Cardinal ;
     NumSectionsPerStack : Integer ;    // No. Z sections per stack
     i,j : Integer ;           // General counter
     FrameNum : Integer ;    // Source frame counter
     OutFrame : Integer ;    // Output frame number
     SectionNum : Integer ;  // Z section counter
     StartSection : Integer ;  // Start section
     EndSection : Integer ;  // End section included in Z proj
     StartStack : Integer ;  // Start stack #
     EndStack : Integer ;    // End stack #
     Stack : Integer ;
     ThresholdLo : Integer ;
     ThresholdHi : Integer ;
     OK : Boolean ;

     PSrcBuf : PIntArray ;       // Source frame buffer pointer
     POutBuf : Pointer ;         // O/P frame buffer pointer

     Runs : PRuns ;              // Pointer to runs array
     MaxRunsPerFrame : Integer ;         // Max. no. of runs allowed
     NumRuns : Integer ;         // No. of runs in Runs
     iRun : Integer ;            // Run index
     jRun : Integer ;            // Run index
     iStart : Integer ;          // Start element of run
     iEnd : Integer ;
     jStart : Integer ;
     jEnd  : Integer ;
     ObjectID : Word ;
     OldID : Word ;
     NumObjects : Word ;
     NumPixels : Integer ;
     RunActive : Boolean ;       // Run active flag
     x : Word ;                  // x (horizontal) pixel position
     y : Word ;                  // y (vertical) pixel position
     z : Word ;                  // z (Depth in stack) pixel position

     ResultsFrmNum : Integer ;   // Index # of results form
     ResultsCol : Integer ;      // Results table column counter
     ResultsRow : Integer ;      // Results table row counter
     StackVar : Integer ;
     TimeVar : Integer ;
     ObjectIDVar : Integer ;
     XVar : Integer ;
     YVar : Integer ;
     ZVar : Integer ;
     VolVar : Integer ;
     NumRecords : Integer ;
     SumX : Single ;
     SumY : Single ;
     SumZ : Single ;
begin


     // Source of stack to be Z projected
     Src := Integer(cbSourceFile.Items.Objects[cbSourceFile.ItemIndex]) ;
     SrcFileName := MainFrm.ViewFrms[Src].FileName ;
     NumComponentsPerPixel := MainFrm.ViewFrms[Src].NumComponentsPerPixel ;
     FrameWidth := MainFrm.ViewFrms[Src].FrameWidth ;
     FrameHeight := MainFrm.ViewFrms[Src].FrameHeight ;
     NumComponentsPerFrame := MainFrm.ViewFrms[Src].FrameWidth
                              *MainFrm.ViewFrms[Src].FrameHeight
                              *NumComponentsPerPixel ;
     NumSectionsPerStack := MainFrm.ViewFrms[Src].NumSectionsPerStack ;

     // Get range of stacks/frames to include
     StartSection := 1 ;
     EndSection := NumSectionsPerStack ;
     if rbAllFrames.Checked then begin
        StartStack := 1 ;
        EndStack := Round( edRange.HiLimit ) ;
        end
     else begin
        StartStack := Round( edRange.LoValue ) ;
        EndStack := Round( edRange.HiValue ) ;
        end ;

     ThresholdLo := Round( edThresholdLo.Value ) ;
     ThresholdHi := Round( edThresholdHi.Value ) ;

     // Create file name
     ParticleFileName := '' ;
     i := 1 ;
     While (SrcFileName[i] <> '.') and (i <= Length(SrcFileName)) do begin
         ParticleFileName := ParticleFileName + SrcFileName[i] ;
         Inc(i) ;
         end ;
     ParticleFileName := ParticleFileName + format('[Objects %d-%d].pic',
                                            [StartStack,EndStack]) ;

     // Let user quit if merge file already exists
     if FileExists(ParticleFileName) then begin
        if (MessageDlg( 'Particle Detection: ' + ParticleFileName + ' already exists! Overwrite it?',
           mtWarning, [mbYes,mbNo], 0 ) = mrNo) then Exit ;
        end ;

     // Close form (if output file is on display)
     MainFrm.CloseViewFrm(ParticleFileName);

     // Create binary particle file
     OK := SaveFile.CreateFile( ParticleFileName,
                                MainFrm.ViewFrms[Src].FrameWidth,
                                MainFrm.ViewFrms[Src].FrameHeight,
                                16,
                                1,
                                True ) ;
     SaveFile.XResolution := MainFrm.ViewFrms[Src].XResolution ;
     SaveFile.YResolution := MainFrm.ViewFrms[Src].YResolution ;
     SaveFile.ZResolution := MainFrm.ViewFrms[Src].ZResolution ;
     SaveFile.TResolution := MainFrm.ViewFrms[Src].TResolution ;
     SaveFile.ResolutionUnit := MainFrm.ViewFrms[Src].PixelUnits ;

     // Allocate buffers
     GetMem( PSrcBuf, NumComponentsPerFrame*4 ) ;
     GetMem( POutBuf, NumComponentsPerFrame*2 ) ;
     MaxRunsPerFrame := (FrameWidth*FrameHeight) div 2 ;
     GetMem( Runs, MaxRunsPerFrame*SizeOf(TRunElement)) ;

     // Prevent multiple activation
     bOK.Enabled := False ;
     Application.ProcessMessages ;

     // Create object measurement results table
     ResultsFrmNum := MainFrm.CreateNewResultsFrm( ChangeFileExt( ParticleFileName, '.res' ));

     StackVar := MainFrm.ResultsFrms[ResultsFrmNum].AddVariable( 'Timepoint' ) ;
     TimeVar := MainFrm.ResultsFrms[ResultsFrmNum].AddVariable( 'Time (s)' ) ;
     ObjectIDVar := MainFrm.ResultsFrms[ResultsFrmNum].AddVariable( 'Tracking ID' ) ;
     XVar := MainFrm.ResultsFrms[ResultsFrmNum].AddVariable( 'Centroid X (µm)' ) ;
     YVar := MainFrm.ResultsFrms[ResultsFrmNum].AddVariable( 'Centroid Y (µm)' ) ;
     ZVar := MainFrm.ResultsFrms[ResultsFrmNum].AddVariable( 'Centroid Z (µm)' ) ;
     VolVar := MainFrm.ResultsFrms[ResultsFrmNum].AddVariable( 'Volume (µm3)' ) ;
     MainFrm.ResultsFrms[ResultsFrmNum].MaxRecords := (EndStack - StartStack + 1)*1000 ;
     NumRecords := 0 ;

     try

     // Process all stacks in file
     OutFrame := 0 ;
     for Stack := StartStack to EndStack do begin

        // Process selected frames from stack
        NumRuns := 0 ;
        for SectionNum := StartSection to EndSection do begin

            // Report progress
            MainFrm.StatusBar.SimpleText := format(
            'Detect Objects: Section %.4d/%.4d of Stack %.4d to %s (%d objects)',
            [SectionNum,EndSection,Stack,ExtractFileName(ParticleFileName),NumRecords] ) ;

            // Load frame from source file
            FrameNum := (Stack-1)*NumSectionsPerStack + SectionNum ;
            MainFrm.ViewFrms[Src].LoadFrame( FrameNum, PSrcBuf ) ;

            // Create binary run-length compressed image of threshold selected points

            for y := 0 to FrameHeight-1 do begin
                RunActive := False ;
                for x := 0 to FrameWidth-1 do begin
                    i := x + y*FrameWidth ;
                    if (PSrcBuf^[i] >= ThresholdLo) and
                       (PSrcBuf^[i] <= ThresholdHi) then begin
                       if not RunActive then begin
                          Inc(NumRuns) ;
                          Runs^[NumRuns-1].x := x ;
                          Runs^[NumRuns-1].y := y ;
                          Runs^[NumRuns-1].z := SectionNum ;
                          Runs^[NumRuns-1].Count := 0 ;
                          Runs^[NumRuns-1].ObjectID := NumRuns ;
                          RunActive := True ;
                          end ;
                       Inc(Runs^[NumRuns-1].Count) ;
                       end
                    else RunActive := False ;
                    end ;
                end ;

            if bOK.Enabled then Break ;
            Application.ProcessMessages ;

            end ;

        // Report progress
        MainFrm.StatusBar.SimpleText := format(
        'Detect Objects: Detecting objects in Stack %.4d to %s (%d objects)',
        [Stack,ExtractFileName(ParticleFileName),NumRecords] ) ;

        // Join together neighboring runs into objects
        for iRun := 0 to NumRuns-1 do begin
            iStart := Runs^[iRun].x ;
            y := Runs^[iRun].y ;
            z := Runs^[iRun].z ;
            iEnd := Runs^[iRun].x + Runs^[iRun].Count - 1 ;
            ObjectID := Runs^[iRun].ObjectID ;
            for jRun := 0 to iRun-1 do begin
                if ((Runs^[jRun].y = (y-1)) and (Runs^[iRun].z = z)) or
                   ((Runs^[jRun].y = y) and (Runs^[iRun].z = (z-1))) then begin
                   jStart := Runs^[jRun].x ;
                   jEnd := Runs^[jRun].x + Runs^[jRun].Count - 1 ;
                   if ((jStart >= iStart) and (jStart <= iEnd)) or
                      ((jEnd >= iStart) and (jEnd <= iEnd)) or
                      ((iStart >= jStart) and (iStart <= jEnd)) or
                      ((iEnd >= jStart) and (iEnd <= jEnd)) then begin
                      OldID := Runs^[jRun].ObjectID ;
                      for i := 0 to iRun do if Runs^[i].ObjectID = OldID then
                          Runs^[i].ObjectID := ObjectID ;
                      end ;
                   end ;
                end ;

            if bOK.Enabled then Break ;
            Application.ProcessMessages ;

            end ;

        // Report progress
        MainFrm.StatusBar.SimpleText := format(
        'Detect Objects: Removing small objects in Stack %.4d to %s (%d objects)',
        [Stack,ExtractFileName(ParticleFileName),NumRecords] ) ;

        // Re-order object ID numbers into contigous 1,2...n sequence
        // and determine number of objects
        NumObjects := 0 ;
        for i := 0 to NumRuns-1 do begin
            if Runs^[i].ObjectID > NumObjects then begin
               ObjectID :=  Runs^[i].ObjectID ;
               Inc(NumObjects) ;
               for j := 0 to NumRuns-1 do if Runs^[j].ObjectID = ObjectID then
                   Runs^[j].ObjectID := NumObjects ;
               end ;
            end ;

        // Remove objects smaller than a certain size
        for i := 1 to NumObjects do begin

            // Determine size of object
            NumPixels := 0 ;
            for iRun := 0 to NumRuns-1 do if Runs^[iRun].ObjectID = i then
                NumPixels := NumPixels + Runs^[iRun].Count ;

            // Erase object if smaller than threshold
            if NumPixels < Round(edObjectSizeThreshold.Value) then begin
               for iRun := 0 to NumRuns-1 do if Runs^[iRun].ObjectID = i then
               Runs^[iRun].ObjectID := 0 ;
               end ;
            end ;

        // Re-order object ID numbers into contigous 1,2...n sequence
        // and determine number of objects
        NumObjects := 0 ;
        for i := 0 to NumRuns-1 do begin
            if Runs^[i].ObjectID > NumObjects then begin
               ObjectID :=  Runs^[i].ObjectID ;
               Inc(NumObjects) ;
               for j := 0 to NumRuns-1 do if Runs^[j].ObjectID = ObjectID then
                   Runs^[j].ObjectID := NumObjects ;
               end ;
            end ;

        for ObjectID := 1 to NumObjects do begin
            NumPixels := 0 ;
            SumX := 0.0 ;
            SumY := 0.0 ;
            SumZ := 0.0 ;
            for iRun := 0 to NumRuns-1 do if Runs^[iRun].ObjectID = ObjectID then begin
                SumY := SumY + Runs^[iRun].y*Runs^[iRun].Count ;
                SumZ := SumZ + Runs^[iRun].z*Runs^[iRun].Count ;
                NumPixels := NumPixels + Runs^[iRun].Count ;
                for i := Runs^[iRun].x to Runs^[iRun].x + Runs^[iRun].Count -1 do
                    SumX := SumX + i ;
                end ;

            // Add results to table
            MainFrm.ResultsFrms[ResultsFrmNum].RecordNum := NumRecords ;
            MainFrm.ResultsFrms[ResultsFrmNum].SetVariable( StackVar, Stack ) ;
            MainFrm.ResultsFrms[ResultsFrmNum].SetVariable( ObjectIDVar, ObjectID ) ;
            MainFrm.ResultsFrms[ResultsFrmNum].SetVariable( TImeVar,
                                               Stack*NumSectionsPerStack*SaveFile.TResolution ) ;
            MainFrm.ResultsFrms[ResultsFrmNum].SetVariable( XVar,
                                               (SumX / NumPixels)*SaveFile.XResolution ) ;
            MainFrm.ResultsFrms[ResultsFrmNum].SetVariable( YVar,
                                               (FrameHeight - (SumY/NumPixels))*SaveFile.YResolution ) ;
            MainFrm.ResultsFrms[ResultsFrmNum].SetVariable( ZVar,
                                               (SumZ / NumPixels)*SaveFile.ZResolution ) ;
            MainFrm.ResultsFrms[ResultsFrmNum].SetVariable( VolVar,
                                               Abs( NumPixels*SaveFile.ZResolution
                                                    *SaveFile.XResolution*SaveFile.YResolution)) ;
            Inc(NumRecords) ;

            if bOK.Enabled then Break ;
            Application.ProcessMessages ;

            end ;

        // Create intensity-coded binary image of objects
        if rbZProjection.Checked then begin
            // Display as 2D Z projections
            // ---------------------------
            for SectionNum := StartSection to EndSection do begin
                for i := 0 to NumComponentsPerFrame-1 do PWordArray(POutBuf)^[i] := 0 ;
                for iRun := 0 to NumRuns-1 do begin
                    iStart := Runs^[iRun].x + Runs^[iRun].y*FrameWidth ;
                    for i := iStart to iStart + Runs^[iRun].Count -1 do
                        PWordArray(POutBuf)^[i] := Runs^[iRun].ObjectID ;
                    end ;

               if bOK.Enabled then Break ;
               Application.ProcessMessages ;

               end ;
            // Save image to destination
            Inc(OutFrame) ;
            SaveFile.SaveFrame( OutFrame, POutBuf ) ;
            // Report progress
            MainFrm.StatusBar.SimpleText := format(
            'Detect Objects: Stack %.4d to %s (%d objects)',
            [Stack,ExtractFileName(ParticleFileName),NumRecords] ) ;
            end
        else begin
            // Display as 3D stacks
            // --------------------
            for SectionNum := StartSection to EndSection do begin
                for i := 0 to NumComponentsPerFrame-1 do PWordArray(POutBuf)^[i] := 0 ;
                for iRun := 0 to NumRuns-1 do if Runs^[iRun].z = SectionNum then begin
                    iStart := Runs^[iRun].x + Runs^[iRun].y*FrameWidth ;
                    for i := iStart to iStart + Runs^[iRun].Count -1 do
                        PWordArray(POutBuf)^[i] := Runs^[iRun].ObjectID ;
                    end ;
                // Save image to destination
                Inc(OutFrame) ;
                SaveFile.SaveFrame( OutFrame, POutBuf ) ;

               if bOK.Enabled then Break ;
               Application.ProcessMessages ;

               end ;
            end ;

        if bOK.Enabled then Break ;
        Application.ProcessMessages ;

        end ;

     // Close file
     SaveFile.CloseFile ;

     MainFrm.StatusBar.SimpleText := format(
                                     'Detect Objects: File %s created (%d bits per pixel) (%d objects)',
                                            [ExtractFileName(ParticleFileName),
                                             MainFrm.ViewFrms[Src].PixelDepth,
                                             NumRecords] ) ;

     // Display objects image
     MainFrm.ResultsFrms[ResultsFrmNum].ViewFrmNum := MainFrm.CreateNewViewFrm(ParticleFileName) ;
     MainFrm.ViewFrms[MainFrm.ResultsFrms[ResultsFrmNum].ViewFrmNum].SetPalette(palFalseColor);

     // Display results
     MainFrm.ResultsFrms[ResultsFrmNum].ShowTable ;
     // Save to file
     MainFrm.ResultsFrms[ResultsFrmNum].SaveTableToFile ;

     finally
        FreeMem( PSrcBuf ) ;
        FreeMem( POutBuf ) ;
        FreeMem( Runs ) ;
        bOK.Enabled := True ;
        end ;

     end;




procedure TParticleDetectionFrm.cbSourceFileChange(Sender: TObject);
// ----------------------
// File selection changed
// ----------------------
//
begin
     NewFile ;
     end;

procedure TParticleDetectionFrm.bCancelClick(Sender: TObject);
begin
     bOK.Enabled := True ;
     end;

end.
