Celestin Apprentice 4/Apprentice-Release4.iso/Source Code/Pascal/Applications/NIH Image 1.59/1.59 Source/fft.p

Text File | 1995-10-10 | 16.8 KB | 714 lines | [TEXT/MPPS]

unit fft;

{
This file contains a Pascal language implementation of the
Fast Hartley Transform algorithm which is covered under
United States Patent Number 4,646,256.

This code may therefore be freely used and distributed only
under the following conditions:

        1)  This header is included in every copy of the code; and
        2)  The code is used for noncommercial research purposes only.

 Firms using this code for commercial purposes will be infringing a United
 States patent and should contact the

                    Office of Technology Licensing
                    Stanford University
                    857 Serra Street, 2nd Floor
                    Stanford, CA   94305-6225
                    [phone]

This implementation is based on Pascal code contibuted
by Arlo Reeves ([e-mail]).
}

interface
    uses
        Memory, QuickDraw, Packages, Menus, Events, Fonts, ToolUtils, globals, Utilities, Graphics, Filters, Analysis;

    procedure doFFT(fftKind: fftType);
    procedure DisplayPowerSpectrum(x: rImagePtr);
    function   isPowerOf2 (x: LongInt): boolean;
    procedure SetFFTWindowName(doInverse: boolean);
    procedure RedisplayPowerSpectrum;
    procedure doSwapQuadrants;


