Implementando algoritmos de papel técnico em C ++ ou MATLAB

Sou graduado em Engenharia Elétrica. Tenho lido muitos artigos técnicos sobre algoritmos de processamento de sinal e imagem (reconstrução, segmentação, filtragem, etc.). A maioria dos algoritmos mostrados nesses trabalhos é definida em tempo e frequência contínuos e geralmente fornece as soluções em termos de equações complicadas. Como você implementaria um documento técnico do zero em C ++ ou MATLAB para replicar os resultados obtidos no referido documento?

Mais especificamente, eu estava olhando o artigo "Um algoritmo geral de reconstrução de feixe de cone", de Wang et al. ( IEEE Trans Med Imaging. 1993; 12 (3): 486-96 ), e fiquei pensando: como eu começo implementando seu algoritmo? A equação 10 fornece a fórmula da imagem reconstruída em. Como você codificaria isso? Você teria um loop for passando por cada voxel e computando a fórmula correspondente? Como você codificaria funções de funções nessa fórmula? Como você avaliaria as funções em pontos arbitrários?

Eu li o livro "Digital Image Processing" de Gonzalez e Woods, mas ainda estou perdido. Também li sobre a série de livros de Receitas Numéricas. Essa seria a maneira correta?

Quais são suas experiências de programação de algoritmos de trabalhos de pesquisa? Alguma dica ou sugestão?

Os algoritmos de processamento de sinais definidos em tempo / espaço / frequência contínuos são tipicamente implementados amostrando o sinal em uma grade discreta e convertendo integrais em somas (e derivadas em diferenças). Os filtros espaciais são implementados por convolução com um núcleo de convolução (ou seja, soma ponderada de vizinhos).

Existe um enorme conhecimento sobre a filtragem de sinais amostrados no domínio do tempo; os filtros no domínio do tempo são implementados como filtros de resposta de impulso finito , em que a amostra de saída atual é calculada como uma soma ponderada das N amostras de entrada anteriores; ou filtros de resposta a impulso infinito, em que a saída atual é uma soma ponderada das entradas e saídas anteriores . Formalmente, os filtros de tempo discretos são descritos usando a transformação z , que é o analógico de tempo discreto da transformação Laplace . A transformação bilinear mapeia uma para a outra ( c2de d2cno Matlab).

Como você avaliaria as funções em pontos arbitrários?

Quando você precisa do valor de um sinal em um ponto que não esteja diretamente na sua grade de amostragem, você interpola o valor de pontos próximos. A interpolação pode ser tão simples quanto escolher a amostra mais próxima, calcular uma média ponderada das amostras mais próximas ou ajustar uma função analítica arbitrariamente complicada aos dados amostrados e avaliar essa função nas coordenadas necessárias. Interpolar em uma grade mais fina e uniforme é aumentar a amostragem . Se o seu sinal (contínuo) original não contiver detalhes (ou seja, frequências) mais finos que metade da grade de amostragem, a função contínua poderá ser perfeitamente reconstruída a partir da versão amostrada (o teorema da amostragem de Nyquist-Shannon ). Para um exemplo de como você pode interpolar em 2D, consulteinterpolação bilinear .

No Matlab, você pode usar interp1ou interp2para interpolar dados 2D 1D ou amostrados regularmente (respectivamente) ou griddatapara interpolar dados 2D amostrados irregularmente.

Você teria um loop for passando por cada voxel e calculando a fórmula correspondente?

Sim, exatamente.

O Matlab evita que você precise fazer isso através de loops explícitos, porque foi projetado para operar em matrizes e vetores (ou seja, matrizes multidimensionais). No Matlab, isso é chamado de "vetorização". Integrais definidas pode ser aproximada com sum, cumsum, trapz, cumtrapz, etc.

Eu li o livro "Digital Image Processing" de Gonzalez e Woods, mas ainda estou perdido. Também li sobre a série de livros de Receitas Numéricas. Essa seria a maneira correta?

Sim, receitas numéricas seria um ótimo começo. É muito prático e abrange a maioria dos métodos numéricos de que você precisará. (Você descobrirá que o Matlab já implementa tudo o que precisa, mas as Receitas numéricas fornecerão um excelente histórico.)

Eu peguei uma classe "algoritmos e estruturas de dados", mas não vejo a relação entre o material apresentado lá e a implementação de algoritmos científicos.

