How to use the pymedphys._imports.numpy.shape function in pymedphys

def create_axes(image, dpcm=100):
    shape = np.shape(image)
    x_span = (np.arange(0, shape[0]) - shape[0] // 2) * 10 / dpcm
    y_span = (np.arange(0, shape[1]) - shape[1] // 2) * 10 / dpcm

    return x_span, y_span

data_samp, int(np.shape(amplitude_vert)[0]) * 10
        )
        data_samp = amplitude_horz[:, slice]
        ampl_resamp_x1[:, slice] = signal.resample(
            data_samp, int(np.shape(amplitude_horz)[0]) * 10
        )

    for slice in tqdm(range(int(volume.shape[2] / 2), volume.shape[2])):
        merge_vol = merge_vol + volume[:, :, slice]
        data_samp = amplitude_vert[:, slice]
        ampl_resamp_y2[:, slice - int(volume.shape[2] / 2)] = signal.resample(
            data_samp, int(np.shape(amplitude_vert)[0]) * 10
        )
        data_samp = amplitude_horz[:, slice]
        ampl_resamp_x2[:, slice - int(volume.shape[2] / 2)] = signal.resample(
            data_samp, int(np.shape(amplitude_horz)[0]) * 10
        )

    fig, ax = plt.subplots(ncols=1, nrows=1, squeeze=True, figsize=(6, 8))

    extent = (0, 0 + (volume.shape[1] * dx), 0, 0 + (volume.shape[0] * dy))

    ax.imshow(merge_vol, extent=extent, aspect="auto")
    ax.set_aspect("equal", "box")
    ax.set_xlabel("x distance [mm]")
    ax.set_ylabel("y distance [mm]")

    ax.hlines(dy * int(volume_resort.shape[0] / 3.25), 0, dx * volume_resort.shape[1])
    ax.text(
        dx * int(volume_resort.shape[1] / 2.25),
        dy * int(volume_resort.shape[0] / 3),
        "Profile 2",

volume
    )  # this will hold the resorted volume 0 to 3 clockwise
    junctions_comb = []
    peaks_figs_comb = []

    # we need to create 4 matches

    # 0,1,2,3 will be tagged top left, top right, bottom right, bottom left
    for i in range(0, int(np.shape(volume)[2])):
        diag_stack = [
            0,
            0,
            0,
            0,
        ]  # we will sum along one direction whichever is biggest will tag the file
        for j in range(0, int(min([np.shape(volume)[0], np.shape(volume)[1]]) / 2)):
            # print('j=',j,int(np.shape(volume)[0] / 2)+j, int(np.shape(volume)[1] / 2)+j)
            diag_stack[0] = (
                diag_stack[0]
                + volume[
                    int(np.shape(volume)[0] / 2) - j,
                    int(np.shape(volume)[1] / 2) - j,
                    i,
                ]
            )
            diag_stack[1] = (
                diag_stack[1]
                + volume[
                    int(np.shape(volume)[0] / 2) - j,
                    int(np.shape(volume)[1] / 2) + j,
                    i,
                ]

pixels_per_mm = 4
    elif np.shape(img) == (512, 512):
        pixels_per_mm = 2
    else:
        raise ValueError(
            "Expect iView images to be either 1024x1024 or 512x512 "
            f"pixels\nShape = {np.shape(img)}"
        )

    img = img[:, 1:-1]

    if img.dtype != np.dtype("uint16"):
        raise ValueError("Expect iView images to have a pixel type of unsigned 16 bit")
    img = 1 - img[::-1, :] / 2 ** 16

    shape = np.shape(img)
    x = np.arange(-shape[1] / 2, shape[1] / 2) / pixels_per_mm
    y = np.arange(-shape[0] / 2, shape[0] / 2) / pixels_per_mm

    return x, y, img

consecutive_keys = find_consecutive_logfiles(
        field_id_key_map, field_id, filehash, index
    )

    logfilepaths = [get_filepath(index, config, key) for key in consecutive_keys]

    logfile_results = calc_and_merge_logfile_mudensity(
        logfilepaths, grid_resolution=grid_resolution
    )

    try:
        assert np.all(logfile_results[0] == mosaiq_results[0])
        assert np.all(logfile_results[1] == mosaiq_results[1])
    except AssertionError:
        print(np.shape(logfile_results[0]))
        print(np.shape(mosaiq_results[0]))
        raise

    grid_xx = logfile_results[0]
    grid_yy = logfile_results[1]

    logfile_mu_density = logfile_results[2]
    mosaiq_mu_density = mosaiq_results[2]

    return grid_xx, grid_yy, logfile_mu_density, mosaiq_mu_density

],
                axis=0,
            )
            maxstack1 = np.amax(stack1)

            stack2 = np.sum(
                volume[
                    int(
                        np.shape(volume)[0] / 2
                        - np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ) : int(
                        np.shape(volume)[0] / 2
                        + np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ),
                    int(
                        np.shape(volume)[1] / 2
                        - np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ) : int(
                        np.shape(volume)[1] / 2
                        + np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ),
                    slice,
                ],
                axis=1,
            )
            maxstack2 = np.amax(stack2)

            if maxstack2 / maxstack1 > 1.1:  # It is a Y field folder
                if ky == 0:
                    yfield = volume[:, :, slice]
                    yfield = yfield[:, :, np.newaxis]
                else:

def get_interpolated_dose(coords_grid, dose_interpolation):
    coords_grid_ij_indexing = np.array(
        [
            np.ravel(coords_grid[:, :, 1]),
            np.ravel(coords_grid[:, :, 0]),
            np.ravel(coords_grid[:, :, 2]),
        ]
    ).T

    interpolated_dose = dose_interpolation(coords_grid_ij_indexing)
    coords_dim = np.shape(coords_grid)
    interpolated_dose = np.reshape(interpolated_dose, (coords_dim[0], coords_dim[1]))

    return interpolated_dose

volume = 1 - volume  # inverting the range

        min_val = np.amin(volume)  # normalizing
        volume = volume - min_val
        volume = volume / (np.amax(volume))

        # print('Volume shape=', np.shape(volume), np.amin([np.shape(volume)[0], np.shape(volume)[1]]))
        kx = 0
        ky = 0
        krot = 0
        for slice in range(0, volume.shape[2]):
            stack1 = np.sum(
                volume[
                    int(
                        np.shape(volume)[0] / 2
                        - np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ) : int(
                        np.shape(volume)[0] / 2
                        + np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ),
                    int(
                        np.shape(volume)[1] / 2
                        - np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ) : int(
                        np.shape(volume)[1] / 2
                        + np.amin([np.shape(volume)[0], np.shape(volume)[1]]) / 2
                    ),
                    slice,
                ],
                axis=0,
            )
            maxstack1 = np.amax(stack1)

]  # the numbers here are reversed because we are going to slide the second graph (the overlay) to minimize the error #for profile 3
    amplitude_vert[:, 1] = volume_resort[
        :, int(volume_resort.shape[1]) - int(volume_resort.shape[1] / 2.8), 2
    ]
    amplitude_vert[:, 3] = volume_resort[
        :, int(volume_resort.shape[1] / 2.8), 3
    ]  # for profile 4
    amplitude_vert[:, 2] = volume_resort[:, int(volume_resort.shape[1] / 2.8), 0]

    plt.figure()
    for slice in tqdm(range(0, int(volume.shape[2] / 2))):
        merge_vol = merge_vol + volume[:, :, slice]

        data_samp = amplitude_vert[:, slice]
        ampl_resamp_y1[:, slice] = signal.resample(
            data_samp, int(np.shape(amplitude_vert)[0]) * 10
        )
        data_samp = amplitude_horz[:, slice]
        ampl_resamp_x1[:, slice] = signal.resample(
            data_samp, int(np.shape(amplitude_horz)[0]) * 10
        )

    for slice in tqdm(range(int(volume.shape[2] / 2), volume.shape[2])):
        merge_vol = merge_vol + volume[:, :, slice]
        data_samp = amplitude_vert[:, slice]
        ampl_resamp_y2[:, slice - int(volume.shape[2] / 2)] = signal.resample(
            data_samp, int(np.shape(amplitude_vert)[0]) * 10
        )
        data_samp = amplitude_horz[:, slice]
        ampl_resamp_x2[:, slice - int(volume.shape[2] / 2)] = signal.resample(
            data_samp, int(np.shape(amplitude_horz)[0]) * 10
        )