implementation

    const
        UpdateTicks = 10; {Show progress 6 times/sec}

    var
        AbortFFT: boolean;
        C, S: rLinePtr;


    function log2 (x: LongInt): LongInt;
        var
            count: LongInt;
    begin
        count := 15;
        while not BTST(x, count) do
            count := count - 1;
        log2 := count;
    end;


    function BitRevX (x, bitlen: LongInt): LongInt;
        var
            i: LongInt;
            temp: longint;
    begin
        temp := 0;
        for i := 0 to bitlen do
            if BTST(x, i) then
                BSET(temp, bitlen - i - 1);
        BitRevX := LoWord(temp);
    end;


    procedure BitRevRArr (x: rLinePtr; bitlen, maxN: LongInt);
        var
            i: LongInt;
            tempArr: rLineType;
    begin
        for i := 0 to maxN - 1 do
            tempArr[i] := x^[BitRevX(i, bitlen)];
        BlockMove(@tempArr, x, maxN * SizeOf(real));
    end;


    procedure transposeR (x: rImagePtr; maxN: LongInt);
        var
            r, c: LongInt;
            rTemp: real;
    begin
        for r := 0 to maxN - 1 do
            for c := r to maxN - 1 do
                if r <> c then
                    begin
                        rTemp := x^[ r * maxN + c];
                        x^[ r * maxN + c] := x^[c * maxN + r];
                        x^[c * maxN + r] := rTemp;
                    end;
    end;


    procedure FHTps (r, maxN: LongInt; inMat: rImagePtr;  var outArr: rLineType);
{ Power Spectrum of one row from 2D Hartley Transform }
        var
            c, base: LongInt;
    begin
        base := r * maxN;
        for c := 0 to maxN - 1 do
            outArr[c] := (sqr(inMat^[base + c]) + sqr(inMat^[((maxN - r) mod maxN) * maxN + (maxN - c) mod maxN])) / 2;
    end;


    procedure dfht3 (x: rLinePtr; inverse: boolean; maxN: LongInt);
    { an optimized real FHT }
        var
            i, stage, gpNum, gpIndex, gpSize, numGps, Nlog2: LongInt;
            bfNum, numBfs: LongInt;
            Ad0, Ad1, Ad2, Ad3, Ad4, CSAd: LongInt;
            rt1, rt2, rt3, rt4: real;
    begin
        Nlog2 := log2(maxN);
        BitRevRArr(x, Nlog2, maxN);            { bitReverse the input array }

        gpSize := 2;                        { first & second stages - do radix 4 butterflies once thru }
        numGps := maxN div 4;
        for gpNum := 0 to numGps - 1 do
            begin
                Ad1 := gpNum * 4;
                Ad2 := Ad1 + 1;
                Ad3 := Ad1 + gpSize;
                Ad4 := Ad2 + gpSize;
                rt1 := x^[Ad1] + x^[Ad2];            { a + b }
                rt2 := x^[Ad1] - x^[Ad2];            { a - b }
                rt3 := x^[Ad3] + x^[Ad4];            { c + d }
                rt4 := x^[Ad3] - x^[Ad4];            { c - d }
                x^[Ad1] := rt1 + rt3;                { a + b + (c + d ) }
                x^[Ad2] := rt2 + rt4;                { a - b + (c - d ) }
                x^[Ad3] := rt1 - rt3;                { a + b - (c + d ) }
                x^[Ad4] := rt2 - rt4;                { a - b - (c - d ) }
            end;

        if Nlog2 > 2 then
            begin                        { third + stages computed here }
                gpSize := 4;
                numBfs := 2;
                numGps := numGps div 2;
                for stage := 2 to Nlog2 - 1 do
                    begin
                        for gpNum := 0 to numGps - 1 do
                            begin
                                Ad0 := gpNum * gpSize * 2;
                                Ad1 := Ad0;                            { 1st butterfly is different from others - no mults needed }
                                Ad2 := Ad1 + gpSize;
                                Ad3 := Ad1 + gpSize div 2;
                                Ad4 := Ad3 + gpSize;
                                rt1 := x^[Ad1];
                                x^[Ad1] := x^[Ad1] + x^[Ad2];
                                x^[Ad2] := rt1 - x^[Ad2];
                                rt1 := x^[Ad3];
                                x^[Ad3] := x^[Ad3] + x^[Ad4];
                                x^[Ad4] := rt1 - x^[Ad4];
                                for bfNum := 1 to numBfs - 1 do
                                    begin        { subsequent BF's dealt with together }
                                        Ad1 := bfNum + Ad0;
                                        Ad2 := Ad1 + gpSize;
                                        Ad3 := gpSize - bfNum + Ad0;
                                        Ad4 := Ad3 + gpSize;

                                        CSAd := bfNum * numGps;
                                        rt1 := x^[Ad2] * C^[CSAd] + x^[Ad4] * S^[CSAd];
                                        rt2 := x^[Ad4] * C^[CSAd] - x^[Ad2] * S^[CSAd];

                                        x^[Ad2] := x^[Ad1] - rt1;
                                        x^[Ad1] := x^[Ad1] + rt1;
                                        x^[Ad4] := x^[Ad3] + rt2;
                                        x^[Ad3] := x^[Ad3] - rt2;

                                    end;         { for bfNum := 0 to ... }
                            end;            { for gpNum := 0 to ... }
                        gpSize := gpSize * 2;
                        numBfs := numBfs * 2;
                        numGps := numGps div 2;
                    end;
            end;        { if }

        if inverse then
            for i := 0 to maxN - 1 do
                x^[i] := x^[i] / maxN;
    end;


    procedure rc2Dfht (x: rImagePtr; inverse: boolean; maxN: LongInt);
{ Row-column 2D FHT }
        var
            row, col, mRow, mCol, NextTicks: LongInt;
            A, B, C, D, E: real;
            theRow: rLinePtr;
    begin
        NextTicks := TickCount + UpdateTicks;
        for row := 0 to maxN - 1 do begin
            theRow := rLinePtr(ord4(x) + row * maxN * SizeOf(real));
            dfht3(theRow, inverse, maxN);
            if TickCount > nextTicks then begin
                UpdateMeter(round(row/maxN * 50.0), 'Computing FHT');
                nextTicks := TickCount + updateTicks;
                if CommandPeriod then begin
                        beep;
                        AbortFFT := true;
                        exit(rc2Dfht)
                    end;
            end;
        end;
        transposeR(x, maxN);
        for row := 0 to maxN - 1 do begin
            theRow := rLinePtr(ord4(x) + row * maxN * SizeOf(real));
            dfht3(theRow, inverse, maxN);
            if TickCount > nextTicks then begin
                UpdateMeter(round(row/maxN * 50.0) + 50, 'Computing FHT');
                nextTicks := TickCount + updateTicks;
                if CommandPeriod then begin
                        beep;
                        AbortFFT := true;
                        exit(rc2Dfht)
                    end;
            end;
        end;
        transposeR(x, maxN);
        for row := 0 to maxN div 2 do             { Now calculate actual Hartley transform }
            for col := 0 to maxN div 2 do
                begin
                    mRow := (maxN - row) mod maxN;
                    mCol := (maxN - col) mod maxN;
                    A := x^[row * maxN + col];                        { see Bracewell, 'Fast 2D Hartley Transf.' IEEE Procs. 9/86 }
                    B := x^[mRow * maxN + col];
                    C := x^[row * maxN + mCol];
                    D := x^[mRow * maxN + mCol];
                    E := ((A + D) - (B + C)) / 2;
                    x^[row * maxN + col] := A - E;
                    x^[mRow * maxN + col] := B + E;
                    x^[row * maxN + mCol] := C + E;
                    x^[mRow * maxN + mCol] := D - E;
                end;
        UpdateMeter(-1, 'Hide Meter');
    end;


procedure SwapQuadrants;
{Swap quadrants 1 and 3 and quadrants 2 and 4 so the}
{power spectrum origin is at the center of the image.}
{2 1}
{3 4}
var
    row, col, halfMaxN, tmp, maxN: LongInt;
    line1, line2: LineType;
begin
    maxN := info^.PixelsPerLine;
    halfMaxN := maxN div 2;
    for row:= 0 to halfMaxN -1 do begin
        GetLine(0, row, maxN, line1);
        GetLine(0, row + halfMaxN, maxN, line2);
        for col := 0 to halfMaxN -1 do begin
            tmp := line1[col];
            line1[col] := line1[col + halfMaxN];
            line1[col + halfMaxN] := tmp;
            tmp := line2[col];
            line2[col] := line2[col + halfMaxN];
            line2[col + halfMaxN] := tmp;
        end;
        PutLine(0, row, maxN, line2);
        PutLine(0, row + halfMaxN, maxN, line1);
    end;
end;



    procedure DisplayRealImage(x: rImagePtr);
        var
            row, col, i, base, maxN: LongInt;
            r, min, max, scale: real;
            line: lineType;
    begin
        maxN := info^.PixelsPerLine;
        min := 1e99;
        max := -1e99;
        i := 1;
        for row := 0 to maxN - 1 do begin
            base := row * maxN;
            for col := 0 to maxN - 1 do begin
                    r := x^[base + col];
                    if r < min then min := r;
                    if r > max then max := r;
            end;
        end;
        scale := 253.0 / (max - min);
        for row := 0 to maxN - 1 do begin
                base := row * maxN;
                for col := 0 to maxN - 1 do begin
                        r := x^[base + col];
                        line[col] := round((r - min) * scale) + 1;
                end;
                PutLine(0, row, maxN, line);
            end;
    end;


    procedure DisplayPSUsingBuffer(fht, ps: rImagePtr; var rLine: rLineType);
    var
        row, col, base, nextTicks, maxN: LongInt;
        r, min, max, scale: real;
        line: lineType;

    begin
        maxN := info^.PixelsPerLine;
        nextTicks := TickCount + updateTicks;
        min := 1e99;
        max := -1e99;
        for row := 0 to maxN - 1 do begin
            FHTps (row, maxN, fht, rLine);
            base := row * maxN;
            for col := 0 to maxN - 1 do begin
                    r := rLine[col];
                    if r < min then min := r;
                    if r > max then max := r;
                    ps^[base + col] := r;
            end;
            if TickCount > nextTicks then begin
                UpdateMeter(round(row/maxN * 80.0), 'Computing Power Spectrum');
                nextTicks := TickCount + updateTicks;
                if CommandPeriod then begin
                        beep;
                        AbortFFT := true;
                        exit(DisplayPSUsingBuffer)
                    end;
            end;
        end;
        if min < 1.0 then
            min := 0.0
        else
            min := ln(min);
        max := ln(max);
        scale := 253.0 / (max - min);
        for row := 0 to maxN - 1 do begin
                base := row * maxN;
                for col := 0 to maxN - 1 do begin
                        r := ps^[base + col];
                        if r < 1.0 then
                            r := 0.0
                        else
                            r := ln(r);
                        line[col] := round((r - min) * scale) + 1;
                end;
                PutLine(0, row, maxN, line);
                if TickCount > nextTicks then begin
                    UpdateMeter(round(row/maxN * 20.0 ) + 80, 'Computing Power Spectrum');
                    nextTicks := TickCount + updateTicks;
                    if CommandPeriod then begin
                            beep;
                            AbortFFT := true;
                            exit(DisplayPSUsingBuffer)
                        end;
                end;
            end;
        SwapQuadrants;
        UpdateMeter(-1, 'Hide Meter');
    end;


    procedure DisplayPowerSpectrum(fht: rImagePtr);
    var
        row, col, nextTicks, maxN: LongInt;
        r, min, max, scale: real;
        line: lineType;
        rLine: rLineType;
        tempH: handle;
        ps: rImagePtr;
    begin
        maxN := info^.PixelsPerLine;
        tempH := GetBigHandle(maxN * maxN * SizeOf(real));
        if tempH <> nil then begin
            hlock(tempH);
            ps := rImagePtr(tempH^);
            DisplayPSUsingBuffer(fht, ps, rLine);
            hunlock(tempH);
            DisposeHandle(tempH);
            exit(DisplayPowerSpectrum);
        end;
        min := 1e99;
        max := -1e99;
        nextTicks := TickCount + updateTicks;
        for row := 0 to maxN - 1 do begin
            FHTps (row, maxN, fht, rLine);
            for col := 0 to maxN - 1 do begin
                    r := rLine[col];
                    if r < min then min := r;
                    if r > max then max := r;
            end;
            if TickCount > nextTicks then begin
                UpdateMeter(round(row/maxN * 35.0), 'Computing Power Spectrum');
                nextTicks := TickCount + updateTicks;
                if CommandPeriod then begin
                        beep;
                        AbortFFT := true;
                        exit(DisplayPowerSpectrum)
                    end;
            end;
        end;
        if min < 1.0 then
            min := 0.0
        else
            min := ln(min);
        max := ln(max);
        scale := 253.0 / (max - min);
        for row := 0 to maxN - 1 do begin
                FHTps (row, maxN, fht, rLine);
                for col := 0 to maxN - 1 do begin
                        r := rLine[col];
                        if r < 1.0 then
                            r := 0.0
                        else
                            r := ln(r);
                        line[col] := round((r - min) * scale) + 1;
                end;
                PutLine(0, row, maxN, line);
                if TickCount > nextTicks then begin
                    UpdateMeter(round(row/maxN * 65.0 ) + 35, 'Computing Power Spectrum');
                    nextTicks := TickCount + updateTicks;
                    if CommandPeriod then begin
                            beep;
                            AbortFFT := true;
                            exit(DisplayPowerSpectrum)
                        end;
                end;
            end;
        SwapQuadrants;
        UpdateMeter(-1, 'Hide Meter');
    end;


procedure ConvertToReal;
var
    row, col, i, sum, base: LongInt;
    width, height: LongInt;
    mean, pixelCount: real;
    line: LineType;
    DataP: rImagePtr;
begin
    with info^ do begin
        width := pixelsPerLine;
        height := nLines;
        hlock(DataH);
        DataP := rImagePtr(DataH^);
    end;
    UpdateMeter(0, 'Converting to Real');
    {
    GetRectHistogram;
    sum := 0;
    pixelCount := width * height;
    for i := 0 to 255 do
        sum := sum + histogram[i] * i;
    mean := sum / pixelCount;
    }
    for row:= 0 to height - 1 do begin
        GetLine(0, row, width, line);
        base := row * width;
        for col := 0 to width - 1 do
            DataP^[base + col] := line[col] {- mean};
    end;
    hunlock(info^.DataH);
end;


function isPowerOf2 (x: LongInt): boolean;
var
    i, sum: integer;
begin
sum := 0;
x := abs(x);
for i := 0 to 15 do
    sum := sum + ord(BTST(x, i));
IsPowerOf2 := (sum <= 1);
end;


function PowerOf2Size: boolean;
var
    width, height: LongInt;

    procedure fail;
    begin
         PutError('A square, power of two size image or selection (128x128, 256x256, etc.) required.');
         AbortMacro;
         PowerOf2Size := false;
         exit(PowerOf2Size);
    end;

begin
    with info^ do begin
        if info = noInfo then
            fail;
        width := pixelsPerLine;
        height := nLines;
        if RoiShowing and (roiType = rectRoi) then with roiRect do begin
            width := right - left;
            height := bottom - top;
        end;
        if not isPowerOf2(width) or (width <> height) then
            fail;
        if width < 4 then begin
                PutMessage('Sequence Length must be >= 4.');
                PowerOf2Size := true;
                exit(PowerOf2Size);
            end;
        PowerOf2Size := true;
    end;
end;


function MakeSinCosTables(maxN: LongInt): boolean;
var
    i: LongInt;
    theta, dTheta: real;
begin
    C := rLinePtr(NewPtr(SizeOf(rLineType)));
    S := rLinePtr(NewPtr(SizeOf(rLineType)));
    if (C = nil) or (S = nil) then begin
        MakeSinCosTables := false;
        PutError('Out of Memory');
        exit(MakeSinCosTables);
    end;
    theta := 0;
    dTheta := 2 * pi / maxN;
    for i := 0 to maxN div 4 - 1 do
        begin
            C^[i] := cos(theta);
            S^[i] := sin(theta);
            theta := theta + dTheta;
        end;
        MakeSinCosTables := true;
end;


function MakeRealImage: boolean;
var
    TempH: handle;
    maxN: LongInt;
begin
    maxN := info^.PixelsPerLine;
    tempH := GetBigHandle(maxN * maxN * SizeOf(real));
    if TempH = nil then begin
        PutMemoryAlert;
        MakeRealImage := false;
        exit(MakeRealImage);
    end;
    if not Duplicate(StringOf('FFT ', nPics + 1:1), false) then begin
        DisposeHandle(handle(TempH));
        exit(MakeRealImage);
    end;
    UpdatePicWindow;
    info^.DataH := tempH;
    ConvertToReal;
    UpdateTitleBar;
    MakeRealImage := true;
end;


procedure SetFFTWindowName(doInverse: boolean);
begin
    with info^ do begin
        if doInverse then
            title := StringOf('Inverse FFT ', picNum:1)
        else
            title := StringOf('FFT ', picNum:1);
        UpdateWindowsMenuItem;
        UpdateTitleBar;
    end;
end;


procedure ApplyFilter(rData: rImagePtr);
var
    row, col, width, height, base, i: LongInt;
    line: LineType;
    passMode: boolean;
    t:FateTable;
begin
    SwapQuadrants;
    with info^ do begin
        width := pixelsPerLine;
        height := nLines;
    end;
    for row:= 0 to height - 1 do begin
        GetLine(0, row, width, line);
        base := row * width;
        for col := 0 to width - 1 do
            rData^[base + col] := line[col]/255.0 * rData^[base + col];
    end;
end;


procedure doMasking(rData: rImagePtr);
var
    row, col, width, height, base, i: LongInt;
    line: LineType;
    passMode: boolean;
    t:FateTable;
begin
    GetRectHistogram;
    if (histogram[0] = 0) and (histogram[255] = 0) then
        exit(doMasking);
    UpdateMeter(0, 'Masking');
    passMode := histogram[255] <> 0;
    if passMode then
        ChangeValues(0,254,0)
    else
        ChangeValues(1,255,255);
    for i := 1 to 3 do
        Filter(UnweightedAvg, 0, t);
    UpdatePicWindow;
    ApplyFilter(rData);
end;


    procedure doFFT(fftKind: fftType);
    var
        startTicks, maxN: LongInt;
        trect: rect;
        RealData: rImagePtr;
        doInverse: boolean;
    begin
            doInverse := fftKind <> ForewardFFT;
            if not PowerOf2Size then
                exit(doFFT);
            startTicks := tickCount;
            if info^.DataH = nil then begin
                if doInverse then begin
                    PutError('A real image is required to do an inverse transform.');
                     AbortMacro;
                    exit(doFFT);
                end;
                if not MakeRealImage then begin
                     AbortMacro;
                    exit(doFFT);
                end
            end else begin
                KillRoi;
                SetFFTWindowName(doInverse);
            end;
            hlock(info^.DataH);
            RealData := rImagePtr(info^.DataH^);
            ShowWatch;
            maxN := info^.PixelsPerLine;
            if not MakeSinCosTables(maxN) then
                exit(doFFT);
            AbortFFT := false;
            ShowMessage(CmdPeriodToStop);
            if doInverse then begin
                if fftKind = InverseFFTWithMask then
                    doMasking(RealData)
                else if fftKind = InverseFFTWithFilter then
                    ApplyFilter(RealData);
                rc2DFHT(RealData, true, maxN);
                if not AbortFFT then
                    DisplayRealImage(RealData);
            end else begin
                rc2DFHT(RealData, false, maxN);
                if not AbortFFT then
                    DisplayPowerSpectrum(RealData);
            end;
            if AbortFFT then
                UpdateMeter(-1, 'Hide');
            hunlock(info^.dataH);
            UpdatePicWindow;
            SetRect(trect, 0, 0, maxN, maxN);
            ShowTime(startTicks, trect, '');
            UpdateWindowsMenuItem;
            DisposePtr(ptr(C));
            DisposePtr(ptr(S));
    end;


procedure RedisplayPowerSpectrum;
    var
        rData: rImagePtr;
    begin
            if info = noInfo then
                exit(RedisplayPowerSpectrum);
            KillRoi;
            if not PowerOf2Size then
                exit(RedisplayPowerSpectrum);
            if info^.DataH = nil then begin
                    PutError('Real image required.');
                    exit(RedisplayPowerSpectrum);
                end;
            hlock(info^.DataH);
            rData := rImagePtr(info^.DataH^);
            DisplayPowerSpectrum(rData);
            hunlock(info^.dataH);
            UpdatePicWindow;
    end;


    Procedure doSwapQuadrants;
    begin
        if info = noInfo then
            exit(doSwapQuadrants);
        KillRoi;
        if not PowerOf2Size then
            exit(doSwapQuadrants);
        SetupUndo;
        WhatToUndo := UndoEdit;
        SwapQuadrants;
        UpdatePicWindow;
    end;


end. {fft Unit}