pearl-area-profiles.ipf

Go to the documentation of this file.

#pragma rtGlobals=3     // Use modern global access method and strict wave access.
#pragma IgorVersion = 6.1
#pragma ModuleName = PearlAreaProfiles
#pragma version = 1.05



threadsafe function /wave ad_extract_rod(dataset, x1, x2, y1, y2, z1, z2, destname, [noavg, sdev, pscale])
    wave dataset
    variable x1, x2, y1, y2, z1, z2
    string destname
    variable noavg
    variable sdev
    variable pscale

    if (wavedims(dataset) != 3)
        return $""
    endif
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    if (ParamIsDefault(sdev))
        sdev = 0
    endif
    if (ParamIsDefault(pscale))
        pscale = 0
    endif
    
    variable p1, p2, q1, q2, r1, r2
    if (pscale)
        p1 = x1
        p2 = x2
        q1 = y1
        q2 = y2
        r1 = z1
        r2 = z2
    else
        p1 = round((x1 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
        p2 = round((x2 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
        q1 = round((y1 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
        q2 = round((y2 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
        r1 = round((z1 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
        r2 = round((z2 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
    endif
    
    if ((numtype(p1) == 2) || (numtype(p2) == 2))
        return ad_extract_rod_x(dataset, min(q1, q2), max(q1, q2), min(r1, r2), max(r1, r2), destname, noavg=noavg, sdev=sdev)
    elseif ((numtype(q1) == 2) || (numtype(q2) == 2))
        return ad_extract_rod_y(dataset, min(p1, p2), max(p1, p2), min(r1, r2), max(r1, r2), destname, noavg=noavg, sdev=sdev)
    elseif ((numtype(r1) == 2) || (numtype(r2) == 2))
        return ad_extract_rod_z(dataset, min(p1, p2), max(p1, p2), min(q1, q2), max(q1, q2), destname, noavg=noavg, sdev=sdev)
    else
        return $""
    endif
end

threadsafe function /wave ad_extract_rod_x(dataset, q1, q2, r1, r2, destname, [noavg, sdev])
    wave dataset
    variable q1, q2, r1, r2
        // -inf < q1 < q2 < +inf, -inf < r1 < r2 < +inf
    string destname
    variable noavg
    variable sdev
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    if (ParamIsDefault(sdev))
        sdev = 0
    endif
    variable avg = !noavg && !sdev

    q1 = max(q1, 0)
    q2 = min(q2, dimsize(dataset, 1) - 1)
    r1 = max(r1, 0)
    r2 = min(r2, dimsize(dataset, 2) - 1)
    
    if (strlen(destname) > 0)
        duplicate /r=[][q1,q1][r1,r1]/o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[][q1,q1][r1,r1] /free dataset, w_dest
    endif
    redimension /n=(dimsize(w_dest, 0)) w_dest
    setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), w_dest
    
    w_dest = 0
    variable qq, rr
    variable nn = 0
    for (qq = q1; qq <= q2; qq += 1)
        for (rr = r1; rr <= r2; rr += 1)
            w_dest += dataset[p][qq][rr]
            nn += 1
        endfor
    endfor
    
    if (sdev)
        duplicate /free w_dest, w_squares
        w_squares = 0
        for (qq = q1; qq <= q2; qq += 1)
            for (rr = r1; rr <= r2; rr += 1)
                w_squares += dataset[p][qq][rr]^2
            endfor
        endfor
    endif
    
    if (avg)
        w_dest /= nn
    elseif (sdev)
        w_dest = sqrt( (w_squares  - w_dest^2 / nn) / (nn - 1) )
    endif
    
    return w_dest
end

threadsafe function /wave ad_extract_rod_y(dataset, p1, p2, r1, r2, destname, [noavg, sdev])
    wave dataset
    variable p1, p2, r1, r2
        // 0 <= p1 < p2 < dimsize(0), 0 <= r1 < r2 < dimsize(2)
    string destname
    variable noavg
    variable sdev
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    if (ParamIsDefault(sdev))
        sdev = 0
    endif
    variable avg = !noavg && !sdev
    
    p1 = max(p1, 0)
    p2 = min(p2, dimsize(dataset, 0) - 1)
    r1 = max(r1, 0)
    r2 = min(r2, dimsize(dataset, 2) - 1)

    if (strlen(destname) > 0)
        duplicate /r=[p1,p1][][r1,r1]/o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[p1,p1][][r1,r1] /free dataset, w_dest
    endif
    redimension /n=(dimsize(w_dest, 1)) w_dest
    setscale /p x dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), w_dest
    
    w_dest = 0
    variable pp, rr
    variable nn = 0
    for (pp = p1; pp <= p2; pp += 1)
        for (rr = r1; rr <= r2; rr += 1)
            w_dest += dataset[pp][p][rr]
            nn += 1
        endfor
    endfor

    if (sdev)
        duplicate /free w_dest, w_squares
        w_squares = 0
        for (pp = p1; pp <= p2; pp += 1)
            for (rr = r1; rr <= r2; rr += 1)
                w_squares += dataset[pp][p][rr]^2
                nn += 1
            endfor
        endfor
    endif
    
    if (avg)
        w_dest /= nn
    elseif (sdev)
        w_dest = sqrt( (w_squares  - w_dest^2 / nn) / (nn - 1) )
    endif
    
    return w_dest
end

threadsafe function /wave ad_extract_rod_z(dataset, p1, p2, q1, q2, destname, [noavg, sdev])
    wave dataset
    variable p1, p2, q1, q2
        // 0 <= p1 < p2 < dimsize(0), 0 <= q1 < q2 < dimsize(1)
    string destname
    variable noavg
    variable sdev
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    if (ParamIsDefault(sdev))
        sdev = 0
    endif
    variable avg = !noavg && !sdev

    p1 = max(p1, 0)
    p2 = min(p2, dimsize(dataset, 0) - 1)
    q1 = max(q1, 0)
    q2 = min(q2, dimsize(dataset, 1) - 1)

    if (strlen(destname) > 0)
        duplicate /r=[p1,p1][q1,q1][]/o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[p1,p1][q1,q1][] /free dataset, w_dest
    endif
    redimension /n=(dimsize(w_dest, 2)) w_dest
    setscale /p x dimoffset(dataset, 2), dimdelta(dataset, 2), waveunits(dataset, 2), w_dest
    
    w_dest = 0
    variable pp, qq
    variable nn = 0
    for (pp = p1; pp <= p2; pp += 1)
        for (qq = q1; qq <= q2; qq += 1)
            w_dest += dataset[pp][qq][p]
            nn += 1
        endfor
    endfor

    if (sdev)
        duplicate /free w_dest, w_squares
        w_squares = 0
        for (pp = p1; pp <= p2; pp += 1)
            for (qq = q1; qq <= q2; qq += 1)
                w_squares += dataset[pp][qq][p]^2
                nn += 1
            endfor
        endfor
    endif
    
    if (avg)
        w_dest /= nn
    elseif (sdev)
        w_dest = sqrt( (w_squares  - w_dest^2 / nn) / (nn - 1) )
    endif
    
    return w_dest
end

threadsafe function /wave ad_extract_slab(dataset, x1, x2, y1, y2, z1, z2, destname, [noavg, pscale])
    wave dataset
    variable x1, x2, y1, y2, z1, z2
    string destname
    variable noavg
    variable pscale
    
    if (wavedims(dataset) != 3)
        return $""
    endif
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    if (ParamIsDefault(pscale))
        pscale = 0
    endif
    
    variable p1, p2, q1, q2, r1, r2
    if (pscale)
        p1 = x1
        p2 = x2
        q1 = y1
        q2 = y2
        r1 = z1
        r2 = z2
    else
        p1 = round((x1 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
        p2 = round((x2 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
        q1 = round((y1 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
        q2 = round((y2 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
        r1 = round((z1 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
        r2 = round((z2 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
    endif
    
    if ((numtype(p1) < 2) && (numtype(p2) < 2))
        return ad_extract_slab_x(dataset, min(p1, p2), max(p1, p2), destname, noavg=noavg)
    elseif ((numtype(q1) < 2) && (numtype(q2) < 2))
        return ad_extract_slab_y(dataset, min(q1, q2), max(q1, q2), destname, noavg=noavg)
    elseif ((numtype(r1) < 2) && (numtype(r2) < 2))
        return ad_extract_slab_z(dataset, min(r1, r2), max(r1, r2), destname, noavg=noavg)
    else
        return $""
    endif
end

threadsafe function /wave ad_extract_slab_x(dataset, p1, p2, destname, [noavg])
    wave dataset
    variable p1, p2
        // x coordinate range (point scaling) to be integrated
        // -inf <= p1 < p2 <= +inf
    string destname // name of destination wave. to be created in current data folder. overrides existing.
        // if empty, the function returns a free wave
    variable noavg // zero or default = average, non-zero = sum
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    p1 = max(p1, 0)
    p2 = min(p2, dimsize(dataset, 0) - 1)

    if (strlen(destname) > 0)
        duplicate /r=[p1,p1][][]/o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[p1,p1][][] /free dataset, w_dest
    endif
    redimension /n=(dimsize(w_dest, 1), dimsize(w_dest, 2)) w_dest
    setscale /p x dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), w_dest
    setscale /p y dimoffset(dataset, 2), dimdelta(dataset, 2), waveunits(dataset, 2), w_dest
    
    w_dest = 0
    variable pp
    variable nn = 0
    for (pp = p1; pp <= p2; pp += 1)
        w_dest += dataset[pp][p][q]
        nn += 1
    endfor

    if (noavg == 0)
        w_dest /= nn
    endif
    
    return w_dest
end

threadsafe function /wave ad_extract_slab_y(dataset, q1, q2, destname, [noavg])
    wave dataset
    variable q1, q2
        // y coordinate range (point scaling) to be integrated
        // -inf <= q1 < q2 <= +inf
    string destname // name of destination wave. to be created in current data folder. overrides existing.
        // if empty, the function returns a free wave
    variable noavg // zero or default = average, non-zero = sum
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    q1 = max(q1, 0)
    q2 = min(q2, dimsize(dataset, 1) - 1)

    if (strlen(destname) > 0)
        duplicate /r=[][q1,q1][]/o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[][q1,q1][] /free dataset, w_dest
    endif
    redimension /n=(dimsize(w_dest, 0), dimsize(w_dest, 2)) w_dest
    setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), w_dest
    setscale /p y dimoffset(dataset, 2), dimdelta(dataset, 2), waveunits(dataset, 2), w_dest
    
    w_dest = 0
    variable qq
    variable nn = 0
    for (qq = q1; qq <= q2; qq += 1)
        w_dest += dataset[p][qq][q]
        nn += 1
    endfor

    if (noavg == 0)
        w_dest /= nn
    endif
    
    return w_dest
end

threadsafe function /wave ad_extract_slab_z(dataset, r1, r2, destname, [noavg])
    wave dataset
    variable r1, r2
        // z coordinate range (point scaling) to be integrated
        // -inf <= r1 < r2 <= +inf
    string destname // name of destination wave. to be created in current data folder. overrides existing.
        // if empty, the function returns a free wave
    variable noavg // zero or default = average, non-zero = sum
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif
    r1 = max(r1, 0)
    r2 = min(r2, dimsize(dataset, 2) - 1)

    if (strlen(destname) > 0)
        duplicate /r=[][][r1,r1]/o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[][][r1,r1] /free dataset, w_dest
    endif
    redimension /n=(dimsize(w_dest, 0), dimsize(w_dest, 1)) w_dest
    setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), w_dest
    setscale /p y dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), w_dest
    
    w_dest = 0
    variable rr
    variable nn = 0
    for (rr = r1; rr <= r2; rr += 1)
        w_dest += dataset[p][q][rr]
        nn += 1
    endfor

    if (noavg == 0)
        w_dest /= nn
    endif
    
    return w_dest
end

threadsafe function /wave ad_profile_x(dataset, q1, q2, destname, [noavg])
    wave dataset
    variable q1, q2 // -inf <= q1 < q2 <= +inf
        // deprecated: q2 = -1 stands for dimsize(0) - 1
    string destname // name of destination wave. to be created in current data folder. overrides existing.
        // if empty, the function returns a free wave
    variable noavg // zero or default = average, non-zero = sum
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif

    if (strlen(destname) > 0)
        duplicate /r=[0,0][] /o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[0,0][] /free dataset, w_dest
    endif
    
    return ad_profile_x_w(dataset, q1, q2, w_dest, noavg=noavg)
end

threadsafe function /wave ad_profile_x_w(dataset, q1, q2, destwave, [noavg])
    wave dataset
    variable q1, q2 // -inf <= q1 < q2 <= +inf
        // deprecated: q2 = -1 stands for dimsize(0) - 1
    wave destwave // existing destination wave
    variable noavg // zero or default = average, non-zero = sum
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif

    redimension /n=(dimsize(dataset, 0)) destwave
    setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), destwave
    setscale d 0, 0, waveunits(dataset, -1), destwave
    
    q1 = max(q1, 0)
    if (q2 < 0)
        q2 = inf
    endif
    q2 = min(q2, dimsize(dataset, 1) - 1)
    
    destwave = 0
    variable qq
    variable nn = 0
    for (qq = q1; qq <= q2; qq += 1)
        destwave += dataset[p][qq]
        nn += 1
    endfor

    if (noavg == 0)
        destwave /= nn
    endif
    
    return destwave
end

threadsafe function /wave ad_profile_y(dataset, p1, p2, destname, [noavg])
    wave dataset
    variable p1, p2 // -inf <= p1 < p2 < inf
        // deprecated: p2 = -1 stands for dimsize(0) - 1
    string destname // name of destination wave. to be created in current data folder. overrides existing.
        // if empty, the function returns a free wave
    variable noavg // zero or default = average, non-zero = sum
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif

    if (strlen(destname) > 0)
        duplicate /r=[][0,0] /o dataset, $destname
        wave w_dest = $destname
    else
        duplicate /r=[][0,0] /free dataset, w_dest
    endif
    MatrixTranspose w_dest
    
    return ad_profile_y_w(dataset, p1, p2, w_dest, noavg=noavg)
end

threadsafe function /wave ad_profile_y_w(dataset, p1, p2, destwave, [noavg])
    wave dataset
    variable p1, p2 // -inf <= p1 < p2 < inf
        // deprecated: p2 = -1 stands for dimsize(0) - 1
    wave destwave // existing destination wave
    variable noavg // zero or default = average, non-zero = sum
    
    if (ParamIsDefault(noavg))
        noavg = 0
    endif

    redimension /n=(dimsize(dataset, 1)) destwave
    setscale /p x dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), destwave
    setscale d 0, 0, waveunits(dataset, -1), destwave
    
    p1 = max(p1, 0)
    if (p2 < 0)
        p2 = inf
    endif
    p2 = min(p2, dimsize(dataset, 0) - 1)
    
    destwave = 0
    variable pp
    variable nn = 0
    for (pp = p1; pp <= p2; pp += 1)
        destwave += dataset[pp][p]
        nn += 1
    endfor

    if (noavg == 0)
        destwave /= nn
    endif
    
    return destwave
end

threadsafe function calc_y_profile_mins(image)
    // experimental
    wave image
    
    wave yminlocs = ad_profile_x(image, 0, 0, "ymins", noavg=1)
    variable nx = dimsize(image, 0)
    variable ix
    for (ix = 0; ix < nx; ix += 1)
        wave profile = ad_profile_y(image, ix, ix, "", noavg=1)
        wavestats /q/m=1 profile
        yminlocs[ix] = v_minloc
    endfor
end

function ad_collect_multiscan_y(dataset, positions, destwave, [noavg])
    wave dataset
    wave positions
    wave destwave
    variable noavg
    
    variable tol = (wavemax(positions) - wavemin(positions)) / numpnts(positions) / 100
    
    duplicate /free positions, positions_sorted
    sort positions_sorted, positions_sorted
    duplicate /free positions_sorted, positions_diff
    differentiate /p /meth=2 positions_sorted /d=positions_diff
    positions_diff[0] = 1
    extract /free positions_sorted, positions_unique, positions_diff > tol
    variable n_unique = numpnts(positions_unique)
    redimension /n=(dimsize(dataset, 0), n_unique) destwave
    
    variable i
    variable nx, ny
    for (i = 0; i < n_unique; i += 1)
        extract /free dataset, data_extract, abs(positions[q] - positions_unique[i]) < tol
        nx = dimsize(dataset, 0)
        ny = dimsize(data_extract, 0) / nx
        redimension /n=(nx, ny) data_extract
        wave profile = ad_profile_x(data_extract, -inf, inf, "", noavg=noavg)
        destwave[][i] = profile[p]
    endfor
end