"Fossies" - the Fresh Open Source Software Archive

Member "sip-0.12.1/src/morphology_int.c" (12 Mar 2012, 8101 Bytes) of package /linux/privat/sip-0.12.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "morphology_int.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 0.5.6_vs_0.12.1.

    1 /*
    2  * -------------------------------------------------------------------------
    3  * SIP - Scilab Image Processing toolbox
    4  * Copyright (C) 2002-2009  Ricardo Fabbri
    5  *
    6  * This program is free software; you can redistribute it and/or modify
    7  * it under the terms of the GNU General Public License as published by
    8  * the Free Software Foundation; either version 2 of the License, or
    9  * (at your option) any later version.
   10  * 
   11  * This program is distributed in the hope that it will be useful,
   12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
   13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   14  * GNU General Public License for more details.
   15  * 
   16  * You should have received a copy of the GNU General Public License
   17  * along with this program; if not, write to the Free Software
   18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
   19  * -------------------------------------------------------------------------
   20  */ 
   21   
   22 
   23 #include <stack-c.h>
   24 #include <stdio.h>
   25 #include <string.h>
   26 #include <animal/animal.h>
   27 #include "sip_common.h"
   28 
   29 
   30 /*----------------------------------------------------------
   31  * [img]=edilate(img [, radius, opt, pad, method])
   32  *    radius defaults to 5
   33  *    opt is "full" if result is not trimmed 
   34  *    opt is "same" if result is trimmed to original size
   35  * @@@@@@@@@
   36  * TODO
   37  * - use IFT fast euclidean propagation
   38  * - work with integers
   39  * - interface must be similar to matlab 
   40  * - fazer opcao: pad/unpad? como matlab faz isto?
   41  *----------------------------------------------------------*/
   42 #define MAX_OPT 10
   43 SipExport int 
   44 edilate_int(char *fname)
   45 {
   46    /* Interface variables */ 
   47    int   r1, c1, p1,  
   48          r2, c2, p2,  
   49          r3, c3, p3,  
   50          r_alg, c_alg, l_alg,
   51          ro, co,
   52          i, irad,
   53          nv=1,
   54          minlhs=1, maxlhs=1, minrhs=1, maxrhs=4;
   55    ImgPUInt32 *im, *tmp;
   56    Img *im_int, *result, *tmp_int;
   57    double *po, radius=5;
   58    char opt[MAX_OPT]="same", *str;
   59    bool stat;
   60    dt_algorithm alg=DT_CUISENAIRE_PMN_1999;
   61 
   62    CheckRhs(minrhs,maxrhs);
   63    CheckLhs(minlhs,maxlhs);
   64 
   65 
   66    GetRhsVar(nv++, "d", &r1, &c1, &p1);  // img 
   67    if (Rhs == 2) {
   68       GetRhsVar(nv++, "d", &r2, &c2, &p2);  // radius
   69       radius = *stk(p2);
   70    } else if (Rhs >= 3) {
   71       GetRhsVar(nv++, "d", &r2, &c2, &p2);  // radius
   72       radius = *stk(p2);
   73       GetRhsVar(nv++, "c", &r3, &c3, &p3);  // opt
   74       strncpy(opt,cstk(p3),MAX_OPT);
   75       if (Rhs == 4) {
   76         GetRhsVar(nv++, "c", &r_alg, &c_alg, &l_alg);
   77         str = cstk(l_alg);
   78         if ( strncasecmp("exact dilations",str,8) == 0 || 
   79                 strcmp("costa-estrozi",str) == 0)
   80           alg=DT_EXACT_DILATIONS;
   81         else if (strcasecmp("cuisenaire pmn",str) == 0)
   82           alg=DT_CUISENAIRE_PMN_1999;
   83         else
   84                 sip_error("invalid second argument -- unknown method");
   85       }
   86    }
   87 
   88    if (alg == DT_EXACT_DILATIONS) {
   89      /* @@@ maybe there's a better way of passing data */
   90      im_int=new_img(c1, r1);
   91      if (!im_int) sip_error("unable to alloc memory");
   92      sci_2D_double_matrix_to_animal(p1,r1,c1,im_int,pixval,1);
   93 
   94      irad   = (int) ceil(radius);
   95      result = new_img(im_int->rows+2*irad, im_int->cols+2*irad);
   96      if (!result) sip_error("unable to alloc memory");
   97 
   98      edilate_np(result,im_int,radius);    
   99 
  100      if (strcmp(opt,"same") == 0) {
  101         tmp_int = result;
  102         result = imtrim (result, irad, irad);
  103         imfree (&tmp_int);
  104      }
  105      imfree(&im_int);
  106      im_int = result;
  107 
  108      ro=im_int->cols; co=im_int->rows;
  109      stat = animal_grayscale_image_to_double_array(fname, im_int, &po);
  110      if (!stat) return false;
  111      imfree(&im_int);
  112 
  113    } else { /* Fast exact Euclidean propagation algorithm */
  114 
  115      /* @@@ maybe there's a better way of passing data */
  116      im = new_img_puint32(c1, r1);
  117      if (!im) sip_error("unable to alloc memory");
  118      sci_2D_double_matrix_to_animal(p1,r1,c1,im,pixval,1);
  119 
  120      for (i=0; i<r1*c1; ++i)
  121         DATA(im)[i] = (PROUND(pixval,*stk(p1+i)) == 0);
  122 
  123      im->isbinary = true;
  124 
  125      if (strcmp(opt,"same") != 0) {
  126        irad = (int) ceil(radius);
  127        tmp = impad_puint32(im, irad, irad, 1);
  128        if (!tmp) sip_error("unable to alloc memory");
  129        imfree_puint32(&im);
  130        im = tmp;
  131      }
  132 
  133      stat = distance_transform_ip_max_dist(im, alg, (int)ceil(radius*radius), false, NULL);
  134      if (!stat) return false; /* @@@ garbage collection */
  135 
  136      ro=im->cols; co=im->rows;
  137 
  138      for (i = 0; i <ro*co; i++)
  139        im->data[i] = (im->data[i] <= radius*radius);
  140      im->isbinary = true;
  141 
  142      stat = animal_grayscale_imgpuint32_to_double_array(fname, im, &po);
  143      if (!stat) return false;
  144      imfree_puint32(&im);
  145    }
  146 
  147    CreateVarFromPtr(nv, "d", &ro, &co, &po);
  148 
  149    /*  Return variables  */
  150    LhsVar(1) = nv;
  151 
  152    free(po);
  153    return true;
  154 }
  155 
  156 /*----------------------------------------------------------------
  157  * out = watershed(img [, markers, nhood])
  158  *
  159  * Description:  Watershed transform for image segmentation
  160  *    - also known as "Sup Reconstruction".
  161  *    - it is better to calculate gradient before using it without markers.
  162  *    - its core is easily extended for chromatic images
  163  * 
  164  * INPUT
  165  *    - im: grayscale image
  166  *    - nhood: connectivity. Must be 4 or 8.
  167  *    - markers: image with the markers (seeds). Each mark must have a
  168  *    unique label from 1 to N, where N is the number of marks.
  169  *       + if markers is "-1" or is omitted, then the regional minima of 
  170  *       "im" will be taken as the markers.
  171  *
  172  * OUTPUT 
  173  *    - an image with the watershed regions, each with a unique number.
  174  *    If minima != NULL, the regions will have the same label as the
  175  *    corresponding supplied markers.
  176  *
  177  * ALGORITHM
  178  *    - Image Foresting Transform (see ift.h)
  179  *
  180  * TODO
  181  *    - allow arbitrary neighborhood (structuring element)
  182  *
  183  * out = watershed(img [, markers, nhood])
  184  *----------------------------------------------------------------*/
  185 SipExport int
  186 watershed_int(char *fname)
  187 {
  188    int   rim, cim, pim, 
  189          r_markers, c_markers, p_markers,
  190             r_nhood, c_nhood, p_nhood,
  191          i, nv=1,
  192          minlhs=1, maxlhs=1, minrhs=1, maxrhs=3;
  193    double *pt;
  194    Img *im;
  195    ImgPUInt32 *markers=NULL, *regions;
  196    nhood *nh;
  197    bool stat;
  198 
  199    CheckRhs(minrhs,maxrhs); CheckLhs(minlhs,maxlhs);
  200 
  201    GetRhsVar(nv++, "d", &rim, &cim, &pim);  // img 
  202 
  203    if (Rhs > 1) {
  204       GetRhsVar(nv++, "d", &r_markers, &c_markers, &p_markers);
  205       if ( *stk(p_markers) != -1 )  {
  206          if (r_markers != rim || c_markers != cim)
  207             sip_error("image and Marker must be the same size");
  208             
  209          markers = new_img_puint32(r_markers,c_markers);
  210          if (!markers) sip_error("not enough memory");
  211 
  212          /* pass transposed image to internal row-wise storage */
  213          for (i=0; i<r_markers*c_markers; i++)
  214             DATA(markers)[i] = (puint32) *stk(p_markers + i);
  215       } 
  216 
  217       if (Rhs == 3) {
  218          GetRhsVar(nv++, "d", &r_nhood, &c_nhood, &p_nhood);
  219          if (r_nhood * c_nhood != 1)
  220             sip_error("nhood must be a scalar");
  221          switch ( (int)*stk(p_nhood) ) {
  222             case 4:
  223                nh = get_4_nhood();
  224                break;
  225             case 8:
  226                nh = get_8_nhood();
  227                break;
  228             default:
  229                sip_error("invalid 3rd parameter (neighborhood)");
  230          }
  231       } else
  232          nh = get_8_nhood();
  233    } else
  234       nh = get_8_nhood();
  235 
  236    /* pass transposed image to internal row-wise storage */
  237    im=new_img(cim,rim); 
  238       if (!im) sip_error("not enough memory");
  239    sci_2D_double_matrix_to_animal(pim,rim,cim,im,pixval,PIXVAL_MAX);
  240 
  241    /* call Watershed */
  242    regions = watershed(im, markers, nh);
  243       if (!regions) sip_error("problem inside watershed C subroutine");
  244 
  245    imfree(&im);
  246    if (markers != NULL)
  247       imfree_puint32(&markers);
  248    free_nhood(&nh);
  249 
  250    /* return to Scilab */
  251    stat = animal_grayscale_imgpuint32_to_double_array(fname, regions, &pt);
  252    if (!stat) return false;
  253 
  254    CreateVarFromPtr(nv, "d", &rim, &cim, &pt);
  255    LhsVar(1) = nv;
  256    free(pt);
  257 
  258    return true;
  259 }