Mover Média Filtro Matlab Código Exemplo

Capítulo 22 Processamento de imagem com MATLAB e GPU Processamento de imagem com MATLAB e GPU 1 Departamento de Oncologia, Universidade de Cambridge, Cambridge, Reino Unido 1. Introdução MATLAB (The MathWorks, Natick, MA, EUA) é um pacote de software para computação numérica que pode ser Usado em várias disciplinas científicas, como matemática, física, eletrônica, engenharia e biologia. Mais de 40 caixas de ferramentas estão disponíveis na versão atual (R2013b lançado em setembro de 2013), que inclui inúmeras funções integradas aprimoradas pelo acesso a uma linguagem de programação de alto nível. Como as imagens podem ser representadas por matrizes 2D ou 3D e o mecanismo de processamento MATLAB depende da representação matricial de todas as entidades, o MATLAB é particularmente adequado para implementação e teste de fluxos de trabalho de processamento de imagem. O Image Processing Toolbox (IPT) inclui todas as ferramentas necessárias para o processamento de imagens de propósito geral, incorporando mais de 300 funções que foram otimizadas para oferecer uma boa precisão e alta velocidade de processamento. Além disso, o built-in Parallel Computing Toolbox (PCT) foi recentemente expandido e agora aceita a aceleração da unidade de processamento gráfico (GPU) para algumas funções do IPT. No entanto, para muitos aplicativos de processamento de imagens, ainda precisamos escrever nosso próprio código, tanto no MATLAB quanto, no caso de aplicativos acelerados por GPU que exigem controle específico sobre os recursos da GPU, na CUDA (Nvidia Corporation, Santa Clara, CA, EUA). Neste capítulo, a primeira parte é dedicada a algumas ferramentas essenciais do IPT que podem ser utilizadas na análise e avaliação de imagens, bem como na extração de informações úteis para posterior processamento e avaliação. Isso inclui a recuperação de informações sobre imagens digitais, ajustes e processamento de imagens, bem como extração de recursos e gerenciamento de vídeo. A segunda parte é dedicada à aceleração GPU de técnicas de processamento de imagem, quer usando as funções PCT integradas ou através da escrita de nossas próprias funções. Cada seção é acompanhada pelo código de exemplo MATLAB. As funções e o código fornecidos neste capítulo são adotados a partir da documentação MATLAB 1, 2, salvo indicação em contrário. 2. Processamento de imagem na CPU 2.1. Conceitos básicos de imagem 2.1.1. Representação de pixels Uma imagem digital é uma representação visual de uma cena que pode ser obtida usando um dispositivo óptico digital. É composto por uma série de elementos de imagem, pixels, e pode ser bidimensional (2D) ou tridimensional (3D). Podem encontrar-se diferentes imagens de profundidade de bits, mas as mais comuns no processamento científico de imagens são as imagens binárias de 1 bit (valores de pixel 0 ou 1), as imagens de escala de cinza de 8 bits (intervalo de pixels de 0 a 255) e as 16 Imagens em cores de bits (intervalo de pixels 0-65535) 3. A Figura 1 mostra a variação da escala de cinza, de preto a branco, para uma imagem de 8 bits. Variação das intensidades de escala de cinza para uma imagem de 8 bits. 2.1.2. MATLAB pixel convention MATLAB usa indexação baseada em um único, onde o primeiro pixel ao longo de qualquer dimensão tem índice1, enquanto muitas outras plataformas são baseadas em zero e consideram o primeiro índice a ser 0. Por convenção, a contagem de índices de pixels começa a partir da parte superior esquerda Canto de uma imagem com o primeiro e segundo índices aumentando para baixo e para a direita, respectivamente. A Figura 2 visualiza a forma como MATLAB indexa uma imagem de pixel de 512512. Esta informação é particularmente importante quando o usuário pretende aplicar uma transformação a um pixel específico ou a uma vizinhança de pixels. Convenção de indexação de pixel MATLABs. 2.1.3. Formatos de imagem Vários formatos de imagem são suportados pelo MATLAB, incluindo os mais usados, como JPEG, TIFF, BMP, GIF e PNG. As imagens podem ser lidas, processadas e depois guardadas em um formato diferente da inicial. Vários parâmetros como a resolução, a profundidade do bit, o nível de compressão ou o espaço de cores podem ser ajustados de acordo com as preferências dos usuários. 2.1.4. Processamento de imagens digitais O processamento de imagens digitais refere-se à modificação de uma imagem digital usando um computador para enfatizar a informação relevante. Isso pode ser conseguido por ambos revelando detalhes latentes e suprimindo o ruído indesejado. O processo geralmente é projetado de acordo com o resultado final desejado, que pode incluir a detecção de objeto de aprimoramento de imagem simples, segmentação ou estimativa de parâmetro de rastreamento ou classificação de condição. Além disso, ao lidar com imagens destinadas a inspeção de pessoas, a estrutura e a função do sistema visual humano podem ser um fator crítico na concepção de qualquer técnica que determine o que pode ser percebido como uma característica facilmente distinguível. 2.2. Pré-processamento de imagem O pré-processamento de imagem é um procedimento que fornece informações iniciais sobre a condição digital de uma imagem candidata. Para receber essas informações, precisamos carregar a imagem na plataforma do software e examinar seus valores de tipo e pixel. 2.2.1. Entrada e saída de imagem Assim como qualquer outro conjunto de dados em MATLAB, as imagens são representadas por uma variável. Se considerarmos um arquivo de imagem com imagem de nome e formato tiff, usando a função imread. A imagem pode ser carregada como uma matriz 2D ou 3D, dependendo do tipo de imagem. A visualização da imagem é conseguida usando a função para mostrar uma nova figura, uma chamada para mostrar deve ser precedida por uma chamada a figura. Se, em vez disso, imshow é usado sozinho, a nova imagem substitui a última na última janela de figura aberta. Salvar uma imagem pode ser conseguida usando a função imwrite onde a imagem desejada (ou seja, variável), o formato e o nome final devem ser especificados. Embora o nome da imagem salva possa ser escolhido livremente, ao construir uma grande pipeline, sugere-se que o nome de cada imagem resultante seja representativo da etapa de processamento para manter a coerência do processo. O exemplo a seguir representa como esta seqüência pode ser alcançada. A função pode também ser acompanhada por um vetor de dois elementos de baixa altura com o qual o usuário especifica a faixa de exibição da imagem da escala de cinza. Se este vetor for deixado vazio. Os valores mínimos e máximos de escala de cinza da imagem são exibidos como pixels em preto e branco, respectivamente, com os valores intermediários como vários tons de cinza. 2.2.2. Conversões do tipo de imagem A conversão do tipo de imagem é uma ferramenta útil, pois o usuário pode converter a imagem de entrada em qualquer tipo desejado. Uma conversão especial e muitas vezes útil é a transformação de um inteiro não assinado em uma representação de dupla precisão. Qualquer algoritmo de processamento de imagem pode resultar em resultados mais precisos, uma vez que esta conversão aumenta a faixa dinâmica de intensidades. O alcance da imagem resultante é de 0,0 a 1,0 com o MATLAB mantendo até 15 dígitos decimais. Os seguintes comandos são exemplos de conversões de imagens. 2.2.3. Informação de pixels Um histograma é uma representação de intensidade útil, pois revela a distribuição de intensidade de pixels. Pode ser obtido usando a função imhist. Esta informação pode, por exemplo, ser usada para selecionar um valor limiar apropriado. Além disso, um perfil de intensidades também pode revelar informações sobre variações locais de intensidade que podem ser usadas para modelar pequenos detalhes. A função de arquivo melhorável pode ser aplicada em pixels pré-selecionados ou como uma função de prompt de comando para o usuário, a fim de selecionar a área desejada manualmente. Exemplo de código para tais processos segue 1, 3. 2.2.4. Ajuste de contraste Uma das principais técnicas de pré-processamento é o ajuste de contraste, pois este processo pode melhorar os recursos desejados, enquanto suprime outros, indesejados. O MATLAB possui várias ferramentas para variar o contraste da imagem. A função imcontrast fornece uma ferramenta de ajuste manual através da qual o usuário pode experimentar e encontrar o contraste ótimo. Os parâmetros resultantes podem então ser adotados, salvos e aplicados em uma pilha de imagens tiradas nas mesmas condições. A função imadjust pode ser usada para especificar um intervalo de intensidade quando o usuário conhece os valores ótimos ou os encontrou usando a ferramenta imcontrast. A mesma função também fornece entrada para o fator gama do ajuste de contraste não-linear de energia. Além disso, uma transformação logarítmica feita sob medida pode ser aplicada 3. A Figura 3 apresenta uma imagem original de escala de cinza e sua contrapartida ajustada ao contraste usando os parâmetros especificados no exemplo anterior. Imagem original (esquerda, 12961936 pixels) e resultado ajustado de contraste (direita). Parâmetros a. B e c podem ser definidos e ajustados pelo usuário, o que significa que qualquer transformação logarítmica customizada pode ser introduzida de acordo com as necessidades específicas. Outras técnicas que podem afetar o contraste são baseadas em histograma. Um histograma representa uma função de densidade de distribuição de intensidade de nível de cinza ou de imagem. Esse conhecimento pode ajudar no processamento posterior, ajudando o usuário a escolher as ferramentas certas 4. O estiramento do histograma pode ser realizado através da função de imadjust, enquanto a equalização do histograma pode ser realizada através da função histeq. A equalização do histograma adaptativo também pode ser aplicada usando a função adapthisteq que considera pequenas vizinhanças de imagens ao invés de toda a imagem como entrada. O parâmetro NumTilesValue assume a forma de um vetor que especifica o número de telhas em cada direção. Outros parâmetros também podem ser especificados na função adapthisteq, como o alcance dinâmico dos dados de saída ou a forma do histograma. A Figura 4 mostra exemplos de equalização de histograma e equalização de histograma adaptativo, respectivamente. Transformação de equalização de histograma (esquerda) e transformação de equalização de histograma adaptativo (direita) da imagem original na Figura 3. Observe os detalhes aprimorados na imagem à direita. 2.2.5. Operações aritméticas As operações aritméticas referem-se a adição, subtração, multiplicação e divisão de duas imagens ou uma imagem e uma constante. As imagens sujeitas a operações aritméticas precisam ter as mesmas dimensões e a representação da escala de cinza. A imagem resultante terá as mesmas dimensões que a entrada. Quando um valor constante é adicionado ou subtraído (em vez de uma segunda imagem), esta constante é adicionada ou subtraída a cada intensidade de pixels, aumentando ou reduzindo a luminosidade da imagem. Na maioria das vezes, essas operações são usadas para aprimoramento de detalhes ou supressão de informações desnecessárias. No código acima, o segundo parâmetro de entrada (im2) pode ser substituído por uma constante escalar. 2.2.6. Transformações diversas Outras transformações de imagem úteis incluem o cultivo, o redimensionamento e a rotação. O corte pode ser usado se o usuário estiver interessado apenas em uma parte específica da imagem de entrada. O usuário pode definir uma região específica de interesse e aplicar qualquer transformação somente a esta parte. O redimensionamento pode ser aplicado para expandir ou reduzir o tamanho da imagem. A redução do tamanho da imagem pode ser especialmente útil para acelerar um processo no caso de imagens maiores ou grandes conjuntos de dados. A rotação pode ser particularmente útil quando uma imagem inclui características de uma direcionalidade particular. O usuário pode especificar o método de interpolação aplicada fora do vizinho mais próximo, bilinear e bicúbico. A inversão de intensidades de escala de cinza pode ser útil quando os objetos interessantes têm intensidades menores do que o plano de fundo. As seguintes funções executam esses processos. 2.3. Processamento de imagem 2.3.1. Thresholding Thresholding é um dos conceitos mais importantes no processamento de imagens, pois ele encontra aplicação em quase todos os projetos. O limiar pode ser manual ou automático, global ou local. No modo manual, o usuário define um valor limiar, geralmente dependendo da concepção da imagem (vários testes podem ser necessários). No modo automático, é necessária uma compreensão mais detalhada da imagem para selecionar o método correto. O IPT fornece a função graythresh baseada no método Otsus e o caractere bimodal de uma imagem 5. Esse limite global criará uma imagem em preto e branco onde os pixels com intensidades acima deste limite tornar-se-ão brancos (valor 1) e os pixels com intensidades abaixo desse limite tornar-se-ão preto (valor 0). Este método pode ser facilmente estendido para multi-thresholding usando a função IPT multithresh. Usando esta função, o usuário especifica um número adequado de níveis de limiar (k) para a imagem. Se esse parâmetro não for fornecido, ele possui a mesma funcionalidade que a função graythresh original. O IPT pode visualizar o resultado da função multithresh usando a função imquantize. O último rotula as várias áreas da imagem de acordo com o número de limiares previamente especificados. A imagem marcada pode então ser transformada em uma imagem RGB, preservando o tipo (por exemplo, uint8) da entrada original. O código a seguir pode ser usado nesses processos. A Figura 5 fornece um exemplo de aplicação de um único e de vários limiares na imagem original da Figura 3. Imagem de um único limite (esquerda) e imagem de vários limiares usando 5 valores de limiar (direita). 2.3.2. Detecção de borda A detecção de borda é uma parte essencial do processamento de imagens, pois geralmente enfatiza o esboço de objetos ou a estrutura interna. Uma vantagem é uma representação da descontinuidade em uma imagem e pode caracterizar uma superfície que separa dois objetos ou um objeto do fundo da imagem 4. Os limites podem ser caracterizados por pixels individuais ou seqüências de pixels conectadas. Tal recurso pode auxiliar no reconhecimento adicional de objetos e, portanto, a detecção de borda é aplicada em muitas seqüências de processamento de imagem. O resultado da detecção de borda é uma imagem binária com bordas apresentadas por pixels brancos. 2.3.3. Operadores de detecção de borda de primeiro orden O IPT inclui os detectores padrão de borda de primeiro orden, como Roberts, Sobel e Prewitt. O detector de borda Roberts conta com 22 máscaras, enquanto as outras duas dependem de 33 máscaras. Um limite opcional pode ser especificado para definir a magnitude do gradiente mínimo. É o código útil para esses detectores. Outros detectores de borda de primeiro orden também podem ser projetados. Exemplos são as máscaras Kirsch e Robinson que não estão incluídas no IPT, mas são fáceis de projetar. São exemplos de detectores de borda direcionais que digitalizam a imagem de diferentes direções para detectar bordas com várias orientações. É utilizado um único kernel que, através de rotações de 0 a 315 em passos de 45, cria oito máscaras diferentes. A imagem é convolvida com cada máscara e os pixels na imagem final são atribuídos a maior magnitude de detecção de borda obtida de qualquer uma das máscaras 4 O código a seguir apresenta esses dois detectores de borda, respectivamente, 6, 4. O detector de pontos, outro exemplo de um detector de borda, detecta pontos brilhantes com base na diferença de intensidade entre um pixel central e seus vizinhos. Um detector de ponto pode ser especificado pelo seguinte código 7. 2.3.4. Operadores de detecção de borda de segunda ordem Além dos detectores de borda de primeira ordem, os detectores de borda de segunda ordem podem encontrar uma ampla aplicação. Tais detectores são, por exemplo, o Canny, zero-cross e Laplacian-of-Gaussian (LoG também chamado Marr-Hildreth). O método Canny usa a derivada de um filtro gaussiano para encontrar o gradiente da imagem original, após o que se baseia nos máximos locais da imagem resultante. 3 O método de cruzamento zero busca cruzamentos de zero depois de um filtro arbitrário ter sido aplicado. Finalmente, o método LoG busca zero cruzamentos após a transformação LoG ter sido aplicada. É o código útil para esses detectores. Neste caso, o limiar refere-se à força de uma borda sigma refere-se ao desvio padrão do filtro gaussiano enquanto o filtro se refere a qualquer filtro que o usuário aplica antes da detecção de borda. Nos métodos LoG e Canny, o limiar e o sigma podem ser deixados não especificados, mas, no caso do método de cruzamento zero, o usuário deve definir um filtro. A Figura 6 apresenta as imagens resultantes após a aplicação dos detectores de borda Roberts e Canny, respectivamente. Aplicação do detector de borda Roberts (esquerda) e do detector de borda Canny com um valor limiar de 0,05 (à direita). 2.3.5. Filtragem de imagens A filtragem espacial é um dos processos mais importantes no processamento de imagens, pois ele pode extrair e processar freqüências específicas a partir de uma imagem enquanto outras freqüências podem ser removidas ou transformadas. Normalmente, a filtragem é usada para aprimoramento de imagem ou remoção de ruído. O IPT inclui as ferramentas padrão necessárias para o filtro de imagem. A função fspecial pode ser usada para o design do filtro. Podem ser introduzidos filtros médios, médios, gaussianos, laplacianos, laplacianos de gaussianos, de movimento, de pré-borda e de borda Sobel. O filtro projetado é aplicado à imagem usando a função imfilter. Seguem exemplos típicos de código. O parâmetro hsize é um vetor que representa o número de linhas e colunas do bairro que é usado ao aplicar o filtro. O parâmetro sigma é o desvio padrão do filtro gaussiano aplicado. Os detectores de borda também podem ser aplicados através da filtragem da imagem com o operador de borda. Um exemplo segue com a aplicação do detector de borda de ponto mencionado anteriormente. Além dos filtros projetados pelo usuário, o IPT inclui filtros que podem ser diretamente aplicados à imagem. Esses exemplos são o filtro médio (medfilt2), o filtro Wiener (wiener2) ou o filtro de estatísticas de ordem 2D (ordfilt2). A vizinhança na função medfilt2 especifica as dimensões da área em que o valor médio do pixel será encontrado. A função ordfilt2 é uma versão generalizada do filtro mediano. Um bairro é definido pelos pixels de domínio não-zero. E cada pixel na imagem é substituído pela ordem - o menor dos seus vizinhos neste domínio 1. Um exemplo pode ser o seguinte comando, onde cada pixel é substituído pelo 6º valor mais pequeno encontrado em sua vizinhança 33. A Figura 7 mostra exemplos de imagens filtradas estatísticas gaussianas e de ordem. Imagens filtradas usando um filtro gaussiano (hsize 9 9 e sigma 1) e um filtro de estatísticas de 6ª ordem (direita). 2.3.6. Processamento de imagem morfológica O processamento de imagens morfológicas refere-se à extração de descritores que descrevem objetos de interesse e, portanto, sua morfologia determina as ferramentas que são utilizadas 8. Os elementos de estruturação são usados ​​para investigar a imagem 3. A função bwmorph executa todas as operações morfológicas com a adição de parâmetros adequados. Uma vez que o tempo de processamento desta função pode aumentar significativamente com a complexidade da imagem, é suportado pelo PCT para maior velocidade de processamento. O processamento morfológico inclui a dilatação, a erosão, a abertura, o fechamento, a transformação do chapéu alto e do chapéu de baixo, a transformação de sucesso ou perda, bem como outros processos que realizam mudanças específicas de pixels. A operação do parâmetro explica o tipo de operador morfológico enquanto n é o número de vezes que este processo deve ser repetido. Se n não está definido, o processo é aplicado apenas uma vez. Processos como dilatação e erosão também podem ser aplicados usando funções individuais quando um elemento estruturante feito sob medida deve ser usado. Exemplos de processos individuais seguem. A Figura 8 apresenta exemplos de dilatação e transformação do chapéu alto com um elemento estruturante de disco do raio 10. Imagem dilatada (esquerda) e uma imagem transformada de topo com um elemento de estrutura de disco do raio 10 (à direita). A transformação de distância geralmente é aplicada a imagens binárias e representa a distância entre um pixel branco e seu pixel zero mais próximo. Os pixels na nova imagem obtêm valores maiores com maior distância de um pixel zero. Essa transformação pode atuar como uma função de segmentação e é freqüentemente usada na segmentação de discos sobrepostos 1, 8. 2.3.7. Processamento de imagem em cores As imagens em cores estão sujeitas a processamento em muitos campos científicos, pois diferentes cores podem representar diferentes recursos. A representação de cores mais usada é RGB (Vermelho-Verde-Azul). A transformação de imagens RGB em intensidade de escala de cinza ou extração de uma cor específica pode ser feita usando o seguinte código: 2.4. Extração de recursos A extração de recursos é o processo através do qual os objetos reconhecidos são avaliados de acordo com alguns critérios geométricos. O primeiro passo deste processo é rotular os objetos de uma imagem binária imbin usando a função bwlabel. A imagem resultante é referida como imagem marcada (imlab). A função varre a imagem de cima para baixo e da esquerda para a direita atribuindo um número a cada pixel indicando a qual objeto pertence. Além disso, o IPT tem a função regionprops que mede os recursos de objetos rotulados, como área, diâmetro equivalente, excentricidade, perímetro e comprimentos de eixo maior e menor. Esta função opera em imagens rotuladas que contém n objetos rotulados. Uma lista completa dos recursos que podem ser calculados usando regionprops pode ser encontrada na documentação 1 do MATLAB IPT. Além dos recursos padrão incluídos no IPT, os recursos personalizados podem ser medidos usando recursos já calculados ou introduzindo medidas completamente novas. Um exemplo poderia ser a medição de uma característica geométrica padrão de objetos, bem como a relação de thinness e compacidade (ou irregularidade) usando um loop for para avaliação de todos os objetos n. Uma vez que o usuário pode ter que lidar com numerosas medidas para muitos objetos, geralmente é útil pré-alocar a memória para reduzir o tempo de processamento. O código a seguir pode ser usado para rotular objetos de imagem, medir os recursos e armazená-los como uma tabela das variáveis ​​1, 3. Os recursos medidos são armazenados em matrizes de estrutura. Normalmente, o processamento adicional de recursos requer transformações de matrizes de estrutura em matrizes. O MATLAB não pode realizar tais transformações sem a aplicação de um passo intermediário: os arrays da estrutura devem primeiro ser transformados em matrizes de células que, por sua vez, podem ser convertidas em matrizes. Observe que a transposição das propriedades das células da matriz celular foi usada para alocar recursos em colunas da matriz em vez de linhas. Por razões de desempenho, geralmente é melhor orientar os dados de forma que sejam processados ​​por coluna, em lugar de linha por linha, neste caso, esperamos passar pela característica de dados por recurso em vez de imagem por imagem. 2.5. Processando séries de imagens Em muitos casos, o gerenciamento de imagens múltiplas pode ser uma tarefa árdua, a menos que um processo automatizado possa ser estabelecido. Supondo que tenhamos um lote de 100 imagens que gostaríamos de processar, usando um loop for e definindo o caminho para o diretório de imagens, podemos carregar, processar e salvar as imagens um de cada vez. Depois de salvar a primeira imagem, o próximo no diretório é automaticamente carregado, processado e salvo. O procedimento continua até a última imagem ter sido salva. O código a seguir executa essa operação. 2.6. Manipulação de vídeo 2.6.1. Vídeo para quadros Uma aplicação interessante para processamento de imagem é o gerenciamento de dados de vídeo. Neste caso, o arquivo de vídeo deve ser dividido em quadros únicos. A função VideoReader pode ser usada para inserir o arquivo como uma variável. Para n quadros, cada quadro é salvo como uma imagem separada em qualquer formato. O código a seguir lê um vídeo (chamado filme) para uma estrutura MATLAB e salva os quadros um a um no formato tiff. 9 2.6.2. Quadros para o vídeo Uma vez que cada quadro é armazenado como uma única imagem, pode ser processado em conformidade, um a um ou em modo de lote. Um possível próximo passo neste processo pode ser combinar as imagens processadas em um único vídeo novamente. Se um novo filme (chamado movienew) for criado a partir de quadros (chamado frame), o código a seguir fornece o backbone para esse processo 9. 3. Processamento de imagem em GPU em MATLAB Grandes quantidades de dados de imagem são produzidas em muitas situações técnicas e experimentais, em particular onde as imagens são adquiridas repetidamente ao longo do tempo ou quando se trata de imagens de maior dimensionalidade do que duas. A imagem de lapso de tempo ea gravação de vídeo podem ser mencionadas como exemplos do primeiro, enquanto que o último pode ser representado por qualquer uma das muitas modalidades de imagens tomográficas presentes. A tomografia computadorizada 4D (CT), onde as imagens de TC 3D são adquiridas em intervalos regulares para monitorar o movimento interno do paciente, é um exemplo de uma aplicação pertencente a ambas as categorias. Muitas vezes é desejável ou mesmo crítico para acelerar a análise e o processamento de tais grandes conjuntos de dados de imagem, especialmente para aplicativos que se deslocam em tempo real ou em tempo real. Devido à natureza inerentemente paralela de muitos algoritmos de processamento de imagem, eles são adequados para implementação em uma unidade de processamento de gráficos (GPU) e, consequentemente, podemos esperar uma aceleração substancial de tal implementação por código executado em uma CPU. No entanto, apesar do fato de as GPUs serem hoje em dia omnipresentes nos computadores de mesa, apenas 34 das várias centenas de funções do IPT são habilitadas pela GPU pelo PCT na versão atual do MATLAB (2013b). Neste sub-capítulo, exploraremos as possibilidades disponíveis para alguém que quer aproveitar o poder de computação da GPU diretamente do MATLAB ou incorporar código GPU externo em programas MATLAB. O foco será em aplicativos de processamento de imagem, mas as técnicas apresentadas podem ser, com pouco ou nenhum esforço, adaptadas a outras aplicações. Na primeira parte desta parte, analisamos como usar as funções de processamento de imagem incorporadas e baseadas em GPU do PCT. Seguindo isso, explicamos como as manipulações em pixels podem ser realizadas usando a versão habilitada pela GPU do arrayfun e como podemos escrever nossas próprias funções de processamento de imagem, fazendo uso de mais de cem funções MATLAB elementares que foram implementadas para serem executadas em GPUs. Na segunda parte desta seção, mostramos como o PCT pode ser usado para chamar as funções do kernel escritas na linguagem de programação CUDA diretamente do MATLAB. Isso nos permite fazer uso das funções do kernel existentes em nossas aplicações MATLAB. Além disso, para quem tem conhecimento da CUDA, torna possível o maior potencial de GPUs do MATLAB e também fornece um modo fácil de testar as funções do kernel em desenvolvimento. A terceira e última parte é dedicada a usuários mais avançados que possam querer usar uma das bibliotecas CUDA fornecidas pela NVIDIA, que preferem escrever seu código em um idioma diferente da CUDA, que não tem acesso à Parallel Computing Toolbox, ou Que têm acesso a uma biblioteca de código GPU existente que eles gostariam de chamar do MATLAB. Observamos como o código CUDA pode ser compilado diretamente em funções MEX usando o PCT, seguido de uma descrição de como o código GPU escrito em CUDA ou OpenCL pode ser acessado a partir do MATLAB compilando-o em uma biblioteca e criando uma função de wrapper MEX em torno dele . Finalmente, mostramos como o código de uma função de wrapper MEX pode ser construído diretamente em nosso compilador externo e, por exemplo, incluído em uma solução existente do Visual Studio (Microsoft Corporation, Redmond, WA, EUA) para que isso seja feito automaticamente ao criar a solução. 3.1. GpuArray e funções habilitadas para GPU incorporadas Para que os exemplos desta peça funcionem, precisamos de um computador equipado com uma GPU NVIDIA de capacidade de cálculo CUDA 1.3 ou superior, que seja configurada e reconhecida corretamente pelo PCT 10. As funções gpuDeviceCount e gpuDevice podem ser usadas para identificar e selecionar uma GPU como descrito na documentação do PCT 11. Para poder processar uma imagem na GPU, os dados correspondentes primeiro devem ser copiados da memória da CPU principal através do barramento PCI para a memória GPU. No MATLAB, os dados da GPU são acessados ​​através de objetos do tipo gpuArray. O comando cria um novo objeto gpuArray chamado imGpu e atribui-lhe uma cópia dos dados da imagem em mim. ImGpu será do mesmo tipo que eu (por exemplo, double. Single. Int32, etc.), o que pode afetar o desempenho da computação GPU conforme discutido abaixo. Deixe-nos por agora assumir que eu sou uma matriz 30723072 de números de ponto flutuante de precisão única (único). Correspondentemente, quando executamos todos os nossos cálculos de GPU, chamamos a função reunir para recuperar o resultado. Por exemplo, copia os dados no resultado do objeto gpuArrayGuem novamente para resultar na memória da CPU. Em geral, a cópia de dados no barramento PCI é relativamente lenta, o que significa que, ao trabalhar com grandes conjuntos de dados, devemos tentar evitar cópias desnecessárias entre CPU e memória GPU. Quando possível, pode ser mais rápido criar filtros, máscaras ou intermediários diretamente na GPU. Para fazer isso, o gpuArray possui vários construtores estáticos correspondentes às funções MATLAB padrão, atualmente no olho. Falso. Inf. Nan. uns . verdade . Zeros. Linspace. Espaço de logs. Rand. Randi e randn. Para pré-alocar e inicializar a memória GPU. Estes podem ser invocados chamando gpuArray. Construtor onde o construtor é substituído pela chamada de função. Por exemplo, cria uma matriz de 3072 3072 de números de pseudorandom normalmente distribuídos com média zero e desvio padrão. Tal como acontece com as funções MATLAB padrão correspondentes, o último argumento de função especifica o tipo de elemento de matriz (neste caso, único) e, se omitido, o padrão é duplicado. Embora isso normalmente não seja um problema ao trabalhar em CPUs modernas, vale a pena ter em mente que as GPUs do consumidor da NVIDIA são muitas vezes várias vezes mais rápidas no processamento de números de ponto flutuante de precisão única (único), em comparação com a dupla precisão (duplo) ou inteiros (int32 ). Isso significa que, quando a dupla precisão não é crucial, é um bom hábito declarar arrays na GPU como única precisão. Como alternativa, os primeiros sete construtores estáticos listados acima podem ser chamados através de sua função MATLAB correspondente anexando a lista de argumentos com gpuArr ay. Por exemplo. Cria um objeto gpuArray contendo 30723072 inteiros de 32 bits (int32) inicializados para zero. Ao chamar essas funções, uma alternativa para especificar explicitamente o tipo está usando o argumento similar. Isso cria uma matriz do mesmo tipo que o argumento que segue o argumento similar, ou seja, cria um objeto gpuArray de valores int32 inicializados para um, enquanto cria uma matriz MATLAB padrão de valores individuais inicializados para um. Isso pode ser útil ao criar novas variáveis ​​em funções que devem ser executadas tanto na CPU quanto na GPU onde não temos conhecimento a priori do tipo de entrada. Para que um objeto gpuArray possa manter números complexos, isso deve ser explicitamente especificado após a construção, usando o argumento complexo ao criá-lo diretamente na GPU ou por moldagem explícita ao copiar dados não complexos, e. GpuArray (complexo (im)). Para inspecionar as propriedades de uma matriz GPU, podemos usar as funções MATLAB padrão, como tamanho. comprimento . Isreal etc. Além disso, a função classUnderlying retorna a classe subjacente à matriz GPU (uma vez que a classe retornará apenas gpuArray) enquanto existeOnGPU retorna true se o argumento for uma matriz de GPU existente na GPU e estiver acessível. Uma vez que nossos dados de imagem estão na memória GPU, temos duas opções para manipulá-los: podemos usar o subconjunto das funções MATLAB internas (incluindo algumas funções disponíveis nas caixas de ferramentas) que são executadas no GPU, ou podemos escrever Nossas próprias funções usando apenas operações em elementos e lançá-las através de arrayfun ou bsxfun. In the first case, we use normal MATLAB syntax with the knowledge that any of these functions are automatically executed on the GPU as long as at least one argument is of type gpuArray . Using imGpu and randNoiseGpu defined above we can create a new, noisy image on the GPU by typing: A list of the GPU-enabled MATLAB functions available on the current system, together with all static constructors of gpuArray . can be displayed by typing methods(gpuArray) . For MATLAB 2013b, the list comprises around 200 standard functions plus any additional functions in installed toolboxes 2 . For example, using the GPU-enabled function imnoise from the IPT, the same result as above can be obtained through: (where a variance of 0.09 equals a standard deviation of 0.3). Another, potentially more useful, GPU-enabled function from the IPT is imfilter . Using imGpu from earlier filters the image using a horizontal Sobel filter. Note that sobelFilter is an ordinary 2D MATLAB array, 1 2 1 0 0 0 -1 -2 -1 . but since imGpu is a GPU array, the GPU-enabled version of imfilter is automatically called and the output, filteredImGpu . will be a GPU array. The second option for manipulating GPU arrays directly from MATLAB is by calling our own functions through the built-in bsxfun or arrayfun functions. As before, if any of the input arguments to the function is a GPU array, the calculation will automatically be carried out on the selected GPU. Thus, a third way of producing a noisy version of imGpu would be to first create the function addAndOffset that performs an element-wise addition of two images and adds a constant offset: and then calling The benefit of writing functions and launching them through bsxfun or arrayfun compared to calling MATLAB functions directly on GPU arrays is a potential speedup for large functions. This is because in the former case, the whole function can automatically be compiled to a GPU function, called CUDA kernel, resulting in a single launch of GPU code (although the overhead for compiling the function will be added to the first call). In contrast, in the latter case, each operation has to be launched in sequence using the precompiled kernel functions available in MATLAB. However, when running on a GPU, arrayfun and bsxfun are limited to element-wise operations. In a typical image processing application, this means that each pixel is unaware of its neighbours, which limits the use to functions where pixels are processed independently of one another. As a result, many image filters cannot be implemented in this way, in which case we are left either to use built-in functions as described earlier, or to write our own kernel functions as described in the next part. Further, since we are constrained to element-wise manipulations, the number of built-in functions at our disposal inside our function is somewhat limited. For a complete list of the available built-in functions, as well as some further limitations when using bsxfun and arrayfun with GPU arrays, see the PCT documentation 12 . Before moving on to the next part we should stress that since GPUs are built to process large amounts of data in parallel, there is no guarantee that running code on a GPU instead of a CPU will always result in a speedup. Although image processing algorithms provide good candidates for substantial speedups, this characteristic of the GPU means that vectorisation of code and simultaneous processing of large amounts of data (i. e. avoiding loops wherever possible) becomes even more crucial than in ordinary MATLAB programs. Further, GPU memory latency and bandwidth often limit the performance of GPU code. This can be alleviated by ensuring that, as far as possible, data that are operated on at the same time are stored nearby in memory. Since arrays are stored in a sequential column-major order in MATLAB, this means avoiding random memory-access patterns where possible and organising our data so that we mostly operate on columns rather than on rows. Finally, when evaluating the performance of our GPU code we should use the function gputimeit . It is called in the same manner as the regular MATLAB function timeit . i. e. it takes as argument a function, which itself does not take any arguments, and times it, but is guaranteed to return the accurate execution time for GPU code (which timeit is not). If this is not feasible, the code section we want to time can be sandwiched between a tic and a toc . as long as we add a call to wait(gpuDevice) just before the toc . This ensures that the time is measured only after execution on the currently selected GPU has finished. (Otherwise MATLAB continues to execute ensuing lines of GPU-independent code, like toc . asynchronously without waiting for the GPU calculation to finish). Since the MATLAB profiler only measures CPU time, we need to employ a similar trick to get accurate GPU timings when profiling: if we sandwich the lines or blocks of GPU code we want to time between two calls to wait(gpuDevice) . the execution time reported for the desired line or block plus the time taken by the second call to wait gives the correct timing. 3.2. Calling CUDA kernel functions from MATLAB By using the techniques described in the previous part we can use MATLAB to perform many of our standard image processing routines on the GPU. However, it does not allow us to test our own CUDA-implemented algorithms or use existing ones in our MATLAB programs, nor does it provide a means to explicitly control GPU resources, such as global and shared memory. In this part we demonstrate how this can be remedied by creating a CUDAKernel object from a kernel function written in CUDA. Instructions on how to write CUDA code is well beyond the scope of this chapter, and therefore this part assumes some previous basic knowledge of this language. A CUDAKernel object required to launch a CUDA kernel function from MATLAB is constructed by calling the static constructor parallel. gpu. CUDAKernel . The constructor takes three MATLAB string arguments: the path to a . ptx file containing the kernel code, the interface of the kernel function, and the name of the desired entry point. For the first argument, a . ptx file with the same name as the source file, containing the corresponding code compiled into pseudo-assembly language, is automatically generated when using the NVIDIA CUDA Compiler (NVCC) with the flag - ptx to compile a . cu file (if it contains one or more kernel functions). (Note that if using an integrated development environment, you might have to instruct it not to delete . ptx files when finishing the build for example Visual Studio 2010 requires that the flag - - keep is used.) The second argument can be either the . cu file corresponding to the . ptx file specified in the first argument, from which the argument list of the desired kernel function can be derived, or the explicit argument list itself. The latter is useful when the . cu file is not available, or when the argument list contains standard data types that have been renamed through the use of the typedef keyword. The third argument specifies which kernel function in the . ptx file to use, and although NVCC mangles function names similar to a C compiler, the names in the . ptx file are guaranteed to contain the original function name. MATLAB uses substring matching when searching for the entry point and it is therefore often enough to provide the name of the original kernel function (see exceptions below). Let us assume that we have access to myFilters. cu containing several kernel functions named myFilter1 . myFilter2 . etc. and its corresponding myFilters. ptx . Then creates a CUDAKernel object called gpuFilter1 that can be used to launch the CUDA kernel function myFilter1 from MATLAB. If we further assume that myFilter1 is declared as the second argument above, myFilters. cu . could equally be replaced by the string const float , float , float . In some cases, care has to be taken when specifying the entry point. For example, if myFilter2 is a templated function instantiated for more than one template argument, each instance of the resulting function will have a name containing the string myFilter2 . Likewise, if another kernel function called myFilter1v2 is declared in the same . cu file, specifying myFilter1 as the third argument of parallel. gpu. CUDAKernel becomes ambiguous. In these cases, we should provide the mangled function names, which are given during compilation with NVCC in verbose mode, i. e. with - - ptxas-options-v specified. The full mangled name of the kernel used by a CUDAKernel object is stored in the object property EntryPoint . and can be obtained by e. g. gpuFilter1.EntryPoint. Once a CUDAKernel object has been created, we need to specify its launch parameters which is done through the ThreadBlockSize . GridSize and SharedMemorySize object properties. Thus, sets the block size to 328 threads, the grid size to 96384 blocks and the shared memory size per block to 43281024 bytes, enough to hold 256 single or int32 values, or 128 double values. In total this will launch 30723072 threads, one per pixel of our sample image. A fourth, read-only property called MaxThreadsPerBlock holds the upper limit for the total number of threads per block that we can use for the kernel function. If the kernel function is dependent on constant GPU memory, this can be set by calling setConstantMemeory . taking as the first parameter the CUDAKerner object, as the second parameter the name of the constant memory symbol and as the third parameter a MATLAB array containing the desired data. For example, we can set the constant memory declared as constant float myConstData 128 in myFilters. cu and needed in myFilter1 To execute our kernel function we call feval with our GPUKernel object as the first parameter, followed by the input parameters for our kernel. For input parameters corresponding to arguments passed by value in the CUDA kernel (here: scalars), MATLAB scalars are normally used (although single element GPU arrays also work), whereas for pointer type arguments, either GPU arrays or MATLAB arrays can be used. In the latter case, these are automatically copied to the GPU. A list of supported data types for the kernel arguments can be found in the PCT documentation 13 . In general these are the CC standard types (along with their corresponding pointers) that have MATLAB counterparts, with the addition of float2 and double2 that map to MATLABs complex single and double types, respectively. CUDAKernel objects have three additional, read-only properties related to input and output. NumRHSArguments and MaxNumLHSArguments respectively hold the expected number of input arguments and the maximum number of output arguments that the objects accepts, and ArgumentTypes holds a cell of strings naming the expected MATLAB input types. Each type is prefixed with either in or inout to signal whether it is input only (corresponding to an argument passed by value or a constant pointer in CUDA) or combined inputoutput (corresponding to a non-constant pointer in CUDA). To function in a MATLAB context, a call to a CUDA kernel through a GPUKernel object must support output variables. Therefore, pointer arguments that are not declared const in the kernel declaration are seen as output variables and, if there is more than one, they are numbered according to their order of appearance in the declaration. This means that calling gpuFilter1 above produces one output, whereas gpuFilter2 created from produces two outputs. With this information we can now call the myFilter1 kernel function through Similarly, we can call myFilter2 through The output from a CUDAKernel object is always a GPU array, which means that if the corresponding input is a MATLAB array it will be created automatically. Consider the kernel with its corresponding CUDAKernel object gpuFilter3 . Since im from our previous examples a MATLAB array, calling automatically copies im to the GPU and creates a new gpuArray object called gpuRes3 to hold the result. 3.3. MEX functions and GPU programming In the previous part we saw how, by running our own CUDA kernels directly from MATLAB, we can overcome some of the limitations present when working only with built-in MATLAB functionality. However, we are still (in release 2013b) limited to using kernel functions that take only standard type arguments, and we can access neither external libraries, such as the NVIDIA Performance Primitives, the NVIDIA CUDA Fast Fourier Transform or the NVIDIA CUDA Random Number Generation libraries, nor the GPUs texture memory with its spatial optimised layout and hardware interpolation features. Further, we need to have an NVIDIA GPU, be writing our code in CUDA and have access to the PCT to use the GPU in our MATLAB code. In this section we look at how we can use MEX functions to partly or fully circumnavigate these limitations. This again assumes previous experience of GPU programming and some knowledge of how to write and use MEX functions. A good overview of the latter can be found in the MATLAB documentation 14 . The first option, which removes the technical restrictions but still requires access to the PCT (running on a 64-bit platform) and a GPU from NVIDIA, is to write MEX functions directly containing CUDA code. The CUDA code itself is written exactly as it would be in any other application, and can either be included in the file containing the MEX function entry point or in a separate file. Although this process is well documented in the PCT documentation 15 and through the PCT examples 16 , we briefly describe it here for consistency. The main advantage of this approach over the one described later is that it enables us to write MEX functions that use GPU arrays as input and output through the underlying CC object mxGPUArray . As all MEX input and output, GPU arrays are passed to MEX functions as pointers to mxArray objects. The first step is therefore to call either mxGPUCreateFromMxArray or mxGPUCopyFromMxArray on the pointer to the mxArray containing the GPU data, in order to obtain a pointer to an mxGPUArray . In the former case, the mxGPUArray becomes read-only, whereas in the latter case the data is copied so that the returned mxGPUArray can be modified. ( mxGPUCreateFromMxArray and mxGPUCopyFromMxArray also accept pointers to mxArray objects containing CPU data, in which case the data is copied to the GPU regardless of the function used.) We can now obtain a raw pointer to device memory that can be passed to CUDA kernels by calling mxGPUGetDataReadOnly (in the case of read-only data) or mxGPUGetData (otherwise) and explicitly casting the returned pointer to the correct type. Information about the number of dimensions, dimension sizes, total number of elements, type and complexity of the underlying data of an mxGPUArray object can be further obtained through the functions mxGPUGetNumberOfDimensions . mxGPUGetDimensions . mxGPUGetNumberOfElements . mxGPUGetClassID . and mxGPUGetComplexity . We can also create a new mxGPUArray object through mxGPUCreateGPUArray . which allocates and, if we want, initialises memory on the GPU for our return data. With this we are in a position where we can treat input from MATLAB just as any other data on the GPU and perform our desired calculations. Once we are ready to return our results to MATLAB we call either mxGPUCreateMxArrayOnCPU or mxGPUCreateMxArrayOnGPU to obtain a pointer to an mxArray object that can be returned to MATLAB through plhs . The former copies the data back to the CPU making the MEX function return a standard MATLAB array whereas in the latter case the data stays on the GPU and a GPU array is returned. Finally, we should call mxGPUDestroyGPUArray on any mxGPUArray objects that we have created. This deletes them from the CPU and, unless they are referenced by another mxArray object (as will be the case if they are returned through plhs ), frees the corresponding GPU memory. Note that there are several other functions in the mxGPU family to examine, duplicate, create and copy mxGPUArray objects, and in particular for working with complex arrays, that work in a similar way and which are described in the PCT documentation 17 . For the above to work, we need to include mxGPUArray. h . in addition to mex. h . in our source file. The source file has to have the extension . cu and it should contain a call to the function mxInitGPU before launching any GPU code in order to initialise the MATLAB GPU library. Provided that the environment variable MWNVCCPATH is set to the NVCC folder path and that a copy of the PCT version of mexopts. bat or mexopts. sh ( matlabroot toolboxdistcompgpuexternsrcmex xxx 64 ) is located in the same folder as the source code, we can compile our CUDA containing MEX functions from MATLAB in the usual way using the mex command. If using external libraries, these also have to be provided, which can normally be done by passing the full library path, including file name, to the mex command after the. c ..cpp or. cu file path. A bare-bone MEX function calling the kernel function myFilter1 from earlier, which takes into account the considerations above but is stripped of any boundary or argument checking, follows: After compilation, assuming the above code is found in mexFilter1.cu . we can call the MEX function like this: Note the explicit type conversion necessary to convert the second argument to single to match the input type. Note also that, depending on the compiler and system settings, in order to have mwSize correctly identified as a 64-bit type, we might need to use the largeArraysDims flag when compiling using the mex command. The rest of this part will be dedicated to describing how we can call GPU code from MATLAB even if we do not have access to the PCT or write our code in CUDA. Since without the PCT, the mex command is not set up to compile anything except standard CC code, we have two alternatives to achieve what we want. The only drawback with these, compared to when using mxGPUArray from the PCT is that it is harder to have data in GPU memory persist when returning to MATLAB between calls to MEX functions. (This can still be achieved by explicitly casting a GPU memory pointer to an integer type of correct length which is returned to MATLAB. The integer is then passed to the next MEX function which casts it back to the correct pointer type. However, the problem with this approach lies in eliminating the risk of memory leaks although solutions for this exist, they are beyond the scope of this chapter.) This means that unless we go through any extra effort, we are left either to perform all our GPU calculations from the same MEX function before returning control to MATLAB or suffer the penalty associated with copying our data back and forth between each call. In many cases, however, this may provide less of a limitation than one might initially think, especially when using MATLAB to test GPU code under development or using MATLAB as a front-end to existing code libraries. The first alternative is to compile our GPU code, regardless of in which language it is written, into a static library using our external compiler, and then to call this library from a MEX function that we compile using the mex command. Since our MEX function cannot call GPU functions directly, a small CC wrapper function has to be written around our GPU code. A wrapper for the myFilter1 CUDA kernel, which we can place either in the same file as the CUDA code or in a separate . cu file, could look like this (again, error checking has been omitted for brevity): Once the static library, normally a . lib file under Windows or a . a file under Linux and Mac OS, has been created, we can write a MEX function calling the wrapper: Note that it is normally better, especially if using a pre-existing library, to use include to include the corresponding header files rather than, as in the example above, to manually enter the function prototype. We should now be able to compile this wrapper using the mex command, remembering to pass the path to our own library as well as to the necessary GPU runtime library. For CUDA, this is cudart. lib on Windows, libcuda. so on Linux, or libcuda. dylib on Mac OS. Thus, assuming that the code above is found in myFilter1Mex. cpp and that the library is called myLibrary. lib . on Windows the call would look like: The second alternative is to build the MEX function directly in our external compiler, without going through the mex command. By doing this, the whole process can be carried out in one step and, if we wish, we are free to keep all our code in a single file. The main advantage of this approach occurs if we are developing or using an existing GPU library which we would like to call from MATLAB, for example for testing purposes. In such a case we can add the MEX file to our normal build routine so that every time we rebuild our library we automatically get the MEX function to call it from MATLAB. A MEX function is a dynamically linked shared library which exports a single entry point called mexFunction . Hence, by mimicking the steps carried out by the mex command (found in mexopts. bat or mexopts. sh ) when calling the external compiler, we can replicate this from outside MATLAB. We can view the exact steps carried out on a particular system by calling the mex command in verbose mode, i. e. using the v flag. Detailed information on how to build MEX functions for Windows and UNIX-like systems (Linux and Mac OS) can be also found in the MATLAB documentation 18 . Here we look at a minimal example from a Windows-centric point of view, although the procedure should be similar on other systems. First, we need to specify that we want to build a dynamically linked shared library, called a dynamic-link library on Windows, and give it the correct file extension, e. g. mexw64 second, we have to provide the compiler with the path to our MEX header files, normally mex. h third, we have to pass the libraries needed to build our MEX function, called libmx. lib . libmex. lib and libmat. lib on Windows, to the linker together with their path and finally, as mentioned above, we need to export the function named mexFunction . In Visual Studio 2010, all of the above steps are done in the Project Properties page of the project in question. Under Configuration properties-gt General . we set the Target Extension to . mexw64 and the Configutation Type to Dynamic Library (.dll) . For the header files we add (MATLABROOT)externinclude to Include Directories under Configuration Properties-gt VCDirectories . For the libraries, we add libmx. liblibmex. liblibmat. lib to Additional Dependencies under Configuration Properties-gt Linker-gt Input and (MATLABROOT)externlibwin64microsoft to Additional Library Directories under Configuration properties-gt Linker-gt General . Finally, we add export:mexFunction to Additional Options under Configuration Properties-gt Linker-gt Command Line . In the above steps it is assumed that the variable MATLABROOT is set to the path of the MATLAB root, otherwise we have to change (MATLABROOT) to, for example, C :Program FilesMATLAB2013b . The code for the MEX example, finally, would look something like this: In the case of an existing code library, we can add a new component (such as a new project to an existing solution in Visual Studio) containing a MEX function calling our desired GPU functions so that it is built directly with the original library without any additional steps. 4. Conclusion In conclusion, MATLAB is a useful tool for prototyping, developing and testing image processing algorithms and pipelines. It provides the user with the option of either using the functions of the IPT or leveraging the capabilities of a high-level programming language combined with many built-in standard functions to create their own algorithms. When high performance or processing of large amounts of data is required, the computational power of a GPU can be exploited. At this point, there are three main options for image processing on GPU in MATLAB: i) we can stick entirely to MATALB code, making use of the built-in, GPU-enabled functions in the IPT and the PCT as well as our own GPU functions built from element-wise operations only ii) we can use the framework provided by the PCT to call our own CUDA kernel functions directly from MATLAB iii) we can write a MEX function in CC that can be called from MATLAB and which in turn calls the GPU code of our choice. 5. Acknowledgements The authors would like to acknowledge the European Commission FP7 ENTERVISION programme, grant agreement no. 264552Calculate Moving Average Posted on April 28th, 2009 in Learn Excel - 191 comments Moving average is frequently used to understand underlying trends and helps in forecasting. MACD or moving average convergence divergence is probably the most used technical analysis tools in stock trading. It is fairly common in several businesses to use moving average of 3 month sales to understand how the trend is. Today we will learn how you can calculate moving average and how average of latest 3 months can be calculated using excel formulas. Calculate Moving Average To calculate moving average, all you need is the good old AVERAGE excel function . Assuming your data is in the range B1:B12, Just enter this formula in the cell D3 AVERAGE(B1:B3) And now copy the formula from D3 to the range D4 to D12 (remember, since you are calculating moving average of 3 months, you will only get 10 values 12-31) That is all you need to calculate moving average. Calculate Moving Average of Latest 3 Months Alone Lets say you need to calculate the average of last 3 months at any point of time. That means when you enter the value for the next month, the average should be automatically adjusted. First let us take a look at the formula and then we will understand how it works. So what the heck the above formula is doing anyway It is counting how many months are already entered 8211 COUNT(B4:B33) Then it is offsetting count minus 3 cells from B4 and fetching 3 cells from there 8211 OFFSET(B4,COUNT(B4:B33)-3,0,3,1). These are nothing but the latest 3 months. Finally it is passing this range to AVERAGE function to calculate the moving average of latest 3 months. Your Home Work Now that you have learned how to calculate moving average using Excel, here is your home work. Lets say you want the number of months used to calculate moving average to be configurable in the cell E1 . ie when E1 is changed from 3 to 6, the moving average table should calculate moving average for 6 months at a time. How do you write the formulas then Don8217t look at the comments, go and figure this out for yourself. If you cant find the answer, come back here and read the comments. Go This post is part of our Spreadcheats series. a 30 day online excel training program for office goers and spreadsheet users. Join today . Share this tip with your friends Hello, just recently found your website and Im loving all the tips. Thank you for all your tutorials. Its exactly I needed however, I ran into a bit a problem as I am also using Vlookup with Offset. For instance, in your example, I would use Vlookup in my template so that as I put in new data every month, it would automatically update the sales data each month. My problem is in my OFFSET formula, I have COUNTA which obviously counts any cells with formulas, even . Any ideas how to incorporate these two functions better, especially when I am trying to graph and average that last 12 months I would appreciate any ideas you or your readers my have. Thanks, again, for the the awesome site Twee. welcome to PHD and thanks for asking a question. I am not sure if I understood it correctly though. Have you tried using count instead of counta You havent shown us the offset formula, without looking that fixing it would be difficult. I need to calculate a 12-month rolling average that will encompass a 24 month period when completed. Can you point me in the right direction as too how to get started My data is vehivle miles and starts on B2 and ends on B25. Help Chandoo, this is a great formula for what I am using except I am trying unsuccessfully to make the formula conditional. I have a spreadsheet, see links below, that tracks all rounds of disc golf played by friends and myself. Ive already got it setup to calculate each of our overall averages and each of our averages on specific courses. What I am trying to do now however is also setup a moving average based off our 5 most recent rounds. Once more data has been entered I will change it to 10, but for now 5 will be just fine. I can get the moving average to work, but I cannot figure out how to add conditional restrictions. IE I want for example just the last 5 rounds that were played by Kevin. After that I will want just the last 5 rounds played by Kevin at the Oshtemo course. The code Im using is below. Code for Cell C9 is listed below. IF(B90,,IF(B9lt6,AVERAGEIF(DiscRoundsA2:A20000,A9,DiscRoundsM2:M20000),AVERAGE(OF FSET(DiscRoundsM2,IF(DiscRoundsA2:A20000A9,COUNT(DiscRoundsM2:M20000),quotquot)-5,0,5,1)))) Essentially if there are 0 rounds it leaves the cell blank. If there are 5 or fewer rounds it just uses the average of all rounds. Finally, if there are 6 or more rounds the code then uses your AVERAGE function from this post. After trying many things however I am uncertain how to conditionally pull the last 5 rounds so that it only pulls the last 5 rounds of the individual named in cell A9. The formula I am referencing is NOT currently in cell C9 on my spreadsheet that is linked. I just have been testing it there. DND: use the following formula in cell C13 onwards AVERAGE(B2:B13) and drag down. Hi, Im sure there is something listed above that is suppose to help, but Im still new to excel and am feeling overwhelmed. I just got a new job and Im tryin to make a good impression, so any help woud be great I have data for each month in 2009, 2010 and 2011 going across and multiple rows of this. Every month at the beginning of the month I need to calculate the sales of the previous year. Currently my formula is SUM(AG4:AR4)SUM(U4:AF4). Example: Current month is March. Info I need is sales total from March 2010-February 2011 divided by March 2009- February 2010 and it works great, but its too time consuming to have to change it every month. Is there a way I can get the formula to automatically change at the beginning of the month I dont know if I did a very good job explaining this or not. Congratulations on your new job. You can drag your formula sideways (to right for eg.) and it shows the s for next month automatically. No, what I need is for the formula to change each month. I have January 2009 through December 2011 boxes going across with data in them. IFERROR(SUM(AG4:AR4)SUM(U4:AF4),0) Next month I need for it go from calculating the sum of 0310 data to 0211 data divided by 0309 data to 0210 data and change to 0410 to 0311 data divided by 0409 data to 0311 data. IFERROR(SUM(AH4:AS4)SUM(V4:AG4),0) What I need is a formula that can refer to the current date and know that on the 1st of each month, it needs to switch the formulas over for the next previous 1-12 months divided by the previous 13-24 months. Im not sure if that makes sense. Basically I use this formula about 8 times on one sheet and I have about 200 sheets. Sorry for the double posting and thank you on the congrats What I need: If the current date is greater than the 1st of the month then the entire cell references to calculate the sales of prev year needs to move to the right by one column This is what Ive come up with. IF(P1gtN1,(SUM(AH4:AS4)SUM(V4:AG4))) p1 is current date n1 is 1st day of month AH4:AS4 is data from 0310-0211 V4:AG4 is data from 0309-0210 Part Im having issues with: How do i make it so that the formula knows exactly what 12 sections to grab and how to get to automatically change at the 1st of the month. Julie. You can use OFFSET formula to solve this. Assuming each column has one month, and first month is in C4 and current date is in P1 The above formula assumes that each column has months in Excel date format. You may want to tweak it until it produces right result. This is probably extremely simple and I am making it more complicated than I need to, but you wrote, The above formula assumes that each column has months in Excel date format. Ive been struggling to do this without having it turn my data into dates. Julie. What I meant is, the row number 4, where you have month names, should contain this data - 1-jan-2009 1-feb-2009 1-mar-2009 Also, I notice few errors in my formula. The correct formula should be, SUM(offset(C5,,datedif(C4,P1,m)1-12,1,12)) SUM(offset(C5,,datedif(C4,P1,m)1-24,1,12)) The above formula assumes dates are in row 4 and values are in row 5. I think that is exactly what I needed. Thank you thank you thank you so much My problem is very similar jasmins (61) and Azrold (74). I have disgusting amounts of data, from D:2 to D:61400 (and correspondingly in E and F, Ill have to do the same thing for these columns as well). Im trying to find the average for batches, such that D2:19, D20:37, D38:55 and so on - clumping 18 rows together and then finding the next average without re-using any previous row. Id also have to likely do this for every 19 and 20 clumps as well, but an example using 18 is fine. Could you annotate the formula you post Im a little confused on what the last 4 numbers mean in the COUNTA part. Thank you so much, this is going to make my life so much easier Laura This is easily done with Average and Offset . Assuming you are doing this in Col J and are averaging Col D J2: AVERAGE(OFFSET(D1,(ROW()-2)J11,,J1)) Where J1 will have the number 18 for a moving total of 18 numbers Copy down Row 2 will average Rows 2-19 Row 3 will average Rows 20-37 etc . You can also add labels in say Col H H2: Rows amp(ROW()-2)J12amp - amp(ROW()-1)J11 Copy down . I have mocked this up at: rapidsharefiles1923874899Averages. xlsx I am beginner trying to: 1. structure a spreadsheet that will then be used to 2. determine the optimal period for my moving average, within the range of a 5 day moving average to a 60 day moving average. Each cell represents the number of sales for that day, ranging from 0 to 100. I would prefer that each month of daily sales be in a new column. Currently I have 3 months of data, but obviously that will grow. So can you please tell me how to set up the spreadsheet and then the appropriate formulas (and their locations) Thank you very much, Hello again Hui, I am struggling yet again with the same spreadsheet you helped me with earlier. As beore, I have the following rows of monthly manually entered data: Volume of Calls Calls Answered age of calls abandoned Average handling time My line manager would now like 2 rows beneath these showing (by using formula): Average speed of answer Average abandoned time And as if that wasnt enough, she would like, for both rows, a summary cell at the end of the 12 months showing the yearly figure :( Many thanks again for any help you are able to give, I am using the vertical version for calculating a moving average. I am stumped when I need to calculate a 6-period moving average. My data starts in column c and the 6-period and 3-period averages are two columns to the right of the last period of data. I add a column for each month, so I currently adjust the formula manually each month: AVERAGE(EC8:EH8) My most recent attempt (that failed) is: AVERAGE(C6,COUNT(C6:EH6),-6,6,1) Please provide an explanation of why this didnt work when responding so I can understand how to create future f ormulas. Thank you so much, Kimber Kimber. Welcome to Chandoo. org and thanks for commenting. I think it is not a good idea to place averages in right most column as it keeps moving. Instead you could modify your sheet so that moving average is placed at left most column (and this will stay there even if you add extra columns to the right). No matter where the average cell is, you can use this formula to calculate the moving average. Afyter having read the whole of this thread I can see Im going to need a combination offset, match, count and averageif but Im not sure where. My problem is as follows: Each month there are over 100 people reporting activity - Column A is their name, Column B is the month, Column C is the year and Columns D through M is their activity in several categories. I need to find their 3 month and six month averages and display that in another worksheet although I could have them displayed in Columns N and O if needed. I use a pivot table to produce sums and total averages but it wont handle moving averages. Any pointers would be greatly appreciated. Thanks, Ben This will average the last MovAvg number of rows including itself (take out the -1 if you want it to not include itself). D75 is the cell that this formula is referencing (my data was very long) MovAvg is how big you want the moving average to be (I assigned this as a named cell (select the cell, Formulas --gt Defined Names --gt Define Name) You can make variable names in a spreadsheet to avoid always having to use rowcolumn.) This starts from the current cell (D75 in this case), goes up MovAvg-1 rows, over 0 columns, selects MovAvg nuber of rows, with 1 column. Passes this to the average function. Hi I read through every post, but havent been able to get this working correctly. How do we calculate the moving average of a percentage This is calculated weekly. Column A - accts met Column B - accts sold Column K - closing Column D - 2 week moving average of the closing Example of week 1 and week 2 Column A, row 7 is 25 and row 8 is 1 Column B, row 7 is 1 and row 8 is 1 Column K, row 7 formula is 125 (4) and row 8 is 11 (100) Column D - The formula in a prior post gives me an answer of 52 2 week avg, but thats not correct. it should be 226 (7) IF(ISERROR(AVERAGE(OFFSET(K7,COUNT(K7:K26)-2,0,2,1))),,AVERAGE(OFFSET(K7,COUNT(K7:K26)-2,0,2,1))) What do i need to change in that formula to use columns A amp B instead of the column K You are trying to average averages, which doesnt work. Try this simple formula beginning in D8: IF(ISBLANK(B8),,(B7B8)(A7A8)) Copy and paste the formula down to D26. This should give you a moving 2 week average. Remember to format column D as a percentage with how ever many decimal points you want. Im pretty much an excel neophyte. I just stumbled across your site amp am looking forward to perusing it at length in the months ahead. Im trying to calculate a 3 month moving average of expenses amp cannot figure out what I am doing wrong. Even after reading this article and the post on offset Im not sure I understand the formula. In my sandbox, I have: Column A - Months A2:A17Sept 2012 - Dec 2013 Column B - Total monthly expenses B2:B8 (B8 because March is the last completed month) - Those totals are 362599,372800,427317,346660,359864,451183,469681 Colum C - 3 Month Moving Average. I put the following formula in C4 (To start calculating in Nov of last year, just for grins). Since there are only three months in the data set at that point, I would assume it calculates the moving average of the first three months. The formula comes up with 469,681. When I average the first three months, I come up with 387,572. What am I doing wrong or misunderstanding Thanks for the help and for putting this website together. Hi Chandoo You have one really useful project here, tons of thanks In the very beginning of this thread Shamsuddin asked something similar to what I need, reverse calculation of values from the moving average. Maybe its stupid, but I cant come up with any ideas except for figure-by-figure lookup. If possible - please advice with this articles data, to get the concept. Actually, Id be happy to get anything, as google was of no use ) Once again - thank you so much for this site Im not really sure what you mean by reverse calculating a moving average Can you explain what your trying to doachieve Posting a sample file might help also Refer: chandoo. orgforumstopicposting-a-sample-workbook Hi Hui, I mean, I have a column of figures (e. g. monthly shipments), which are calculated as moving average based on another data set (e. g. monthly manufacturing output). Smth like this: (A1) Jan Feb Mar Apr May Jun Mfg Ship 100 500 450 600 600 700 Where Ship average(B2:C2) I know only shipments volumes, and have to find out respective mfg volumes. Generally speaking, the question is how we can find initial data with only MA on hand Suppose, this thread may not be the one for asking this (if you agree - maybe you know where to ask). Its just that Shamsuddins question was the most relevant result out of 10 google pages Mey To calculate the original data from a Moving Average (MA) you need two MAs eg a 9 and a 10 day MA or 1 MA and 1 piece of data From these you can recalculate the previous result But if you have a formula Average(B2:C2) you should have access to the data If it is a 2 day MA like your formula above MAAverage(B2:C2) MA(B2C2)2 if you know B2 C2(2MA)-B2 If you have a set of data you can share I can give a better solution Refer: chandoo. orgforumstopicposting-a-sample-workbook Great website. Forgive this question. I used to be an Expert in Lotus 123 decades ago, but I find Excel somewhat backwards in its progressions to Lotus 123, so I am starting over with Excel 2010. I am a logical person and I try to understand what the formulas do when I use them. I notice that there are not but 14 sales figures in column B, yet somehow we are counting from B4 to B33. I tested the formula out using: AVERAGE(OFFSET(B4,COUNT(B4:B14)-3,0,3,1)) and I get the same result as if I used AVERAGE(OFFSET(B4,COUNT(B4:B33)-3,0,3,1)). My first rule of old school spreadsheet creation is never to build a data table larger than the data provided if it is static (that is, not expanding in data). As a result, I have no real clue as to how OFFSET works. Is there a clear explanation of OFFSET with a singular example of it being used outside of the average and all by itself The reason I came here is to build a spreadsheet model that would use iterative calculations to find the best fit for profit data (that is maximizing profit) when the a short moving average of the cumulative profit curve (or equity curve) crosses OVER the longer term moving average of the equity curve. I find nothing that allows expansion of moving averages from 3 periods to say 100 periods (for both averages). By using the MA cross over to determine which trades to take, one can find an optimal level of profit to run the model from (which could be tweaked when the model is reoptimized). I can find nothing in most Excel books that cover this, and this kind of calculations should be relatively simple to pull off. Where could I find such information Thanks again for the wonderful website. Just in case you havent found it yet, heres a link for the OFFSET function: I have a question. I already have a 3 day moving average that I was given in my problem. Is it related to the average of stocks. The questions says that you have 1 stock that you PLAN on selling on day 10. My 3 day moving average is an integration from a, b where at and bt3 at any time. If you want to find the price you expect to sell the share for, do you integrate from 6,9 9,11 7,10. Do you want the far end of day 10, the middle of day 10, or leave day 10 out I am not sure what time frame to put this 3 day average between. Again, my function represents up to day 14, but I need the price at day 10. ivan Santos says: Im looking to see the moving average for a call center. im trying to find the index for every month for a full year. i only have 2 years worth of data and im wanting forecast out for 2014 in quarters. can i use this method for this I have a problem in average, I want to calculate the average of highlighted rows only in coloumn F on colomn G which also has highlighted blank cells Hi, I am working on a spreadsheet that has the past four years of weekly data but the current years data is incomplete as it only gets entered each week. Is there a way of setting up a formula that will calculate an average based on the number of weeks that have data in them For eg. in the middle of the year it will create an average based on cells 2-27 26 but the next week it would be cells 2-28 27. Its doing my head in and I dont want to have to manually adjust the average every week. Great site by the way Very helpful. ) Rosie Yes this can be done Can you please ask the question at the Forums and attach a sample file chandoo. orgforum Ok here is my question that has been plaguing me for the last 2 12 months and I havent found a solution anywhere on the web: I have a sales team and I need a moving avg but with a fix format and a shifting date rage that is fixed as well. i. e. Sales person 1115 2115 3115 12114 11114 10114 ME 1 2 0 4 5 6 What I am trying to do is this: Lets say today date is 3115 I need a way to go back 3 (6 and 12 as well) months from the current date and avg the sales numbers. The hard part is I would like to just change the year of the dates so I dont have to mess with the format or if I hire(fire) someone. So in the above example I would have the formula take the 6 1 2 (9)3 3 but then as time would go on this would keep going but once the new year began in JAN 2016 it would have to use the figures from the past 2015 data (3,6 and 12 Month rolling avgs). I hope that this clear and I would love to get some help with this. Thank you in advance. Can you please ask the question in the Chandoo. org Forums at: forum. chandoo. org Attach a sample file to simplify the process Ok I have posted to the forums and uploaded a sample file. 8230 Calculate Moving Average Chandoo. org 8211 Learn Moving average is frequently used to understand underlying trends and helps in forecasting. MACD or moving average convergence divergence is probably the 8230 Amelia McCabe says: Looking for a little help. I have tried what I think is a modified version of this formula that is not really working. I have a row of data (one number per month) that I need a continuous average for based on the number of months of entered data not on 12 months. Data are in cells b53 to m53. So I tried to modify this formula as follow (it did not work) and I wonder if I can use this formula this way at all since my data is in a row not a column. AVERAGE(OFFSET(B53COUNT(B53:M53)-12,0,1,12)). Have also tried the arguments as 0,0,1,12 and -1,0,1,12. Please help me understand if I am up the totally wrong tree or just on the wrong branch. Amelia Without seeing the data id suggest that AVERAGE(OFFSET(B53,COUNT(B53:M53)-12,0,1,12)) should be: AVERAGE(OFFSET(B53. 1,COUNT(B53:M53))) One issue with the original formula is that there are 12 cells between B53:M53, If only 5 have data in them, then you take 12 away, the offset is trying to offset B53, a negative 7 columns, which will force an error You may also be able to use the Averageifs function Possibly: Averageifs(B53:M53,B53:M53,0) Are you able to post a sample file in the Chandoo. org Forums forum. chandoo. org