O material tratado nos cursos "Algoritmos e estruturas de dados" tende a se concentrar em estruturas como listas, matrizes, árvores e gráficos contendo números inteiros ou seqüências de caracteres e operações como classificação e seleção: problemas para os quais normalmente existe um único resultado correto. Quando se trata de algoritmos científicos, isso é apenas metade da história. A outra metade diz respeito a métodos para estimar números reais e funções analíticas. Você encontrará isso em um curso sobre "Métodos Numéricos" (ou "Análise Numérica"; como este- role para baixo nos slides): como estimar funções especiais, como estimar integrais e derivadas, etc. Aqui uma das principais tarefas é estimar a precisão do seu resultado, e um padrão comum é repetir uma rotina que melhore uma estimar até que seja suficientemente preciso. (Você pode se perguntar como o Matlab sabe fazer algo tão simples quanto estimar um valor sin(x)para alguns x.)

Como um exemplo simples, aqui está um pequeno script que calcula uma transformação de radônio de uma imagem no Matlab. A transformação de radônio faz projeções de uma imagem sobre um conjunto de ângulos de projeção. Em vez de tentar calcular a projeção em um ângulo arbitrário, eu giro a imagem inteira usando imrotate, para que a projeção seja sempre vertical. Em seguida, podemos fazer a projeção simplesmente usando sum, uma vez que a summatriz retorna um vetor contendo a soma sobre cada coluna.

Você pode escrever o seu próprio, imrotatese preferir, usando interp2.

%%# Home-made Radon Tranform

%# load a density map (image).  
A = phantom;

n_pixels = size(A, 1);  %# image width (assume square)

%# At what rotation angles do we want to take projections?
n_thetas = 101;
thetas = linspace(0, 180, n_thetas);

result = zeros(n_thetas, n_pixels);

%# Loop over angles
for ii=1:length(thetas)
    theta = thetas(ii);
    rotated_image = imrotate(A, theta, 'crop');
    result(ii, :) = sum(rotated_image);
end

%# display the result
imagesc(thetas, 1:n_pixels, result.');
xlabel('projection angle [degrees]');

O que antes era parte integrante da densidade ao longo de um raio é agora uma soma sobre uma coluna de uma imagem com amostragem discreta, que por sua vez foi encontrada interpolando a imagem original sobre um sistema de coordenadas transformado.

Adicionando à excelente explicação de nibot , apenas mais alguns pontos.

• Ambientes de computação numérica como MATLAB, Octave ou SciPy / NumPy economizarão muito esforço em comparação com fazer tudo sozinho em uma linguagem de programação genérica como C ++. Malabarismo com doublematrizes e loops simplesmente não se compara a ter tipos de dados como números complexos e operações como integrais na ponta dos dedos. (É factível, com certeza, e um bom código C ++ pode ser uma ordem de magnitude mais rápida, com boas abstrações e modelos de biblioteca pode até ser razoavelmente limpo e claro, mas é definitivamente mais fácil começar, por exemplo, com MATLAB.)

• O MATLAB também possui "kits de ferramentas" para, por exemplo, Processamento de Imagem e Processamento de Sinal Digital , que podem ajudar bastante, dependendo do que você faz.

• O Processamento de sinais digitais da Mitra é um bom livro para aprender (no MATLAB!) Os conceitos básicos de tempo discreto, filtros, transformações etc., que é um conhecimento praticamente obrigatório para a execução de algoritmos técnicos decentes.

Métodos numéricos. Geralmente, é um curso universitário e um livro da divisão superior.

O DSP geralmente fica próximo à interseção de métodos numéricos e implementação eficiente. Se você ignorar a eficiência, o que poderá estar procurando é qualquer método de aproximação numérica que possa produzir um resultado "suficientemente preciso" para as equações de interesse do documento técnico. Às vezes, pode-se lidar com dados amostrados, em que os teoremas da amostragem colocam alguns limites no método de aquisição de dados (pré-filtragem) e no intervalo ou na qualidade dos resultados que você pode obter com esses dados.

Às vezes, o Matlab, receitas numéricas ou várias bibliotecas de processamento de imagem / sinal terão algoritmos ou códigos eficientes para a solução numérica desejada. Mas, às vezes, você pode ter que fazer o seu próprio, por isso, ajuda a conhecer a matemática por trás de vários métodos de solução numérica. E esse é um grande assunto por si só.