"Fossies" - the Fresh Open Source Software Archive

Member "sip-0.12.1/src/analysis_int.c" (26 Mar 2012, 10442 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 "analysis_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  * [skl,dt,lbl]=skel(img [, opt, algorithm])
   32  *  - opt = "interior" or "full"(default)
   33  * - algorithm = "exact euclidean" or "fast euclidean"(default)
   34  *--------------------------------------------------------------*/
   35 SipExport int 
   36 skel_int(char *fname)
   37 {
   38    /* Interface variables */ 
   39    int   r1, c1, p1,  // img
   40          r2, c2, // skl
   41          r3, c3, // dt
   42          r4, c4, // lbl
   43          r_opt, c_opt, l_opt, opt=INTERIOR,
   44             r_alg, c_alg, l_alg, alg=SKL_IFT,
   45          i,
   46          minlhs=1, maxlhs=3, minrhs=1, maxrhs=3;
   47    Img *im; //*aux;
   48    annotated_skl *skl;
   49 
   50    double *pt2, *pt3, *pt4;
   51 
   52    CheckRhs(minrhs,maxrhs);
   53    CheckLhs(minlhs,maxlhs);
   54 
   55 
   56    GetRhsVar(1, "d", &r1, &c1, &p1);  // img 
   57    if (Rhs >= 2) {
   58       GetRhsVar(2, "c", &r_opt, &c_opt, &l_opt);
   59       if (strcmp("interior",cstk(l_opt)) == 0)
   60          opt = INTERIOR;
   61       else if (strcmp("both",cstk(l_opt)) == 0)
   62          opt = BOTH;
   63       else if (strcmp("exterior",cstk(l_opt)) == 0)
   64          opt = EXTERIOR;
   65       else
   66          sip_error("invalid second argument");
   67         if (Rhs == 3) {
   68             GetRhsVar(3,"c",&r_alg,&c_alg,&l_alg);
   69             if (strcmp("exact euclidean",cstk(l_alg)) == 0)
   70                 alg=SKL_COSTA_ESTROZI;
   71             else if (strcmp("fast euclidean",cstk(l_alg)) == 0)
   72                 alg=SKL_IFT;
   73             else
   74                 sip_error("invalid third argument");
   75         }
   76    }
   77 
   78 
   79    // pass transposed image to internal row-wise storage
   80    im=new_img(c1,r1); 
   81 
   82    sci_2D_double_matrix_to_animal(p1,r1,c1,im,pixval,1);
   83     
   84    im->isbinary = true;
   85 
   86    // regularize image
   87    //aux = im;
   88    //im = imregularize(im,NULL);
   89    //imfree(&aux);
   90 
   91    skl = msskl(im,opt,alg);
   92    if (!skl)
   93       sip_error("internal msskl routine failed");
   94    imfree(&im);
   95 
   96    /* @@@ change here for integer data */
   97    pt2 = (double *)calloc(r1*c1, sizeof(double));
   98       if (!pt2)
   99          sip_error("unable to alloc memory");
  100    pt3 = (double *)calloc(r1*c1, sizeof(double));
  101       if (!pt3)
  102          sip_error("unable to alloc memory");
  103    pt4 = (double *)calloc(r1*c1, sizeof(double));
  104       if (!pt4)
  105           sip_error("unable to alloc memory");
  106 
  107    for (i=0; i<r1*c1; i++) {
  108       pt2[i] = skl->skl->data[i];
  109       pt3[i] = skl->dt->data[i];
  110       pt4[i] = skl->lbl->data[i];
  111    }
  112 
  113    free_ann_skl(&skl);
  114    r2 = r3 = r4 = r1;
  115    c2 = c3 = c4 = c1;
  116    CreateVarFromPtr(2, "d", &r2, &c2, &pt2);
  117    CreateVarFromPtr(3, "d", &r3, &c3, &pt3);
  118    CreateVarFromPtr(4, "d", &r4, &c4, &pt4);
  119 
  120    free(pt2);
  121    free(pt3);
  122    free(pt4);
  123 
  124    /*  Return variables  */
  125    LhsVar(1) = 2;
  126    LhsVar(2) = 3;
  127    LhsVar(3) = 4;
  128    return true;
  129 }
  130 
  131 /*----------------------------------------------------------------
  132  * out = thin(img)
  133  * Simple thinning algorithms, only Zhang-Suen supported for now.
  134  * The "skel" routine provide better and more complex algorithms.
  135  * In this funcion only "boundary deletion" algorithms shall be
  136  * provided.
  137  *----------------------------------------------------------------*/
  138 SipExport int 
  139 thin_int(char *fname)
  140 {
  141    int   rim, cim, pim,  // img
  142          i,s,
  143          minlhs=1, maxlhs=1, minrhs=1, maxrhs=1;
  144    double *pt;
  145    bool stat;
  146    Img *im;
  147 
  148    CheckRhs(minrhs,maxrhs);
  149    CheckLhs(minlhs,maxlhs);
  150 
  151    GetRhsVar(1, "d", &rim, &cim, &pim);  // img 
  152 
  153    // pass transposed image to internal row-wise storage
  154    im=new_img(cim,rim); 
  155    sci_2D_double_matrix_to_animal(pim,rim,cim,im,pixval,1);
  156     
  157    im->isbinary = true;
  158 
  159    s= thinzs_np(im);
  160       if (!s) sip_error("thin: problem inside thinzs_np C subroutine");
  161 
  162 
  163 
  164    stat = animal_grayscale_image_to_double_array(fname,im,&pt);
  165    if (!stat) return false;
  166    imfree(&im);
  167 
  168    CreateVarFromPtr(2, "d", &rim, &cim, &pt);
  169    LhsVar(1) = 2;
  170 
  171    free(pt);
  172    return true;
  173 }
  174 
  175 /*----------------------------------------------------------------
  176  * i=percol(img [, direction])
  177  *  i is 1 if image is percolated; 0 otherwise. Direction is 1 if
  178  *  percolation is to be tested horizontally; 0 if it is to be
  179  *  tested vertically.
  180  *----------------------------------------------------------------*/
  181 
  182 SipExport int 
  183 percol_int(char *fname)
  184 {
  185    int   r1, c1, p1,  
  186          r2, c2, p2,  
  187          ro, co, po, nv,
  188          i, direction, stat,
  189          minlhs=1, maxlhs=1, minrhs=1, maxrhs=2;
  190    Img *im; 
  191 
  192    CheckRhs(minrhs,maxrhs);
  193    CheckLhs(minlhs,maxlhs);
  194 
  195 
  196    nv = 1;
  197    GetRhsVar(nv++, "d", &r1, &c1, &p1);  // img 
  198    if (Rhs == 1)
  199       direction = 0;
  200    else {
  201       GetRhsVar(nv++, "d", &r2, &c2, &p2);  // direction
  202       direction = ((int)*stk(p2) != 1);
  203    }
  204 
  205    im=new_img(c1, r1);
  206    sci_2D_double_matrix_to_animal(p1,r1,c1,im,pixval,1);
  207 
  208    stat = percol(im,direction);
  209    imfree(&im);
  210 
  211    ro=co=1;
  212    CreateVar(nv, "d", &ro, &co, &po);
  213    *stk(po) = stat;
  214 
  215    LhsVar(1) = nv;
  216    return true;
  217 }
  218 
  219 /*----------------------------------------------------------------
  220  * out = bwdist(img [,method, max_dist])
  221  * Distance transforms. 
  222  *
  223  * "method" may be: 
  224  *    - "euclidean" : default euclidean method (Lotufo-Zampirolli)
  225  *    - "lotufo-zampirolli" (fast exact euclidean)
  226  *    - "costa-estrozi" (exact euclidean)
  227  *    - "IFT" - Image Foresting Transform (fast very accurate euclidean)
  228  *
  229  * In the future:
  230  *    - "chessboard"
  231  *    - "chamfer"
  232  *    - (...)
  233  *
  234  * TODO
  235  *    - add an output Label parameter (discrete Voronoi diagram)
  236  *----------------------------------------------------------------*/
  237 
  238 SipExport int 
  239 bwdist_int(char *fname)
  240 {
  241    int   rim, cim, pim,  // img
  242             r_alg, c_alg, l_alg,
  243             r_d, c_d, l_d,
  244          i, nv=1,
  245          minlhs=1, maxlhs=2, minrhs=1, maxrhs=4;
  246 
  247    dt_algorithm alg=DT_MEIJSTER_2000;
  248    double *pt, *pt_lbl, max_dist = (double)((puint32) -1);
  249    char *str;
  250    bool noexec=false, stat, is1const, use_label = false;
  251    Img *im;
  252    ImgPUInt32 *dt, *imlabel=NULL;
  253 
  254    CheckRhs(minrhs,maxrhs); CheckLhs(minlhs,maxlhs);
  255 
  256    if (Lhs >= 2) {
  257      use_label = true;
  258      if (Rhs == 1)
  259        alg = DT_MAURER2003;
  260    }
  261 
  262    GetRhsVar(nv++, "d", &rim, &cim, &pim);  // img 
  263    if (Rhs >= 2) {
  264       GetRhsVar(nv++, "c", &r_alg, &c_alg, &l_alg);
  265       str=cstk(l_alg);
  266       if (strcasecmp("lotufo-zampirolli",str) == 0)
  267          alg=DT_LOTUFO_ZAMPIROLLI;
  268       else if ( strncasecmp("exact dilations",str,8) == 0 || 
  269               /* backward-compat:*/ strcmp("costa-estrozi",str) == 0)
  270          alg=DT_EXACT_DILATIONS;
  271       else if ( strcasecmp("IFT",str) == 0 || 
  272                 strncasecmp("IFT 8",str,5) == 0)
  273          alg=DT_IFT;
  274       else if ( strncasecmp("IFT 4",str,5) == 0)
  275          alg=DT_IFT_4;
  276       else if ( strncasecmp("meijster",str,3) == 0)
  277          alg=DT_MEIJSTER_2000;
  278       else if ( strncasecmp("maurer",str,3) == 0)
  279          alg=DT_MAURER2003;
  280       else if ( strncasecmp("euclidean",str,6) == 0 || 
  281                 strcasecmp("cuisenaire pmn",str) == 0)
  282          alg=DT_CUISENAIRE_PMN_1999;
  283       else if ( strncasecmp("cuisenaire pmon",str,15) == 0)
  284          alg=DT_CUISENAIRE_PMON_1999;
  285       else if ( strncasecmp("cuisenaire psn4",str,15) == 0)
  286          alg=DT_CUISENAIRE_PSN4_1999;
  287       else if ( strncasecmp("cuisenaire psn8",str,15) == 0)
  288          alg=DT_CUISENAIRE_PSN8_1999;
  289       else if ( strncasecmp("noexec",str,5) == 0) 
  290          noexec = true;
  291          /* undocumented option used to see how much overhead does
  292           * this interface function imposes for a particular image */
  293       else
  294               sip_error("invalid second argument -- unknown or unimplemented method");
  295       if (Rhs >= 3) {
  296         GetRhsVar(nv++, "d", &r_d, &c_d, &l_d);
  297         max_dist = *stk(l_d);
  298       }
  299 
  300       if (use_label && alg!=DT_MAURER2003)
  301         sip_error("such algorithm choice does not currently support label");
  302    }
  303    // pass transposed image to internal row-wise storage
  304    im = new_img(cim,rim); 
  305 
  306    sci_2D_double_matrix_to_animal(pim,rim,cim,im,pixval,1);
  307 
  308    is1const=true;
  309    for (i=0; i<rim*cim; ++i)
  310       if (DATA(im)[i]==0) {
  311          is1const=false;
  312          break;
  313       }
  314    if (is1const) {
  315       sip_warning("the input image is constant and different than 0");
  316       sip_warning("the distance transform is undefined for this case");
  317    }
  318 
  319    im->isbinary = true;
  320 
  321    if (noexec) {
  322       dt = new_img_puint32(im->rows, im->cols);
  323       if (use_label)
  324         imlabel = new_img_puint32(cim,rim); 
  325    } else {
  326      if (max_dist == (double)(puint32)-1) {
  327        if (use_label) {
  328          dt = distance_transform_label(im, alg, true, &imlabel);
  329        } else
  330          dt = distance_transform(im, alg);
  331      } else {
  332        if (use_label)
  333          sip_error("label + max_dist not yet implemented.");
  334        dt = distance_transform_max_dist(im, alg, max_dist*max_dist, false, &imlabel);
  335      }
  336    }
  337 
  338    if (!dt) sip_error("problem inside distance_transform C subroutine");
  339    if (use_label && !imlabel) sip_error("problem with label inside distance_transform C subroutine");
  340    
  341    imfree(&im);  /* FIXME: use better err treatment */
  342 
  343    stat = animal_grayscale_imgpuint32_to_double_array(fname,dt,&pt);
  344    if (!stat) return false;
  345    imfree_puint32(&dt);
  346 
  347    CreateVarFromPtr(nv, "d", &rim, &cim, &pt);
  348    LhsVar(1) = nv++;
  349    free(pt);
  350 
  351    if (use_label) {
  352      stat = animal_grayscale_imgpuint32_to_double_array(fname,imlabel,&pt_lbl);
  353      if (!stat) return false;
  354 
  355      imfree_puint32(&imlabel);
  356      CreateVarFromPtr(nv, "d", &rim, &cim, &pt_lbl);
  357      LhsVar(2) = nv++;
  358      free(pt_lbl);
  359    }
  360 
  361    return true;
  362 }